ap.h   ap.h 
skipping to change at line 41 skipping to change at line 41
#include <list.h> #include <list.h>
#include <vector.h> #include <vector.h>
#else #else
#include <list> #include <list>
#include <vector> #include <vector>
#endif #endif
/******************************************************************** /********************************************************************
Array bounds check Array bounds check
********************************************************************/ ********************************************************************/
#define AP_ASSERT #define NO_AP_ASSERT
#ifndef AP_ASSERT // #ifndef AP_ASSERT //
#define NO_AP_ASSERT // This code avoids definition of the #define NO_AP_ASSERT // This code avoids definition of the
#endif // both AP_ASSERT and NO_AP_ASSERT symbols #endif // both AP_ASSERT and NO_AP_ASSERT symbols
#ifdef NO_AP_ASSERT // #ifdef NO_AP_ASSERT //
#ifdef AP_ASSERT // #ifdef AP_ASSERT //
#undef NO_AP_ASSERT // #undef NO_AP_ASSERT //
#endif // #endif //
#endif // #endif //
skipping to change at line 182 skipping to change at line 182
const double abscomplex(const complex &z); const double abscomplex(const complex &z);
const complex conj(const complex &z); const complex conj(const complex &z);
const complex csqr(const complex &z); const complex csqr(const complex &z);
/******************************************************************** /********************************************************************
Templates for vector operations Templates for vector operations
********************************************************************/ ********************************************************************/
#include "apvt.h" #include "apvt.h"
/******************************************************************** /********************************************************************
BLAS functions Level 1 BLAS functions
********************************************************************/
double vdotproduct(const double *v0, int stride0, const double *v1, int str
ide1, int n);
complex vdotproduct(const complex *v0, int stride0, const char *conj0, cons
t complex *v1, int stride1, const char *conj1, int n);
void vmove(double *vdst, int stride_dst, const double* vsrc, int stride_s
rc, int n);
void vmove(complex *vdst, int stride_dst, const complex* vsrc, int stride_s
rc, const char *conj_src, int n);
void vmoveneg(double *vdst, int stride_dst, const double* vsrc, int strid
e_src, int n);
void vmoveneg(complex *vdst, int stride_dst, const complex* vsrc, int strid
e_src, const char *conj_src, int n);
void vmove(double *vdst, int stride_dst, const double* vsrc, int stride_s
rc, int n, double alpha);
void vmove(complex *vdst, int stride_dst, const complex* vsrc, int stride_s
rc, const char *conj_src, int n, double alpha);
void vmove(complex *vdst, int stride_dst, const complex* vsrc, int stride_s
rc, const char *conj_src, int n, complex alpha);
void vadd(double *vdst, int stride_dst, const double *vsrc, int stride_sr
c, int n);
void vadd(complex *vdst, int stride_dst, const complex *vsrc, int stride_sr
c, const char *conj_src, int n);
void vadd(double *vdst, int stride_dst, const double *vsrc, int stride_sr
c, int n, double alpha);
void vadd(complex *vdst, int stride_dst, const complex *vsrc, int stride_sr
c, const char *conj_src, int n, double alpha);
void vadd(complex *vdst, int stride_dst, const complex *vsrc, int stride_sr
c, const char *conj_src, int n, complex alpha);
void vsub(double *vdst, int stride_dst, const double *vsrc, int stride_sr
c, int n);
void vsub(complex *vdst, int stride_dst, const complex *vsrc, int stride_sr
c, const char *conj_src, int n);
void vsub(double *vdst, int stride_dst, const double *vsrc, int stride_sr
c, int n, double alpha);
void vsub(complex *vdst, int stride_dst, const complex *vsrc, int stride_sr
c, const char *conj_src, int n, double alpha);
void vsub(complex *vdst, int stride_dst, const complex *vsrc, int stride_sr
c, const char *conj_src, int n, complex alpha);
void vmul(double *vdst, int stride_dst, int n, double alpha);
void vmul(complex *vdst, int stride_dst, int n, double alpha);
void vmul(complex *vdst, int stride_dst, int n, complex alpha);
/********************************************************************
Obsolete BLAS functions
********************************************************************/ ********************************************************************/
double vdotproduct(const double *v1, const double *v2, int N); double vdotproduct(const double *v1, const double *v2, int N);
complex vdotproduct(const complex *v1, const complex *v2, int N); complex vdotproduct(const complex *v1, const complex *v2, int N);
void vmove(double *vdst, const double* vsrc, int N); void vmove(double *vdst, const double* vsrc, int N);
void vmove(complex *vdst, const complex* vsrc, int N); void vmove(complex *vdst, const complex* vsrc, int N);
void vmoveneg(double *vdst, const double *vsrc, int N); void vmoveneg(double *vdst, const double *vsrc, int N);
void vmoveneg(complex *vdst, const complex *vsrc, int N); void vmoveneg(complex *vdst, const complex *vsrc, int N);
skipping to change at line 297 skipping to change at line 331
{ {
if( Aligned ) if( Aligned )
ap::afree(m_Vec); ap::afree(m_Vec);
else else
delete[] m_Vec; delete[] m_Vec;
} }
m_iLow = iLow; m_iLow = iLow;
m_iHigh = iHigh; m_iHigh = iHigh;
m_iVecSize = iHigh-iLow+1; m_iVecSize = iHigh-iLow+1;
if( Aligned ) if( Aligned )
m_Vec = (T*)ap::amalloc(m_iVecSize*sizeof(T), 16); m_Vec = (T*)ap::amalloc((size_t)(m_iVecSize*sizeof(T)), 16);
else else
m_Vec = new T[m_iVecSize]; m_Vec = new T[(size_t)m_iVecSize];
}; };
void setlength(int iLen) void setlength(int iLen)
{ {
setbounds(0, iLen-1); setbounds(0, iLen-1);
} }
void setcontent( int iLow, int iHigh, const T *pContent ) void setcontent( int iLow, int iHigh, const T *pContent )
{ {
setbounds(iLow, iHigh); setbounds(iLow, iHigh);
skipping to change at line 472 skipping to change at line 506
int n2 = iHigh2-iLow2+1; int n2 = iHigh2-iLow2+1;
m_iVecSize = n1*n2; m_iVecSize = n1*n2;
if( Aligned ) if( Aligned )
{ {
//if( n2%2!=0 ) //if( n2%2!=0 )
while( (n2*sizeof(T))%16!=0 ) while( (n2*sizeof(T))%16!=0 )
{ {
n2++; n2++;
m_iVecSize += n1; m_iVecSize += n1;
} }
m_Vec = (T*)ap::amalloc(m_iVecSize*sizeof(T), 16); m_Vec = (T*)ap::amalloc((size_t)(m_iVecSize*sizeof(T)), 16);
} }
else else
m_Vec = new T[m_iVecSize]; m_Vec = new T[(size_t)m_iVecSize];
m_iLow1 = iLow1; m_iLow1 = iLow1;
m_iHigh1 = iHigh1; m_iHigh1 = iHigh1;
m_iLow2 = iLow2; m_iLow2 = iLow2;
m_iHigh2 = iHigh2; m_iHigh2 = iHigh2;
m_iConstOffset = -m_iLow2-m_iLow1*n2; m_iConstOffset = -m_iLow2-m_iLow1*n2;
m_iLinearMember = n2; m_iLinearMember = n2;
}; };
void setlength(int iLen1, int iLen2) void setlength(int iLen1, int iLen2)
{ {
skipping to change at line 539 skipping to change at line 573
return const_raw_vector<T>(&((*this)(iRowStart, iColumn)), iRow End-iRowStart+1, m_iLinearMember); return const_raw_vector<T>(&((*this)(iRowStart, iColumn)), iRow End-iRowStart+1, m_iLinearMember);
}; };
const_raw_vector<T> getrow(int iRow, int iColumnStart, int iColumnEnd) const const_raw_vector<T> getrow(int iRow, int iColumnStart, int iColumnEnd) const
{ {
if( (iColumnStart>iColumnEnd) || wrongRow(iRow) || wrongColumn(iCol umnStart) || wrongColumn(iColumnEnd)) if( (iColumnStart>iColumnEnd) || wrongRow(iRow) || wrongColumn(iCol umnStart) || wrongColumn(iColumnEnd))
return const_raw_vector<T>(0, 0, 1); return const_raw_vector<T>(0, 0, 1);
else else
return const_raw_vector<T>(&((*this)(iRow, iColumnStart)), iCol umnEnd-iColumnStart+1, 1); return const_raw_vector<T>(&((*this)(iRow, iColumnStart)), iCol umnEnd-iColumnStart+1, 1);
}; };
int getstride() const
{
return m_iLinearMember;
};
private: private:
bool wrongRow(int i) const { return i<m_iLow1 || i>m_iHigh1; }; bool wrongRow(int i) const { return i<m_iLow1 || i>m_iHigh1; };
bool wrongColumn(int j) const { return j<m_iLow2 || j>m_iHigh2; }; bool wrongColumn(int j) const { return j<m_iLow2 || j>m_iHigh2; };
T *m_Vec; T *m_Vec;
long m_iVecSize; long m_iVecSize;
long m_iLow1, m_iLow2, m_iHigh1, m_iHigh2; long m_iLow1, m_iLow2, m_iHigh1, m_iHigh2;
long m_iConstOffset, m_iLinearMember; long m_iConstOffset, m_iLinearMember;
}; };
 End of changes. 7 change blocks. 
6 lines changed or deleted 64 lines changed or added


 cdet.h   cdet.h 
skipping to change at line 27 skipping to change at line 27
>>> END OF LICENSE >>> >>> END OF LICENSE >>>
*************************************************************************/ *************************************************************************/
#ifndef _cdet_h #ifndef _cdet_h
#define _cdet_h #define _cdet_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "clu.h" #include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
/************************************************************************* /*************************************************************************
Determinant calculation of the matrix given by its LU decomposition. Determinant calculation of the matrix given by its LU decomposition.
Input parameters: Input parameters:
A - LU decomposition of the matrix (output of A - LU decomposition of the matrix (output of
RMatrixLU subroutine). RMatrixLU subroutine).
Pivots - table of permutations which were made during Pivots - table of permutations which were made during
the LU decomposition. the LU decomposition.
Output of RMatrixLU subroutine. Output of RMatrixLU subroutine.
skipping to change at line 63 skipping to change at line 69
A - matrix, array[0..N-1, 0..N-1] A - matrix, array[0..N-1, 0..N-1]
N - size of matrix A. N - size of matrix A.
Result: determinant of matrix A. Result: determinant of matrix A.
-- ALGLIB -- -- ALGLIB --
Copyright 2005 by Bochkanov Sergey Copyright 2005 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
ap::complex cmatrixdet(ap::complex_2d_array a, int n); ap::complex cmatrixdet(ap::complex_2d_array a, int n);
ap::complex complexdeterminantlu(const ap::complex_2d_array& a,
const ap::integer_1d_array& pivots,
int n);
ap::complex complexdeterminant(ap::complex_2d_array a, int n);
#endif #endif
 End of changes. 2 change blocks. 
7 lines changed or deleted 7 lines changed or added


 cinverse.h   cinverse.h 
skipping to change at line 33 skipping to change at line 33
>>> END OF LICENSE >>> >>> END OF LICENSE >>>
*************************************************************************/ *************************************************************************/
#ifndef _cinverse_h #ifndef _cinverse_h
#define _cinverse_h #define _cinverse_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "clu.h" #include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "ctrinverse.h" #include "ctrinverse.h"
/************************************************************************* /*************************************************************************
Inversion of a complex matrix given by its LU decomposition. Inversion of a complex matrix given by its LU decomposition.
Input parameters: Input parameters:
A - LU decomposition of the matrix (output of CMatrixLU subrout ine). A - LU decomposition of the matrix (output of CMatrixLU subrout ine).
Pivots - table of permutations which were made during the LU decompo sition Pivots - table of permutations which were made during the LU decompo sition
(the output of CMatrixLU subroutine). (the output of CMatrixLU subroutine).
N - size of matrix A. N - size of matrix A.
skipping to change at line 82 skipping to change at line 88
Result: Result:
True, if the matrix is not singular. True, if the matrix is not singular.
False, if the matrix is singular. False, if the matrix is singular.
-- ALGLIB -- -- ALGLIB --
Copyright 2005 by Bochkanov Sergey Copyright 2005 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
bool cmatrixinverse(ap::complex_2d_array& a, int n); bool cmatrixinverse(ap::complex_2d_array& a, int n);
bool complexinverselu(ap::complex_2d_array& a,
const ap::integer_1d_array& pivots,
int n);
bool complexinverse(ap::complex_2d_array& a, int n);
#endif #endif
 End of changes. 2 change blocks. 
7 lines changed or deleted 7 lines changed or added


 densesolver.h   densesolver.h 
skipping to change at line 35 skipping to change at line 35
#include "ialglib.h" #include "ialglib.h"
#include "reflections.h" #include "reflections.h"
#include "bidiagonal.h" #include "bidiagonal.h"
#include "qr.h" #include "qr.h"
#include "lq.h" #include "lq.h"
#include "blas.h" #include "blas.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "lu.h" #include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "trlinsolve.h" #include "trlinsolve.h"
#include "safesolve.h"
#include "rcond.h" #include "rcond.h"
#include "tsort.h" #include "tsort.h"
#include "xblas.h" #include "xblas.h"
struct densesolverreport struct densesolverreport
{ {
double r1; double r1;
double rinf; double rinf;
}; };
skipping to change at line 58 skipping to change at line 64
{ {
double r2; double r2;
ap::real_2d_array cx; ap::real_2d_array cx;
int n; int n;
int k; int k;
}; };
/************************************************************************* /*************************************************************************
Dense solver. Dense solver.
This subroutine solves a system A*X=B, where A is NxN non-denegerate This subroutine solves a system A*x=b, where A is NxN non-denegerate
real matrix, X and B are NxM real matrices. real matrix, x and b are vectors.
Additional features include: Algorithm features:
* automatic detection of degenerate cases * automatic detection of degenerate cases
* iterative improvement * condition number estimation
* iterative refinement
* O(N^3) complexity
INPUT PARAMETERS INPUT PARAMETERS
A - array[0..N-1,0..N-1], system matrix A - array[0..N-1,0..N-1], system matrix
N - size of A N - size of A
B - array[0..N-1,0..M-1], right part B - array[0..N-1], right part
M - size of right part
OUTPUT PARAMETERS OUTPUT PARAMETERS
Info - return code: Info - return code:
* -3 if A is singular, or VERY close to singular. * -3 A is singular, or VERY close to singular.
X is filled by zeros in such cases. X is filled by zeros in such cases.
* -1 if N<=0 or M<=0 was passed * -1 N<=0 was passed
* 1 if task is solved (matrix A may be near singular, * 1 task is solved (but matrix A may be ill-conditioned
,
check R1/RInf parameters for condition numbers). check R1/RInf parameters for condition numbers).
Rep - solver report, see below for more info Rep - solver report, see below for more info
X - array[0..N-1,0..M-1], it contains: X - array[0..N-1], it contains:
* solution of A*X=B if A is non-singular (well-conditioned * solution of A*x=b if A is non-singular (well-conditioned
or ill-conditioned, but not very close to singular) or ill-conditioned, but not very close to singular)
* zeros, if A is singular or VERY close to singular * zeros, if A is singular or VERY close to singular
(in this case Info=-3). (in this case Info=-3).
SOLVER REPORT SOLVER REPORT
Subroutine sets following fields of the Rep structure: Subroutine sets following fields of the Rep structure:
* R1 reciprocal of condition number: 1/cond(A), 1-norm. * R1 reciprocal of condition number: 1/cond(A), 1-norm.
* RInf reciprocal of condition number: 1/cond(A), inf-norm. * RInf reciprocal of condition number: 1/cond(A), inf-norm.
SEE ALSO: -- ALGLIB --
DenseSolverR() - solves A*x = b, where x and b are Nx1 matrices. Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void rmatrixsolve(const ap::real_2d_array& a,
int n,
const ap::real_1d_array& b,
int& info,
densesolverreport& rep,
ap::real_1d_array& x);
/*************************************************************************
Dense solver.
Similar to RMatrixSolve() but solves task with multiple right parts (where
b and x are NxM matrices).
Algorithm features:
* automatic detection of degenerate cases
* condition number estimation
* optional iterative refinement
* O(N^3+M*N^2) complexity
INPUT PARAMETERS
A - array[0..N-1,0..N-1], system matrix
N - size of A
B - array[0..N-1,0..M-1], right part
M - right part size
RFS - iterative refinement switch:
* True - refinement is used.
Less performance, more precision.
* False - refinement is not used.
More performance, less precision.
OUTPUT PARAMETERS
Info - same as in RMatrixSolve
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB -- -- ALGLIB --
Copyright 24.08.2009 by Bochkanov Sergey Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void rmatrixsolvem(const ap::real_2d_array& a, void rmatrixsolvem(const ap::real_2d_array& a,
int n, int n,
const ap::real_2d_array& b, const ap::real_2d_array& b,
int m, int m,
bool rfs,
int& info,
densesolverreport& rep,
ap::real_2d_array& x);
/*************************************************************************
Dense solver.
This subroutine solves a system A*X=B, where A is NxN non-denegerate
real matrix given by its LU decomposition, X and B are NxM real matrices.
Algorithm features:
* automatic detection of degenerate cases
* O(N^2) complexity
* condition number estimation
No iterative refinement is provided because exact form of original matrix
is not known to subroutine. Use RMatrixSolve or RMatrixMixedSolve if you
need iterative refinement.
INPUT PARAMETERS
LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
P - array[0..N-1], pivots array, RMatrixLU result
N - size of A
B - array[0..N-1], right part
OUTPUT PARAMETERS
Info - same as in RMatrixSolve
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void rmatrixlusolve(const ap::real_2d_array& lua,
const ap::integer_1d_array& p,
int n,
const ap::real_1d_array& b,
int& info,
densesolverreport& rep,
ap::real_1d_array& x);
/*************************************************************************
Dense solver.
Similar to RMatrixLUSolve() but solves task with multiple right parts
(where b and x are NxM matrices).
Algorithm features:
* automatic detection of degenerate cases
* O(M*N^2) complexity
* condition number estimation
No iterative refinement is provided because exact form of original matrix
is not known to subroutine. Use RMatrixSolve or RMatrixMixedSolve if you
need iterative refinement.
INPUT PARAMETERS
LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
P - array[0..N-1], pivots array, RMatrixLU result
N - size of A
B - array[0..N-1,0..M-1], right part
M - right part size
OUTPUT PARAMETERS
Info - same as in RMatrixSolve
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void rmatrixlusolvem(const ap::real_2d_array& lua,
const ap::integer_1d_array& p,
int n,
const ap::real_2d_array& b,
int m,
int& info, int& info,
densesolverreport& rep, densesolverreport& rep,
ap::real_2d_array& x); ap::real_2d_array& x);
/************************************************************************* /*************************************************************************
Dense solver. Dense solver.
This subroutine solves a system A*x=b, where BOTH ORIGINAL A AND ITS
LU DECOMPOSITION ARE KNOWN. You can use it if for some reasons you have
both A and its LU decomposition.
Algorithm features:
* automatic detection of degenerate cases
* condition number estimation
* iterative refinement
* O(N^2) complexity
INPUT PARAMETERS
A - array[0..N-1,0..N-1], system matrix
LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
P - array[0..N-1], pivots array, RMatrixLU result
N - size of A
B - array[0..N-1], right part
OUTPUT PARAMETERS
Info - same as in RMatrixSolveM
Rep - same as in RMatrixSolveM
X - same as in RMatrixSolveM
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void rmatrixmixedsolve(const ap::real_2d_array& a,
const ap::real_2d_array& lua,
const ap::integer_1d_array& p,
int n,
const ap::real_1d_array& b,
int& info,
densesolverreport& rep,
ap::real_1d_array& x);
/*************************************************************************
Dense solver.
Similar to RMatrixMixedSolve() but solves task with multiple right parts
(where b and x are NxM matrices).
Algorithm features:
* automatic detection of degenerate cases
* condition number estimation
* iterative refinement
* O(M*N^2) complexity
INPUT PARAMETERS
A - array[0..N-1,0..N-1], system matrix
LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
P - array[0..N-1], pivots array, RMatrixLU result
N - size of A
B - array[0..N-1,0..M-1], right part
M - right part size
OUTPUT PARAMETERS
Info - same as in RMatrixSolveM
Rep - same as in RMatrixSolveM
X - same as in RMatrixSolveM
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void rmatrixmixedsolvem(const ap::real_2d_array& a,
const ap::real_2d_array& lua,
const ap::integer_1d_array& p,
int n,
const ap::real_2d_array& b,
int m,
int& info,
densesolverreport& rep,
ap::real_2d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixSolveM(), but for complex matrices.
Algorithm features:
* automatic detection of degenerate cases
* condition number estimation
* iterative refinement
* O(N^3+M*N^2) complexity
INPUT PARAMETERS
A - array[0..N-1,0..N-1], system matrix
N - size of A
B - array[0..N-1,0..M-1], right part
M - right part size
RFS - iterative refinement switch:
* True - refinement is used.
Less performance, more precision.
* False - refinement is not used.
More performance, less precision.
OUTPUT PARAMETERS
Info - same as in RMatrixSolve
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void cmatrixsolvem(const ap::complex_2d_array& a,
int n,
const ap::complex_2d_array& b,
int m,
bool rfs,
int& info,
densesolverreport& rep,
ap::complex_2d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixSolve(), but for complex matrices.
Algorithm features:
* automatic detection of degenerate cases
* condition number estimation
* iterative refinement
* O(N^3) complexity
INPUT PARAMETERS
A - array[0..N-1,0..N-1], system matrix
N - size of A
B - array[0..N-1], right part
OUTPUT PARAMETERS
Info - same as in RMatrixSolve
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void cmatrixsolve(const ap::complex_2d_array& a,
int n,
const ap::complex_1d_array& b,
int& info,
densesolverreport& rep,
ap::complex_1d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixLUSolveM(), but for complex matrices.
Algorithm features:
* automatic detection of degenerate cases
* O(M*N^2) complexity
* condition number estimation
No iterative refinement is provided because exact form of original matrix
is not known to subroutine. Use CMatrixSolve or CMatrixMixedSolve if you
need iterative refinement.
INPUT PARAMETERS
LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
P - array[0..N-1], pivots array, RMatrixLU result
N - size of A
B - array[0..N-1,0..M-1], right part
M - right part size
OUTPUT PARAMETERS
Info - same as in RMatrixSolve
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void cmatrixlusolvem(const ap::complex_2d_array& lua,
const ap::integer_1d_array& p,
int n,
const ap::complex_2d_array& b,
int m,
int& info,
densesolverreport& rep,
ap::complex_2d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixLUSolve(), but for complex matrices.
Algorithm features:
* automatic detection of degenerate cases
* O(N^2) complexity
* condition number estimation
No iterative refinement is provided because exact form of original matrix
is not known to subroutine. Use CMatrixSolve or CMatrixMixedSolve if you
need iterative refinement.
INPUT PARAMETERS
LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
P - array[0..N-1], pivots array, CMatrixLU result
N - size of A
B - array[0..N-1], right part
OUTPUT PARAMETERS
Info - same as in RMatrixSolve
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void cmatrixlusolve(const ap::complex_2d_array& lua,
const ap::integer_1d_array& p,
int n,
const ap::complex_1d_array& b,
int& info,
densesolverreport& rep,
ap::complex_1d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixMixedSolveM(), but for complex matrices.
Algorithm features:
* automatic detection of degenerate cases
* condition number estimation
* iterative refinement
* O(M*N^2) complexity
INPUT PARAMETERS
A - array[0..N-1,0..N-1], system matrix
LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
P - array[0..N-1], pivots array, CMatrixLU result
N - size of A
B - array[0..N-1,0..M-1], right part
M - right part size
OUTPUT PARAMETERS
Info - same as in RMatrixSolveM
Rep - same as in RMatrixSolveM
X - same as in RMatrixSolveM
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void cmatrixmixedsolvem(const ap::complex_2d_array& a,
const ap::complex_2d_array& lua,
const ap::integer_1d_array& p,
int n,
const ap::complex_2d_array& b,
int m,
int& info,
densesolverreport& rep,
ap::complex_2d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixMixedSolve(), but for complex matrices.
Algorithm features:
* automatic detection of degenerate cases
* condition number estimation
* iterative refinement
* O(N^2) complexity
INPUT PARAMETERS
A - array[0..N-1,0..N-1], system matrix
LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
P - array[0..N-1], pivots array, CMatrixLU result
N - size of A
B - array[0..N-1], right part
OUTPUT PARAMETERS
Info - same as in RMatrixSolveM
Rep - same as in RMatrixSolveM
X - same as in RMatrixSolveM
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void cmatrixmixedsolve(const ap::complex_2d_array& a,
const ap::complex_2d_array& lua,
const ap::integer_1d_array& p,
int n,
const ap::complex_1d_array& b,
int& info,
densesolverreport& rep,
ap::complex_1d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixSolveM(), but for symmetric positive definite
matrices.
Algorithm features:
* automatic detection of degenerate cases
* condition number estimation
* O(N^3+M*N^2) complexity
* matrix is represented by its upper or lower triangle
No iterative refinement is provided because such partial representation of
matrix does not allow efficient calculation of extra-precise matrix-vector
products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
need iterative refinement.
INPUT PARAMETERS
A - array[0..N-1,0..N-1], system matrix
N - size of A
IsUpper - what half of A is provided
B - array[0..N-1,0..M-1], right part
M - right part size
OUTPUT PARAMETERS
Info - same as in RMatrixSolve.
Returns -3 for non-SPD matrices.
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void spdmatrixsolvem(const ap::real_2d_array& a,
int n,
bool isupper,
const ap::real_2d_array& b,
int m,
int& info,
densesolverreport& rep,
ap::real_2d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixSolve(), but for SPD matrices.
Algorithm features:
* automatic detection of degenerate cases
* condition number estimation
* O(N^3) complexity
* matrix is represented by its upper or lower triangle
No iterative refinement is provided because such partial representation of
matrix does not allow efficient calculation of extra-precise matrix-vector
products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
need iterative refinement.
INPUT PARAMETERS
A - array[0..N-1,0..N-1], system matrix
N - size of A
IsUpper - what half of A is provided
B - array[0..N-1], right part
OUTPUT PARAMETERS
Info - same as in RMatrixSolve
Returns -3 for non-SPD matrices.
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void spdmatrixsolve(const ap::real_2d_array& a,
int n,
bool isupper,
const ap::real_1d_array& b,
int& info,
densesolverreport& rep,
ap::real_1d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixLUSolveM(), but for SPD matrices represented
by their Cholesky decomposition.
Algorithm features:
* automatic detection of degenerate cases
* O(M*N^2) complexity
* condition number estimation
* matrix is represented by its upper or lower triangle
No iterative refinement is provided because such partial representation of
matrix does not allow efficient calculation of extra-precise matrix-vector
products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
need iterative refinement.
INPUT PARAMETERS
CHA - array[0..N-1,0..N-1], Cholesky decomposition,
SPDMatrixCholesky result
N - size of CHA
IsUpper - what half of CHA is provided
B - array[0..N-1,0..M-1], right part
M - right part size
OUTPUT PARAMETERS
Info - same as in RMatrixSolve
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void spdmatrixcholeskysolvem(const ap::real_2d_array& cha,
int n,
bool isupper,
const ap::real_2d_array& b,
int m,
int& info,
densesolverreport& rep,
ap::real_2d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixLUSolve(), but for SPD matrices represented
by their Cholesky decomposition.
Algorithm features:
* automatic detection of degenerate cases
* O(N^2) complexity
* condition number estimation
* matrix is represented by its upper or lower triangle
No iterative refinement is provided because such partial representation of
matrix does not allow efficient calculation of extra-precise matrix-vector
products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
need iterative refinement.
INPUT PARAMETERS
CHA - array[0..N-1,0..N-1], Cholesky decomposition,
SPDMatrixCholesky result
N - size of A
IsUpper - what half of CHA is provided
B - array[0..N-1], right part
OUTPUT PARAMETERS
Info - same as in RMatrixSolve
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void spdmatrixcholeskysolve(const ap::real_2d_array& cha,
int n,
bool isupper,
const ap::real_1d_array& b,
int& info,
densesolverreport& rep,
ap::real_1d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixSolveM(), but for Hermitian positive definite
matrices.
Algorithm features:
* automatic detection of degenerate cases
* condition number estimation
* O(N^3+M*N^2) complexity
* matrix is represented by its upper or lower triangle
No iterative refinement is provided because such partial representation of
matrix does not allow efficient calculation of extra-precise matrix-vector
products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
need iterative refinement.
INPUT PARAMETERS
A - array[0..N-1,0..N-1], system matrix
N - size of A
IsUpper - what half of A is provided
B - array[0..N-1,0..M-1], right part
M - right part size
OUTPUT PARAMETERS
Info - same as in RMatrixSolve.
Returns -3 for non-HPD matrices.
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void hpdmatrixsolvem(const ap::complex_2d_array& a,
int n,
bool isupper,
const ap::complex_2d_array& b,
int m,
int& info,
densesolverreport& rep,
ap::complex_2d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixSolve(), but for Hermitian positive definite
matrices.
Algorithm features:
* automatic detection of degenerate cases
* condition number estimation
* O(N^3) complexity
* matrix is represented by its upper or lower triangle
No iterative refinement is provided because such partial representation of
matrix does not allow efficient calculation of extra-precise matrix-vector
products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
need iterative refinement.
INPUT PARAMETERS
A - array[0..N-1,0..N-1], system matrix
N - size of A
IsUpper - what half of A is provided
B - array[0..N-1], right part
OUTPUT PARAMETERS
Info - same as in RMatrixSolve
Returns -3 for non-HPD matrices.
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void hpdmatrixsolve(const ap::complex_2d_array& a,
int n,
bool isupper,
const ap::complex_1d_array& b,
int& info,
densesolverreport& rep,
ap::complex_1d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixLUSolveM(), but for HPD matrices represented
by their Cholesky decomposition.
Algorithm features:
* automatic detection of degenerate cases
* O(M*N^2) complexity
* condition number estimation
* matrix is represented by its upper or lower triangle
No iterative refinement is provided because such partial representation of
matrix does not allow efficient calculation of extra-precise matrix-vector
products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
need iterative refinement.
INPUT PARAMETERS
CHA - array[0..N-1,0..N-1], Cholesky decomposition,
HPDMatrixCholesky result
N - size of CHA
IsUpper - what half of CHA is provided
B - array[0..N-1,0..M-1], right part
M - right part size
OUTPUT PARAMETERS
Info - same as in RMatrixSolve
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void hpdmatrixcholeskysolvem(const ap::complex_2d_array& cha,
int n,
bool isupper,
const ap::complex_2d_array& b,
int m,
int& info,
densesolverreport& rep,
ap::complex_2d_array& x);
/*************************************************************************
Dense solver. Same as RMatrixLUSolve(), but for HPD matrices represented
by their Cholesky decomposition.
Algorithm features:
* automatic detection of degenerate cases
* O(N^2) complexity
* condition number estimation
* matrix is represented by its upper or lower triangle
No iterative refinement is provided because such partial representation of
matrix does not allow efficient calculation of extra-precise matrix-vector
products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
need iterative refinement.
INPUT PARAMETERS
CHA - array[0..N-1,0..N-1], Cholesky decomposition,
SPDMatrixCholesky result
N - size of A
IsUpper - what half of CHA is provided
B - array[0..N-1], right part
OUTPUT PARAMETERS
Info - same as in RMatrixSolve
Rep - same as in RMatrixSolve
X - same as in RMatrixSolve
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void hpdmatrixcholeskysolve(const ap::complex_2d_array& cha,
int n,
bool isupper,
const ap::complex_1d_array& b,
int& info,
densesolverreport& rep,
ap::complex_1d_array& x);
/*************************************************************************
Dense solver.
This subroutine finds solution of the linear system A*X=B with non-square, This subroutine finds solution of the linear system A*X=B with non-square,
possibly degenerate A. System is solved in the least squares sense, and possibly degenerate A. System is solved in the least squares sense, and
general least squares solution X = X0 + CX*y which minimizes |A*X-B| is general least squares solution X = X0 + CX*y which minimizes |A*X-B| is
returned. If A is non-degenerate, solution in the usual sense is returned returned. If A is non-degenerate, solution in the usual sense is returned
Additional features include: Algorithm features:
* iterative improvement * automatic detection of degenerate cases
* iterative refinement
* O(N^3) complexity
INPUT PARAMETERS INPUT PARAMETERS
A - array[0..NRows-1,0..NCols-1], system matrix A - array[0..NRows-1,0..NCols-1], system matrix
NRows - vertical size of A NRows - vertical size of A
NCols - horizontal size of A NCols - horizontal size of A
B - array[0..NCols-1], right part B - array[0..NCols-1], right part
Threshold- a number in [0,1]. Singular values beyond Threshold are Threshold- a number in [0,1]. Singular values beyond Threshold are
considered zero. Set it to 0.0, if you don't understand considered zero. Set it to 0.0, if you don't understand
what it means, so the solver will choose good value on its what it means, so the solver will choose good value on its
own. own.
skipping to change at line 159 skipping to change at line 870
*************************************************************************/ *************************************************************************/
void rmatrixsolvels(const ap::real_2d_array& a, void rmatrixsolvels(const ap::real_2d_array& a,
int nrows, int nrows,
int ncols, int ncols,
const ap::real_1d_array& b, const ap::real_1d_array& b,
double threshold, double threshold,
int& info, int& info,
densesolverlsreport& rep, densesolverlsreport& rep,
ap::real_1d_array& x); ap::real_1d_array& x);
/*************************************************************************
Dense solver.
Similar to RMatrixSolveM() but solves task with one right part (where b/x
are vectors, not matrices).
See RMatrixSolveM() description for more information about subroutine
parameters.
-- ALGLIB --
Copyright 24.08.2009 by Bochkanov Sergey
*************************************************************************/
void rmatrixsolve(const ap::real_2d_array& a,
int n,
const ap::real_1d_array& b,
int& info,
densesolverreport& rep,
ap::real_1d_array& x);
#endif #endif
 End of changes. 15 change blocks. 
36 lines changed or deleted 729 lines changed or added


 det.h   det.h 
skipping to change at line 27 skipping to change at line 27
>>> END OF LICENSE >>> >>> END OF LICENSE >>>
*************************************************************************/ *************************************************************************/
#ifndef _det_h #ifndef _det_h
#define _det_h #define _det_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "lu.h" #include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
/************************************************************************* /*************************************************************************
Determinant calculation of the matrix given by its LU decomposition. Determinant calculation of the matrix given by its LU decomposition.
Input parameters: Input parameters:
A - LU decomposition of the matrix (output of A - LU decomposition of the matrix (output of
RMatrixLU subroutine). RMatrixLU subroutine).
Pivots - table of permutations which were made during Pivots - table of permutations which were made during
the LU decomposition. the LU decomposition.
Output of RMatrixLU subroutine. Output of RMatrixLU subroutine.
skipping to change at line 63 skipping to change at line 69
A - matrix, array[0..N-1, 0..N-1] A - matrix, array[0..N-1, 0..N-1]
N - size of matrix A. N - size of matrix A.
Result: determinant of matrix A. Result: determinant of matrix A.
-- ALGLIB -- -- ALGLIB --
Copyright 2005 by Bochkanov Sergey Copyright 2005 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
double rmatrixdet(ap::real_2d_array a, int n); double rmatrixdet(ap::real_2d_array a, int n);
double determinantlu(const ap::real_2d_array& a,
const ap::integer_1d_array& pivots,
int n);
double determinant(ap::real_2d_array a, int n);
#endif #endif
 End of changes. 2 change blocks. 
7 lines changed or deleted 7 lines changed or added


 hqrnd.h   hqrnd.h 
skipping to change at line 107 skipping to change at line 107
double hqrndnormal(hqrndstate& state); double hqrndnormal(hqrndstate& state);
/************************************************************************* /*************************************************************************
Random number generator: random X and Y such that X^2+Y^2=1 Random number generator: random X and Y such that X^2+Y^2=1
State structure must be initialized with HQRNDRandomize() or HQRNDSeed(). State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
-- ALGLIB -- -- ALGLIB --
Copyright 02.12.2009 by Bochkanov Sergey Copyright 02.12.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
double hqrndunit2(hqrndstate& state, double& x, double& y); void hqrndunit2(hqrndstate& state, double& x, double& y);
/************************************************************************* /*************************************************************************
Random number generator: normal numbers Random number generator: normal numbers
This function generates two independent random numbers from normal This function generates two independent random numbers from normal
distribution. Its performance is equal to that of HQRNDNormal() distribution. Its performance is equal to that of HQRNDNormal()
State structure must be initialized with HQRNDRandomize() or HQRNDSeed(). State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
-- ALGLIB -- -- ALGLIB --
 End of changes. 1 change blocks. 
1 lines changed or deleted 1 lines changed or added


 ialglib.h   ialglib.h 
/******************************************************************** /********************************************************************
Stub file for assembly optimized ALGLIB subroutines. optimized ALGLIB subroutines.
********************************************************************/ ********************************************************************/
#ifndef IALGLIB_H #ifndef IALGLIB_H
#define IALGLIB_H #define IALGLIB_H
#include "ap.h" #include "ap.h"
#define ALGLIB_INTERCEPTS_ABLAS
namespace ialglib
{
void mv_32(const double *a, const double *x, double *y, int stride, double
alpha, double beta);
void mv(int m, int n, const double *a, const double *x, double *y, int stri
de, double alpha, double beta);
void mv_generic(int m, int n, const double *a, const double *x, double *y,
int stride, double alpha, double beta);
void mv_complex(int m, int n, const double *a, const double *x, ap::complex
*cy, double *dy, int stride, ap::complex alpha, ap::complex beta);
void mv_complex_generic(int m, int n, const double *a, const double *x, ap:
:complex *cy, double *dy, int stride, ap::complex alpha, ap::complex beta);
void vzero(int n, double *p, int stride);
void vzero_complex(int n, ap::complex *p, int stride);
void vcopy(int n, const double *a, int stridea, double *b, int strideb);
void vcopy_complex(int n, const ap::complex *a, int stridea, double *b, int
strideb, char *conj);
void vcopy_complex(int n, const double *a, int stridea, double *b, int stri
deb, char *conj);
void mcopyblock(int m, int n, const double *a, int op, int stride, double *
b);
void mcopyunblock(int m, int n, const double *a, int op, double *b, int str
ide);
void mcopyblock_complex(int m, int n, const ap::complex *a, int op, int str
ide, double *b);
void mcopyunblock_complex(int m, int n, const double *a, int op, ap::comple
x* b, int stride);
bool _i_rmatrixgemmf(int m,
int n,
int k,
double alpha,
const ap::real_2d_array& a,
int ia,
int ja,
int optypea,
const ap::real_2d_array& b,
int ib,
int jb,
int optypeb,
double beta,
ap::real_2d_array& c,
int ic,
int jc);
bool _i_cmatrixgemmf(int m,
int n,
int k,
ap::complex alpha,
const ap::complex_2d_array& a,
int ia,
int ja,
int optypea,
const ap::complex_2d_array& b,
int ib,
int jb,
int optypeb,
ap::complex beta,
ap::complex_2d_array& c,
int ic,
int jc);
bool _i_cmatrixrighttrsmf(int m,
int n,
const ap::complex_2d_array& a,
int i1,
int j1,
bool isupper,
bool isunit,
int optype,
ap::complex_2d_array& x,
int i2,
int j2);
bool _i_rmatrixrighttrsmf(int m,
int n,
const ap::real_2d_array& a,
int i1,
int j1,
bool isupper,
bool isunit,
int optype,
ap::real_2d_array& x,
int i2,
int j2);
bool _i_cmatrixlefttrsmf(int m,
int n,
const ap::complex_2d_array& a,
int i1,
int j1,
bool isupper,
bool isunit,
int optype,
ap::complex_2d_array& x,
int i2,
int j2);
bool _i_rmatrixlefttrsmf(int m,
int n,
const ap::real_2d_array& a,
int i1,
int j1,
bool isupper,
bool isunit,
int optype,
ap::real_2d_array& x,
int i2,
int j2);
bool _i_cmatrixsyrkf(int n,
int k,
double alpha,
const ap::complex_2d_array& a,
int ia,
int ja,
int optypea,
double beta,
ap::complex_2d_array& c,
int ic,
int jc,
bool isupper);
bool _i_rmatrixsyrkf(int n,
int k,
double alpha,
const ap::real_2d_array& a,
int ia,
int ja,
int optypea,
double beta,
ap::real_2d_array& c,
int ic,
int jc,
bool isupper);
bool _i_cmatrixrank1f(int m,
int n,
ap::complex_2d_array& a,
int ia,
int ja,
ap::complex_1d_array& u,
int uoffs,
ap::complex_1d_array& v,
int voffs);
bool _i_rmatrixrank1f(int m,
int n,
ap::real_2d_array& a,
int ia,
int ja,
ap::real_1d_array& u,
int uoffs,
ap::real_1d_array& v,
int voffs);
}
#endif #endif
 End of changes. 3 change blocks. 
1 lines changed or deleted 151 lines changed or added


 inv.h   inv.h 
skipping to change at line 33 skipping to change at line 33
>>> END OF LICENSE >>> >>> END OF LICENSE >>>
*************************************************************************/ *************************************************************************/
#ifndef _inv_h #ifndef _inv_h
#define _inv_h #define _inv_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "lu.h" #include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "trinverse.h" #include "trinverse.h"
/************************************************************************* /*************************************************************************
Inversion of a matrix given by its LU decomposition. Inversion of a matrix given by its LU decomposition.
Input parameters: Input parameters:
A - LU decomposition of the matrix (output of RMatrixLU subrout ine). A - LU decomposition of the matrix (output of RMatrixLU subrout ine).
Pivots - table of permutations which were made during the LU decompo sition Pivots - table of permutations which were made during the LU decompo sition
(the output of RMatrixLU subroutine). (the output of RMatrixLU subroutine).
N - size of matrix A. N - size of matrix A.
skipping to change at line 82 skipping to change at line 88
Result: Result:
True, if the matrix is not singular. True, if the matrix is not singular.
False, if the matrix is singular. False, if the matrix is singular.
-- ALGLIB -- -- ALGLIB --
Copyright 2005 by Bochkanov Sergey Copyright 2005 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
bool rmatrixinverse(ap::real_2d_array& a, int n); bool rmatrixinverse(ap::real_2d_array& a, int n);
bool inverselu(ap::real_2d_array& a,
const ap::integer_1d_array& pivots,
int n);
bool inverse(ap::real_2d_array& a, int n);
#endif #endif
 End of changes. 2 change blocks. 
7 lines changed or deleted 7 lines changed or added


 lbfgs.h   lbfgs.h 
skipping to change at line 138 skipping to change at line 138
a many subsequent tasks with same N/M values. a many subsequent tasks with same N/M values.
First call MUST be without this flag bit set, First call MUST be without this flag bit set,
subsequent calls of MinLBFGS with same LBFGSSta te subsequent calls of MinLBFGS with same LBFGSSta te
structure can set Flags to 1. structure can set Flags to 1.
Output parameters: Output parameters:
State - structure used for reverse communication. State - structure used for reverse communication.
See also MinLBFGSIteration, MinLBFGSResults See also MinLBFGSIteration, MinLBFGSResults
NOTE:
Passing EpsG=0, EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to
automatic stopping criterion selection (small EpsX).
-- ALGLIB -- -- ALGLIB --
Copyright 14.11.2007 by Bochkanov Sergey Copyright 14.11.2007 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void minlbfgs(const int& n, void minlbfgs(int n,
const int& m, int m,
ap::real_1d_array& x, const ap::real_1d_array& x,
const double& epsg, double epsg,
const double& epsf, double epsf,
const double& epsx, double epsx,
const int& maxits, int maxits,
int flags, int flags,
lbfgsstate& state); lbfgsstate& state);
/************************************************************************* /*************************************************************************
One L-BFGS iteration One L-BFGS iteration
Called after initialization with MinLBFGS. Called after initialization with MinLBFGS.
See HTML documentation for examples. See HTML documentation for examples.
Input parameters: Input parameters:
 End of changes. 2 change blocks. 
7 lines changed or deleted 12 lines changed or added


 lda.h   lda.h 
skipping to change at line 34 skipping to change at line 34
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "blas.h" #include "blas.h"
#include "rotations.h" #include "rotations.h"
#include "tdevd.h" #include "tdevd.h"
#include "sblas.h" #include "sblas.h"
#include "reflections.h" #include "reflections.h"
#include "tridiagonal.h" #include "tridiagonal.h"
#include "sevd.h" #include "sevd.h"
#include "cholesky.h" #include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "spdinverse.h" #include "spdinverse.h"
/************************************************************************* /*************************************************************************
Multiclass Fisher LDA Multiclass Fisher LDA
Subroutine finds coefficients of linear combination which optimally separat es Subroutine finds coefficients of linear combination which optimally separat es
training set on classes. training set on classes.
INPUT PARAMETERS: INPUT PARAMETERS:
XY - training set, array[0..NPoints-1,0..NVars]. XY - training set, array[0..NPoints-1,0..NVars].
 End of changes. 1 change blocks. 
1 lines changed or deleted 6 lines changed or added


 logit.h   logit.h 
skipping to change at line 29 skipping to change at line 29
*************************************************************************/ *************************************************************************/
#ifndef _logit_h #ifndef _logit_h
#define _logit_h #define _logit_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "descriptivestatistics.h" #include "descriptivestatistics.h"
#include "mlpbase.h" #include "mlpbase.h"
#include "cholesky.h" #include "reflections.h"
#include "spdsolve.h" #include "bidiagonal.h"
#include "qr.h"
#include "lq.h"
#include "blas.h"
#include "rotations.h"
#include "bdsvd.h"
#include "svd.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "trlinsolve.h"
#include "safesolve.h"
#include "rcond.h"
#include "tsort.h" #include "tsort.h"
#include "xblas.h"
#include "densesolver.h"
#include "bdss.h" #include "bdss.h"
struct logitmodel struct logitmodel
{ {
ap::real_1d_array w; ap::real_1d_array w;
}; };
struct logitmcstate struct logitmcstate
{ {
bool brackt; bool brackt;
 End of changes. 2 change blocks. 
2 lines changed or deleted 19 lines changed or added


 lsfit.h   lsfit.h 
skipping to change at line 28 skipping to change at line 28
>>> END OF LICENSE >>> >>> END OF LICENSE >>>
*************************************************************************/ *************************************************************************/
#ifndef _lsfit_h #ifndef _lsfit_h
#define _lsfit_h #define _lsfit_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "blas.h" #include "blas.h"
#include "trinverse.h"
#include "cholesky.h"
#include "spdsolve.h"
#include "lbfgs.h"
#include "minlm.h"
#include "reflections.h" #include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "trinverse.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "bidiagonal.h" #include "bidiagonal.h"
#include "qr.h" #include "qr.h"
#include "lq.h" #include "lq.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "lu.h"
#include "trlinsolve.h" #include "trlinsolve.h"
#include "safesolve.h"
#include "rcond.h" #include "rcond.h"
#include "tsort.h"
#include "xblas.h"
#include "densesolver.h"
#include "lbfgs.h"
#include "minlm.h"
#include "spline3.h" #include "spline3.h"
#include "leastsquares.h" #include "leastsquares.h"
/************************************************************************* /*************************************************************************
Least squares fitting report: Least squares fitting report:
TaskRCond reciprocal of task's condition number TaskRCond reciprocal of task's condition number
RMSError RMS error RMSError RMS error
AvgError average error AvgError average error
AvgRelError average relative error (for non-zero Y[I]) AvgRelError average relative error (for non-zero Y[I])
MaxError maximum error MaxError maximum error
 End of changes. 5 change blocks. 
6 lines changed or deleted 13 lines changed or added


 mannwhitneyu.h   mannwhitneyu.h 
skipping to change at line 27 skipping to change at line 27
>>> END OF LICENSE >>> >>> END OF LICENSE >>>
*************************************************************************/ *************************************************************************/
#ifndef _mannwhitneyu_h #ifndef _mannwhitneyu_h
#define _mannwhitneyu_h #define _mannwhitneyu_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "normaldistr.h"
/************************************************************************* /*************************************************************************
Mann-Whitney U-test Mann-Whitney U-test
This test checks hypotheses about whether X and Y are samples of two This test checks hypotheses about whether X and Y are samples of two
continuous distributions of the same shape and same median or whether continuous distributions of the same shape and same median or whether
their medians are different. their medians are different.
The following tests are performed: The following tests are performed:
* two-tailed test (null hypothesis - the medians are equal) * two-tailed test (null hypothesis - the medians are equal)
* left-tailed test (null hypothesis - the median of the first sample * left-tailed test (null hypothesis - the median of the first sample
 End of changes. 1 change blocks. 
2 lines changed or deleted 0 lines changed or added


 minlm.h   minlm.h 
skipping to change at line 28 skipping to change at line 28
>>> END OF LICENSE >>> >>> END OF LICENSE >>>
*************************************************************************/ *************************************************************************/
#ifndef _minlm_h #ifndef _minlm_h
#define _minlm_h #define _minlm_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "blas.h" #include "blas.h"
#include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "trinverse.h" #include "trinverse.h"
#include "cholesky.h" #include "ablasf.h"
#include "spdsolve.h" #include "ablas.h"
#include "trfac.h"
#include "bidiagonal.h"
#include "qr.h"
#include "lq.h"
#include "rotations.h"
#include "bdsvd.h"
#include "svd.h"
#include "trlinsolve.h"
#include "safesolve.h"
#include "rcond.h"
#include "tsort.h"
#include "xblas.h"
#include "densesolver.h"
#include "lbfgs.h" #include "lbfgs.h"
struct lmstate struct lmstate
{ {
bool wrongparams; bool wrongparams;
int n; int n;
int m; int m;
double epsf; double epsf;
double epsx; double epsx;
int maxits; int maxits;
skipping to change at line 71 skipping to change at line 88
ap::real_2d_array model; ap::real_2d_array model;
ap::real_1d_array work; ap::real_1d_array work;
ap::rcommstate rstate; ap::rcommstate rstate;
int repiterationscount; int repiterationscount;
int repterminationtype; int repterminationtype;
int repnfunc; int repnfunc;
int repnjac; int repnjac;
int repngrad; int repngrad;
int repnhess; int repnhess;
int repncholesky; int repncholesky;
int solverinfo;
densesolverreport solverrep;
}; };
struct lmreport struct lmreport
{ {
int iterationscount; int iterationscount;
int terminationtype; int terminationtype;
int nfunc; int nfunc;
int njac; int njac;
int ngrad; int ngrad;
int nhess; int nhess;
skipping to change at line 280 skipping to change at line 299
State - algorithm state (used by MinLMIteration). State - algorithm state (used by MinLMIteration).
Output parameters: Output parameters:
X - array[0..N-1], solution X - array[0..N-1], solution
Rep - optimization report: Rep - optimization report:
* Rep.TerminationType completetion code: * Rep.TerminationType completetion code:
* -1 incorrect parameters were specified * -1 incorrect parameters were specified
* 1 relative function improvement is no more than * 1 relative function improvement is no more than
EpsF. EpsF.
* 2 relative step is no more than EpsX. * 2 relative step is no more than EpsX.
* 4 gradient norm is no more than EpsG
* 5 MaxIts steps was taken * 5 MaxIts steps was taken
* Rep.IterationsCount contains iterations count * Rep.IterationsCount contains iterations count
* Rep.NFunc - number of function calculations * Rep.NFunc - number of function calculations
* Rep.NJac - number of Jacobi matrix calculations * Rep.NJac - number of Jacobi matrix calculations
* Rep.NGrad - number of gradient calculations * Rep.NGrad - number of gradient calculations
* Rep.NHess - number of Hessian calculations * Rep.NHess - number of Hessian calculations
* Rep.NCholesky - number of Cholesky decomposition calculat ions * Rep.NCholesky - number of Cholesky decomposition calculat ions
-- ALGLIB -- -- ALGLIB --
Copyright 10.03.2009 by Bochkanov Sergey Copyright 10.03.2009 by Bochkanov Sergey
 End of changes. 4 change blocks. 
3 lines changed or deleted 21 lines changed or added


 mlpe.h   mlpe.h 
skipping to change at line 30 skipping to change at line 30
#ifndef _mlpe_h #ifndef _mlpe_h
#define _mlpe_h #define _mlpe_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "mlpbase.h" #include "mlpbase.h"
#include "trinverse.h" #include "trinverse.h"
#include "lbfgs.h" #include "lbfgs.h"
#include "cholesky.h" #include "reflections.h"
#include "spdsolve.h" #include "bidiagonal.h"
#include "mlptrain.h" #include "qr.h"
#include "lq.h"
#include "blas.h"
#include "rotations.h"
#include "bdsvd.h"
#include "svd.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "trlinsolve.h"
#include "safesolve.h"
#include "rcond.h"
#include "tsort.h" #include "tsort.h"
#include "xblas.h"
#include "densesolver.h"
#include "mlptrain.h"
#include "descriptivestatistics.h" #include "descriptivestatistics.h"
#include "bdss.h" #include "bdss.h"
/************************************************************************* /*************************************************************************
Neural networks ensemble Neural networks ensemble
*************************************************************************/ *************************************************************************/
struct mlpensemble struct mlpensemble
{ {
ap::integer_1d_array structinfo; ap::integer_1d_array structinfo;
int ensemblesize; int ensemblesize;
 End of changes. 2 change blocks. 
3 lines changed or deleted 20 lines changed or added


 mlptrain.h   mlptrain.h 
skipping to change at line 30 skipping to change at line 30
#ifndef _mlptrain_h #ifndef _mlptrain_h
#define _mlptrain_h #define _mlptrain_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "mlpbase.h" #include "mlpbase.h"
#include "trinverse.h" #include "trinverse.h"
#include "lbfgs.h" #include "lbfgs.h"
#include "cholesky.h" #include "reflections.h"
#include "spdsolve.h" #include "bidiagonal.h"
#include "qr.h"
#include "lq.h"
#include "blas.h"
#include "rotations.h"
#include "bdsvd.h"
#include "svd.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "trlinsolve.h"
#include "safesolve.h"
#include "rcond.h"
#include "tsort.h"
#include "xblas.h"
#include "densesolver.h"
/************************************************************************* /*************************************************************************
Training report: Training report:
* NGrad - number of gradient calculations * NGrad - number of gradient calculations
* NHess - number of Hessian calculations * NHess - number of Hessian calculations
* NCholesky - number of Cholesky decompositions * NCholesky - number of Cholesky decompositions
*************************************************************************/ *************************************************************************/
struct mlpreport struct mlpreport
{ {
int ngrad; int ngrad;
 End of changes. 1 change blocks. 
2 lines changed or deleted 20 lines changed or added


 polint.h   polint.h 
skipping to change at line 30 skipping to change at line 30
#ifndef _polint_h #ifndef _polint_h
#define _polint_h #define _polint_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "tsort.h" #include "tsort.h"
#include "ratinterpolation.h" #include "ratinterpolation.h"
#include "blas.h" #include "blas.h"
#include "trinverse.h"
#include "cholesky.h"
#include "spdsolve.h"
#include "lbfgs.h"
#include "minlm.h"
#include "reflections.h" #include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "trinverse.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "bidiagonal.h" #include "bidiagonal.h"
#include "qr.h" #include "qr.h"
#include "lq.h" #include "lq.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "lu.h"
#include "trlinsolve.h" #include "trlinsolve.h"
#include "safesolve.h"
#include "rcond.h" #include "rcond.h"
#include "xblas.h"
#include "densesolver.h"
#include "lbfgs.h"
#include "minlm.h"
#include "spline3.h" #include "spline3.h"
#include "leastsquares.h" #include "leastsquares.h"
#include "lsfit.h" #include "lsfit.h"
#include "ratint.h" #include "ratint.h"
#include "taskgen.h" #include "taskgen.h"
/************************************************************************* /*************************************************************************
Polynomial fitting report: Polynomial fitting report:
TaskRCond reciprocal of task's condition number TaskRCond reciprocal of task's condition number
RMSError RMS error RMSError RMS error
 End of changes. 5 change blocks. 
6 lines changed or deleted 12 lines changed or added


 ratint.h   ratint.h 
skipping to change at line 30 skipping to change at line 30
#ifndef _ratint_h #ifndef _ratint_h
#define _ratint_h #define _ratint_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "tsort.h" #include "tsort.h"
#include "ratinterpolation.h" #include "ratinterpolation.h"
#include "blas.h" #include "blas.h"
#include "trinverse.h"
#include "cholesky.h"
#include "spdsolve.h"
#include "lbfgs.h"
#include "minlm.h"
#include "reflections.h" #include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "trinverse.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "bidiagonal.h" #include "bidiagonal.h"
#include "qr.h" #include "qr.h"
#include "lq.h" #include "lq.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "lu.h"
#include "trlinsolve.h" #include "trlinsolve.h"
#include "safesolve.h"
#include "rcond.h" #include "rcond.h"
#include "xblas.h"
#include "densesolver.h"
#include "lbfgs.h"
#include "minlm.h"
#include "spline3.h" #include "spline3.h"
#include "leastsquares.h" #include "leastsquares.h"
#include "lsfit.h" #include "lsfit.h"
/************************************************************************* /*************************************************************************
Barycentric interpolant. Barycentric interpolant.
*************************************************************************/ *************************************************************************/
struct barycentricinterpolant struct barycentricinterpolant
{ {
int n; int n;
skipping to change at line 72 skipping to change at line 78
Barycentric fitting report: Barycentric fitting report:
TaskRCond reciprocal of task's condition number TaskRCond reciprocal of task's condition number
RMSError RMS error RMSError RMS error
AvgError average error AvgError average error
AvgRelError average relative error (for non-zero Y[I]) AvgRelError average relative error (for non-zero Y[I])
MaxError maximum error MaxError maximum error
*************************************************************************/ *************************************************************************/
struct barycentricfitreport struct barycentricfitreport
{ {
double taskrcond; double taskrcond;
double dbest; int dbest;
double rmserror; double rmserror;
double avgerror; double avgerror;
double avgrelerror; double avgrelerror;
double maxerror; double maxerror;
}; };
/************************************************************************* /*************************************************************************
Rational interpolation using barycentric formula Rational interpolation using barycentric formula
F(t) = SUM(i=0,n-1,w[i]*f[i]/(t-x[i])) / SUM(i=0,n-1,w[i]/(t-x[i])) F(t) = SUM(i=0,n-1,w[i]*f[i]/(t-x[i])) / SUM(i=0,n-1,w[i]/(t-x[i]))
skipping to change at line 310 skipping to change at line 316
const ap::real_1d_array& y, const ap::real_1d_array& y,
int n, int n,
int d, int d,
barycentricinterpolant& b); barycentricinterpolant& b);
/************************************************************************* /*************************************************************************
Weghted rational least squares fitting using Floater-Hormann rational Weghted rational least squares fitting using Floater-Hormann rational
functions with optimal D chosen from [0,9], with constraints and functions with optimal D chosen from [0,9], with constraints and
individual weights. individual weights.
Equidistant grid with M+1 node on [min(x),max(x)] is used to build basis Equidistant grid with M node on [min(x),max(x)] is used to build basis
functions. Different values of D are tried, optimal D (least WEIGHTED root functions. Different values of D are tried, optimal D (least WEIGHTED root
mean square error) is chosen. Task is linear, so linear least squares mean square error) is chosen. Task is linear, so linear least squares
solver is used. Complexity of this computational scheme is O(N*M^2) solver is used. Complexity of this computational scheme is O(N*M^2)
(mostly dominated by the least squares solver). (mostly dominated by the least squares solver).
SEE ALSO SEE ALSO
* BarycentricFitFloaterHormann(), "lightweight" fitting without invididual * BarycentricFitFloaterHormann(), "lightweight" fitting without invididual
weights and constraints. weights and constraints.
INPUT PARAMETERS: INPUT PARAMETERS:
 End of changes. 7 change blocks. 
8 lines changed or deleted 14 lines changed or added


 rcond.h   rcond.h 
skipping to change at line 33 skipping to change at line 33
>>> END OF LICENSE >>> >>> END OF LICENSE >>>
*************************************************************************/ *************************************************************************/
#ifndef _rcond_h #ifndef _rcond_h
#define _rcond_h #define _rcond_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "lu.h" #include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "trlinsolve.h" #include "trlinsolve.h"
#include "safesolve.h"
/************************************************************************* /*************************************************************************
Estimate of a matrix condition number (1-norm) Estimate of a matrix condition number (1-norm)
The algorithm calculates a lower bound of the condition number. In this cas e, The algorithm calculates a lower bound of the condition number. In this cas e,
the algorithm does not return a lower bound of the condition number, but an the algorithm does not return a lower bound of the condition number, but an
inverse number (to avoid an overflow in case of a singular matrix). inverse number (to avoid an overflow in case of a singular matrix).
Input parameters: Input parameters:
A - matrix. Array whose indexes range within [0..N-1, 0..N-1]. A - matrix. Array whose indexes range within [0..N-1, 0..N-1].
N - size of matrix A. N - size of matrix A.
Result: 1/LowerBound(cond(A)) Result: 1/LowerBound(cond(A))
NOTE:
if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
returned in such cases.
*************************************************************************/ *************************************************************************/
double rmatrixrcond1(const ap::real_2d_array& a, int n); double rmatrixrcond1(ap::real_2d_array a, int n);
/************************************************************************* /*************************************************************************
Estimate of the condition number of a matrix given by its LU decomposition (1-norm) Estimate of a matrix condition number (infinity-norm).
The algorithm calculates a lower bound of the condition number. In this cas e, The algorithm calculates a lower bound of the condition number. In this cas e,
the algorithm does not return a lower bound of the condition number, but an the algorithm does not return a lower bound of the condition number, but an
inverse number (to avoid an overflow in case of a singular matrix). inverse number (to avoid an overflow in case of a singular matrix).
Input parameters: Input parameters:
LUDcmp - LU decomposition of a matrix in compact form. Output of A - matrix. Array whose indexes range within [0..N-1, 0..N-1].
the RMatrixLU subroutine. N - size of matrix A.
N - size of matrix A.
Result: 1/LowerBound(cond(A))
NOTE:
if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
returned in such cases.
*************************************************************************/
double rmatrixrcondinf(ap::real_2d_array a, int n);
/*************************************************************************
Condition number estimate of a symmetric positive definite matrix.
The algorithm calculates a lower bound of the condition number. In this cas
e,
the algorithm does not return a lower bound of the condition number, but an
inverse number (to avoid an overflow in case of a singular matrix).
It should be noted that 1-norm and inf-norm of condition numbers of symmetr
ic
matrices are equal, so the algorithm doesn't take into account the
differences between these types of norms.
Input parameters:
A - symmetric positive definite matrix which is given by its
upper or lower triangle depending on the value of
IsUpper. Array with elements [0..N-1, 0..N-1].
N - size of matrix A.
IsUpper - storage format.
Result:
1/LowerBound(cond(A)), if matrix A is positive definite,
-1, if matrix A is not positive definite, and its condition number
could not be found by this algorithm.
NOTE:
if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
returned in such cases.
*************************************************************************/
double spdmatrixrcond(ap::real_2d_array a, int n, bool isupper);
/*************************************************************************
Condition number estimate of a Hermitian positive definite matrix.
The algorithm calculates a lower bound of the condition number. In this cas
e,
the algorithm does not return a lower bound of the condition number, but an
inverse number (to avoid an overflow in case of a singular matrix).
It should be noted that 1-norm and inf-norm of condition numbers of symmetr
ic
matrices are equal, so the algorithm doesn't take into account the
differences between these types of norms.
Input parameters:
A - Hermitian positive definite matrix which is given by its
upper or lower triangle depending on the value of
IsUpper. Array with elements [0..N-1, 0..N-1].
N - size of matrix A.
IsUpper - storage format.
Result:
1/LowerBound(cond(A)), if matrix A is positive definite,
-1, if matrix A is not positive definite, and its condition number
could not be found by this algorithm.
NOTE:
if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
returned in such cases.
*************************************************************************/
double hpdmatrixrcond(ap::complex_2d_array a, int n, bool isupper);
/*************************************************************************
Estimate of a matrix condition number (1-norm)
The algorithm calculates a lower bound of the condition number. In this cas
e,
the algorithm does not return a lower bound of the condition number, but an
inverse number (to avoid an overflow in case of a singular matrix).
Input parameters:
A - matrix. Array whose indexes range within [0..N-1, 0..N-1].
N - size of matrix A.
Result: 1/LowerBound(cond(A)) Result: 1/LowerBound(cond(A))
NOTE:
if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
returned in such cases.
*************************************************************************/ *************************************************************************/
double rmatrixlurcond1(const ap::real_2d_array& ludcmp, int n); double cmatrixrcond1(ap::complex_2d_array a, int n);
/************************************************************************* /*************************************************************************
Estimate of a matrix condition number (infinity-norm). Estimate of a matrix condition number (infinity-norm).
The algorithm calculates a lower bound of the condition number. In this cas e, The algorithm calculates a lower bound of the condition number. In this cas e,
the algorithm does not return a lower bound of the condition number, but an the algorithm does not return a lower bound of the condition number, but an
inverse number (to avoid an overflow in case of a singular matrix). inverse number (to avoid an overflow in case of a singular matrix).
Input parameters: Input parameters:
A - matrix. Array whose indexes range within [0..N-1, 0..N-1]. A - matrix. Array whose indexes range within [0..N-1, 0..N-1].
N - size of matrix A. N - size of matrix A.
Result: 1/LowerBound(cond(A)) Result: 1/LowerBound(cond(A))
NOTE:
if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
returned in such cases.
*************************************************************************/ *************************************************************************/
double rmatrixrcondinf(const ap::real_2d_array& a, int n); double cmatrixrcondinf(ap::complex_2d_array a, int n);
/*************************************************************************
Estimate of the condition number of a matrix given by its LU decomposition
(1-norm)
The algorithm calculates a lower bound of the condition number. In this cas
e,
the algorithm does not return a lower bound of the condition number, but an
inverse number (to avoid an overflow in case of a singular matrix).
Input parameters:
LUA - LU decomposition of a matrix in compact form. Output of
the RMatrixLU subroutine.
N - size of matrix A.
Result: 1/LowerBound(cond(A))
NOTE:
if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
returned in such cases.
*************************************************************************/
double rmatrixlurcond1(const ap::real_2d_array& lua, int n);
/************************************************************************* /*************************************************************************
Estimate of the condition number of a matrix given by its LU decomposition Estimate of the condition number of a matrix given by its LU decomposition
(infinity norm). (infinity norm).
The algorithm calculates a lower bound of the condition number. In this cas e, The algorithm calculates a lower bound of the condition number. In this cas e,
the algorithm does not return a lower bound of the condition number, but an the algorithm does not return a lower bound of the condition number, but an
inverse number (to avoid an overflow in case of a singular matrix). inverse number (to avoid an overflow in case of a singular matrix).
Input parameters: Input parameters:
LUDcmp - LU decomposition of a matrix in compact form. Output of LUA - LU decomposition of a matrix in compact form. Output of
the RMatrixLU subroutine. the RMatrixLU subroutine.
N - size of matrix A. N - size of matrix A.
Result: 1/LowerBound(cond(A)) Result: 1/LowerBound(cond(A))
NOTE:
if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
returned in such cases.
*************************************************************************/
double rmatrixlurcondinf(const ap::real_2d_array& lua, int n);
/*************************************************************************
Condition number estimate of a symmetric positive definite matrix given by
Cholesky decomposition.
The algorithm calculates a lower bound of the condition number. In this
case, the algorithm does not return a lower bound of the condition number,
but an inverse number (to avoid an overflow in case of a singular matrix).
It should be noted that 1-norm and inf-norm condition numbers of symmetric
matrices are equal, so the algorithm doesn't take into account the
differences between these types of norms.
Input parameters:
CD - Cholesky decomposition of matrix A,
output of SMatrixCholesky subroutine.
N - size of matrix A.
Result: 1/LowerBound(cond(A))
NOTE:
if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
returned in such cases.
*************************************************************************/ *************************************************************************/
double rmatrixlurcondinf(const ap::real_2d_array& ludcmp, int n); double spdmatrixcholeskyrcond(const ap::real_2d_array& a, int n, bool isupp
er);
/*************************************************************************
Condition number estimate of a Hermitian positive definite matrix given by
Cholesky decomposition.
The algorithm calculates a lower bound of the condition number. In this
case, the algorithm does not return a lower bound of the condition number,
but an inverse number (to avoid an overflow in case of a singular matrix).
It should be noted that 1-norm and inf-norm condition numbers of symmetric
matrices are equal, so the algorithm doesn't take into account the
differences between these types of norms.
Input parameters:
CD - Cholesky decomposition of matrix A,
output of SMatrixCholesky subroutine.
N - size of matrix A.
double rcond1(ap::real_2d_array a, int n); Result: 1/LowerBound(cond(A))
double rcond1lu(const ap::real_2d_array& ludcmp, int n); NOTE:
if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
returned in such cases.
*************************************************************************/
double hpdmatrixcholeskyrcond(const ap::complex_2d_array& a,
int n,
bool isupper);
double rcondinf(ap::real_2d_array a, int n); /*************************************************************************
Estimate of the condition number of a matrix given by its LU decomposition
(1-norm)
double rcondinflu(const ap::real_2d_array& ludcmp, int n); The algorithm calculates a lower bound of the condition number. In this cas
e,
the algorithm does not return a lower bound of the condition number, but an
inverse number (to avoid an overflow in case of a singular matrix).
Input parameters:
LUA - LU decomposition of a matrix in compact form. Output of
the CMatrixLU subroutine.
N - size of matrix A.
Result: 1/LowerBound(cond(A))
NOTE:
if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
returned in such cases.
*************************************************************************/
double cmatrixlurcond1(const ap::complex_2d_array& lua, int n);
/*************************************************************************
Estimate of the condition number of a matrix given by its LU decomposition
(infinity norm).
The algorithm calculates a lower bound of the condition number. In this cas
e,
the algorithm does not return a lower bound of the condition number, but an
inverse number (to avoid an overflow in case of a singular matrix).
Input parameters:
LUA - LU decomposition of a matrix in compact form. Output of
the CMatrixLU subroutine.
N - size of matrix A.
Result: 1/LowerBound(cond(A))
NOTE:
if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
returned in such cases.
*************************************************************************/
double cmatrixlurcondinf(const ap::complex_2d_array& lua, int n);
/*************************************************************************
Threshold for rcond: matrices with condition number beyond this threshold
are considered singular.
Threshold must be far enough from underflow, at least Sqr(Threshold) must
be greater than underflow.
*************************************************************************/
double rcondthreshold();
#endif #endif
 End of changes. 17 change blocks. 
14 lines changed or deleted 238 lines changed or added


 spddet.h   spddet.h 
skipping to change at line 27 skipping to change at line 27
>>> END OF LICENSE >>> >>> END OF LICENSE >>>
*************************************************************************/ *************************************************************************/
#ifndef _spddet_h #ifndef _spddet_h
#define _spddet_h #define _spddet_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "cholesky.h" #include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
/************************************************************************* /*************************************************************************
Determinant calculation of the matrix given by the Cholesky decomposition. Determinant calculation of the matrix given by the Cholesky decomposition.
Input parameters: Input parameters:
A - Cholesky decomposition, A - Cholesky decomposition,
output of SMatrixCholesky subroutine. output of SMatrixCholesky subroutine.
N - size of matrix A. N - size of matrix A.
As the determinant is equal to the product of squares of diagonal elements, As the determinant is equal to the product of squares of diagonal elements,
skipping to change at line 69 skipping to change at line 75
Result: Result:
determinant of matrix A. determinant of matrix A.
If matrix A is not positive definite, then subroutine returns -1. If matrix A is not positive definite, then subroutine returns -1.
-- ALGLIB -- -- ALGLIB --
Copyright 2005-2008 by Bochkanov Sergey Copyright 2005-2008 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
double spdmatrixdet(ap::real_2d_array a, int n, bool isupper); double spdmatrixdet(ap::real_2d_array a, int n, bool isupper);
double determinantcholesky(const ap::real_2d_array& a, int n);
double determinantspd(ap::real_2d_array a, int n, bool isupper);
#endif #endif
 End of changes. 2 change blocks. 
5 lines changed or deleted 7 lines changed or added


 spdgevd.h   spdgevd.h 
skipping to change at line 27 skipping to change at line 27
>>> END OF LICENSE >>> >>> END OF LICENSE >>>
*************************************************************************/ *************************************************************************/
#ifndef _spdgevd_h #ifndef _spdgevd_h
#define _spdgevd_h #define _spdgevd_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "cholesky.h" #include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "sblas.h" #include "sblas.h"
#include "blas.h" #include "blas.h"
#include "trinverse.h" #include "trinverse.h"
#include "rotations.h" #include "rotations.h"
#include "tdevd.h" #include "tdevd.h"
#include "reflections.h"
#include "tridiagonal.h" #include "tridiagonal.h"
#include "sevd.h" #include "sevd.h"
/************************************************************************* /*************************************************************************
Algorithm for solving the following generalized symmetric positive-definite Algorithm for solving the following generalized symmetric positive-definite
eigenproblem: eigenproblem:
A*x = lambda*B*x (1) or A*x = lambda*B*x (1) or
A*B*x = lambda*x (2) or A*B*x = lambda*x (2) or
B*A*x = lambda*x (3). B*A*x = lambda*x (3).
where A is a symmetric matrix, B - symmetric positive-definite matrix. where A is a symmetric matrix, B - symmetric positive-definite matrix.
skipping to change at line 158 skipping to change at line 163
*************************************************************************/ *************************************************************************/
bool smatrixgevdreduce(ap::real_2d_array& a, bool smatrixgevdreduce(ap::real_2d_array& a,
int n, int n,
bool isuppera, bool isuppera,
const ap::real_2d_array& b, const ap::real_2d_array& b,
bool isupperb, bool isupperb,
int problemtype, int problemtype,
ap::real_2d_array& r, ap::real_2d_array& r,
bool& isupperr); bool& isupperr);
bool generalizedsymmetricdefiniteevd(ap::real_2d_array a,
int n,
bool isuppera,
const ap::real_2d_array& b,
bool isupperb,
int zneeded,
int problemtype,
ap::real_1d_array& d,
ap::real_2d_array& z);
bool generalizedsymmetricdefiniteevdreduce(ap::real_2d_array& a,
int n,
bool isuppera,
const ap::real_2d_array& b,
bool isupperb,
int problemtype,
ap::real_2d_array& r,
bool& isupperr);
#endif #endif
 End of changes. 3 change blocks. 
21 lines changed or deleted 7 lines changed or added


 spdinverse.h   spdinverse.h 
skipping to change at line 33 skipping to change at line 33
>>> END OF LICENSE >>> >>> END OF LICENSE >>>
*************************************************************************/ *************************************************************************/
#ifndef _spdinverse_h #ifndef _spdinverse_h
#define _spdinverse_h #define _spdinverse_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "cholesky.h" #include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
/************************************************************************* /*************************************************************************
Inversion of a symmetric positive definite matrix which is given Inversion of a symmetric positive definite matrix which is given
by Cholesky decomposition. by Cholesky decomposition.
Input parameters: Input parameters:
A - Cholesky decomposition of the matrix to be inverted: A - Cholesky decomposition of the matrix to be inverted:
A=U A=U
Output of CholeskyDecomposition subroutine. Output of CholeskyDecomposition subroutine.
Array with elements [0..N-1, 0..N-1]. Array with elements [0..N-1, 0..N-1].
skipping to change at line 93 skipping to change at line 99
is used, and the elements below the main diagonal are not is used, and the elements below the main diagonal are not
used nor changed. The same applies if IsUpper = False. used nor changed. The same applies if IsUpper = False.
Result: Result:
True, if the matrix is positive definite. True, if the matrix is positive definite.
False, if the matrix is not positive definite (and it could not be False, if the matrix is not positive definite (and it could not be
inverted by this algorithm). inverted by this algorithm).
*************************************************************************/ *************************************************************************/
bool spdmatrixinverse(ap::real_2d_array& a, int n, bool isupper); bool spdmatrixinverse(ap::real_2d_array& a, int n, bool isupper);
bool inversecholesky(ap::real_2d_array& a, int n, bool isupper);
bool inversesymmetricpositivedefinite(ap::real_2d_array& a,
int n,
bool isupper);
#endif #endif
 End of changes. 2 change blocks. 
7 lines changed or deleted 7 lines changed or added


 spline1d.h   spline1d.h 
skipping to change at line 29 skipping to change at line 29
*************************************************************************/ *************************************************************************/
#ifndef _spline1d_h #ifndef _spline1d_h
#define _spline1d_h #define _spline1d_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "spline3.h" #include "spline3.h"
#include "blas.h" #include "blas.h"
#include "trinverse.h"
#include "cholesky.h"
#include "spdsolve.h"
#include "lbfgs.h"
#include "minlm.h"
#include "reflections.h" #include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "trinverse.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "bidiagonal.h" #include "bidiagonal.h"
#include "qr.h" #include "qr.h"
#include "lq.h" #include "lq.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "lu.h"
#include "trlinsolve.h" #include "trlinsolve.h"
#include "safesolve.h"
#include "rcond.h" #include "rcond.h"
#include "tsort.h"
#include "xblas.h"
#include "densesolver.h"
#include "lbfgs.h"
#include "minlm.h"
#include "leastsquares.h" #include "leastsquares.h"
#include "lsfit.h" #include "lsfit.h"
/************************************************************************* /*************************************************************************
1-dimensional spline inteprolant 1-dimensional spline inteprolant
*************************************************************************/ *************************************************************************/
struct spline1dinterpolant struct spline1dinterpolant
{ {
int n; int n;
int k; int k;
 End of changes. 5 change blocks. 
6 lines changed or deleted 13 lines changed or added


 spline2d.h   spline2d.h 
skipping to change at line 29 skipping to change at line 29
*************************************************************************/ *************************************************************************/
#ifndef _spline2d_h #ifndef _spline2d_h
#define _spline2d_h #define _spline2d_h
#include "ap.h" #include "ap.h"
#include "ialglib.h" #include "ialglib.h"
#include "spline3.h" #include "spline3.h"
#include "blas.h" #include "blas.h"
#include "trinverse.h"
#include "cholesky.h"
#include "spdsolve.h"
#include "lbfgs.h"
#include "minlm.h"
#include "reflections.h" #include "reflections.h"
#include "creflections.h"
#include "hqrnd.h"
#include "matgen.h"
#include "trinverse.h"
#include "ablasf.h"
#include "ablas.h"
#include "trfac.h"
#include "bidiagonal.h" #include "bidiagonal.h"
#include "qr.h" #include "qr.h"
#include "lq.h" #include "lq.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "lu.h"
#include "trlinsolve.h" #include "trlinsolve.h"
#include "safesolve.h"
#include "rcond.h" #include "rcond.h"
#include "tsort.h"
#include "xblas.h"
#include "densesolver.h"
#include "lbfgs.h"
#include "minlm.h"
#include "leastsquares.h" #include "leastsquares.h"
#include "lsfit.h" #include "lsfit.h"
#include "spline1d.h" #include "spline1d.h"
/************************************************************************* /*************************************************************************
2-dimensional spline inteprolant 2-dimensional spline inteprolant
*************************************************************************/ *************************************************************************/
struct spline2dinterpolant struct spline2dinterpolant
{ {
int k; int k;
 End of changes. 5 change blocks. 
6 lines changed or deleted 13 lines changed or added


 xblas.h   xblas.h 
skipping to change at line 26 skipping to change at line 26
INPUT PARAMETERS INPUT PARAMETERS
A - array[0..N-1], vector 1 A - array[0..N-1], vector 1
B - array[0..N-1], vector 2 B - array[0..N-1], vector 2
N - vectors length, N<2^29. N - vectors length, N<2^29.
Temp - array[0..N-1], pre-allocated temporary storage Temp - array[0..N-1], pre-allocated temporary storage
OUTPUT PARAMETERS OUTPUT PARAMETERS
R - (A,B) R - (A,B)
RErr - estimate of error. This estimate accounts for both errors RErr - estimate of error. This estimate accounts for both errors
during calculation of (A,B) and errors introduced by during calculation of (A,B) and errors introduced by
rounding of A/B to fit in double (about 1 ulp). rounding of A and B to fit in double (about 1 ulp).
-- ALGLIB -- -- ALGLIB --
Copyright 24.08.2009 by Bochkanov Sergey Copyright 24.08.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void xdot(const ap::real_1d_array& a, void xdot(const ap::real_1d_array& a,
const ap::real_1d_array& b, const ap::real_1d_array& b,
int n, int n,
ap::real_1d_array& temp, ap::real_1d_array& temp,
double& r, double& r,
double& rerr); double& rerr);
/*************************************************************************
More precise complex dot-product. Absolute error of subroutine result is
about 1 ulp of max(MX,V), where:
MX = max( |a[i]*b[i]| )
V = |(a,b)|
INPUT PARAMETERS
A - array[0..N-1], vector 1
B - array[0..N-1], vector 2
N - vectors length, N<2^29.
Temp - array[0..2*N-1], pre-allocated temporary storage
OUTPUT PARAMETERS
R - (A,B)
RErr - estimate of error. This estimate accounts for both errors
during calculation of (A,B) and errors introduced by
rounding of A and B to fit in double (about 1 ulp).
-- ALGLIB --
Copyright 27.01.2010 by Bochkanov Sergey
*************************************************************************/
void xcdot(const ap::complex_1d_array& a,
const ap::complex_1d_array& b,
int n,
ap::real_1d_array& temp,
ap::complex& r,
double& rerr);
#endif #endif
 End of changes. 2 change blocks. 
1 lines changed or deleted 29 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/