dng/AlgLib_ver3.19
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
AlgLib_ver3.19/solvers.mqh
super.admin 9e263d779c convert
2025-05-30 14:39:48 +02:00

7154 lines
295 KiB
MQL5
Raw Permalink Blame History

//+------------------------------------------------------------------+
//| solvers.mqh |
//| Copyright 2003-2022 Sergey Bochkanov (ALGLIB project) |
//| Copyright 2012-2023, MetaQuotes Ltd. |
//| https://www.mql5.com |
//+------------------------------------------------------------------+
//| Implementation of ALGLIB library in MetaQuotes Language 5 |
//| |
//| The features of the library include: |
//| - Linear algebra (direct algorithms, EVD, SVD) |
//| - Solving systems of linear and non-linear equations |
//| - Interpolation |
//| - Optimization |
//| - FFT (Fast Fourier Transform) |
//| - Numerical integration |
//| - Linear and nonlinear least-squares fitting |
//| - Ordinary differential equations |
//| - Computation of special functions |
//| - Descriptive statistics and hypothesis testing |
//| - Data analysis - classification, regression |
//| - Implementing linear algebra algorithms, interpolation, etc. |
//| in high-precision arithmetic (using MPFR) |
//| |
//| This file is free software; you can redistribute it and/or |
//| modify it under the terms of the GNU General Public License as |
//| published by the Free Software Foundation (www.fsf.org); either |
//| version 2 of the License, or (at your option) any later version. |
//| |
//| This program is distributed in the hope that it will be useful, |
//| but WITHOUT ANY WARRANTY; without even the implied warranty of |
//| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
//| GNU General Public License for more details. |
//+------------------------------------------------------------------+
//#include "matrix.mqh"
#include "ap.mqh"
#include "alglibinternal.mqh"
#include "linalg.mqh"
//+------------------------------------------------------------------+
//| Auxiliary class for CDenseSolver |
//+------------------------------------------------------------------+
class CDenseSolverReport
{
public:
double m_r1;
double m_rinf;
CDenseSolverReport(void) { ZeroMemory(this); }
~CDenseSolverReport(void) {}
void Copy(CDenseSolverReport &obj);
};
//+------------------------------------------------------------------+
//| Copy |
//+------------------------------------------------------------------+
void CDenseSolverReport::Copy(CDenseSolverReport &obj)
{
//--- copy variables
m_r1=obj.m_r1;
m_rinf=obj.m_rinf;
}
//+------------------------------------------------------------------+
//| This class is a shell for class CDenseSolverReport |
//+------------------------------------------------------------------+
class CDenseSolverReportShell
{
private:
CDenseSolverReport m_innerobj;
public:
//--- constructors, destructor
CDenseSolverReportShell(void) {}
CDenseSolverReportShell(CDenseSolverReport &obj) { m_innerobj.Copy(obj); }
~CDenseSolverReportShell(void) {}
//--- methods
double GetR1(void);
void SetR1(const double d);
double GetRInf(void);
void SetRInf(const double d);
CDenseSolverReport *GetInnerObj(void);
};
//+------------------------------------------------------------------+
//| Returns the value of the variable r1 |
//+------------------------------------------------------------------+
double CDenseSolverReportShell::GetR1(void)
{
return(m_innerobj.m_r1);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable r1 |
//+------------------------------------------------------------------+
void CDenseSolverReportShell::SetR1(const double d)
{
m_innerobj.m_r1=d;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable rinf |
//+------------------------------------------------------------------+
double CDenseSolverReportShell::GetRInf(void)
{
return(m_innerobj.m_rinf);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable rinf |
//+------------------------------------------------------------------+
void CDenseSolverReportShell::SetRInf(const double d)
{
m_innerobj.m_rinf=d;
}
//+------------------------------------------------------------------+
//| Return object of class |
//+------------------------------------------------------------------+
CDenseSolverReport *CDenseSolverReportShell::GetInnerObj(void)
{
return(GetPointer(m_innerobj));
}
//+------------------------------------------------------------------+
//| Auxiliary class for CDenseSolver |
//+------------------------------------------------------------------+
class CDenseSolverLSReport
{
public:
double m_r2;
CMatrixDouble m_cx;
int m_n;
int m_k;
//--- constructor, destructor
CDenseSolverLSReport(void) { m_r2=0; m_n=0; m_k=0; }
~CDenseSolverLSReport(void) {}
//--- copy
void Copy(CDenseSolverLSReport &obj);
};
//+------------------------------------------------------------------+
//| Copy |
//+------------------------------------------------------------------+
void CDenseSolverLSReport::Copy(CDenseSolverLSReport &obj)
{
//--- copy variables
m_r2=obj.m_r2;
m_n=obj.m_n;
m_k=obj.m_k;
//--- copy matrix
m_cx=obj.m_cx;
}
//+------------------------------------------------------------------+
//| This class is a shell for class CDenseSolverLSReport |
//+------------------------------------------------------------------+
class CDenseSolverLSReportShell
{
private:
CDenseSolverLSReport m_innerobj;
public:
//--- constructors, destructor
CDenseSolverLSReportShell(void) {}
CDenseSolverLSReportShell(CDenseSolverLSReport &obj) { m_innerobj.Copy(obj); }
~CDenseSolverLSReportShell(void) {}
//--- methods
double GetR2(void);
void SetR2(const double d);
int GetN(void);
void SetN(const int i);
int GetK(void);
void SetK(const int i);
CDenseSolverLSReport *GetInnerObj(void);
};
//+------------------------------------------------------------------+
//| Returns the value of the variable r2 |
//+------------------------------------------------------------------+
double CDenseSolverLSReportShell::GetR2(void)
{
//--- return result
return(m_innerobj.m_r2);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable r2 |
//+------------------------------------------------------------------+
void CDenseSolverLSReportShell::SetR2(const double d)
{
m_innerobj.m_r2=d;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable n |
//+------------------------------------------------------------------+
int CDenseSolverLSReportShell::GetN(void)
{
return(m_innerobj.m_n);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable n |
//+------------------------------------------------------------------+
void CDenseSolverLSReportShell::SetN(const int i)
{
m_innerobj.m_n=i;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable k |
//+------------------------------------------------------------------+
int CDenseSolverLSReportShell::GetK(void)
{
return(m_innerobj.m_k);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable k |
//+------------------------------------------------------------------+
void CDenseSolverLSReportShell::SetK(const int i)
{
m_innerobj.m_k=i;
}
//+------------------------------------------------------------------+
//| Return object of class |
//+------------------------------------------------------------------+
CDenseSolverLSReport *CDenseSolverLSReportShell::GetInnerObj(void)
{
return(GetPointer(m_innerobj));
}
//+------------------------------------------------------------------+
//| Dense solver |
//+------------------------------------------------------------------+
class CDenseSolver
{
private:
static void RMatrixLUSolveInternal(CMatrixDouble &lua,int &p[],const double scalea,const int n,CMatrixDouble &a,const bool havea,CMatrixDouble &b,const int m,int &info,CDenseSolverReport &rep,CMatrixDouble &x);
static void SPDMatrixCholeskySolveInternal(CMatrixDouble &cha,const double sqrtscalea,const int n,const bool IsUpper,CMatrixDouble &a,const bool havea,CMatrixDouble &b,const int m,int &info,CDenseSolverReport &rep,CMatrixDouble &x);
static void CMatrixLUSolveInternal(CMatrixComplex &lua,int &p[],const double scalea,const int n,CMatrixComplex &a,const bool havea,CMatrixComplex &b,const int m,int &info,CDenseSolverReport &rep,CMatrixComplex &x);
static void HPDMatrixCholeskySolveInternal(CMatrixComplex &cha,const double sqrtscalea,const int n,const bool IsUpper,CMatrixComplex &a,const bool havea,CMatrixComplex &b,const int m,int &info,CDenseSolverReport &rep,CMatrixComplex &x);
static int CDenseSolverRFSMax(const int n,const double r1,const double rinf);
static int CDenseSolverRFSMaxV2(const int n,const double r2);
static void RBasicLUSolve(CMatrixDouble &lua,int &p[],const double scalea,const int n,double &xb[],double &tmp[]);
static void SPDBasicCholeskySolve(CMatrixDouble &cha,const double sqrtscalea,const int n,const bool IsUpper,double &xb[],double &tmp[]);
static void CBasicLUSolve(CMatrixComplex &lua,int &p[],const double scalea,const int n,complex &xb[],complex &tmp[]);
static void HPDBasicCholeskySolve(CMatrixComplex &cha,const double sqrtscalea,const int n,const bool IsUpper,complex &xb[],complex &tmp[]);
public:
static void RMatrixSolve(CMatrixDouble &a,const int n,double &b[],int &info,CDenseSolverReport &rep,double &x[]);
static void RMatrixSolve(CMatrixDouble &a,const int n,CRowDouble &b,int &info,CDenseSolverReport &rep,CRowDouble &x);
static void RMatrixSolveM(CMatrixDouble &a,const int n,CMatrixDouble &b,const int m,const bool rfs,int &info,CDenseSolverReport &rep,CMatrixDouble &x);
static void RMatrixLUSolve(CMatrixDouble &lua,int &p[],const int n,double &b[],int &info,CDenseSolverReport &rep,double &x[]);
static void RMatrixLUSolveM(CMatrixDouble &lua,int &p[],const int n,CMatrixDouble &b,const int m,int &info,CDenseSolverReport &rep,CMatrixDouble &x);
static void RMatrixMixedSolve(CMatrixDouble &a,CMatrixDouble &lua,int &p[],const int n,double &b[],int &info,CDenseSolverReport &rep,double &x[]);
static void RMatrixMixedSolveM(CMatrixDouble &a,CMatrixDouble &lua,int &p[],const int n,CMatrixDouble &b,const int m,int &info,CDenseSolverReport &rep,CMatrixDouble &x);
static void CMatrixSolveM(CMatrixComplex &a,const int n,CMatrixComplex &b,const int m,const bool rfs,int &info,CDenseSolverReport &rep,CMatrixComplex &x);
static void CMatrixSolve(CMatrixComplex &a,const int n,complex &b[],int &info,CDenseSolverReport &rep,complex &x[]);
static void CMatrixLUSolveM(CMatrixComplex &lua,int &p[],const int n,CMatrixComplex &b,const int m,int &info,CDenseSolverReport &rep,CMatrixComplex &x);
static void CMatrixLUSolve(CMatrixComplex &lua,int &p[],const int n,complex &b[],int &info,CDenseSolverReport &rep,complex &x[]);
static void CMatrixMixedSolveM(CMatrixComplex &a,CMatrixComplex &lua,int &p[],const int n,CMatrixComplex &b,const int m,int &info,CDenseSolverReport &rep,CMatrixComplex &x);
static void CMatrixMixedSolve(CMatrixComplex &a,CMatrixComplex &lua,int &p[],const int n,complex &b[],int &info,CDenseSolverReport &rep,complex &x[]);
static void SPDMatrixSolveM(CMatrixDouble &a,const int n,const bool IsUpper,CMatrixDouble &b,const int m,int &info,CDenseSolverReport &rep,CMatrixDouble &x);
static void SPDMatrixSolve(CMatrixDouble &a,const int n,const bool IsUpper,double &b[],int &info,CDenseSolverReport &rep,double &x[]);
static void SPDMatrixCholeskySolveM(CMatrixDouble &cha,const int n,const bool IsUpper,CMatrixDouble &b,const int m,int &info,CDenseSolverReport &rep,CMatrixDouble &x);
static void SPDMatrixCholeskySolve(CMatrixDouble &cha,const int n,const bool IsUpper,double &b[],int &info,CDenseSolverReport &rep,double &x[]);
static void SPDMatrixCholeskySolve(CMatrixDouble &cha,const int n,const bool IsUpper,CRowDouble &b,int &info,CDenseSolverReport &rep,CRowDouble &x);
static void HPDMatrixSolveM(CMatrixComplex &a,const int n,const bool IsUpper,CMatrixComplex &b,const int m,int &info,CDenseSolverReport &rep,CMatrixComplex &x);
static void HPDMatrixSolve(CMatrixComplex &a,const int n,const bool IsUpper,complex &b[],int &info,CDenseSolverReport &rep,complex &x[]);
static void HPDMatrixCholeskySolveM(CMatrixComplex &cha,const int n,const bool IsUpper,CMatrixComplex &b,const int m,int &info,CDenseSolverReport &rep,CMatrixComplex &x);
static void HPDMatrixCholeskySolve(CMatrixComplex &cha,const int n,const bool IsUpper,complex &b[],int &info,CDenseSolverReport &rep,complex &x[]);
static void RMatrixSolveLS(CMatrixDouble &a,const int nrows,const int ncols,double &b[],double threshold,int &info,CDenseSolverLSReport &rep,double &x[]);
};
//+------------------------------------------------------------------+
//| Dense solver. |
//| This subroutine solves a system A*x=b, where A is NxN |
//| non-denegerate real matrix, x and b are vectors. |
//| 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 - return code: |
//| * -3 A is singular, or VERY close to singular.|
//| X is filled by zeros in such cases. |
//| * -1 N<=0 was passed |
//| * 1 task is solved (but matrix A may be |
//| ill-conditioned, check R1/RInf parameters|
//| for condition numbers). |
//| Rep - solver report, see below for more info |
//| X - array[0..N-1], it contains: |
//| * solution of A*x=b if A is non-singular |
//| (well-conditioned or ill-conditioned, but not |
//| very close to singular) |
//| * zeros, if A is singular or VERY close to |
//| singular (in this case Info=-3). |
//| SOLVER REPORT |
//| Subroutine sets following fields of the Rep structure: |
//| * R1 reciprocal of condition number: 1/cond(A), 1-norm. |
//| * RInf reciprocal of condition number: 1/cond(A), inf-norm. |
//+------------------------------------------------------------------+
void CDenseSolver::RMatrixSolve(CMatrixDouble &a,const int n,double &b[],
int &info,CDenseSolverReport &rep,
double &x[])
{
//--- create matrix
CMatrixDouble bm;
CMatrixDouble xm;
//--- initialization
info=0;
//--- check
if(n<=0)
{
info=-1;
return;
}
//--- allocation
bm.Resize(n,1);
//--- filling
for(int i_=0; i_<n; i_++)
bm.Set(i_,0,b[i_]);
//--- function call
RMatrixSolveM(a,n,bm,1,true,info,rep,xm);
//--- allocation
ArrayResize(x,n);
//--- copy
for(int i_=0; i_<n; i_++)
x[i_]=xm[i_][0];
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
void CDenseSolver::RMatrixSolve(CMatrixDouble &a,const int n,CRowDouble &b,
int &info,CDenseSolverReport &rep,
CRowDouble &x)
{
//--- create matrix
CMatrixDouble bm;
CMatrixDouble xm;
//--- initialization
info=0;
//--- check
if(n<=0)
{
info=-1;
return;
}
//--- allocation
bm.Resize(n,1);
//--- filling
bm.Col(0,b);
//--- function call
RMatrixSolveM(a,n,bm,1,true,info,rep,xm);
//--- copy
x=xm.Col(0)+0;
}
//+------------------------------------------------------------------+
//| 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 |
//+------------------------------------------------------------------+
void CDenseSolver::RMatrixSolveM(CMatrixDouble &a,const int n,
CMatrixDouble &b,const int m,
const bool rfs,int &info,
CDenseSolverReport &rep,
CMatrixDouble &x)
{
double scalea=0;
//--- create matrix
CMatrixDouble da;
CMatrixDouble emptya;
//--- create array
int p[];
//--- initialization
info=0;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
//--- allocation
da.Resize(n,n);
//--- 1. scale matrix,max(|A[i][j]|)
//--- 2. factorize scaled matrix
//--- 3. solve
scalea=0;
for(int i=0; i<n; i++)
{
for(int j=0; j<n; j++)
scalea=MathMax(scalea,MathAbs(a[i][j]));
}
//--- check
if(scalea==0.0)
scalea=1;
//--- change values
scalea=1/scalea;
for(int i=0; i<n; i++)
{
for(int i_=0; i_<n; i_++)
da.Set(i,i_,a[i][i_]);
}
//--- function call
CTrFac::RMatrixLU(da,n,n,p);
//--- check
if(rfs)
RMatrixLUSolveInternal(da,p,scalea,n,a,true,b,m,info,rep,x);
else
RMatrixLUSolveInternal(da,p,scalea,n,emptya,false,b,m,info,rep,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 |
//+------------------------------------------------------------------+
void CDenseSolver::RMatrixLUSolve(CMatrixDouble &lua,int &p[],
const int n,double &b[],
int &info,CDenseSolverReport &rep,
double &x[])
{
//--- create matrix
CMatrixDouble bm;
CMatrixDouble xm;
//--- initialization
info=0;
//--- check
if(n<=0)
{
info=-1;
return;
}
//--- allocation
bm.Resize(n,1);
//--- filling
for(int i_=0; i_<n; i_++)
bm.Set(i_,0,b[i_]);
//--- function call
RMatrixLUSolveM(lua,p,n,bm,1,info,rep,xm);
//--- allocation
ArrayResize(x,n);
//--- copy
for(int i_=0; i_<n; i_++)
x[i_]=xm[i_][0];
}
//+------------------------------------------------------------------+
//| 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 |
//+------------------------------------------------------------------+
void CDenseSolver::RMatrixLUSolveM(CMatrixDouble &lua,int &p[],
const int n,CMatrixDouble &b,
const int m,int &info,
CDenseSolverReport &rep,
CMatrixDouble &x)
{
double scalea=0;
//--- create matrix
CMatrixDouble emptya;
//--- initialization
info=0;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
//--- 1. scale matrix,max(|U[i][j]|)
//--- we assume that LU is in its normal form,i.e. |L[i][j]|<=1
//--- 2. solve
scalea=0;
for(int i=0; i<n; i++)
{
for(int j=i; j<n; j++)
scalea=MathMax(scalea,MathAbs(lua[i][j]));
}
//--- check
if(scalea==0.0)
scalea=1;
//--- change values
scalea=1/scalea;
//--- function call
RMatrixLUSolveInternal(lua,p,scalea,n,emptya,false,b,m,info,rep,x);
}
//+------------------------------------------------------------------+
//| 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 |
//+------------------------------------------------------------------+
void CDenseSolver::RMatrixMixedSolve(CMatrixDouble &a,CMatrixDouble &lua,
int &p[],const int n,double &b[],
int &info,CDenseSolverReport &rep,
double &x[])
{
//--- create matrix
CMatrixDouble bm;
CMatrixDouble xm;
//--- initialization
info=0;
//--- check
if(n<=0)
{
info=-1;
return;
}
//--- allocation
bm.Resize(n,1);
//--- filling
for(int i_=0; i_<n; i_++)
bm.Set(i_,0,b[i_]);
//--- function call
RMatrixMixedSolveM(a,lua,p,n,bm,1,info,rep,xm);
//--- allocation
ArrayResize(x,n);
//--- copy
for(int i_=0; i_<n; i_++)
x[i_]=xm[i_][0];
}
//+------------------------------------------------------------------+
//| 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 |
//+------------------------------------------------------------------+
void CDenseSolver::RMatrixMixedSolveM(CMatrixDouble &a,CMatrixDouble &lua,
int &p[],const int n,CMatrixDouble &b,
const int m,int &info,
CDenseSolverReport &rep,
CMatrixDouble &x)
{
double scalea=0;
//--- initialization
info=0;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
//--- 1. scale matrix,max(|A[i][j]|)
//--- 2. factorize scaled matrix
//--- 3. solve
scalea=0;
for(int i=0; i<n; i++)
{
for(int j=0; j<n; j++)
scalea=MathMax(scalea,MathAbs(a[i][j]));
}
//--- check
if(scalea==0.0)
scalea=1;
//--- change values
scalea=1/scalea;
//--- function call
RMatrixLUSolveInternal(lua,p,scalea,n,a,true,b,m,info,rep,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 |
//+------------------------------------------------------------------+
void CDenseSolver::CMatrixSolveM(CMatrixComplex &a,const int n,
CMatrixComplex &b,const int m,
const bool rfs,int &info,
CDenseSolverReport &rep,
CMatrixComplex &x)
{
double scalea=0;
//--- create array
int p[];
//--- create matrix
CMatrixComplex da;
CMatrixComplex emptya;
//--- initialization
info=0;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
//--- allocation
da.Resize(n,n);
//--- 1. scale matrix,max(|A[i][j]|)
//--- 2. factorize scaled matrix
//--- 3. solve
scalea=0;
for(int i=0; i<n; i++)
{
for(int j=0; j<n; j++)
scalea=MathMax(scalea,CMath::AbsComplex(a[i][j]));
}
//--- check
if(scalea==0.0)
scalea=1;
//--- change values
scalea=1/scalea;
for(int i=0; i<n; i++)
{
for(int i_=0; i_<n; i_++)
da.Set(i,i_,a[i][i_]);
}
//--- function call
CTrFac::CMatrixLU(da,n,n,p);
//--- check
if(rfs)
CMatrixLUSolveInternal(da,p,scalea,n,a,true,b,m,info,rep,x);
else
CMatrixLUSolveInternal(da,p,scalea,n,emptya,false,b,m,info,rep,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 |
//+------------------------------------------------------------------+
void CDenseSolver::CMatrixSolve(CMatrixComplex &a,const int n,
complex &b[],int &info,
CDenseSolverReport &rep,complex &x[])
{
//--- create matrix
CMatrixComplex bm;
CMatrixComplex xm;
//--- initialization
info=0;
//--- check
if(n<=0)
{
info=-1;
return;
}
//--- allocation
bm.Resize(n,1);
//--- filling
for(int i_=0; i_<n; i_++)
bm.Set(i_,0,b[i_]);
//--- function call
CMatrixSolveM(a,n,bm,1,true,info,rep,xm);
//--- allocation
ArrayResize(x,n);
//--- copy
for(int i_=0; i_<n; i_++)
x[i_]=xm[i_][0];
}
//+------------------------------------------------------------------+
//| 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 |
//+------------------------------------------------------------------+
void CDenseSolver::CMatrixLUSolveM(CMatrixComplex &lua,int &p[],
const int n,CMatrixComplex &b,
const int m,int &info,
CDenseSolverReport &rep,
CMatrixComplex &x)
{
double scalea=0;
//--- create matrix
CMatrixComplex emptya;
//--- initialization
info=0;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
//--- 1. scale matrix,max(|U[i][j]|)
//--- we assume that LU is in its normal form,i.e. |L[i][j]|<=1
//--- 2. solve
scalea=0;
for(int i=0; i<n; i++)
{
for(int j=i; j<n; j++)
scalea=MathMax(scalea,CMath::AbsComplex(lua[i][j]));
}
//--- check
if(scalea==0.0)
scalea=1;
//--- change values
scalea=1/scalea;
//--- function call
CMatrixLUSolveInternal(lua,p,scalea,n,emptya,false,b,m,info,rep,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 |
//+------------------------------------------------------------------+
void CDenseSolver::CMatrixLUSolve(CMatrixComplex &lua,int &p[],
const int n,complex &b[],int &info,
CDenseSolverReport &rep,complex &x[])
{
//--- create matrix
CMatrixComplex bm;
CMatrixComplex xm;
//--- initialization
info=0;
//--- check
if(n<=0)
{
info=-1;
return;
}
//--- allocation
bm.Resize(n,1);
//--- filling
for(int i_=0; i_<n; i_++)
bm.Set(i_,0,b[i_]);
//--- function call
CMatrixLUSolveM(lua,p,n,bm,1,info,rep,xm);
//--- allocation
ArrayResize(x,n);
//--- copy
for(int i_=0; i_<n; i_++)
x[i_]=xm[i_][0];
}
//+------------------------------------------------------------------+
//| 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 |
//+------------------------------------------------------------------+
void CDenseSolver::CMatrixMixedSolveM(CMatrixComplex &a,CMatrixComplex &lua,
int &p[],const int n,CMatrixComplex &b,
const int m,int &info,CDenseSolverReport &rep,
CMatrixComplex &x)
{
double scalea=0;
//--- initialization
info=0;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
//--- 1. scale matrix,max(|A[i][j]|)
//--- 2. factorize scaled matrix
//--- 3. solve
scalea=0;
for(int i=0; i<n; i++)
{
for(int j=0; j<n; j++)
scalea=MathMax(scalea,CMath::AbsComplex(a[i][j]));
}
//--- check
if(scalea==0.0)
scalea=1;
//--- change values
scalea=1/scalea;
//--- function call
CMatrixLUSolveInternal(lua,p,scalea,n,a,true,b,m,info,rep,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 |
//+------------------------------------------------------------------+
void CDenseSolver::CMatrixMixedSolve(CMatrixComplex &a,CMatrixComplex &lua,
int &p[],const int n,complex &b[],
int &info,CDenseSolverReport &rep,
complex &x[])
{
//--- create matrix
CMatrixComplex bm;
CMatrixComplex xm;
//--- initialization
info=0;
//--- check
if(n<=0)
{
info=-1;
return;
}
//--- allocation
bm.Resize(n,1);
//--- filling
for(int i_=0; i_<n; i_++)
bm.Set(i_,0,b[i_]);
//--- function call
CMatrixMixedSolveM(a,lua,p,n,bm,1,info,rep,xm);
//--- allocation
ArrayResize(x,n);
//--- copy
for(int i_=0; i_<n; i_++)
x[i_]=xm[i_][0];
}
//+------------------------------------------------------------------+
//| 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 |
//+------------------------------------------------------------------+
void CDenseSolver::SPDMatrixSolveM(CMatrixDouble &a,const int n,
const bool IsUpper,CMatrixDouble &b,
const int m,int &info,
CDenseSolverReport &rep,
CMatrixDouble &x)
{
//--- create variables
double sqrtscalea=0;
int i=0;
int j=0;
int j1=0;
int j2=0;
int i_=0;
//--- create matrix
CMatrixDouble da;
//--- initialization
info=0;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
//--- allocation
da.Resize(n,n);
//--- 1. scale matrix,max(|A[i][j]|)
//--- 2. factorize scaled matrix
//--- 3. solve
sqrtscalea=0;
for(i=0; i<n; i++)
{
//--- check
if(IsUpper)
{
j1=i;
j2=n-1;
}
else
{
j1=0;
j2=i;
}
//--- calculation
for(j=j1; j<=j2; j++)
sqrtscalea=MathMax(sqrtscalea,MathAbs(a[i][j]));
}
//--- check
if(sqrtscalea==0.0)
sqrtscalea=1;
//--- change values
sqrtscalea=1/sqrtscalea;
sqrtscalea=MathSqrt(sqrtscalea);
for(i=0; i<n; i++)
{
//--- check
if(IsUpper)
{
j1=i;
j2=n-1;
}
else
{
j1=0;
j2=i;
}
//--- calculation
for(i_=j1; i_<=j2; i_++)
da.Set(i,i_,a[i][i_]);
}
//--- check
if(!CTrFac::SPDMatrixCholesky(da,n,IsUpper))
{
//--- allocation
x.Resize(n,m);
for(i=0; i<n; i++)
{
for(j=0; j<=m-1; j++)
x.Set(i,j,0);
}
//--- change values
rep.m_r1=0;
rep.m_rinf=0;
info=-3;
//--- exit the function
return;
}
info=1;
//--- function call
SPDMatrixCholeskySolveInternal(da,sqrtscalea,n,IsUpper,a,true,b,m,info,rep,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 |
//+------------------------------------------------------------------+
void CDenseSolver::SPDMatrixSolve(CMatrixDouble &a,const int n,
const bool IsUpper,double &b[],
int &info,CDenseSolverReport &rep,
double &x[])
{
//--- create matrix
CMatrixDouble bm;
CMatrixDouble xm;
//--- initialization
info=0;
//--- check
if(n<=0)
{
info=-1;
return;
}
//--- allocation
bm.Resize(n,1);
//--- filling
for(int i_=0; i_<n; i_++)
bm.Set(i_,0,b[i_]);
//--- function call
SPDMatrixSolveM(a,n,IsUpper,bm,1,info,rep,xm);
//--- allocation
ArrayResize(x,n);
//--- copy
for(int i_=0; i_<n; i_++)
x[i_]=xm[i_][0];
}
//+------------------------------------------------------------------+
//| 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 |
//+------------------------------------------------------------------+
void CDenseSolver::SPDMatrixCholeskySolveM(CMatrixDouble &cha,const int n,
const bool IsUpper,CMatrixDouble &b,
const int m,int &info,
CDenseSolverReport &rep,
CMatrixDouble &x)
{
//--- create variables
double sqrtscalea=0;
int i=0;
int j=0;
int j1=0;
int j2=0;
//--- create matrix
CMatrixDouble emptya;
//--- initialization
info=0;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
//--- 1. scale matrix,max(|U[i][j]|)
//--- 2. factorize scaled matrix
//--- 3. solve
sqrtscalea=0;
for(i=0; i<n; i++)
{
//--- check
if(IsUpper)
{
j1=i;
j2=n-1;
}
else
{
j1=0;
j2=i;
}
//--- calculation
for(j=j1; j<=j2; j++)
sqrtscalea=MathMax(sqrtscalea,MathAbs(cha[i][j]));
}
//--- check
if(sqrtscalea==0.0)
sqrtscalea=1;
//--- change values
sqrtscalea=1/sqrtscalea;
//--- function call
SPDMatrixCholeskySolveInternal(cha,sqrtscalea,n,IsUpper,emptya,false,b,m,info,rep,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 |
//+------------------------------------------------------------------+
void CDenseSolver::SPDMatrixCholeskySolve(CMatrixDouble &cha,const int n,
const bool IsUpper,double &b[],
int &info,CDenseSolverReport &rep,
double &x[])
{
CRowDouble X;
CRowDouble B=b;
SPDMatrixCholeskySolve(cha,n,IsUpper,B,info,rep,X);
X.ToArray(x);
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
void CDenseSolver::SPDMatrixCholeskySolve(CMatrixDouble &cha,const int n,
const bool IsUpper,CRowDouble &b,
int &info,CDenseSolverReport &rep,
CRowDouble &x)
{
//--- create matrix
CMatrixDouble bm;
CMatrixDouble xm;
//--- initialization
info=0;
//--- check
if(n<=0)
{
info=-1;
return;
}
//--- allocation
bm.Resize(n,1);
//--- filling
bm.Col(0,b);
//--- function call
SPDMatrixCholeskySolveM(cha,n,IsUpper,bm,1,info,rep,xm);
//--- copy
x=xm.Col(0)+0;
x.Resize(n);
}
//+------------------------------------------------------------------+
//| 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 |
//+------------------------------------------------------------------+
void CDenseSolver::HPDMatrixSolveM(CMatrixComplex &a,const int n,
const bool IsUpper,CMatrixComplex &b,
const int m,int &info,
CDenseSolverReport &rep,CMatrixComplex &x)
{
//--- create variables
double sqrtscalea=0;
int i=0;
int j=0;
int j1=0;
int j2=0;
int i_=0;
//--- create matrix
CMatrixComplex da;
//--- initialization
info=0;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
//--- allocation
da.Resize(n,n);
//--- 1. scale matrix,max(|A[i][j]|)
//--- 2. factorize scaled matrix
//--- 3. solve
sqrtscalea=0;
for(i=0; i<n; i++)
{
//--- check
if(IsUpper)
{
j1=i;
j2=n-1;
}
else
{
j1=0;
j2=i;
}
//--- calculation
for(j=j1; j<=j2; j++)
sqrtscalea=MathMax(sqrtscalea,CMath::AbsComplex(a[i][j]));
}
//--- check
if(sqrtscalea==0.0)
sqrtscalea=1;
//--- change values
sqrtscalea=1/sqrtscalea;
sqrtscalea=MathSqrt(sqrtscalea);
for(i=0; i<n; i++)
{
//--- check
if(IsUpper)
{
j1=i;
j2=n-1;
}
else
{
j1=0;
j2=i;
}
//--- calculation
for(i_=j1; i_<=j2; i_++)
da.Set(i,i_,a[i][i_]);
}
//--- check
if(!CTrFac::HPDMatrixCholesky(da,n,IsUpper))
{
//--- allocation
x.Resize(n,m);
for(i=0; i<n; i++)
{
for(j=0; j<=m-1; j++)
x.Set(i,j,0.0);
}
//--- change values
rep.m_r1=0;
rep.m_rinf=0;
info=-3;
//--- exit the function
return;
}
info=1;
//--- function call
HPDMatrixCholeskySolveInternal(da,sqrtscalea,n,IsUpper,a,true,b,m,info,rep,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 |
//+------------------------------------------------------------------+
void CDenseSolver::HPDMatrixSolve(CMatrixComplex &a,const int n,
const bool IsUpper,complex &b[],
int &info,CDenseSolverReport &rep,
complex &x[])
{
//--- create matrix
CMatrixComplex bm;
CMatrixComplex xm;
//--- initialization
info=0;
//--- check
if(n<=0)
{
info=-1;
return;
}
//--- allocation
bm.Resize(n,1);
//--- filling
for(int i_=0; i_<n; i_++)
bm.Set(i_,0,b[i_]);
//--- function call
HPDMatrixSolveM(a,n,IsUpper,bm,1,info,rep,xm);
//--- allocation
ArrayResize(x,n);
//--- copy
for(int i_=0; i_<n; i_++)
x[i_]=xm[i_][0];
}
//+------------------------------------------------------------------+
//| 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 |
//+------------------------------------------------------------------+
void CDenseSolver::HPDMatrixCholeskySolveM(CMatrixComplex &cha,
const int n,const bool IsUpper,
CMatrixComplex &b,const int m,
int &info,CDenseSolverReport &rep,
CMatrixComplex &x)
{
//--- create variables
double sqrtscalea=0;
int j1=0;
int j2=0;
//--- create matrix
CMatrixComplex emptya;
//--- initialization
info=0;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
//--- 1. scale matrix,max(|U[i][j]|)
//--- 2. factorize scaled matrix
//--- 3. solve
sqrtscalea=0;
for(int i=0; i<n; i++)
{
//--- check
if(IsUpper)
{
j1=i;
j2=n-1;
}
else
{
j1=0;
j2=i;
}
//--- calculation
for(int j=j1; j<=j2; j++)
{
sqrtscalea=MathMax(sqrtscalea,CMath::AbsComplex(cha[i][j]));
}
}
//--- check
if(sqrtscalea==0.0)
{
sqrtscalea=1;
}
//--- change values
sqrtscalea=1/sqrtscalea;
//--- function call
HPDMatrixCholeskySolveInternal(cha,sqrtscalea,n,IsUpper,emptya,false,b,m,info,rep,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 |
//+------------------------------------------------------------------+
void CDenseSolver::HPDMatrixCholeskySolve(CMatrixComplex &cha,
const int n,const bool IsUpper,
complex &b[],int &info,
CDenseSolverReport &rep,
complex &x[])
{
//--- create matrix
CMatrixComplex bm;
CMatrixComplex xm;
//--- initialization
info=0;
//--- check
if(n<=0)
{
info=-1;
return;
}
//--- allocation
bm.Resize(n,1);
//--- filling
for(int i_=0; i_<n; i_++)
bm.Set(i_,0,b[i_]);
//--- function call
HPDMatrixCholeskySolveM(cha,n,IsUpper,bm,1,info,rep,xm);
//--- allocation
ArrayResize(x,n);
//--- copy
for(int i_=0; i_<n; i_++)
x[i_]=xm[i_][0];
}
//+------------------------------------------------------------------+
//| Dense solver. |
//| 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 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 |
//| Algorithm features: |
//| * automatic detection of degenerate cases |
//| * iterative refinement |
//| * O(N^3) complexity |
//| INPUT PARAMETERS |
//| A - array[0..NRows-1,0..NCols-1], system matrix |
//| NRows - vertical size of A |
//| NCols - horizontal size of A |
//| B - array[0..NCols-1], right part |
//| Threshold- a number in [0,1]. Singular values beyond |
//| Threshold are considered zero. Set it to 0.0, |
//| if you don't understand what it means, so the |
//| solver will choose good value on its own. |
//| OUTPUT PARAMETERS |
//| Info - return code: |
//| * -4 SVD subroutine failed |
//| * -1 if NRows<=0 or NCols<=0 or Threshold<0 |
//| was passed |
//| * 1 if task is solved |
//| Rep - solver report, see below for more info |
//| X - array[0..N-1,0..M-1], it contains: |
//| * solution of A*X=B if A is non-singular |
//| (well-conditioned or ill-conditioned, but not |
//| very close to singular) |
//| * zeros, if A is singular or VERY close to |
//| singular (in this case Info=-3). |
//| SOLVER REPORT |
//| Subroutine sets following fields of the Rep structure: |
//| * R2 reciprocal of condition number: 1/cond(A), 2-norm. |
//| * N = NCols |
//| * K dim(Null(A)) |
//| * CX array[0..N-1,0..K-1], kernel of A. |
//| Columns of CX store such vectors that A*CX[i]=0. |
//+------------------------------------------------------------------+
void CDenseSolver::RMatrixSolveLS(CMatrixDouble &a,const int nrows,
const int ncols,double &b[],
double threshold,int &info,
CDenseSolverLSReport &rep,
double &x[])
{
//--- create matrix
CMatrixDouble u;
CMatrixDouble vt;
//--- create arrays
double sv[];
double rp[];
double utb[];
double sutb[];
double tmp[];
double ta[];
double tx[];
double buf[];
double w[];
//--- create variables
int i=0;
int j=0;
int nsv=0;
int kernelidx=0;
double v=0;
double verr=0;
bool svdfailed;
bool zeroa;
int rfs=0;
int nrfs=0;
bool terminatenexttime;
bool smallerr;
int i_=0;
//--- initialization
info=0;
//--- check
if((nrows<=0 || ncols<=0) || threshold<0.0)
{
info=-1;
return;
}
//--- check
if(threshold==0.0)
threshold=1000*CMath::m_machineepsilon;
//--- Factorize A first
svdfailed=!CSingValueDecompose::RMatrixSVD(a,nrows,ncols,1,2,2,sv,u,vt);
//--- check
if(sv[0]==0.0)
zeroa=true;
else
zeroa=false;
//--- check
if(svdfailed || zeroa)
{
//--- check
if(svdfailed)
info=-4;
else
info=1;
//--- allocation
ArrayResize(x,ncols);
for(i=0; i<=ncols-1; i++)
x[i]=0;
//--- change values
rep.m_n=ncols;
rep.m_k=ncols;
rep.m_cx.Resize(ncols,ncols);
for(i=0; i<=ncols-1; i++)
{
for(j=0; j<=ncols-1; j++)
{
//--- check
if(i==j)
rep.m_cx.Set(i,j,1);
else
rep.m_cx.Set(i,j,0);
}
}
rep.m_r2=0;
//--- exit the function
return;
}
nsv=MathMin(ncols,nrows);
//--- check
if(nsv==ncols)
rep.m_r2=sv[nsv-1]/sv[0];
else
rep.m_r2=0;
//--- change values
rep.m_n=ncols;
info=1;
//--- Iterative refinement of xc combined with solution:
//--- 1. xc=0
//--- 2. calculate r=bc-A*xc using extra-precise dot product
//--- 3. solve A*y=r
//--- 4. update x:=x+r
//--- 5. goto 2
//--- This cycle is executed until one of two things happens:
//--- 1. maximum number of iterations reached
//--- 2. last iteration decreased error to the lower limit
ArrayResize(utb,nsv);
ArrayResize(sutb,nsv);
ArrayResize(x,ncols);
ArrayResize(tmp,ncols);
ArrayResize(ta,ncols+1);
ArrayResize(tx,ncols+1);
ArrayResize(buf,ncols+1);
//--- initialization
for(i=0; i<=ncols-1; i++)
x[i]=0;
kernelidx=nsv;
for(i=0; i<=nsv-1; i++)
{
//--- check
if(sv[i]<=threshold*sv[0])
{
kernelidx=i;
break;
}
}
//--- change values
rep.m_k=ncols-kernelidx;
nrfs=CDenseSolverRFSMaxV2(ncols,rep.m_r2);
terminatenexttime=false;
//--- allocation
ArrayResize(rp,nrows);
for(rfs=0; rfs<=nrfs; rfs++)
{
//--- check
if(terminatenexttime)
break;
//--- calculate right part
if(rfs==0)
{
for(i_=0; i_<=nrows-1; i_++)
rp[i_]=b[i_];
}
else
{
smallerr=true;
for(i=0; i<=nrows-1; i++)
{
//--- copy
for(i_=0; i_<=ncols-1; i_++)
ta[i_]=a[i][i_];
ta[ncols]=-1;
//--- copy
for(i_=0; i_<=ncols-1; i_++)
tx[i_]=x[i_];
tx[ncols]=b[i];
//--- function call
CXblas::XDot(ta,tx,ncols+1,buf,v,verr);
rp[i]=-v;
smallerr=smallerr && MathAbs(v)<4*verr;
}
//--- check
if(smallerr)
terminatenexttime=true;
}
//--- solve A*dx=rp
for(i=0; i<=ncols-1; i++)
tmp[i]=0;
for(i=0; i<=nsv-1; i++)
utb[i]=0;
//--- change values
for(i=0; i<=nrows-1; i++)
{
v=rp[i];
for(i_=0; i_<=nsv-1; i_++)
utb[i_]=utb[i_]+v*u[i][i_];
}
for(i=0; i<=nsv-1; i++)
{
//--- check
if(i<kernelidx)
sutb[i]=utb[i]/sv[i];
else
sutb[i]=0;
}
//--- change values
for(i=0; i<=nsv-1; i++)
{
v=sutb[i];
for(i_=0; i_<=ncols-1; i_++)
tmp[i_]=tmp[i_]+v*vt[i][i_];
}
//--- update x: x:=x+dx
for(i_=0; i_<=ncols-1; i_++)
x[i_]=x[i_]+tmp[i_];
}
//--- fill CX
if(rep.m_k>0)
{
//--- allocation
rep.m_cx.Resize(ncols,rep.m_k);
for(i=0; i<=rep.m_k-1; i++)
{
for(i_=0; i_<=ncols-1; i_++)
rep.m_cx.Set(i_,i,vt[kernelidx+i][i_]);
}
}
}
//+------------------------------------------------------------------+
//| Internal LU solver |
//+------------------------------------------------------------------+
void CDenseSolver::RMatrixLUSolveInternal(CMatrixDouble &lua,int &p[],
const double scalea,const int n,
CMatrixDouble &a,const bool havea,
CMatrixDouble &b,const int m,
int &info,CDenseSolverReport &rep,
CMatrixDouble &x)
{
//--- create variables
int i=0;
int j=0;
int k=0;
int rfs=0;
int nrfs=0;
double v=0;
double verr=0;
double mxb=0;
double scaleright=0;
bool smallerr;
bool terminatenexttime;
int i_=0;
//--- create arrays
double xc[];
double y[];
double bc[];
double xa[];
double xb[];
double tx[];
//--- initialization
info=0;
//--- check
if(!CAp::Assert(scalea>0.0))
return;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
for(i=0; i<n; i++)
{
//--- check
if(p[i]>n-1 || p[i]<i)
{
info=-1;
return;
}
}
//--- allocation
x.Resize(n,m);
ArrayResize(y,n);
ArrayResize(xc,n);
ArrayResize(bc,n);
ArrayResize(tx,n+1);
ArrayResize(xa,n+1);
ArrayResize(xb,n+1);
//--- estimate condition number,test for near singularity
rep.m_r1=CRCond::RMatrixLURCond1(lua,n);
rep.m_rinf=CRCond::RMatrixLURCondInf(lua,n);
//--- check
if(rep.m_r1<CRCond::RCondThreshold() || rep.m_rinf<CRCond::RCondThreshold())
{
for(i=0; i<n; i++)
{
for(j=0; j<=m-1; j++)
x.Set(i,j,0);
}
//--- change values
rep.m_r1=0;
rep.m_rinf=0;
info=-3;
//--- exit the function
return;
}
info=1;
//--- solve
for(k=0; k<=m-1; k++)
{
//--- copy B to contiguous storage
for(i_=0; i_<n; i_++)
bc[i_]=b[i_][k];
//--- Scale right part:
//--- * MX stores max(|Bi|)
//--- * ScaleRight stores actual scaling applied to B when solving systems
//--- it is chosen to make |scaleRight*b| close to 1.
mxb=0;
for(i=0; i<n; i++)
mxb=MathMax(mxb,MathAbs(bc[i]));
//--- check
if(mxb==0.0)
mxb=1;
//--- change values
scaleright=1/mxb;
//--- First,non-iterative part of solution process.
//--- We use separate code for this task because
//--- XDot is quite slow and we want to save time.
for(i_=0; i_<n; i_++)
xc[i_]=bc[i_]*scaleright;
//--- function call
RBasicLUSolve(lua,p,scalea,n,xc,tx);
//--- Iterative refinement of xc:
//--- * calculate r=bc-A*xc using extra-precise dot product
//--- * solve A*y=r
//--- * update x:=x+r
//--- This cycle is executed until one of two things happens:
//--- 1. maximum number of iterations reached
//--- 2. last iteration decreased error to the lower limit
if(havea)
{
//--- calculation
nrfs=CDenseSolverRFSMax(n,rep.m_r1,rep.m_rinf);
terminatenexttime=false;
for(rfs=0; rfs<=nrfs-1; rfs++)
{
//--- check
if(terminatenexttime)
break;
//--- generate right part
smallerr=true;
for(i_=0; i_<n; i_++)
xb[i_]=xc[i_];
for(i=0; i<n; i++)
{
//--- change values
for(i_=0; i_<n; i_++)
xa[i_]=a[i][i_]*scalea;
xa[n]=-1;
xb[n]=bc[i]*scaleright;
//--- function call
CXblas::XDot(xa,xb,n+1,tx,v,verr);
y[i]=-v;
smallerr=smallerr && MathAbs(v)<4*verr;
}
//--- check
if(smallerr)
terminatenexttime=true;
//--- solve and update
RBasicLUSolve(lua,p,scalea,n,y,tx);
for(i_=0; i_<n; i_++)
xc[i_]=xc[i_]+y[i_];
}
}
//--- Store xc.
//--- Post-scale result.
v=scalea*mxb;
for(i_=0; i_<n; i_++)
x.Set(i_,k,xc[i_]*v);
}
}
//+------------------------------------------------------------------+
//| Internal Cholesky solver |
//+------------------------------------------------------------------+
void CDenseSolver::SPDMatrixCholeskySolveInternal(CMatrixDouble &cha,
const double sqrtscalea,
const int n,const bool IsUpper,
CMatrixDouble &a,
const bool havea,
CMatrixDouble &b,
const int m,int &info,
CDenseSolverReport &rep,
CMatrixDouble &x)
{
//--- create variables
int i=0;
int j=0;
int k=0;
double v=0;
double mxb=0;
double scaleright=0;
int i_=0;
//--- create arrays
double xc[];
double y[];
double bc[];
double xa[];
double xb[];
double tx[];
//--- initialization
info=0;
//--- check
if(!CAp::Assert(sqrtscalea>0.0))
return;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
//--- allocation
x.Resize(n,m);
ArrayResize(y,n);
ArrayResize(xc,n);
ArrayResize(bc,n);
ArrayResize(tx,n+1);
ArrayResize(xa,n+1);
ArrayResize(xb,n+1);
//--- estimate condition number,test for near singularity
rep.m_r1=CRCond::SPDMatrixCholeskyRCond(cha,n,IsUpper);
rep.m_rinf=rep.m_r1;
//--- check
if(rep.m_r1<CRCond::RCondThreshold())
{
for(i=0; i<n; i++)
{
for(j=0; j<=m-1; j++)
x.Set(i,j,0);
}
//--- change values
rep.m_r1=0;
rep.m_rinf=0;
info=-3;
//--- exit the function
return;
}
info=1;
//--- solve
for(k=0; k<=m-1; k++)
{
//--- copy B to contiguous storage
for(i_=0; i_<n; i_++)
bc[i_]=b[i_][k];
//--- Scale right part:
//--- * MX stores max(|Bi|)
//--- * ScaleRight stores actual scaling applied to B when solving systems
//--- it is chosen to make |scaleRight*b| close to 1.
mxb=0;
for(i=0; i<n; i++)
mxb=MathMax(mxb,MathAbs(bc[i]));
//--- check
if(mxb==0.0)
mxb=1;
scaleright=1/mxb;
//--- First,non-iterative part of solution process.
//--- We use separate code for this task because
//--- XDot is quite slow and we want to save time.
for(i_=0; i_<n; i_++)
xc[i_]=bc[i_]*scaleright;
//--- function call
SPDBasicCholeskySolve(cha,sqrtscalea,n,IsUpper,xc,tx);
//--- Store xc.
//--- Post-scale result.
v=CMath::Sqr(sqrtscalea)*mxb;
for(i_=0; i_<n; i_++)
x.Set(i_,k,xc[i_]*v);
}
}
//+------------------------------------------------------------------+
//| Internal LU solver |
//+------------------------------------------------------------------+
void CDenseSolver::CMatrixLUSolveInternal(CMatrixComplex &lua,int &p[],
const double scalea,const int n,
CMatrixComplex &a,const bool havea,
CMatrixComplex &b,const int m,
int &info,CDenseSolverReport &rep,
CMatrixComplex &x)
{
//--- create variables
int i=0;
int j=0;
int k=0;
int rfs=0;
int nrfs=0;
complex v=0;
double verr=0;
double mxb=0;
double scaleright=0;
bool smallerr;
bool terminatenexttime;
int i_=0;
//--- create arrays
complex xc[];
complex y[];
complex bc[];
complex xa[];
complex xb[];
complex tx[];
double tmpbuf[];
//--- initialization
info=0;
//--- check
if(!CAp::Assert(scalea>0.0))
return;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
for(i=0; i<n; i++)
{
//--- check
if(p[i]>n-1 || p[i]<i)
{
info=-1;
return;
}
}
//--- allocation
x.Resize(n,m);
ArrayResize(y,n);
ArrayResize(xc,n);
ArrayResize(bc,n);
ArrayResize(tx,n);
ArrayResize(xa,n+1);
ArrayResize(xb,n+1);
ArrayResize(tmpbuf,2*n+2);
//--- estimate condition number,test for near singularity
rep.m_r1=CRCond::CMatrixLURCond1(lua,n);
rep.m_rinf=CRCond::CMatrixLURCondInf(lua,n);
//--- check
if(rep.m_r1<CRCond::RCondThreshold() || rep.m_rinf<CRCond::RCondThreshold())
{
for(i=0; i<n; i++)
{
for(j=0; j<=m-1; j++)
x.Set(i,j,0.0);
}
//--- change values
rep.m_r1=0;
rep.m_rinf=0;
info=-3;
//--- exit the function
return;
}
info=1;
//--- solve
for(k=0; k<=m-1; k++)
{
//--- copy B to contiguous storage
for(i_=0; i_<n; i_++)
bc[i_]=b[i_][k];
//--- Scale right part:
//--- * MX stores max(|Bi|)
//--- * ScaleRight stores actual scaling applied to B when solving systems
//--- it is chosen to make |scaleRight*b| close to 1.
mxb=0;
for(i=0; i<n; i++)
mxb=MathMax(mxb,CMath::AbsComplex(bc[i]));
//--- check
if(mxb==0.0)
mxb=1;
scaleright=1/mxb;
//--- First,non-iterative part of solution process.
//--- We use separate code for this task because
//--- XDot is quite slow and we want to save time.
for(i_=0; i_<n; i_++)
xc[i_]=bc[i_]*scaleright;
//--- function call
CBasicLUSolve(lua,p,scalea,n,xc,tx);
//--- Iterative refinement of xc:
//--- * calculate r=bc-A*xc using extra-precise dot product
//--- * solve A*y=r
//--- * update x:=x+r
//--- This cycle is executed until one of two things happens:
//--- 1. maximum number of iterations reached
//--- 2. last iteration decreased error to the lower limit
if(havea)
{
//--- calculation
nrfs=CDenseSolverRFSMax(n,rep.m_r1,rep.m_rinf);
terminatenexttime=false;
for(rfs=0; rfs<=nrfs-1; rfs++)
{
//--- check
if(terminatenexttime)
break;
//--- generate right part
smallerr=true;
//--- copy
for(i_=0; i_<n; i_++)
xb[i_]=xc[i_];
for(i=0; i<n; i++)
{
//--- change values
for(i_=0; i_<n; i_++)
xa[i_]=a[i][i_]*scalea;
xa[n]=-1;
xb[n]=bc[i]*scaleright;
//--- function call
CXblas::XCDot(xa,xb,n+1,tmpbuf,v,verr);
y[i]=-v;
smallerr=smallerr && CMath::AbsComplex(v)<4*verr;
}
//--- check
if(smallerr)
terminatenexttime=true;
//--- solve and update
CBasicLUSolve(lua,p,scalea,n,y,tx);
for(i_=0; i_<n; i_++)
xc[i_]=xc[i_]+y[i_];
}
}
//--- Store xc.
//--- Post-scale result.
v=scalea*mxb;
for(i_=0; i_<n; i_++)
x.Set(i_,k,xc[i_]*v);
}
}
//+------------------------------------------------------------------+
//| Internal Cholesky solver |
//+------------------------------------------------------------------+
void CDenseSolver::HPDMatrixCholeskySolveInternal(CMatrixComplex &cha,
const double sqrtscalea,
const int n,
const bool IsUpper,
CMatrixComplex &a,
const bool havea,
CMatrixComplex &b,
const int m,int &info,
CDenseSolverReport &rep,
CMatrixComplex &x)
{
//--- create variables
int i=0;
int j=0;
int k=0;
double v=0;
double mxb=0;
double scaleright=0;
int i_=0;
//--- create arrays
complex xc[];
complex y[];
complex bc[];
complex xa[];
complex xb[];
complex tx[];
//--- initialization
info=0;
//--- check
if(!CAp::Assert(sqrtscalea>0.0))
return;
//--- prepare: check inputs,allocate space...
if(n<=0 || m<=0)
{
info=-1;
return;
}
//--- allocation
x.Resize(n,m);
ArrayResize(y,n);
ArrayResize(xc,n);
ArrayResize(bc,n);
ArrayResize(tx,n+1);
ArrayResize(xa,n+1);
ArrayResize(xb,n+1);
//--- estimate condition number,test for near singularity
rep.m_r1=CRCond::HPDMatrixCholeskyRCond(cha,n,IsUpper);
rep.m_rinf=rep.m_r1;
//--- check
if(rep.m_r1<CRCond::RCondThreshold())
{
for(i=0; i<n; i++)
{
for(j=0; j<=m-1; j++)
x.Set(i,j,0.0);
}
//--- change values
rep.m_r1=0;
rep.m_rinf=0;
info=-3;
//--- exit the function
return;
}
info=1;
//--- solve
for(k=0; k<=m-1; k++)
{
//--- copy B to contiguous storage
for(i_=0; i_<n; i_++)
bc[i_]=b[i_][k];
//--- Scale right part:
//--- * MX stores max(|Bi|)
//--- * ScaleRight stores actual scaling applied to B when solving systems
//--- it is chosen to make |scaleRight*b| close to 1.
mxb=0;
for(i=0; i<n; i++)
mxb=MathMax(mxb,CMath::AbsComplex(bc[i]));
//--- check
if(mxb==0.0)
mxb=1;
scaleright=1/mxb;
//--- First,non-iterative part of solution process.
//--- We use separate code for this task because
//--- XDot is quite slow and we want to save time.
for(i_=0; i_<n; i_++)
xc[i_]=bc[i_]*scaleright;
//--- function call
HPDBasicCholeskySolve(cha,sqrtscalea,n,IsUpper,xc,tx);
//--- Store xc.
//--- Post-scale result.
v=CMath::Sqr(sqrtscalea)*mxb;
for(i_=0; i_<n; i_++)
x.Set(i_,k,xc[i_]*v);
}
}
//+------------------------------------------------------------------+
//| Internal subroutine. |
//| Returns maximum count of RFS iterations as function of: |
//| 1. machine epsilon |
//| 2. task size. |
//| 3. condition number |
//+------------------------------------------------------------------+
int CDenseSolver::CDenseSolverRFSMax(const int n,const double r1,
const double rinf)
{
//--- return result
return(5);
}
//+------------------------------------------------------------------+
//| Internal subroutine. |
//| Returns maximum count of RFS iterations as function of: |
//| 1. machine epsilon |
//| 2. task size. |
//| 3. norm-2 condition number |
//+------------------------------------------------------------------+
int CDenseSolver::CDenseSolverRFSMaxV2(const int n,const double r2)
{
//--- return result
return(CDenseSolverRFSMax(n,0,0));
}
//+------------------------------------------------------------------+
//| Basic LU solver for ScaleA*PLU*x = y. |
//| This subroutine assumes that: |
//| * L is well-scaled, and it is U which needs scaling by ScaleA. |
//| * A=PLU is well-conditioned, so no zero divisions or overflow may|
//| occur |
//+------------------------------------------------------------------+
void CDenseSolver::RBasicLUSolve(CMatrixDouble &lua,int &p[],
const double scalea,const int n,
double &xb[],double &tmp[])
{
//--- create variables
int i=0;
double v=0;
int i_=0;
//--- swap
for(i=0; i<n; i++)
{
//--- check
if(p[i]!=i)
{
v=xb[i];
xb[i]=xb[p[i]];
xb[p[i]]=v;
}
}
//--- calculation
for(i=1; i<n; i++)
{
v=0.0;
//--- change value
for(i_=0; i_<=i-1; i_++)
v+=lua[i][i_]*xb[i_];
//--- shift
xb[i]=xb[i]-v;
}
//--- change value
xb[n-1]=xb[n-1]/(lua[n-1][n-1]*scalea);
for(i=n-2; i>=0; i--)
{
//--- calculation
for(i_=i+1; i_<n; i_++)
tmp[i_]=lua[i][i_]*scalea;
v=0.0;
//--- change value
for(i_=i+1; i_<n; i_++)
v+=tmp[i_]*xb[i_];
//--- get result
xb[i]=(xb[i]-v)/(lua[i][i]*scalea);
}
}
//+------------------------------------------------------------------+
//| Basic Cholesky solver for ScaleA*Cholesky(A)'*x = y. |
//| This subroutine assumes that: |
//| * A*ScaleA is well scaled |
//| * A is well-conditioned, so no zero divisions or overflow may |
//| occur |
//+------------------------------------------------------------------+
void CDenseSolver::SPDBasicCholeskySolve(CMatrixDouble &cha,
const double sqrtscalea,
const int n,const bool IsUpper,
double &xb[],double &tmp[])
{
//--- create variables
int i=0;
double v=0;
int i_=0;
//--- A=L*L' or A=U'*U
if(IsUpper)
{
//--- Solve U'*y=b first.
for(i=0; i<n; i++)
{
xb[i]=xb[i]/(sqrtscalea*cha[i][i]);
//--- check
if(i<n-1)
{
v=xb[i];
//--- calculation
for(i_=i+1; i_<n; i_++)
tmp[i_]=sqrtscalea*cha[i][i_];
for(i_=i+1; i_<n; i_++)
xb[i_]=xb[i_]-v*tmp[i_];
}
}
//--- Solve U*x=y then.
for(i=n-1; i>=0; i--)
{
//--- check
if(i<n-1)
{
for(i_=i+1; i_<n; i_++)
tmp[i_]=sqrtscalea*cha[i][i_];
v=0.0;
//--- change value
for(i_=i+1; i_<n; i_++)
v+=tmp[i_]*xb[i_];
//--- shift
xb[i]=xb[i]-v;
}
xb[i]=xb[i]/(sqrtscalea*cha[i][i]);
}
}
else
{
//--- Solve L*y=b first
for(i=0; i<n; i++)
{
//--- check
if(i>0)
{
for(i_=0; i_<=i-1; i_++)
tmp[i_]=sqrtscalea*cha[i][i_];
//--- change value
v=0.0;
for(i_=0; i_<=i-1; i_++)
v+=tmp[i_]*xb[i_];
//--- shift
xb[i]=xb[i]-v;
}
xb[i]=xb[i]/(sqrtscalea*cha[i][i]);
}
//--- Solve L'*x=y then.
for(i=n-1; i>=0; i--)
{
xb[i]=xb[i]/(sqrtscalea*cha[i][i]);
//--- check
if(i>0)
{
v=xb[i];
//--- calculation
for(i_=0; i_<=i-1; i_++)
tmp[i_]=sqrtscalea*cha[i][i_];
for(i_=0; i_<=i-1; i_++)
xb[i_]=xb[i_]-v*tmp[i_];
}
}
}
}
//+------------------------------------------------------------------+
//| Basic LU solver for ScaleA*PLU*x = y. |
//| This subroutine assumes that: |
//| * L is well-scaled, and it is U which needs scaling by ScaleA. |
//| * A=PLU is well-conditioned, so no zero divisions or overflow may|
//| occur |
//+------------------------------------------------------------------+
void CDenseSolver::CBasicLUSolve(CMatrixComplex &lua,int &p[],
const double scalea,const int n,
complex &xb[],complex &tmp[])
{
//--- create variables
int i=0;
complex v=0;
int i_=0;
//--- swap
for(i=0; i<n; i++)
{
//--- check
if(p[i]!=i)
{
v=xb[i];
xb[i]=xb[p[i]];
xb[p[i]]=v;
}
}
for(i=1; i<n; i++)
{
v=0.0;
//--- calculation
for(i_=0; i_<=i-1; i_++)
v+=lua[i][i_]*xb[i_];
//--- shift
xb[i]=xb[i]-v;
}
//--- change values
xb[n-1]=xb[n-1]/(lua[n-1][n-1]*scalea);
for(i=n-2; i>=0; i--)
{
for(i_=i+1; i_<n; i_++)
tmp[i_]=lua[i][i_]*scalea;
//--- calculation
v=0.0;
for(i_=i+1; i_<n; i_++)
v+=tmp[i_]*xb[i_];
//--- get result
xb[i]=(xb[i]-v)/(lua[i][i]*scalea);
}
}
//+------------------------------------------------------------------+
//| Basic Cholesky solver for ScaleA*Cholesky(A)'*x = y. |
//| This subroutine assumes that: |
//| * A*ScaleA is well scaled |
//| * A is well-conditioned, so no zero divisions or overflow may |
//| occur |
//+------------------------------------------------------------------+
void CDenseSolver::HPDBasicCholeskySolve(CMatrixComplex &cha,
const double sqrtscalea,
const int n,const bool IsUpper,
complex &xb[],complex &tmp[])
{
//--- create variables
int i=0;
complex v=0;
int i_=0;
//--- A=L*L' or A=U'*U
if(IsUpper)
{
//--- Solve U'*y=b first.
for(i=0; i<n; i++)
{
xb[i]=xb[i]/(CMath::Conj(cha[i][i])*sqrtscalea);
//--- check
if(i<n-1)
{
v=xb[i];
//--- calculation
for(i_=i+1; i_<n; i_++)
tmp[i_]=CMath::Conj(cha[i][i_])*sqrtscalea;
for(i_=i+1; i_<n; i_++)
xb[i_]=xb[i_]-v*tmp[i_];
}
}
//--- Solve U*x=y then.
for(i=n-1; i>=0; i--)
{
//--- check
if(i<n-1)
{
for(i_=i+1; i_<n; i_++)
tmp[i_]=cha[i][i_]*sqrtscalea;
//--- change value
v=0.0;
for(i_=i+1; i_<n; i_++)
v+=tmp[i_]*xb[i_];
//--- shift
xb[i]=xb[i]-v;
}
xb[i]=xb[i]/(cha[i][i]*sqrtscalea);
}
}
else
{
//--- Solve L*y=b first
for(i=0; i<n; i++)
{
//--- check
if(i>0)
{
for(i_=0; i_<=i-1; i_++)
tmp[i_]=cha[i][i_]*sqrtscalea;
//--- change value
v=0.0;
for(i_=0; i_<=i-1; i_++)
v+=tmp[i_]*xb[i_];
//--- shift
xb[i]=xb[i]-v;
}
xb[i]=xb[i]/(cha[i][i]*sqrtscalea);
}
//--- Solve L'*x=y then.
for(i=n-1; i>=0; i--)
{
xb[i]=xb[i]/(CMath::Conj(cha[i][i])*sqrtscalea);
//--- check
if(i>0)
{
v=xb[i];
//--- calculation
for(i_=0; i_<=i-1; i_++)
tmp[i_]=CMath::Conj(cha[i][i_])*sqrtscalea;
for(i_=0; i_<=i-1; i_++)
xb[i_]=xb[i_]-v*tmp[i_];
}
}
}
}
//+------------------------------------------------------------------+
//| Auxiliary class for CNlEq |
//+------------------------------------------------------------------+
class CNlEqState
{
public:
//--- variables
int m_n;
int m_m;
double m_epsf;
int m_maxits;
bool m_xrep;
double m_stpmax;
double m_f;
bool m_needf;
bool m_needfij;
bool m_xupdated;
RCommState m_rstate;
int m_repiterationscount;
int m_repnfunc;
int m_repnjac;
int m_repterminationtype;
double m_fbase;
double m_fprev;
//--- arrays
double m_x[];
double m_fi[];
double m_xbase[];
double m_candstep[];
double m_rightpart[];
double m_cgbuf[];
//--- matrix
CMatrixDouble m_j;
//--- constructor, destructor
CNlEqState(void);
~CNlEqState(void) {}
//--- copy
void Copy(CNlEqState &obj);
};
//+------------------------------------------------------------------+
//| Constructor |
//+------------------------------------------------------------------+
CNlEqState::CNlEqState(void)
{
m_n=0;
m_m=0;
m_epsf=0;
m_maxits=0;
m_xrep=false;
m_stpmax=0;
m_f=0;
m_needf=false;
m_needfij=false;
m_xupdated=false;
m_repiterationscount=0;
m_repnfunc=0;
m_repnjac=0;
m_repterminationtype=0;
m_fbase=0;
m_fprev=0;
}
//+------------------------------------------------------------------+
//| Copy |
//+------------------------------------------------------------------+
void CNlEqState::Copy(CNlEqState &obj)
{
//--- copy variables
m_n=obj.m_n;
m_m=obj.m_m;
m_epsf=obj.m_epsf;
m_maxits=obj.m_maxits;
m_xrep=obj.m_xrep;
m_stpmax=obj.m_stpmax;
m_f=obj.m_f;
m_needf=obj.m_needf;
m_needfij=obj.m_needfij;
m_xupdated=obj.m_xupdated;
m_repiterationscount=obj.m_repiterationscount;
m_repnfunc=obj.m_repnfunc;
m_repnjac=obj.m_repnjac;
m_repterminationtype=obj.m_repterminationtype;
m_fbase=obj.m_fbase;
m_fprev=obj.m_fprev;
m_rstate.Copy(obj.m_rstate);
//--- copy arrays
ArrayCopy(m_x,obj.m_x);
ArrayCopy(m_fi,obj.m_fi);
ArrayCopy(m_xbase,obj.m_xbase);
ArrayCopy(m_candstep,obj.m_candstep);
ArrayCopy(m_rightpart,obj.m_rightpart);
ArrayCopy(m_cgbuf,obj.m_cgbuf);
//--- copy matrix
m_j=obj.m_j;
}
//+------------------------------------------------------------------+
//| This class is a shell for class CNlEqState |
//+------------------------------------------------------------------+
class CNlEqStateShell
{
private:
CNlEqState m_innerobj;
public:
//--- constructors, destructor
CNlEqStateShell(void) {}
CNlEqStateShell(CNlEqState &obj) { m_innerobj.Copy(obj); }
~CNlEqStateShell(void) {}
//--- methods
bool GetNeedF(void);
void SetNeedF(const bool b);
bool GetNeedFIJ(void);
void SetNeedFIJ(const bool b);
bool GetXUpdated(void);
void SetXUpdated(const bool b);
double GetF(void);
void SetF(const double d);
CNlEqState *GetInnerObj(void);
};
//+------------------------------------------------------------------+
//| Returns the value of the variable needf |
//+------------------------------------------------------------------+
bool CNlEqStateShell::GetNeedF(void)
{
return(m_innerobj.m_needf);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable needf |
//+------------------------------------------------------------------+
void CNlEqStateShell::SetNeedF(const bool b)
{
m_innerobj.m_needf=b;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable needfij |
//+------------------------------------------------------------------+
bool CNlEqStateShell::GetNeedFIJ(void)
{
return(m_innerobj.m_needfij);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable needfij |
//+------------------------------------------------------------------+
void CNlEqStateShell::SetNeedFIJ(const bool b)
{
m_innerobj.m_needfij=b;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable xupdated |
//+------------------------------------------------------------------+
bool CNlEqStateShell::GetXUpdated(void)
{
return(m_innerobj.m_xupdated);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable xupdated |
//+------------------------------------------------------------------+
void CNlEqStateShell::SetXUpdated(const bool b)
{
m_innerobj.m_xupdated=b;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable f |
//+------------------------------------------------------------------+
double CNlEqStateShell::GetF(void)
{
return(m_innerobj.m_f);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable f |
//+------------------------------------------------------------------+
void CNlEqStateShell::SetF(const double d)
{
m_innerobj.m_f=d;
}
//+------------------------------------------------------------------+
//| Return object of class |
//+------------------------------------------------------------------+
CNlEqState *CNlEqStateShell::GetInnerObj(void)
{
return(GetPointer(m_innerobj));
}
//+------------------------------------------------------------------+
//| Auxiliary class for CNlEq |
//+------------------------------------------------------------------+
class CNlEqReport
{
public:
//--- variables
int m_iterationscount;
int m_nfunc;
int m_njac;
int m_terminationtype;
//--- constructor, destructor
CNlEqReport(void) { ZeroMemory(this); }
~CNlEqReport(void) {}
//--- copy
void Copy(CNlEqReport &obj);
};
//+------------------------------------------------------------------+
//| Copy |
//+------------------------------------------------------------------+
void CNlEqReport::Copy(CNlEqReport &obj)
{
//--- copy variables
m_iterationscount=obj.m_iterationscount;
m_nfunc=obj.m_nfunc;
m_njac=obj.m_njac;
m_terminationtype=obj.m_terminationtype;
}
//+------------------------------------------------------------------+
//| This class is a shell for class CNlEqReport |
//+------------------------------------------------------------------+
class CNlEqReportShell
{
private:
CNlEqReport m_innerobj;
public:
//--- constructors, destructor
CNlEqReportShell(void) {}
CNlEqReportShell(CNlEqReport &obj) { m_innerobj.Copy(obj); }
~CNlEqReportShell(void) {}
//--- methods
int GetIterationsCount(void);
void SetIterationsCount(const int i);
int GetNFunc(void);
void SetNFunc(const int i);
int GetNJac(void);
void SetNJac(const int i);
int GetTerminationType(void);
void SetTerminationType(const int i);
CNlEqReport *GetInnerObj(void);
};
//+------------------------------------------------------------------+
//| Returns the value of the variable iterationscount |
//+------------------------------------------------------------------+
int CNlEqReportShell::GetIterationsCount(void)
{
return(m_innerobj.m_iterationscount);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable iterationscount |
//+------------------------------------------------------------------+
void CNlEqReportShell::SetIterationsCount(const int i)
{
m_innerobj.m_iterationscount=i;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable nfunc |
//+------------------------------------------------------------------+
int CNlEqReportShell::GetNFunc(void)
{
return(m_innerobj.m_nfunc);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable nfunc |
//+------------------------------------------------------------------+
void CNlEqReportShell::SetNFunc(const int i)
{
m_innerobj.m_nfunc=i;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable njac |
//+------------------------------------------------------------------+
int CNlEqReportShell::GetNJac(void)
{
return(m_innerobj.m_njac);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable njac |
//+------------------------------------------------------------------+
void CNlEqReportShell::SetNJac(const int i)
{
m_innerobj.m_njac=i;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable terminationtype |
//+------------------------------------------------------------------+
int CNlEqReportShell::GetTerminationType(void)
{
return(m_innerobj.m_terminationtype);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable terminationtype |
//+------------------------------------------------------------------+
void CNlEqReportShell::SetTerminationType(const int i)
{
m_innerobj.m_terminationtype=i;
}
//+------------------------------------------------------------------+
//| Return object of class |
//+------------------------------------------------------------------+
CNlEqReport *CNlEqReportShell::GetInnerObj(void)
{
return(GetPointer(m_innerobj));
}
//+------------------------------------------------------------------+
//| Solving systems of nonlinear equations |
//+------------------------------------------------------------------+
class CNlEq
{
private:
//--- private methods
static void ClearRequestFields(CNlEqState &state);
static bool IncreaseLambda(double &lambdav,double &nu,const double lambdaup);
static void DecreaseLambda(double &lambdav,double &nu,const double lambdadown);
//--- auxiliary functions forNlEqiteration
static void Func_case_rcomm(CNlEqState &state,const int n,const int m,const int i,const bool b,const double lambdaup,const double lambdadown,const double lambdav,const double rho,const double mu,const double stepnorm);
static void Func_case_7(CNlEqState &state,const int n);
static bool Func_case_5(CNlEqState &state,double &lambdaup,double &lambdadown,double &lambdav,double &rho);
static bool Func_case_11(CNlEqState &state,const double stepnorm);
static int Func_case_10(CNlEqState &state,const int n,const int m,const int i,const bool b,const double lambdaup,const double lambdadown,const double lambdav,const double rho,const double mu,const double stepnorm);
static int Func_case_9(CNlEqState &state,int &n,int &m,int &i,bool &b,const double lambdaup,const double lambdadown,double &lambdav,const double rho,const double mu,double &stepnorm);
public:
//--- constant
static const int m_armijomaxfev;
//--- public methods
static void NlEqCreateLM(const int n,const int m,double &x[],CNlEqState &state);
static void NlEqSetCond(CNlEqState &state,double epsf,const int maxits);
static void NlEqSetXRep(CNlEqState &state,const bool needxrep);
static void NlEqSetStpMax(CNlEqState &state,const double stpmax);
static void NlEqResults(CNlEqState &state,double &x[],CNlEqReport &rep);
static void NlEqResultsBuf(CNlEqState &state,double &x[],CNlEqReport &rep);
static void NlEqRestartFrom(CNlEqState &state,double &x[]);
static bool NlEqIteration(CNlEqState &state);
};
//+------------------------------------------------------------------+
//| Initialize constant |
//+------------------------------------------------------------------+
const int CNlEq::m_armijomaxfev=20;
//+------------------------------------------------------------------+
//| LEVENBERG-MARQUARDT-LIKE NONLINEAR SOLVER |
//| DESCRIPTION: |
//| This algorithm solves system of nonlinear equations |
//| F[0](x[0], ..., x[n-1]) = 0 |
//| F[1](x[0], ..., x[n-1]) = 0 |
//| ... |
//| F[M-1](x[0], ..., x[n-1]) = 0 |
//| with M/N do not necessarily coincide. Algorithm converges |
//| quadratically under following conditions: |
//| * the solution set XS is nonempty |
//| * for some xs in XS there exist such neighbourhood N(xs) |
//| that: |
//| * vector function F(x) and its Jacobian J(x) are |
//| continuously differentiable on N |
//| * ||F(x)|| provides local error bound on N, i.e. there |
//| exists such c1, that ||F(x)||>c1*distance(x,XS) |
//| Note that these conditions are much more weaker than usual |
//| non-singularity conditions. For example, algorithm will converge |
//| for any affine function F (whether its Jacobian singular or not).|
//| REQUIREMENTS: |
//| Algorithm will request following information during its |
//| operation: |
//| * function vector F[] and Jacobian matrix at given point X |
//| * value of merit function f(x)=F[0]^2(x)+...+F[M-1]^2(x) at given|
//| point X |
//| USAGE: |
//| 1. User initializes algorithm state with NLEQCreateLM() call |
//| 2. User tunes solver parameters with NLEQSetCond(), |
//| NLEQSetStpMax() and other functions |
//| 3. User calls NLEQSolve() function which takes algorithm state |
//| and pointers (delegates, etc.) to callback functions which |
//| calculate merit function value and Jacobian. |
//| 4. User calls NLEQResults() to get solution |
//| 5. Optionally, user may call NLEQRestartFrom() to solve another |
//| problem with same parameters (N/M) but another starting point |
//| and/or another function vector. NLEQRestartFrom() allows to |
//| reuse already initialized structure. |
//| INPUT PARAMETERS: |
//| N - space dimension, N>1: |
//| * if provided, only leading N elements of X are |
//| used |
//| * if not provided, determined automatically from |
//| size of X |
//| M - system size |
//| X - starting point |
//| OUTPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//| NOTES: |
//| 1. you may tune stopping conditions with NLEQSetCond() function |
//| 2. if target function contains exp() or other fast growing |
//| functions, and optimization algorithm makes too large steps |
//| which leads to overflow, use NLEQSetStpMax() function to bound|
//| algorithm's steps. |
//| 3. this algorithm is a slightly modified implementation of the |
//| method described in 'Levenberg-Marquardt method for |
//| constrained nonlinear equations with strong local convergence |
//| properties' by Christian Kanzow Nobuo Yamashita and Masao |
//| Fukushima and further developed in 'On the convergence of a |
//| New Levenberg-Marquardt Method' by Jin-yan Fan and Ya-Xiang |
//| Yuan. |
//+------------------------------------------------------------------+
void CNlEq::NlEqCreateLM(const int n,const int m,double &x[],
CNlEqState &state)
{
//--- check
if(!CAp::Assert(n>=1,__FUNCTION__+": N<1!"))
return;
//--- check
if(!CAp::Assert(m>=1,__FUNCTION__+": M<1!"))
return;
//--- check
if(!CAp::Assert(CAp::Len(x)>=n,__FUNCTION__+": Length(X)<N!"))
return;
//--- check
if(!CAp::Assert(CApServ::IsFiniteVector(x,n),__FUNCTION__+": X contains infinite or NaN values!"))
return;
//--- Initialize
state.m_n=n;
state.m_m=m;
NlEqSetCond(state,0,0);
NlEqSetXRep(state,false);
NlEqSetStpMax(state,0);
//--- allocation
ArrayResize(state.m_x,n);
ArrayResize(state.m_xbase,n);
state.m_j.Resize(m,n);
ArrayResize(state.m_fi,m);
ArrayResize(state.m_rightpart,n);
ArrayResize(state.m_candstep,n);
//--- function call
NlEqRestartFrom(state,x);
}
//+------------------------------------------------------------------+
//| This function sets stopping conditions for the nonlinear solver |
//| INPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//| EpsF - >=0 |
//| The subroutine finishes its work if on k+1-th |
//| iteration the condition ||F||<=EpsF is satisfied |
//| MaxIts - maximum number of iterations. If MaxIts=0, the |
//| number of iterations is unlimited. |
//| Passing EpsF=0 and MaxIts=0 simultaneously will lead to |
//| automatic stopping criterion selection (small EpsF). |
//| NOTES: |
//+------------------------------------------------------------------+
void CNlEq::NlEqSetCond(CNlEqState &state,double epsf,const int maxits)
{
//--- check
if(!CAp::Assert(CMath::IsFinite(epsf),__FUNCTION__+": EpsF is not finite number!"))
return;
//--- check
if(!CAp::Assert(epsf>=0.0,__FUNCTION__+": negative EpsF!"))
return;
//--- check
if(!CAp::Assert(maxits>=0,__FUNCTION__+": negative MaxIts!"))
return;
//--- check
if(epsf==0.0 && maxits==0)
epsf=1.0E-6;
//--- change values
state.m_epsf=epsf;
state.m_maxits=maxits;
}
//+------------------------------------------------------------------+
//| This function turns on/off reporting. |
//| INPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//| NeedXRep- whether iteration reports are needed or not |
//| If NeedXRep is True, algorithm will call rep() callback function |
//| if it is provided to NLEQSolve(). |
//+------------------------------------------------------------------+
void CNlEq::NlEqSetXRep(CNlEqState &state,const bool needxrep)
{
//--- change value
state.m_xrep=needxrep;
}
//+------------------------------------------------------------------+
//| This function sets maximum step length |
//| INPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//| StpMax - maximum step length, >=0. Set StpMax to 0.0, if |
//| you don't want to limit step length. |
//| Use this subroutine when target function contains exp() or other |
//| fast growing functions, and algorithm makes too large steps which|
//| lead to overflow. This function allows us to reject steps that |
//| are too large (and therefore expose us to the possible overflow) |
//| without actually calculating function value at the x+stp*d. |
//+------------------------------------------------------------------+
void CNlEq::NlEqSetStpMax(CNlEqState &state,const double stpmax)
{
//--- check
if(!CAp::Assert(CMath::IsFinite(stpmax),__FUNCTION__+": StpMax is not finite!"))
return;
//--- check
if(!CAp::Assert(stpmax>=0.0,__FUNCTION__+": StpMax<0!"))
return;
//--- change value
state.m_stpmax=stpmax;
}
//+------------------------------------------------------------------+
//| NLEQ solver results |
//| INPUT PARAMETERS: |
//| State - algorithm state. |
//| OUTPUT PARAMETERS: |
//| X - array[0..N-1], solution |
//| Rep - optimization report: |
//| * Rep.TerminationType completetion code: |
//| * -4 ERROR: algorithm has converged to the|
//| stationary point Xf which is local |
//| minimum of f=F[0]^2+...+F[m-1]^2, |
//| but is not solution of nonlinear |
//| system. |
//| * 1 sqrt(f)<=EpsF. |
//| * 5 MaxIts steps was taken |
//| * 7 stopping conditions are too |
//| stringent, further improvement is |
//| impossible |
//| * Rep.IterationsCount contains iterations count |
//| * NFEV countains number of function calculations |
//| * ActiveConstraints contains number of active |
//| constraints |
//+------------------------------------------------------------------+
void CNlEq::NlEqResults(CNlEqState &state,double &x[],CNlEqReport &rep)
{
ArrayResize(x,0);
//--- function call
NlEqResultsBuf(state,x,rep);
}
//+------------------------------------------------------------------+
//| NLEQ solver results |
//| Buffered implementation of NLEQResults(), which uses |
//| pre-allocated buffer to store X[]. If buffer size is too small, |
//| it resizes buffer. It is intended to be used in the inner cycles |
//| of performance critical algorithms where array reallocation |
//| penalty is too large to be ignored. |
//+------------------------------------------------------------------+
void CNlEq::NlEqResultsBuf(CNlEqState &state,double &x[],CNlEqReport &rep)
{
//--- check
if(CAp::Len(x)<state.m_n)
ArrayResize(x,state.m_n);
//--- copy
for(int i_=0; i_<state.m_n; i_++)
x[i_]=state.m_xbase[i_];
//--- change values
rep.m_iterationscount=state.m_repiterationscount;
rep.m_nfunc=state.m_repnfunc;
rep.m_njac=state.m_repnjac;
rep.m_terminationtype=state.m_repterminationtype;
}
//+------------------------------------------------------------------+
//| This subroutine restarts CG algorithm from new point. All |
//| optimization parameters are left unchanged. |
//| This function allows to solve multiple optimization problems |
//| (which must have same number of dimensions) without object |
//| reallocation penalty. |
//| INPUT PARAMETERS: |
//| State - structure used for reverse communication |
//| previously allocated with MinCGCreate call. |
//| X - new starting point. |
//| BndL - new lower bounds |
//| BndU - new upper bounds |
//+------------------------------------------------------------------+
void CNlEq::NlEqRestartFrom(CNlEqState &state,double &x[])
{
//--- check
if(!CAp::Assert(CAp::Len(x)>=state.m_n,__FUNCTION__+": Length(X)<N!"))
return;
//--- check
if(!CAp::Assert(CApServ::IsFiniteVector(x,state.m_n),__FUNCTION__+": X contains infinite or NaN values!"))
return;
//--- copy
for(int i_=0; i_<state.m_n; i_++)
state.m_x[i_]=x[i_];
//--- allocation
state.m_rstate.ia.Resize(3);
ArrayResizeAL(state.m_rstate.ba,1);
state.m_rstate.ra.Resize(6);
state.m_rstate.stage=-1;
//--- function call
ClearRequestFields(state);
}
//+------------------------------------------------------------------+
//| Clears request fileds (to be sure that we don't forgot to clear |
//| something) |
//+------------------------------------------------------------------+
void CNlEq::ClearRequestFields(CNlEqState &state)
{
//--- change values
state.m_needf=false;
state.m_needfij=false;
state.m_xupdated=false;
}
//+------------------------------------------------------------------+
//| Increases lambda,returns False when there is a danger of |
//| overflow |
//+------------------------------------------------------------------+
bool CNlEq::IncreaseLambda(double &lambdav,double &nu,const double lambdaup)
{
//--- create variables
double lnlambda=MathLog(lambdav);
double lnlambdaup=MathLog(lambdaup);
double lnnu=MathLog(nu);
double lnmax=0.5*MathLog(CMath::m_maxrealnumber);
//--- check
if(lnlambda+lnlambdaup+lnnu>lnmax)
return(false);
//--- check
if(lnnu+MathLog(2)>lnmax)
return(false);
//--- change values
lambdav=lambdav*lambdaup*nu;
nu=nu*2;
//--- return result
return(true);
}
//+------------------------------------------------------------------+
//| Decreases lambda, but leaves it unchanged when there is danger of|
//| underflow. |
//+------------------------------------------------------------------+
void CNlEq::DecreaseLambda(double &lambdav,double &nu,const double lambdadown)
{
//--- initialization
nu=1;
//--- check
if(MathLog(lambdav)+MathLog(lambdadown)<MathLog(CMath::m_minrealnumber))
lambdav=CMath::m_minrealnumber;
else
lambdav=lambdav*lambdadown;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
bool CNlEq::NlEqIteration(CNlEqState &state)
{
//--- create variables
int n=0;
int m=0;
int i=0;
double lambdaup=0;
double lambdadown=0;
double lambdav=0;
double rho=0;
double mu=0;
double stepnorm=0;
bool b;
int i_=0;
int temp;
//--- This code initializes locals by:
//--- * random values determined during code
//--- generation - on first subroutine call
//--- * values from previous call - on subsequent calls
if(state.m_rstate.stage>=0)
{
//--- initialization
n=state.m_rstate.ia[0];
m=state.m_rstate.ia[1];
i=state.m_rstate.ia[2];
b=state.m_rstate.ba[0];
lambdaup=state.m_rstate.ra[0];
lambdadown=state.m_rstate.ra[1];
lambdav=state.m_rstate.ra[2];
rho=state.m_rstate.ra[3];
mu=state.m_rstate.ra[4];
stepnorm=state.m_rstate.ra[5];
}
else
{
//--- initialization
n=-983;
m=-989;
i=-834;
b=false;
lambdaup=-287;
lambdadown=364;
lambdav=214;
rho=-338;
mu=-686;
stepnorm=912;
}
//--- check
if(state.m_rstate.stage==0)
{
//--- change values
state.m_needf=false;
state.m_repnfunc=state.m_repnfunc+1;
//--- copy
for(i_=0; i_<n; i_++)
state.m_xbase[i_]=state.m_x[i_];
//--- change values
state.m_fbase=state.m_f;
state.m_fprev=CMath::m_maxrealnumber;
//--- check
if(!state.m_xrep)
{
//--- check
if(!Func_case_5(state,lambdaup,lambdadown,lambdav,rho))
return(false);
//--- function call
Func_case_7(state,n);
//--- Saving state
Func_case_rcomm(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- return result
return(true);
}
//--- progress report
ClearRequestFields(state);
state.m_xupdated=true;
state.m_rstate.stage=1;
//--- Saving state
Func_case_rcomm(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- return result
return(true);
}
//--- check
if(state.m_rstate.stage==1)
{
//--- change value
state.m_xupdated=false;
//--- check
if(!Func_case_5(state,lambdaup,lambdadown,lambdav,rho))
return(false);
//--- function call
Func_case_7(state,n);
//--- Saving state
Func_case_rcomm(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- return result
return(true);
}
//--- check
if(state.m_rstate.stage==2)
{
//--- change values
state.m_needfij=false;
state.m_repnfunc=state.m_repnfunc+1;
state.m_repnjac=state.m_repnjac+1;
//--- function call
CAblas::RMatrixMVect(n,m,state.m_j,0,0,1,state.m_fi,0,state.m_rightpart,0);
for(i_=0; i_<n; i_++)
state.m_rightpart[i_]=-1*state.m_rightpart[i_];
//--- Inner cycle: find good lambda
temp=Func_case_9(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- check
if(temp==-1)
return(false);
//--- check
if(temp==1)
return(true);
//--- Saving state
Func_case_rcomm(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- return result
return(true);
}
//--- check
if(state.m_rstate.stage==3)
{
//--- change values
state.m_needf=false;
state.m_repnfunc=state.m_repnfunc+1;
//--- check
if(state.m_f<state.m_fbase)
{
//--- function value decreased,move on
DecreaseLambda(lambdav,rho,lambdadown);
temp=Func_case_10(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- check
if(temp==-1)
return(false);
//--- check
if(temp==1)
return(true);
//--- Saving state
Func_case_rcomm(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- return result
return(true);
}
//--- check
if(!IncreaseLambda(lambdav,rho,lambdaup))
{
//--- Lambda is too large (near overflow),force zero step and break
stepnorm=0;
for(i_=0; i_<n; i_++)
state.m_x[i_]=state.m_xbase[i_];
state.m_f=state.m_fbase;
temp=Func_case_10(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- check
if(temp==-1)
return(false);
//--- check
if(temp==1)
return(true);
//--- Saving state
Func_case_rcomm(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- return result
return(true);
}
temp=Func_case_9(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- check
if(temp==-1)
return(false);
//--- check
if(temp==1)
return(true);
//--- Saving state
Func_case_rcomm(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- return result
return(true);
}
//--- check
if(state.m_rstate.stage==4)
{
//--- change value
state.m_xupdated=false;
//--- check
if(!Func_case_11(state,stepnorm))
return(false);
//--- Now,iteration is finally over
Func_case_7(state,n);
//--- Saving state
Func_case_rcomm(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- return result
return(true);
}
//--- Routine body
//--- Prepare
n=state.m_n;
m=state.m_m;
state.m_repterminationtype=0;
state.m_repiterationscount=0;
state.m_repnfunc=0;
state.m_repnjac=0;
//--- Calculate F/G,initialize algorithm
ClearRequestFields(state);
state.m_needf=true;
state.m_rstate.stage=0;
//--- Saving state
Func_case_rcomm(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- return result
return(true);
}
//+------------------------------------------------------------------+
//| Auxiliary function for NlEqiteration. Is a product to get rid of |
//| the operator unconditional jump goto. |
//+------------------------------------------------------------------+
void CNlEq::Func_case_rcomm(CNlEqState &state,const int n,const int m,
const int i,const bool b,const double lambdaup,
const double lambdadown,const double lambdav,
const double rho,const double mu,const double stepnorm)
{
//--- save
state.m_rstate.ia.Set(0,n);
state.m_rstate.ia.Set(1,m);
state.m_rstate.ia.Set(2,i);
state.m_rstate.ba[0]=b;
state.m_rstate.ra.Set(0,lambdaup);
state.m_rstate.ra.Set(1,lambdadown);
state.m_rstate.ra.Set(2,lambdav);
state.m_rstate.ra.Set(3,rho);
state.m_rstate.ra.Set(4,mu);
state.m_rstate.ra.Set(5,stepnorm);
}
//+------------------------------------------------------------------+
//| Auxiliary function for NlEqiteration. Is a product to get rid of |
//| the operator unconditional jump goto. |
//+------------------------------------------------------------------+
void CNlEq::Func_case_7(CNlEqState &state,const int n)
{
//--- Get Jacobian;
//--- before we get to this point we already have State.XBase filled
//--- with current point and State.FBase filled with function value
//--- at XBase
ClearRequestFields(state);
state.m_needfij=true;
//--- copy
for(int i_=0; i_<n; i_++)
state.m_x[i_]=state.m_xbase[i_];
state.m_rstate.stage=2;
}
//+------------------------------------------------------------------+
//| Auxiliary function for NlEqiteration. Is a product to get rid of |
//| the operator unconditional jump goto. |
//+------------------------------------------------------------------+
bool CNlEq::Func_case_5(CNlEqState &state,double &lambdaup,
double &lambdadown,double &lambdav,
double &rho)
{
//--- check
if(state.m_f<=CMath::Sqr(state.m_epsf))
{
state.m_repterminationtype=1;
return(false);
}
//--- change values
lambdaup=10;
lambdadown=0.3;
lambdav=0.001;
rho=1;
//--- return result
return(true);
}
//+------------------------------------------------------------------+
//| Auxiliary function for NlEqiteration. Is a product to get rid of |
//| the operator unconditional jump goto. |
//+------------------------------------------------------------------+
bool CNlEq::Func_case_11(CNlEqState &state,const double stepnorm)
{
//--- Test stopping conditions on F,step (zero/non-zero) and MaxIts;
//--- If one of the conditions is met,RepTerminationType is changed.
if(MathSqrt(state.m_f)<=state.m_epsf)
state.m_repterminationtype=1;
//--- check
if(stepnorm==0.0 && state.m_repterminationtype==0)
state.m_repterminationtype=-4;
//--- check
if(state.m_repiterationscount>=state.m_maxits && state.m_maxits>0)
state.m_repterminationtype=5;
//--- check
if(state.m_repterminationtype!=0)
return(false);
//--- return result
return(true);
}
//+------------------------------------------------------------------+
//| Auxiliary function for NlEqiteration. Is a product to get rid of |
//| the operator unconditional jump goto. |
//+------------------------------------------------------------------+
int CNlEq::Func_case_10(CNlEqState &state,const int n,const int m,
const int i,const bool b,const double lambdaup,
const double lambdadown,const double lambdav,
const double rho,const double mu,
const double stepnorm)
{
//--- Accept step:
//--- * new position
//--- * new function value
state.m_fbase=state.m_f;
for(int i_=0; i_<n; i_++)
state.m_xbase[i_]=state.m_xbase[i_]+stepnorm*state.m_candstep[i_];
state.m_repiterationscount=state.m_repiterationscount+1;
//--- Report new iteration
if(!state.m_xrep)
{
//--- check
if(!Func_case_11(state,stepnorm))
return(-1);
//--- Now,iteration is finally over
Func_case_7(state,n);
//--- Saving state
Func_case_rcomm(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- return result
return(1);
}
//--- function call
ClearRequestFields(state);
//--- change values
state.m_xupdated=true;
state.m_f=state.m_fbase;
//--- copy
for(int i_=0; i_<n; i_++)
state.m_x[i_]=state.m_xbase[i_];
state.m_rstate.stage=4;
//--- return result
return(0);
}
//+------------------------------------------------------------------+
//| Auxiliary function for NlEqiteration. Is a product to get rid of |
//| the operator unconditional jump goto. |
//+------------------------------------------------------------------+
int CNlEq::Func_case_9(CNlEqState &state,int &n,int &m,int &i,bool &b,
const double lambdaup,const double lambdadown,
double &lambdav,const double rho,const double mu,
double &stepnorm)
{
//--- Solve (J^T*J + (Lambda+Mu)*I)*y=J^T*F
//--- to get step d=-y where:
//--- * Mu=||F|| - is damping parameter for nonlinear system
//--- * Lambda - is additional Levenberg-Marquardt parameter
//--- for better convergence when far away from minimum
for(i=0; i<n; i++)
state.m_candstep[i]=0;
//--- function call
CFbls::FblsSolveCGx(state.m_j,m,n,lambdav,state.m_rightpart,state.m_candstep,state.m_cgbuf);
//--- Normalize step (it must be no more than StpMax)
stepnorm=0;
for(i=0; i<n; i++)
{
//--- check
if(state.m_candstep[i]!=0.0)
{
stepnorm=1;
break;
}
}
CLinMin::LinMinNormalized(state.m_candstep,stepnorm,n);
//--- check
if(state.m_stpmax!=0.0)
stepnorm=MathMin(stepnorm,state.m_stpmax);
//--- Test new step - is it good enough?
//--- * if not,Lambda is increased and we try again.
//--- * if step is good,we decrease Lambda and move on.
//--- We can break this cycle on two occasions:
//--- * step is so small that x+step==x (in floating point arithmetics)
//--- * lambda is so large
for(int i_=0; i_<n; i_++)
state.m_x[i_]=state.m_xbase[i_];
for(int i_=0; i_<n; i_++)
state.m_x[i_]=state.m_x[i_]+stepnorm*state.m_candstep[i_];
b=true;
for(i=0; i<n; i++)
{
//--- check
if(state.m_x[i]!=state.m_xbase[i])
{
b=false;
break;
}
}
//--- check
if(b)
{
//--- Step is too small,force zero step and break
stepnorm=0;
for(int i_=0; i_<n; i_++)
state.m_x[i_]=state.m_xbase[i_];
state.m_f=state.m_fbase;
//--- function call
int temp=Func_case_10(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- check
if(temp!=0)
return(temp);
//--- Saving state
Func_case_rcomm(state,n,m,i,b,lambdaup,lambdadown,lambdav,rho,mu,stepnorm);
//--- return result
return(1);
}
//--- function call
ClearRequestFields(state);
//--- change values
state.m_needf=true;
state.m_rstate.stage=3;
//--- return result
return(0);
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
struct CPolynomialSolverReport
{
double m_maxerr;
//--- constructor / destructor
CPolynomialSolverReport(void) { m_maxerr=0; }
~CPolynomialSolverReport(void) {}
void Copy(const CPolynomialSolverReport &obj) { m_maxerr=obj.m_maxerr; }
//--- overloading
void operator=(const CPolynomialSolverReport &obj) { Copy(obj); }
};
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
class CPolynomialSolver
{
public:
static void PolynomialSolve(CRowDouble &a,int n,CRowComplex &x,CPolynomialSolverReport &rep);
};
//+------------------------------------------------------------------+
//| Polynomial root finding. |
//| This function returns all roots of the polynomial |
//| P(x) = a0 + a1*x + a2*x^2 + ... + an*x^n |
//| Both real and complex roots are returned (see below). |
//| INPUT PARAMETERS: |
//| A - array[N+1], polynomial coefficients: |
//| * A[0] is constant term |
//| * A[N] is a coefficient of X^N |
//| N - polynomial degree |
//| OUTPUT PARAMETERS: |
//| X - array of complex roots: |
//| * for isolated real root, X[I] is strictly real: |
//| IMAGE(X[I])=0 |
//| * complex roots are always returned in pairs-roots |
//| occupy positions I and I+1, with: |
//| * X[I+1]=Conj(X[I]) |
//| * IMAGE(X[I]) > 0 |
//| * IMAGE(X[I+1]) = -IMAGE(X[I]) < 0 |
//| * multiple real roots may have non-zero imaginary |
//| part due to roundoff errors. There is no reliable|
//| way to distinguish real root of multiplicity 2 |
//| from two complex roots in the presence of |
//| roundoff errors. |
//| Rep - report, additional information, following fields |
//| are set: |
//| * Rep.MaxErr - max( |P(xi)| ) for i=0..N-1. This |
//| field allows to quickly estimate "quality" of the|
//| roots being returned. |
//| NOTE: this function uses companion matrix method to find roots. |
//| In case internal EVD solver fails do find eigenvalues, |
//| exception is generated. |
//| NOTE: roots are not "polished" and no matrix balancing is |
//| performed for them. |
//+------------------------------------------------------------------+
void CPolynomialSolver::PolynomialSolve(CRowDouble &A,int n,
CRowComplex &x,
CPolynomialSolverReport &rep)
{
//--- create variables
CMatrixDouble c;
CMatrixDouble vl;
CMatrixDouble vr;
CRowDouble wr;
CRowDouble wi;
int i=0;
int j=0;
bool status=false;
int nz=0;
int ne=0;
complex v=0;
complex vv=0;
CRowDouble a=A;
x.Resize(0);
//--- check
if(!CAp::Assert(n>0,__FUNCTION__+": N<=0"))
return;
if(!CAp::Assert(CAp::Len(a)>n,__FUNCTION__+": Length(A)<N+1"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(a,n+1),__FUNCTION__+": A contains infitite numbers"))
return;
if(!CAp::Assert(a[n]!=0.0,__FUNCTION__+": A[N]=0"))
return;
//--- Prepare
x.Resize(n);
//--- Normalize A:
//--- * analytically determine NZ zero roots
//--- * quick exit for NZ=N
//--- * make residual NE-th degree polynomial monic
//--- (here NE=N-NZ)
nz=0;
while(nz<n && a[nz]==0.0)
nz++;
ne=n-nz;
for(i=nz; i<=n; i++)
a.Set(i-nz,a[i]/a[n]);
//--- For NZ<N, build companion matrix and find NE non-zero roots
if(ne>0)
{
c=matrix<double>::Zeros(ne,ne);
c.Set(0,ne-1,-a[0]);
for(i=1; i<ne; i++)
{
c.Set(i,i-1,1);
c.Set(i,ne-1,-a[i]);
}
status=CEigenVDetect::RMatrixEVD(c,ne,0,wr,wi,vl,vr);
//--- check
if(!CAp::Assert(status,__FUNCTION__+": inernal error - EVD solver failed"))
return;
for(i=0; i<ne; i++)
{
x.SetRe(i,wr[i]);
x.SetIm(i,wi[i]);
}
}
//--- Remaining NZ zero roots
for(i=ne; i<n; i++)
x.Set(i,0.0);
//--- Rep
rep.m_maxerr=0;
for(i=0; i<ne; i++)
{
v=0;
vv=1;
for(j=0; j<=ne; j++)
{
v=v+a[j]*vv;
vv=vv*x[i];
}
rep.m_maxerr=MathMax(rep.m_maxerr,CMath::AbsComplex(v));
}
}
//+------------------------------------------------------------------+
//| This structure is a sparse solver report (both direct and |
//| iterative solvers use this structure). |
//| Following fields can be accessed by users: |
//| * TerminationType (specific error codes depend on the solver |
//| being used, with positive values ALWAYS signaling that |
//| something useful is returned in X, and negative values |
//| ALWAYS meaning critical failures. |
//| * NMV - number of matrix - vector products performed (0 for |
//| direct solvers) |
//| * IterationsCount - inner iterations count (0 for direct |
//| solvers) |
//| * R2 - squared residual |
//+------------------------------------------------------------------+
struct CSparseSolverReport
{
int m_iterationscount;
int m_nmv;
int m_terminationtype;
double m_r2;
//--- constructor / destructor
CSparseSolverReport(void) { ZeroMemory(this); }
~CSparseSolverReport(void) {}
//---
void Copy(const CSparseSolverReport &obj)
{
m_iterationscount=obj.m_iterationscount;
m_nmv=obj.m_nmv;
m_terminationtype=obj.m_terminationtype;
m_r2=obj.m_r2;
}
//--- overloading
void operator=(const CSparseSolverReport &obj) { Copy(obj); }
};
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
class CDirectSparseSolvers
{
public:
static void SparseSPDSolveSKS(CSparseMatrix &a,bool IsUpper,CRowDouble &b,CRowDouble &x,CSparseSolverReport &rep);
static void SparseSPDSolve(CSparseMatrix &a,bool IsUpper,CRowDouble &b,CRowDouble &x,CSparseSolverReport &rep);
static void SparseSPDCholeskySolve(CSparseMatrix &a,bool IsUpper,CRowDouble &b,CRowDouble &x,CSparseSolverReport &rep);
static void SparseSolve(CSparseMatrix &a,CRowDouble &b,CRowDouble &x,CSparseSolverReport &rep);
static void SparseLUSolve(CSparseMatrix &a,CRowInt &p,CRowInt &q,CRowDouble &b,CRowDouble &x,CSparseSolverReport &rep);
static void InitSparseSolverReport(CSparseSolverReport &rep);
};
//+------------------------------------------------------------------+
//| Sparse linear solver for A*x = b with N*N sparse real symmetric |
//| positive definite matrix A, N * 1 vectors x and b. |
//| This solver converts input matrix to SKS format, performs |
//| Cholesky factorization using SKS Cholesky subroutine (works well |
//| for limited bandwidth matrices) and uses sparse triangular |
//| solvers to get solution of the original system. |
//| INPUT PARAMETERS: |
//| A - sparse matrix, must be NxN exactly |
//| IsUpper - which half of A is provided (another half is |
//| ignored) |
//| B - array[0..N - 1], right part |
//| OUTPUT PARAMETERS: |
//| X - array[N], it contains: |
//| * rep.m_terminationtype > 0 => solution |
//| * rep.m_terminationtype = -3 => filled by zeros |
//| Rep - solver report, following fields are set: |
//| * rep.m_terminationtype - solver status; > 0 for |
//| success, set to - 3 on |
//| failure (degenerate or |
//| non-SPD system). |
//+------------------------------------------------------------------+
void CDirectSparseSolvers::SparseSPDSolveSKS(CSparseMatrix &a,
bool IsUpper,
CRowDouble &b,
CRowDouble &x,
CSparseSolverReport &rep)
{
//--- create variables
CSparseMatrix a2 ;
int n=CSparse::SparseGetNRows(a);
x.Resize(0);
//--- check
if(!CAp::Assert(n>0,__FUNCTION__+": N<=0"))
return;
if(!CAp::Assert(CSparse::SparseGetNRows(a)==n,__FUNCTION__+": rows(A)!=N"))
return;
if(!CAp::Assert(CSparse::SparseGetNCols(a)==n,__FUNCTION__+": cols(A)!=N"))
return;
if(!CAp::Assert(CAp::Len(b)>=n,__FUNCTION__+": length(B)<N"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(b,n),__FUNCTION__+": B contains infinities or NANs"))
return;
InitSparseSolverReport(rep);
CSparse::SparseCopyToSKS(a,a2);
if(!CTrFac::SparseCholeskySkyLine(a2,n,IsUpper))
{
rep.m_terminationtype=-3;
x=vector<double>::Zeros(n);
return;
}
x=b;
x.Resize(n);
if(IsUpper)
{
CSparse::SparseTRSV(a2,IsUpper,false,1,x);
CSparse::SparseTRSV(a2,IsUpper,false,0,x);
}
else
{
CSparse::SparseTRSV(a2,IsUpper,false,0,x);
CSparse::SparseTRSV(a2,IsUpper,false,1,x);
}
rep.m_terminationtype=1;
}
//+------------------------------------------------------------------+
//| Sparse linear solver for A*x = b with N*N sparse real symmetric |
//| positive definite matrix A, N * 1 vectors x and b. |
//| This solver converts input matrix to CRS format, performs |
//| Cholesky factorization using supernodal Cholesky decomposition |
//| with permutation-reducing ordering and uses sparse triangular |
//| solver to get solution of the original system. |
//| INPUT PARAMETERS: |
//| A - sparse matrix, must be NxN exactly |
//| IsUpper - which half of A is provided (another half is |
//| ignored) |
//| B - array[N], right part |
//| OUTPUT PARAMETERS: |
//| X - array[N], it contains: |
//| * rep.m_terminationtype > 0 => solution |
//| * rep.m_terminationtype = -3 => filled by zeros |
//| Rep - solver report, following fields are set: |
//| * rep.m_terminationtype - solver status; > 0 for |
//| success, set to - 3 on |
//| failure (degenerate or |
//| non-SPD system). |
//+------------------------------------------------------------------+
void CDirectSparseSolvers::SparseSPDSolve(CSparseMatrix &a,
bool IsUpper,
CRowDouble &b,
CRowDouble &x,
CSparseSolverReport &rep)
{
//--- create variables
int n=CSparse::SparseGetNRows(a);
double v=0;
CSparseMatrix a2;
CRowInt p;
x.Resize(0);
//--- check
if(!CAp::Assert(n>0,__FUNCTION__+": N<=0"))
return;
if(!CAp::Assert(CSparse::SparseGetNRows(a)==n,__FUNCTION__+": rows(A)!=N"))
return;
if(!CAp::Assert(CSparse::SparseGetNCols(a)==n,__FUNCTION__+": cols(A)!=N"))
return;
if(!CAp::Assert(CAp::Len(b)>=n,__FUNCTION__+": length(B)<N"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(b,n),__FUNCTION__+": B contains infinities or NANs"))
return;
InitSparseSolverReport(rep);
CSparse::SparseCopyToCRS(a,a2);
if(!CTrFac::SparseCholeskyP(a2,IsUpper,p))
{
rep.m_terminationtype=-3;
x=vector<double>::Zeros(n);
return;
}
CAblasF::RCopyAllocV(n,b,x);
for(int i=0; i<n; i++)
x.Swap(i,p[i]);
if(IsUpper)
{
CSparse::SparseTRSV(a2,IsUpper,false,1,x);
CSparse::SparseTRSV(a2,IsUpper,false,0,x);
}
else
{
CSparse::SparseTRSV(a2,IsUpper,false,0,x);
CSparse::SparseTRSV(a2,IsUpper,false,1,x);
}
for(int i=n-1; i>=0; i--)
x.Swap(i,p[i]);
rep.m_terminationtype=1;
}
//+------------------------------------------------------------------+
//| Sparse linear solver for A*x = b with N*N real symmetric positive|
//| definite matrix A given by its Cholesky decomposition, and N * 1 |
//| vectors x and b. |
//| IMPORTANT: this solver requires input matrix to be in the SKS |
//| (Skyline) or CRS (compressed row storage) format. An |
//| exception will be generated if you pass matrix in some|
//| other format. |
//| INPUT PARAMETERS: |
//| A - sparse NxN matrix stored in CRs or SKS format, must|
//| be NxN exactly |
//| IsUpper - which half of A is provided (another half is |
//| ignored) |
//| B - array[N], right part |
//| OUTPUT PARAMETERS: |
//| X - array[N], it contains: |
//| * rep.m_terminationtype > 0 => solution |
//| * rep.m_terminationtype = -3 => filled by zeros |
//| Rep - solver report, following fields are set: |
//| * rep.m_terminationtype - solver status; > 0 for |
//| success, set to - 3 on |
//| failure (degenerate or |
//| non-SPD system). |
//+------------------------------------------------------------------+
void CDirectSparseSolvers::SparseSPDCholeskySolve(CSparseMatrix &a,
bool IsUpper,
CRowDouble &b,
CRowDouble &x,
CSparseSolverReport &rep)
{
int n=CSparse::SparseGetNRows(a);
x.Resize(0);
//--- check
if(!CAp::Assert(n>0,__FUNCTION__+": N<=0"))
return;
if(!CAp::Assert(CSparse::SparseGetNRows(a)==n,__FUNCTION__+": rows(A)!=N"))
return;
if(!CAp::Assert(CSparse::SparseGetNCols(a)==n,__FUNCTION__+": cols(A)!=N"))
return;
if(!CAp::Assert(CSparse::SparseIsSKS(a) || CSparse::SparseIsCRS(a),__FUNCTION__+": A is not an SKS/CRS matrix"))
return;
if(!CAp::Assert(CAp::Len(b)>=n,__FUNCTION__+": length(B)<N"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(b,n),__FUNCTION__+": B contains infinities or NANs"))
return;
InitSparseSolverReport(rep);
for(int i=0; i<n; i++)
{
if(CSparse::SparseGet(a,i,i)==0.0)
{
rep.m_terminationtype=-3;
x=vector<double>::Zeros(n);
return;
}
}
x=b;
x.Resize(n);
if(IsUpper)
{
CSparse::SparseTRSV(a,IsUpper,false,1,x);
CSparse::SparseTRSV(a,IsUpper,false,0,x);
}
else
{
CSparse::SparseTRSV(a,IsUpper,false,0,x);
CSparse::SparseTRSV(a,IsUpper,false,1,x);
}
rep.m_terminationtype=1;
}
//+------------------------------------------------------------------+
//| Sparse linear solver for A*x = b with general (nonsymmetric) N*N |
//| sparse real matrix A, N * 1 vectors x and b. |
//| This solver converts input matrix to CRS format, performs LU |
//| factorization and uses sparse triangular solvers to get solution |
//| of the original system. |
//| INPUT PARAMETERS: |
//| A - sparse matrix, must be NxN exactly, any storage |
//| format |
//| N - size of A, N > 0 |
//| B - array[0..N - 1], right part |
//| OUTPUT PARAMETERS: |
//| X - array[N], it contains: |
//| * rep.m_terminationtype > 0 => solution |
//| * rep.m_terminationtype = -3 => filled by zeros |
//| Rep - solver report, following fields are set: |
//| * rep.m_terminationtype - solver status; > 0 for |
//| success, set to - 3 on failure (degenerate |
//| system). |
//+------------------------------------------------------------------+
void CDirectSparseSolvers::SparseSolve(CSparseMatrix &a,
CRowDouble &b,
CRowDouble &x,
CSparseSolverReport &rep)
{
//--- create variables
int n=CSparse::SparseGetNRows(a);
double v=0;
CSparseMatrix a2;
CRowInt pivp;
CRowInt pivq;
x.Resize(0);
//--- check
if(!CAp::Assert(n>0,__FUNCTION__+": N<=0"))
return;
if(!CAp::Assert(CSparse::SparseGetNRows(a)==n,__FUNCTION__+": rows(A)!=N"))
return;
if(!CAp::Assert(CSparse::SparseGetNCols(a)==n,__FUNCTION__+": cols(A)!=N"))
return;
if(!CAp::Assert(CAp::Len(b)>=n,__FUNCTION__+": length(B)<N"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(b,n),__FUNCTION__+": B contains infinities or NANs"))
return;
InitSparseSolverReport(rep);
CSparse::SparseCopyToCRS(a,a2);
if(!CTrFac::SparseLU(a2,0,pivp,pivq))
{
rep.m_terminationtype=-3;
x=vector<double>::Zeros(n);
return;
}
x=b;
x.Resize(n);
for(int i=0; i<n; i++)
x.Swap(i,pivp[i]);
CSparse::SparseTRSV(a2,false,true,0,x);
CSparse::SparseTRSV(a2,true,false,0,x);
for(int i=n-1; i>=0; i--)
x.Swap(i,pivq[i]);
rep.m_terminationtype=1;
}
//+------------------------------------------------------------------+
//| Sparse linear solver for A*x = b with general (nonsymmetric) N*N |
//| sparse real matrix A given by its LU factorization, N*1 vectors x|
//| and b. |
//| IMPORTANT: this solver requires input matrix to be in the CRS |
//| sparse storage format. An exception will be generated |
//| if you pass matrix in some other format (HASH or SKS).|
//| INPUT PARAMETERS: |
//| A - LU factorization of the sparse matrix, must be NxN |
//| exactly in CRS storage format |
//| P, Q - pivot indexes from LU factorization |
//| N - size of A, N > 0 |
//| B - array[0..N - 1], right part |
//| OUTPUT PARAMETERS: |
//| X - array[N], it contains: |
//| * rep.m_terminationtype > 0 => solution |
//| * rep.m_terminationtype = -3 => filled by zeros |
//| Rep - solver report, following fields are set: |
//| * rep.m_terminationtype - solver status; > 0 for |
//| success, set to - 3 on |
//| failure (degenerate |
//| system). |
//+------------------------------------------------------------------+
void CDirectSparseSolvers::SparseLUSolve(CSparseMatrix &a,
CRowInt &p,
CRowInt &q,
CRowDouble &b,
CRowDouble &x,
CSparseSolverReport &rep)
{
//--- create variables
int n=CSparse::SparseGetNRows(a);
int i=0;
int j=0;
double v=0;
x.Resize(0);
//--- check
if(!CAp::Assert(n>0,__FUNCTION__+": N<=0"))
return;
if(!CAp::Assert(CSparse::SparseGetNRows(a)==n,__FUNCTION__+": rows(A)!=N"))
return;
if(!CAp::Assert(CSparse::SparseGetNCols(a)==n,__FUNCTION__+": cols(A)!=N"))
return;
if(!CAp::Assert(CSparse::SparseIsCRS(a),__FUNCTION__+": A is not an SKS matrix"))
return;
if(!CAp::Assert(CAp::Len(b)>=n,__FUNCTION__+": length(B)<N"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(b,n),__FUNCTION__+": B contains infinities or NANs"))
return;
if(!CAp::Assert(CAp::Len(p)>=n,__FUNCTION__+": length(P)<N"))
return;
if(!CAp::Assert(CAp::Len(q)>=n,__FUNCTION__+": length(Q)<N"))
return;
for(i=0; i<n; i++)
{
//--- check
if(!CAp::Assert(p[i]>=i && p[i]<n,__FUNCTION__+": P is corrupted"))
return;
if(!CAp::Assert(q[i]>=i && q[i]<n,__FUNCTION__+": Q is corrupted"))
return;
}
InitSparseSolverReport(rep);
x.Resize(n);
for(i=0; i<n; i++)
{
if(a.m_DIdx[i]==a.m_UIdx[i] || a.m_Vals[a.m_DIdx[i]]==0.0)
{
rep.m_terminationtype=-3;
x=vector<double>::Zeros(n);
return;
}
}
x=b;
x.Resize(n);
for(i=0; i<n; i++)
x.Swap(i,p[i]);
CSparse::SparseTRSV(a,false,true,0,x);
CSparse::SparseTRSV(a,true,false,0,x);
for(i=n-1; i>=0; i--)
x.Swap(i,q[i]);
rep.m_terminationtype=1;
}
//+------------------------------------------------------------------+
//| Reset report fields |
//+------------------------------------------------------------------+
void CDirectSparseSolvers::InitSparseSolverReport(CSparseSolverReport &rep)
{
rep.m_terminationtype=0;
rep.m_nmv=0;
rep.m_iterationscount=0;
rep.m_r2=0;
}
//+------------------------------------------------------------------+
//| This object stores state of the sparse linear solver object. |
//| You should use ALGLIB functions to work with this object. |
//| Never try to access its fields directly! |
//+------------------------------------------------------------------+
struct CSparseSolverState
{
int m_algotype;
int m_gmresk;
int m_maxits;
int m_n;
int m_repiterationscount;
int m_repnmv;
int m_repterminationtype;
int m_requesttype;
double m_epsf;
double m_reply1;
double m_repr2;
bool m_running;
bool m_userterminationneeded;
bool m_xrep;
RCommState m_rstate;
CSparseMatrix m_convbuf;
CRowDouble m_ax;
CRowDouble m_b;
CRowDouble m_wrkb;
CRowDouble m_x0;
CRowDouble m_x;
CRowDouble m_xf;
CFblsGMRESState m_gmressolver;
//--- constructor / destructor
CSparseSolverState(void);
~CSparseSolverState(void) {}
//---
void Copy(const CSparseSolverState &obj);
//--- overloading
void operator=(const CSparseSolverState &obj) { Copy(obj); }
};
//+------------------------------------------------------------------+
//| Constructor |
//+------------------------------------------------------------------+
CSparseSolverState::CSparseSolverState(void)
{
m_algotype=0;
m_gmresk=0;
m_maxits=0;
m_n=0;
m_repiterationscount=0;
m_repnmv=0;
m_repterminationtype=0;
m_requesttype=0;
m_epsf=0;
m_reply1=0;
m_repr2=0;
m_running=false;
m_userterminationneeded=false;
m_xrep=false;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
void CSparseSolverState::Copy(const CSparseSolverState &obj)
{
m_rstate=obj.m_rstate;
m_algotype=obj.m_algotype;
m_gmresk=obj.m_gmresk;
m_maxits=obj.m_maxits;
m_n=obj.m_n;
m_repiterationscount=obj.m_repiterationscount;
m_repnmv=obj.m_repnmv;
m_repterminationtype=obj.m_repterminationtype;
m_requesttype=obj.m_requesttype;
m_epsf=obj.m_epsf;
m_reply1=obj.m_reply1;
m_repr2=obj.m_repr2;
m_running=obj.m_running;
m_userterminationneeded=obj.m_userterminationneeded;
m_xrep=obj.m_xrep;
m_convbuf=obj.m_convbuf;
m_ax=obj.m_ax;
m_b=obj.m_b;
m_wrkb=obj.m_wrkb;
m_x0=obj.m_x0;
m_x=obj.m_x;
m_xf=obj.m_xf;
m_gmressolver=obj.m_gmressolver;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
class CIterativeSparse
{
public:
static void SparseSolveSymmetricGMRES(CSparseMatrix &a,bool IsUpper,CRowDouble &b,int k,double epsf,int maxits,CRowDouble &x,CSparseSolverReport &rep);
static void SparseSolveGMRES(CSparseMatrix &a,CRowDouble &b,int k,double epsf,int maxits,CRowDouble &x,CSparseSolverReport &rep);
static void SparseSolverCreate(int n,CSparseSolverState &state);
static void SparseSolverSetAlgoGMRES(CSparseSolverState &state,int k);
static void SparseSolverSetStartingPoint(CSparseSolverState &state,CRowDouble &x);
static void SparseSolverSetCond(CSparseSolverState &state,double epsf,int maxits);
static void SparseSolverSolveSymmetric(CSparseSolverState &state,CSparseMatrix &a,bool IsUpper,CRowDouble &b);
static void SparseSolverSolve(CSparseSolverState &state,CSparseMatrix &a,CRowDouble &b);
static void SparseSolverResults(CSparseSolverState &state,CRowDouble &x,CSparseSolverReport &rep);
static void SparseSolverSetXRep(CSparseSolverState &state,bool needxrep);
static void SparseSolverOOCStart(CSparseSolverState &state,CRowDouble &b);
static bool SparseSolverOOCContinue(CSparseSolverState &state);
static void SparseSolverOOCGetRequestInfo(CSparseSolverState &state,int &requesttype);
static void SparseSolverOOCGetRequestData(CSparseSolverState &state,CRowDouble &x);
static void SparseSolverOOCGetRequestData1(CSparseSolverState &state,double &v);
static void SparseSolverOOCSendResult(CSparseSolverState &state,CRowDouble &ax);
static void SparseSolverOOCStop(CSparseSolverState &state,CRowDouble &x,CSparseSolverReport &rep);
static void SparseSolverRequestTermination(CSparseSolverState &state);
static bool SparseSolverIteration(CSparseSolverState &state);
private:
static void ClearRequestFields(CSparseSolverState &state);
static void ClearReportFields(CSparseSolverState &state);
};
//+------------------------------------------------------------------+
//| Solving sparse symmetric linear system A*x = b using GMRES(k) |
//| method. Sparse symmetric A is given by its lower or upper |
//| triangle. |
//| NOTE: use SparseSolveGMRES() to solve system with nonsymmetric A.|
//| This function provides convenience API for an 'expert' interface |
//| provided by SparseSolverState class. Use SparseSolver API if you |
//| need advanced functions like providing initial point, using |
//| out-of-core API and so on. |
//| INPUT PARAMETERS: |
//| A - sparse symmetric NxN matrix in any sparse storage |
//| format. Using CRS format is recommended because it |
//| avoids internal conversion. An exception will be |
//| generated if A is not NxN matrix (where N is a size|
//| specified during solver object creation). |
//| IsUpper - whether upper or lower triangle of A is used: |
//| * IsUpper = True => only upper triangle is used and|
//| lower triangle is not referenced at all |
//| * IsUpper = False => only lower triangle is used |
//| and upper triangle is not referenced at all |
//| B - right part, array[N] |
//| K - k parameter for GMRES(k), k >= 0. Zero value means |
//| that algorithm will choose it automatically. |
//| EpsF - stopping condition, EpsF >= 0. The algorithm will |
//| stop when residual will decrease below EpsF* | B |.|
//| Having EpsF = 0 means that this stopping condition |
//| is ignored. |
//| MaxIts - stopping condition, MaxIts >= 0. The algorithm will|
//| stop after performing MaxIts iterations. Zero value|
//| means no limit. |
//| NOTE: having both EpsF = 0 and MaxIts = 0 means that stopping |
//| criteria will be chosen automatically. |
//| OUTPUT PARAMETERS: |
//| X - array[N], the solution |
//| Rep - solution report: |
//| * Rep.TerminationType completion code: |
//| * -5 CG method was used for a matrix which is |
//| not positive definite |
//| * -4 overflow / underflow during solution (ill |
//| conditioned problem) |
//| * 1 || residual || <= EpsF* || b || |
//| * 5 MaxIts steps was taken |
//| * 7 rounding errors prevent further progress, |
//| best point found is returned |
//| * 8 the algorithm was terminated early with |
//| SparseSolverRequestTermination() being |
//| called from other thread. |
//| * Rep.IterationsCount contains iterations count |
//| * Rep.NMV contains number of matrix - vector |
//| calculations |
//| * Rep.R2 contains squared residual |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolveSymmetricGMRES(CSparseMatrix &a,
bool IsUpper,
CRowDouble &b,
int k,
double epsf,
int maxits,
CRowDouble &x,
CSparseSolverReport &rep)
{
//--- create variables
int n=CSparse::SparseGetNRows(a);
CSparseMatrix convbuf;
CSparseSolverState solver;
x.Resize(0);
//--- Test inputs
if(!CAp::Assert(n>=1,__FUNCTION__+": tried to automatically detect N from sizeof(A),got nonpositive size"))
return;
if(!CAp::Assert(CSparse::SparseGetNRows(a)==n,__FUNCTION__+": rows(A)!=N"))
return;
if(!CAp::Assert(CSparse::SparseGetNCols(a)==n,__FUNCTION__+": cols(A)!=N"))
return;
if(!CAp::Assert(CAp::Len(b)>=n,__FUNCTION__+": length(B)<N"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(b,n),__FUNCTION__+": B contains NAN/INF"))
return;
if(!CAp::Assert(MathIsValidNumber(epsf) && epsf>=0.0,__FUNCTION__+": EpsF<0 or infinite"))
return;
if(!CAp::Assert(maxits>=0,__FUNCTION__+": MaxIts<0"))
return;
if(epsf==0.0 && maxits==0)
epsf=1.0E-6;
//--- If A is non-CRS, perform conversion
if(!CSparse::SparseIsCRS(a))
{
CSparse::SparseCopyToCRSBuf(a,convbuf);
SparseSolveSymmetricGMRES(convbuf,IsUpper,b,k,epsf,maxits,x,rep);
return;
}
//--- Solve using temporary solver object
SparseSolverCreate(n,solver);
SparseSolverSetAlgoGMRES(solver,k);
SparseSolverSetCond(solver,epsf,maxits);
SparseSolverSolveSymmetric(solver,a,IsUpper,b);
SparseSolverResults(solver,x,rep);
}
//+------------------------------------------------------------------+
//| Solving sparse linear system A*x = b using GMRES(k) method. |
//| This function provides convenience API for an 'expert' interface |
//| provided by SparseSolverState class. Use SparseSolver API if you |
//| need advanced functions like providing initial point, using |
//| out-of-core API and so on. |
//| INPUT PARAMETERS: |
//| A - sparse NxN matrix in any sparse storage format. |
//| Using CRS format is recommended because it avoids |
//| internal conversion. An exception will be generated|
//| if A is not NxN matrix (where N is a size specified|
//| during solver object creation). |
//| B - right part, array[N] |
//| K - k parameter for GMRES(k), k >= 0. Zero value means |
//| that algorithm will choose it automatically. |
//| EpsF - stopping condition, EpsF >= 0. The algorithm will |
//| stop when residual will decrease below EpsF*| B |. |
//| Having EpsF = 0 means that this stopping condition |
//| is ignored. |
//| MaxIts - stopping condition, MaxIts >= 0. The algorithm will|
//| stop after performing MaxIts iterations. Zero value|
//| means no limit. |
//| NOTE: having both EpsF = 0 and MaxIts = 0 means that stopping |
//| criteria will be chosen automatically. |
//| OUTPUT PARAMETERS: |
//| X - array[N], the solution |
//| Rep - solution report: |
//| * Rep.TerminationType completion code: |
//| * -5 CG method was used for a matrix which is |
//| not positive definite |
//| * -4 overflow / underflow during solution (ill |
//| conditioned problem) |
//| * 1 || residual || <= EpsF* || b || |
//| * 5 MaxIts steps was taken |
//| * 7 rounding errors prevent further progress, |
//| best point found is returned |
//| * 8 the algorithm was terminated early with |
//| SparseSolverRequestTermination() being |
//| called from other thread. |
//| * Rep.IterationsCount contains iterations count |
//| * Rep.NMV contains number of matrix - vector |
//| calculations |
//| * Rep.R2 contains squared residual |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolveGMRES(CSparseMatrix &a,
CRowDouble &b,
int k,
double epsf,
int maxits,
CRowDouble &x,
CSparseSolverReport &rep)
{
//--- create variables
int n=CSparse::SparseGetNRows(a);
CSparseMatrix convbuf;
CSparseSolverState solver;
x.Resize(0);
//--- Test inputs
if(!CAp::Assert(n>=1,__FUNCTION__+": tried to automatically detect N from sizeof(A),got nonpositive size"))
return;
if(!CAp::Assert(CSparse::SparseGetNRows(a)==n,__FUNCTION__+": rows(A)!=N"))
return;
if(!CAp::Assert(CSparse::SparseGetNCols(a)==n,__FUNCTION__+": cols(A)!=N"))
return;
if(!CAp::Assert(CAp::Len(b)>=n,__FUNCTION__+": length(B)<N"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(b,n),__FUNCTION__+": B contains NAN/INF"))
return;
if(!CAp::Assert(MathIsValidNumber(epsf) && epsf>=0.0,__FUNCTION__+": EpsF<0 or infinite"))
return;
if(!CAp::Assert(maxits>=0,__FUNCTION__+": MaxIts<0"))
return;
if(epsf==0.0 && maxits==0)
epsf=1.0E-6;
//--- If A is non-CRS, perform conversion
if(!CSparse::SparseIsCRS(a))
{
CSparse::SparseCopyToCRSBuf(a,convbuf);
SparseSolveGMRES(convbuf,b,k,epsf,maxits,x,rep);
return;
}
//--- Solve using temporary solver object
SparseSolverCreate(n,solver);
SparseSolverSetAlgoGMRES(solver,k);
SparseSolverSetCond(solver,epsf,maxits);
SparseSolverSolve(solver,a,b);
SparseSolverResults(solver,x,rep);
}
//+------------------------------------------------------------------+
//| This function initializes sparse linear iterative solver object. |
//| This solver can be used to solve nonsymmetric and symmetric |
//| positive definite NxN(square) linear systems. |
//| The solver provides 'expert' API which allows advanced control |
//| over algorithms being used, including ability to get progress |
//| report, terminate long-running solver from other thread, |
//| out-of-core solution and so on. |
//| NOTE: there are also convenience functions that allows quick one-|
//| line to the solvers: |
//| * SparseSolveCG() to solve SPD linear systems |
//| * SparseSolveGMRES() to solve unsymmetric linear systems. |
//| NOTE: if you want to solve MxN(rectangular) linear problem you |
//| may use LinLSQR solver provided by ALGLIB. |
//| USAGE (A is given by the SparseMatrix structure): |
//| 1. User initializes algorithm state with SparseSolverCreate() |
//| call |
//| 2. User selects algorithm with one of the SparseSolverSetAlgo??|
//| functions. By default, GMRES(k) is used with automatically |
//| chosen k |
//| 3. Optionally, user tunes solver parameters, sets starting |
//| point, etc. |
//| 4. Depending on whether system is symmetric or not, user calls:|
//| * SparseSolverSolveSymmetric() for a symmetric system given |
//| by its lower or upper triangle |
//| * SparseSolverSolve() for a nonsymmetric system or a |
//| symmetric one given by the full matrix |
//| 5. User calls SparseSolverResults() to get the solution |
//| It is possible to call SparseSolverSolve???() again to solve |
//| another task with same dimensionality but different matrix and/or|
//| right part without reinitializing SparseSolverState structure. |
//| USAGE(out-of-core mode): |
//| 1. User initializes algorithm state with SparseSolverCreate() |
//| call |
//| 2. User selects algorithm with one of the SparseSolverSetAlgo??|
//| functions. By default, GMRES(k) is used with automatically |
//| chosen k |
//| 3. Optionally, user tunes solver parameters, sets starting |
//| point, etc. |
//| 4. After that user should work with out-of-core interface in a |
//| loop like one given below: |
//| > CAlgLib::SparseSolverOOCStart(state) |
//| > while CAlgLib::SparseSolverOOCContinue(state) do |
//| > CAlgLib::SparseSolver))CGetRequestInfo(state, |
//| RequestType) |
//| > CAlgLib::SparseSolverOOCGetRequeStdata(state, X) |
//| > if RequestType = 0 then |
//| > [calculate Y = A * X, with X = R ^ N] |
//| > CAlgLib:SparseSolverOOCSendResult(state, Y) |
//| > CAlgLib::SparseSolverOOCStop(state, X, Report) |
//| INPUT PARAMETERS: |
//| N - problem dimensionality (fixed at start-up) |
//| OUTPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverCreate(int n,
CSparseSolverState &state)
{
//--- check
if(!CAp::Assert(n>=1,__FUNCTION__+": N<=0"))
return;
state.m_n=n;
state.m_running=false;
state.m_userterminationneeded=false;
CAblasF::RSetAllocV(state.m_n,0.0,state.m_x0);
CAblasF::RSetAllocV(state.m_n,0.0,state.m_x);
CAblasF::RSetAllocV(state.m_n,0.0,state.m_ax);
CAblasF::RSetAllocV(state.m_n,0.0,state.m_xf);
CAblasF::RSetAllocV(state.m_n,0.0,state.m_b);
CAblasF::RSetAllocV(state.m_n,0.0,state.m_wrkb);
state.m_reply1=0.0;
SparseSolverSetXRep(state,false);
SparseSolverSetCond(state,0.0,0);
SparseSolverSetAlgoGMRES(state,0);
ClearRequestFields(state);
ClearReportFields(state);
}
//+------------------------------------------------------------------+
//| This function sets the solver algorithm to GMRES(k). |
//| NOTE: if you do not need advanced functionality of the |
//| SparseSolver API, you may use convenience functions |
//| SparseSolveGMRES() and SparseSolveSymmetricGMRES(). |
//| INPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//| K - GMRES parameter, K >= 0: |
//| * recommended values are in 10..100 range |
//| * larger values up to N are possible but have |
//| little sense - the algorithm will be slower than |
//| any dense solver. |
//| * values above N are truncated down to N |
//| * zero value means that default value is chosen. |
//| This value is 50 in the current version, but it |
//| may change in future ALGLIB releases. |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverSetAlgoGMRES(CSparseSolverState &state,
int k)
{
//--- check
if(!CAp::Assert(k>=0,__FUNCTION__+": K<0"))
return;
state.m_algotype=0;
if(k==0)
k=50;
state.m_gmresk=MathMin(k,state.m_n);
}
//+------------------------------------------------------------------+
//| This function sets starting point. |
//| By default, zero starting point is used. |
//| INPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//| X - starting point, array[N] |
//| OUTPUT PARAMETERS: |
//| State - new starting point was set |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverSetStartingPoint(CSparseSolverState &state,
CRowDouble &x)
{
//--- check
if(!CAp::Assert(state.m_n<=CAp::Len(x),__FUNCTION__+": Length(X)<N"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(x,state.m_n),__FUNCTION__+": X contains infinite or NaN values!"))
return;
state.m_x0=x;
state.m_x0.Resize(state.m_n);
}
//+------------------------------------------------------------------+
//| This function sets stopping criteria. |
//| INPUT PARAMETERS: |
//| EpsF - algorithm will be stopped if norm of residual is |
//| less than EpsF* || b ||. |
//| MaxIts - algorithm will be stopped if number of iterations |
//| is more than MaxIts. |
//| OUTPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//| NOTES: If both EpsF and MaxIts are zero then small EpsF will be |
//| set to small value. |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverSetCond(CSparseSolverState &state,
double epsf,
int maxits)
{
//--- check
if(!CAp::Assert(MathIsValidNumber(epsf) && epsf>=0.0,__FUNCTION__+": EpsF is negative or contains infinite or NaN values"))
return;
if(!CAp::Assert(maxits>=0,__FUNCTION__+": MaxIts is negative"))
return;
if(epsf==0.0 && maxits==0)
{
state.m_epsf=1.0E-6;
state.m_maxits=0;
}
else
{
state.m_epsf=epsf;
state.m_maxits=maxits;
}
}
//+------------------------------------------------------------------+
//| Procedure for the solution of A*x = b with sparse symmetric A |
//| given by its lower or upper triangle. |
//| This function will work with any solver algorithm being used, |
//| SPD one (like CG) or not(like GMRES). Using unsymmetric solvers |
//| (like GMRES) on SPD problems is suboptimal, but still possible. |
//| NOTE: the solver behavior is ill-defined for a situation when a |
//| SPD solver is used on indefinite matrix. It may solve the |
//| problem up to desired precision (sometimes, rarely) or |
//| return with error code signalling violation of underlying |
//| assumptions. |
//| INPUT PARAMETERS: |
//| State - algorithm state |
//| A - sparse symmetric NxN matrix in any sparse storage |
//| format. Using CRS format is recommended because it |
//| avoids internal conversion. An exception will be |
//| generated if A is not NxN matrix (where N is a size|
//| specified during solver object creation). |
//| IsUpper - whether upper or lower triangle of A is used: |
//| * IsUpper = True => only upper triangle is used and|
//| lower triangle is not referenced at all |
//| * IsUpper = False => only lower triangle is used |
//| and upper triangle is not referenced at all |
//| B - right part, array[N] |
//| RESULT: |
//| This function returns no result. You can get the solution by |
//| calling SparseSolverResults() |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverSolveSymmetric(CSparseSolverState &state,
CSparseMatrix &a,
bool IsUpper,
CRowDouble &b)
{
int n=state.m_n;
//--- Test inputs
if(!CAp::Assert(CSparse::SparseGetNRows(a)==n,__FUNCTION__+": rows(A)!=N"))
return;
if(!CAp::Assert(CSparse::SparseGetNCols(a)==n,__FUNCTION__+": cols(A)!=N"))
return;
if(!CAp::Assert(CAp::Len(b)>=n,__FUNCTION__+": length(B)<N"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(b,n),__FUNCTION__+": B contains NAN/INF"))
return;
//--- If A is non-CRS, perform conversion
if(!CSparse::SparseIsCRS(a))
{
CSparse::SparseCopyToCRSBuf(a,state.m_convbuf);
SparseSolverSolveSymmetric(state,state.m_convbuf,IsUpper,b);
return;
}
//--- Solve using out-of-core API
SparseSolverOOCStart(state,b);
while(SparseSolverOOCContinue(state))
{
if(state.m_requesttype==-1)
{
//--- Skip location reports
continue;
}
//--- check
if(!CAp::Assert(state.m_requesttype==0,__FUNCTION__+": integrity check 7372 failed"))
return;
CSparse::SparseSMV(a,IsUpper,state.m_x,state.m_ax);
}
}
//+------------------------------------------------------------------+
//| Procedure for the solution of A*x = b with sparse nonsymmetric A |
//| IMPORTANT: this function will work with any solver algorithm |
//| being used, symmetric solver like CG, or not. However,|
//| using symmetric solvers on nonsymmetric problems is |
//| dangerous. It may solve the problem up to desired |
//| precision (sometimes, rarely) or terminate with error |
//| code signalling violation of underlying assumptions. |
//| INPUT PARAMETERS: |
//| State - algorithm state |
//| A - sparse NxN matrix in any sparse storage format. |
//| Using CRS format is recommended because it avoids |
//| internal conversion. An exception will be generated|
//| if A is not NxN matrix (where N is a size specified|
//| during solver object creation). |
//| B - right part, array[N] |
//| RESULT: |
//| This function returns no result. |
//| You can get the solution by calling SparseSolverResults() |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverSolve(CSparseSolverState &state,
CSparseMatrix &a,
CRowDouble &b)
{
int n=state.m_n;
//--- Test inputs
if(!CAp::Assert(CSparse::SparseGetNRows(a)==n,__FUNCTION__+": rows(A)!=N"))
return;
if(!CAp::Assert(CSparse::SparseGetNCols(a)==n,__FUNCTION__+": cols(A)!=N"))
return;
if(!CAp::Assert(CAp::Len(b)>=n,__FUNCTION__+": length(B)<N"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(b,n),__FUNCTION__+": B contains NAN/INF"))
return;
//--- If A is non-CRS, perform conversion
if(!CSparse::SparseIsCRS(a))
{
CSparse::SparseCopyToCRSBuf(a,state.m_convbuf);
SparseSolverSolve(state,state.m_convbuf,b);
return;
}
//--- Solve using out-of-core API
SparseSolverOOCStart(state,b);
while(SparseSolverOOCContinue(state))
{
if(state.m_requesttype==-1)
{
//--- Skip location reports
continue;
}
if(!CAp::Assert(state.m_requesttype==0,__FUNCTION__+": integrity check 7372 failed"))
return;
CSparse::SparseMV(a,state.m_x,state.m_ax);
}
}
//+------------------------------------------------------------------+
//| Sparse solver results. |
//| This function must be called after calling one of the |
//| SparseSolverSolve() functions. |
//| INPUT PARAMETERS: |
//| State - algorithm state |
//| OUTPUT PARAMETERS: |
//| X - array[N], solution |
//| Rep - solution report: |
//| * Rep.TerminationType completion code: |
//| * -5 CG method was used for a matrix which is |
//| not positive definite |
//| * -4 overflow/underflow during solution (ill |
//| conditioned problem) |
//| * 1 || residual || <= EpsF* || b || |
//| * 5 MaxIts steps was taken |
//| * 7 rounding errors prevent further progress, |
//| best point found is returned |
//| * 8 the algorithm was terminated early with |
//| SparseSolverRequestTermination() being |
//| called from other thread. |
//| * Rep.IterationsCount contains iterations count |
//| * Rep.NMV contains number of matrix - vector |
//| calculations |
//| * Rep.R2 contains squared residual |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverResults(CSparseSolverState &state,
CRowDouble &x,
CSparseSolverReport &rep)
{
x.Resize(0);
SparseSolverOOCStop(state,x,rep);
}
//+------------------------------------------------------------------+
//|This function turns on/off reporting during out-of-core processing|
//| When the solver works in the out-of-core mode, it can be |
//| configured to report its progress by returning current location. |
//| These location reports implemented as a special kind of the |
//| out-of-core request: |
//| * SparseSolverOOCGetRequestInfo() returns - 1 |
//| * SparseSolverOOCGetRequestData() returns current location |
//| * SparseSolverOOCGetRequestData1() returns squared norm of |
//| the residual |
//| * SparseSolverOOCSendResult() shall NOT be called |
//| This function has no effect when SparseSolverSolve() is used |
//| because this function has no method of reporting its progress. |
//| NOTE: when used with GMRES(k), this function reports progress |
//| every k-th iteration. |
//| INPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//| NeedXRep - whether iteration reports are needed or not |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverSetXRep(CSparseSolverState &state,
bool needxrep)
{
state.m_xrep=needxrep;
}
//+------------------------------------------------------------------+
//| This function initiates out-of-core mode of the sparse solver. It|
//| should be used in conjunction with other out-of-core - related |
//| functions of this subspackage in a loop like one given below: |
//| > CAlgLib::SparseSolverOOCStart(state) |
//| > while CAlgLib::SparseSolverOOCContinue(state) do |
//| > CAlgLib::SparseSolverOOCGetRequestInfo(state, RequestType) |
//| > CAlgLib::SparseSolverOOCGetRequestData(state, out X) |
//| > if RequestType = 0 then |
//| > [calculate Y = A * X, with X = R ^ N] |
//| > CAlgLib::SparseSolverOOCSendResult(state, in Y) |
//| > CAlgLib::SparseSolverOOCStop(state, out X, out Report) |
//| INPUT PARAMETERS: |
//| State - solver object |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverOOCStart(CSparseSolverState &state,
CRowDouble &b)
{
state.m_rstate.ia.Resize(1);
state.m_rstate.ra.Resize(3);
state.m_rstate.stage=-1;
ClearRequestFields(state);
ClearReportFields(state);
state.m_running=true;
state.m_userterminationneeded=false;
CAblasF::RCopyV(state.m_n,b,state.m_b);
}
//+------------------------------------------------------------------+
//| This function performs iterative solution of the linear system in|
//| the out-of-core mode. It should be used in conjunction with other|
//| out-of-core - related functions of this subspackage in a loop |
//| like one given below: |
//| > CAlgLib::SparseSolverOOCStart(state) |
//| > while CAlgLib::SparseSolverOOCContinue(state) do |
//| > CAlgLib::SparseSolverOOCGetRequestInfo(state, RequestType) |
//| > CAlgLib::SparseSolverOOCGetRequestData(state, out X) |
//| > if RequestType = 0 then |
//| > [calculate Y = A * X, with X = R ^ N] |
//| > CAlgLib::SparseSolverOOCSendResult(state, in Y) |
//| > CAlgLib::SparseSolverOOCStop(state, out X, out Report) |
//| INPUT PARAMETERS: |
//| State - solver object |
//+------------------------------------------------------------------+
bool CIterativeSparse::SparseSolverOOCContinue(CSparseSolverState &state)
{
//--- check
if(!CAp::Assert(state.m_running,__FUNCTION__+": the solver is not running"))
return(false);
bool result=SparseSolverIteration(state);
state.m_running=result;
//--- return result
return(result);
}
//+------------------------------------------------------------------+
//| This function is used to retrieve information about out-of-core |
//| request sent by the solver: |
//| * RequestType = 0 means that matrix - vector products A * x is |
//| requested |
//| * RequestType = -1 means that solver reports its progress; this|
//| request is returned only when reports are |
//| activated wit SparseSolverSetXRep(). |
//| This function returns just request type; in order to get contents|
//| of the trial vector, use SparseSolverOOCGetRequestData(). |
//| It should be used in conjunction with other out-of-core - related|
//| functions of this subspackage in a loop like one given below: |
//| > CAlgLib::SparseSolverOOCStart(state) |
//| > while CAlgLib::SparseSolverOOCContinue(state) do |
//| > CAlgLib::SparseSolverOOCGetRequestInfo(state, RequestType) |
//| > CAlgLib::SparseSolverOOCGetRequestData(state, out X) |
//| > if RequestType = 0 then |
//| > [calculate Y = A * X, with X = R ^ N] |
//| > CAlgLib::SparseSolverOOCSendResult(state, in Y) |
//| > CAlgLib::SparseSolverOOCStop(state, out X, out Report) |
//| INPUT PARAMETERS: |
//| State - solver running in out-of-core mode |
//| OUTPUT PARAMETERS: |
//| RequestType - type of the request to process: |
//| * 0 for matrix - vector product A * x, with A |
//| being NxN system matrix and X being |
//| N-dimensional vector |
//| *-1 for location and residual report |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverOOCGetRequestInfo(CSparseSolverState &state,
int &requesttype)
{
requesttype=0;
//--- check
if(!CAp::Assert(state.m_running,__FUNCTION__+": the solver is not running"))
return;
requesttype=state.m_requesttype;
}
//+------------------------------------------------------------------+
//| This function is used to retrieve vector associated with |
//| out-of-core request sent by the solver to user code. Depending on|
//| the request type(returned by the SparseSolverOOCGetRequestInfo())|
//| this vector should be multiplied by A or subjected to another |
//| processing. |
//| It should be used in conjunction with other out-of-core - related|
//| functions of this subspackage in a loop like one given below: |
//| > CAlgLib::SparseSolverOOCStart(state) |
//| > while CAlgLib::SparseSolverOOCContinue(state) do |
//| > CAlgLib::SparseSolverOOCGetRequestInfo(state, RequestType) |
//| > CAlgLib::SparseSolverOOCGetRequestData(state, out X) |
//| > if RequestType = 0 then |
//| > [calculate Y = A * X, with X = R ^ N] |
//| > CAlgLib::SparseSolverOOCSendResult(state, in Y) |
//| > CAlgLib::SparseSolverOOCStop(state, out X, out Report) |
//| INPUT PARAMETERS: |
//| State - solver running in out-of-core mode |
//| X - possibly preallocated storage; reallocated if |
//| needed, left unchanged, if large enough to store|
//| request data. |
//| OUTPUT PARAMETERS: |
//| X - array[N] or larger, leading N elements are |
//| filled with vector X. |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverOOCGetRequestData(CSparseSolverState &state,
CRowDouble &x)
{
//--- check
if(!CAp::Assert(state.m_running,__FUNCTION__+": the solver is not running"))
return;
x=state.m_x;
}
//+------------------------------------------------------------------+
//| This function is used to retrieve scalar value associated with |
//| out-of-core request sent by the solver to user code. In the |
//| current ALGLIB version this function is used to retrieve squared |
//| residual norm during progress reports. |
//| INPUT PARAMETERS: |
//| State - solver running in out-of-core mode |
//| OUTPUT PARAMETERS: |
//| V - scalar value associated with the current request|
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverOOCGetRequestData1(CSparseSolverState &state,
double &v)
{
v=0;
//--- check
if(!CAp::Assert(state.m_running,__FUNCTION__+": the solver is not running"))
return;
v=state.m_reply1;
}
//+------------------------------------------------------------------+
//| This function is used to send user reply to out-of-core request |
//| sent by the solver. Usually it is product A*x for vector X |
//| returned by the solver. |
//| It should be used in conjunction with other out-of-core - related|
//| functions of this subspackage in a loop like one given below: |
//| > CAlgLib::SparseSolverOOCStart(state) |
//| > while CAlgLib::SparseSolverOOCContinue(state) do |
//| > CAlgLib::SparseSolverOOCGetRequestInfo(state, RequestType) |
//| > CAlgLib::SparseSolverOOCGetRequestData(state, out X) |
//| > if RequestType = 0 then |
//| > [calculate Y = A * X, with X = R ^ N] |
//| > CAlgLib::SparseSolverOOCSendResult(state, in Y) |
//| > CAlgLib::SparseSolverOOCStop(state, out X, out Report) |
//| INPUT PARAMETERS: |
//| State - solver running in out-of-core mode |
//| AX - array[N] or larger, leading N elements contain A*x |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverOOCSendResult(CSparseSolverState &state,
CRowDouble &ax)
{
//--- check
if(!CAp::Assert(state.m_running,__FUNCTION__+": the solver is not running"))
return;
if(!CAp::Assert(state.m_requesttype==0,__FUNCTION__+": this request type does not accept replies"))
return;
CAblasF::RCopyV(state.m_n,ax,state.m_ax);
}
//+------------------------------------------------------------------+
//| This function finalizes out-of-core mode of the linear solver. It|
//| should be used in conjunction with other out-of-core - related |
//| functions of this subspackage in a loop like one given below: |
//| > CAlgLib::SparseSolverOOCStart(state) |
//| > while CAlgLib::SparseSolverOOCContinue(state) do |
//| > CAlgLib::SparseSolverOOCGetRequestInfo(state, RequestType) |
//| > CAlgLib::SparseSolverOOCGetRequestData(state, out X) |
//| > if RequestType = 0 then |
//| > [calculate Y = A * X, with X = R ^ N] |
//| > CAlgLib::SparseSolverOOCSendResult(state, in Y) |
//| > CAlgLib::SparseSolverOOCStop(state, out X, out Report) |
//| INPUT PARAMETERS: |
//| State - solver state |
//| OUTPUT PARAMETERS: |
//| X - array[N], the solution. |
//| Zero - filled on the failure(Rep.TerminationType < 0). |
//| Rep - report with additional info: |
//| * Rep.TerminationType completion code: |
//| * -5 CG method was used for a matrix which is |
//| not positive definite |
//| * -4 overflow/underflow during solution (ill |
//| conditioned problem) |
//| * 1 || residual || <= EpsF* || b || |
//| * 5 MaxIts steps was taken |
//| * 7 rounding errors prevent further progress, |
//| best point found is returned |
//| * 8 the algorithm was terminated early with |
//| SparseSolverRequestTermination() being |
//| called from other thread. |
//| * Rep.IterationsCount contains iterations count |
//| * Rep.NMV contains number of matrix - vector |
//| calculations |
//| * Rep.R2 contains squared residual |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverOOCStop(CSparseSolverState &state,
CRowDouble &x,
CSparseSolverReport &rep)
{
x.Resize(0);
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": the solver is still running"))
return;
x.Resize(state.m_n);
CAblasF::RCopyV(state.m_n,state.m_xf,x);
CDirectSparseSolvers::InitSparseSolverReport(rep);
rep.m_iterationscount=state.m_repiterationscount;
rep.m_nmv=state.m_repnmv;
rep.m_terminationtype=state.m_repterminationtype;
rep.m_r2=state.m_repr2;
}
//+------------------------------------------------------------------+
//| This subroutine submits request for termination of the running |
//| solver. It can be called from some other thread which wants the |
//| solver to terminate or when processing an out-of-core request. |
//| As result, solver stops at point which was "current accepted" |
//| when the termination request was submitted and returns error code|
//| 8 (successful termination). Such termination is a smooth process |
//| which properly deallocates all temporaries. |
//| INPUT PARAMETERS: |
//| State - solver structure |
//| NOTE: calling this function on solver which is NOT running will |
//| have no effect. |
//| NOTE: multiple calls to this function are possible. First call is|
//| counted, subsequent calls are silently ignored. |
//| NOTE: solver clears termination flag on its start, it means that |
//| if some other thread will request termination too soon, its|
//| request will went unnoticed. |
//+------------------------------------------------------------------+
void CIterativeSparse::SparseSolverRequestTermination(CSparseSolverState &state)
{
state.m_userterminationneeded=true;
}
//+------------------------------------------------------------------+
//| Reverse communication sparse iteration subroutine |
//+------------------------------------------------------------------+
bool CIterativeSparse::SparseSolverIteration(CSparseSolverState &state)
{
//--- create variables
int outeridx=0;
double res=0;
double prevres=0;
double res0=0;
int label=-1;
//--- Reverse communication preparations
//--- This code initializes locals by:
//--- * random values determined during code
//--- generation - on first subroutine call
//--- * values from previous call - on subsequent calls
if(state.m_rstate.stage>=0)
{
outeridx=state.m_rstate.ia[0];
res=state.m_rstate.ra[0];
prevres=state.m_rstate.ra[1];
res0=state.m_rstate.ra[2];
}
else
{
outeridx=359;
res=-58;
prevres=-919;
res0=-909;
}
switch(state.m_rstate.stage)
{
case 0:
label=0;
break;
case 1:
label=1;
break;
case 2:
label=2;
break;
case 3:
label=3;
break;
case 4:
label=4;
break;
default:
//--- Routine body
state.m_running=true;
ClearRequestFields(state);
ClearReportFields(state);
//--- GMRES?
//
if(state.m_algotype!=0)
{
label=5;
break;
}
if(CAblasF::RDotV2(state.m_n,state.m_x0)!=0.0)
{
label=7;
break;
}
//--- Starting point is default one (zero), quick initialization
CAblasF::RSetV(state.m_n,0.0,state.m_xf);
CAblasF::RCopyV(state.m_n,state.m_b,state.m_wrkb);
label=8;
break;
}
//--- main loop
while(label>=0)
{
switch(label)
{
case 7:
//--- Non-zero starting point is provided,
CAblasF::RCopyV(state.m_n,state.m_x0,state.m_xf);
state.m_requesttype=0;
CAblasF::RCopyV(state.m_n,state.m_x0,state.m_x);
state.m_rstate.stage=0;
label=-1;
break;
case 0:
state.m_requesttype=-999;
state.m_repnmv=state.m_repnmv+1;
CAblasF::RCopyV(state.m_n,state.m_b,state.m_wrkb);
CAblasF::RAddV(state.m_n,-1.0,state.m_ax,state.m_wrkb);
case 8:
outeridx=0;
state.m_repterminationtype=5;
state.m_repr2=CAblasF::RDotV2(state.m_n,state.m_wrkb);
res0=MathSqrt(CAblasF::RDotV2(state.m_n,state.m_b));
res=MathSqrt(state.m_repr2);
if(!state.m_xrep)
{
label=9;
break;
}
//--- Report initial point
state.m_requesttype=-1;
state.m_reply1=res*res;
CAblasF::RCopyV(state.m_n,state.m_xf,state.m_x);
state.m_rstate.stage=1;
label=-1;
break;
case 1:
state.m_requesttype=-999;
case 9:
case 11:
if(!(res>0.0 && (state.m_maxits==0 || state.m_repiterationscount<state.m_maxits)))
{
label=12;
break;
}
//--- Solve with GMRES(k) for current residual.
//--- We set EpsF-based stopping condition for GMRES(k). It allows us
//--- to quickly detect sufficient decrease in the residual. We still
//--- have to recompute residual after the GMRES round because residuals
//--- computed by GMRES are different from the true one (due to restarts).
//--- However, checking residual decrease within GMRES still gives us
//--- an opportunity to stop early without waiting for GMRES round to
//--- complete.
CFbls::FblsGMRESCreate(state.m_wrkb,state.m_n,state.m_gmresk,state.m_gmressolver);
state.m_gmressolver.m_epsres=state.m_epsf*res0/res;
case 13:
if(!CFbls::FblsGMRESIteration(state.m_gmressolver))
{
label=14;
break;
}
state.m_requesttype=0;
CAblasF::RCopyV(state.m_n,state.m_gmressolver.m_x,state.m_x);
state.m_rstate.stage=2;
label=-1;
break;
case 2:
state.m_requesttype=-999;
CAblasF::RCopyV(state.m_n,state.m_ax,state.m_gmressolver.m_ax);
state.m_repnmv++;
if(state.m_userterminationneeded)
{
//--- User requested termination
state.m_repterminationtype=8;
return(false);
}
label=13;
break;
case 14:
state.m_repiterationscount+=state.m_gmressolver.m_itsperformed;
CAblasF::RAddV(state.m_n,1.0,state.m_gmressolver.m_xs,state.m_xf);
//--- Update residual, evaluate residual decrease, terminate if needed
state.m_requesttype=0;
CAblasF::RCopyV(state.m_n,state.m_xf,state.m_x);
state.m_rstate.stage=3;
label=-1;
break;
case 3:
state.m_requesttype=-999;
state.m_repnmv++;
CAblasF::RCopyV(state.m_n,state.m_b,state.m_wrkb);
CAblasF::RAddV(state.m_n,-1.0,state.m_ax,state.m_wrkb);
state.m_repr2=CAblasF::RDotV2(state.m_n,state.m_wrkb);
prevres=res;
res=MathSqrt(state.m_repr2);
if(!state.m_xrep)
{
label=15;
break;
}
//--- Report initial point
state.m_requesttype=-1;
state.m_reply1=res*res;
CAblasF::RCopyV(state.m_n,state.m_xf,state.m_x);
state.m_rstate.stage=4;
label=-1;
break;
case 4:
state.m_requesttype=-999;
case 15:
if(res<=(state.m_epsf*res0))
{
//--- Residual decrease condition met, stopping
state.m_repterminationtype=1;
label=12;
break;
}
if(res>=(prevres*(1-MathSqrt(CMath::m_machineepsilon))))
{
//--- The algorithm stagnated
state.m_repterminationtype=7;
label=12;
break;
}
if(state.m_userterminationneeded)
{
//--- User requested termination
state.m_repterminationtype=8;
return(false);
}
outeridx++;
label=11;
break;
case 12:
return(false);
case 5:
CAp::Assert(false,__FUNCTION__+": integrity check failed (unexpected algo)");
return(false);
}
}
//--- Saving state
state.m_rstate.ia.Set(0,outeridx);
state.m_rstate.ra.Set(0,res);
state.m_rstate.ra.Set(1,prevres);
state.m_rstate.ra.Set(2,res0);
//--- return result
return(true);
}
//+------------------------------------------------------------------+
//| Clears request fileds (to be sure that we don't forgot to clear |
//| something) |
//+------------------------------------------------------------------+
void CIterativeSparse::ClearRequestFields(CSparseSolverState &state)
{
state.m_requesttype=-999;
}
//+------------------------------------------------------------------+
//| Clears report fileds (to be sure that we don't forgot to clear |
//| something) |
//+------------------------------------------------------------------+
void CIterativeSparse::ClearReportFields(CSparseSolverState &state)
{
state.m_repiterationscount=0;
state.m_repnmv=0;
state.m_repterminationtype=0;
state.m_repr2=0;
}
//+------------------------------------------------------------------+
//| This object stores state of the linear CG method. |
//| You should use ALGLIB functions to work with this object. |
//| Never try to access its fields directly! |
//+------------------------------------------------------------------+
struct CLinCGState
{
int m_itsbeforerestart;
int m_itsbeforerupdate;
int m_maxits;
int m_n;
int m_prectype;
int m_repiterationscount;
int m_repnmv;
int m_repterminationtype;
double m_alpha;
double m_beta;
double m_epsf;
double m_meritfunction;
double m_r2;
double m_vmv;
bool m_needmtv;
bool m_needmv2;
bool m_needmv;
bool m_needprec;
bool m_needvmv;
bool m_running;
bool m_xrep;
bool m_xupdated;
RCommState m_rstate;
CRowDouble m_b;
CRowDouble m_cr;
CRowDouble m_cx;
CRowDouble m_cz;
CRowDouble m_mv;
CRowDouble m_p;
CRowDouble m_pv;
CRowDouble m_r;
CRowDouble m_rx;
CRowDouble m_startx;
CRowDouble m_tmpd;
CRowDouble m_x;
CRowDouble m_z;
//--- constructor / destructor
CLinCGState(void);
~CLinCGState(void) {}
//---
void Copy(const CLinCGState &obj);
//--- overloading
void operator=(const CLinCGState &obj) { Copy(obj); }
};
//+------------------------------------------------------------------+
//| Constructor |
//+------------------------------------------------------------------+
CLinCGState::CLinCGState(void)
{
m_itsbeforerestart=0;
m_itsbeforerupdate=0;
m_maxits=0;
m_n=0;
m_prectype=0;
m_repiterationscount=0;
m_repnmv=0;
m_repterminationtype=0;
m_alpha=0;
m_beta=0;
m_epsf=0;
m_meritfunction=0;
m_r2=0;
m_vmv=0;
m_needmtv=false;
m_needmv2=false;
m_needmv=false;
m_needprec=false;
m_needvmv=false;
m_running=false;
m_xrep=false;
m_xupdated=false;
}
//+------------------------------------------------------------------+
//| Copy |
//+------------------------------------------------------------------+
void CLinCGState::Copy(const CLinCGState &obj)
{
m_itsbeforerestart=obj.m_itsbeforerestart;
m_itsbeforerupdate=obj.m_itsbeforerupdate;
m_maxits=obj.m_maxits;
m_n=obj.m_n;
m_prectype=obj.m_prectype;
m_repiterationscount=obj.m_repiterationscount;
m_repnmv=obj.m_repnmv;
m_repterminationtype=obj.m_repterminationtype;
m_alpha=obj.m_alpha;
m_beta=obj.m_beta;
m_epsf=obj.m_epsf;
m_meritfunction=obj.m_meritfunction;
m_r2=obj.m_r2;
m_vmv=obj.m_vmv;
m_needmtv=obj.m_needmtv;
m_needmv2=obj.m_needmv2;
m_needmv=obj.m_needmv;
m_needprec=obj.m_needprec;
m_needvmv=obj.m_needvmv;
m_running=obj.m_running;
m_xrep=obj.m_xrep;
m_xupdated=obj.m_xupdated;
m_rstate=obj.m_rstate;
m_b=obj.m_b;
m_cr=obj.m_cr;
m_cx=obj.m_cx;
m_cz=obj.m_cz;
m_mv=obj.m_mv;
m_p=obj.m_p;
m_pv=obj.m_pv;
m_r=obj.m_r;
m_rx=obj.m_rx;
m_startx=obj.m_startx;
m_tmpd=obj.m_tmpd;
m_x=obj.m_x;
m_z=obj.m_z;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
struct CLinCGReport
{
int m_iterationscount;
int m_nmv;
int m_terminationtype;
double m_r2;
//--- constructor / destructor
CLinCGReport(void) { ZeroMemory(this); }
~CLinCGReport(void) {}
//---
void Copy(const CLinCGReport &obj);
//--- overloading
void operator=(const CLinCGReport &obj) { Copy(obj); }
};
//+------------------------------------------------------------------+
//| Copy |
//+------------------------------------------------------------------+
void CLinCGReport::Copy(const CLinCGReport &obj)
{
m_iterationscount=obj.m_iterationscount;
m_nmv=obj.m_nmv;
m_terminationtype=obj.m_terminationtype;
m_r2=obj.m_r2;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
class CLinCG
{
public:
//--- constants
static const double m_defaultprecision;
static void LinCGCreate(int n,CLinCGState &state);
static void LinCGSetStartingPoint(CLinCGState &state,CRowDouble &x);
static void LinCGSetB(CLinCGState &state,CRowDouble &b);
static void LinCGSetPrecUnit(CLinCGState &state);
static void LinCGSetPrecDiag(CLinCGState &state);
static void LinCGSetCond(CLinCGState &state,double epsf,int maxits);
static bool LinCGIteration(CLinCGState &state);
static void LinCGSolveSparse(CLinCGState &state,CSparseMatrix &a,bool IsUpper,CRowDouble &b);
static void LinCGResult(CLinCGState &state,CRowDouble &x,CLinCGReport &rep);
static void LinCGSetRestartFreq(CLinCGState &state,int srf);
static void LinCGSetRUpdateFreq(CLinCGState &state,int freq);
static void LinCGSetXRep(CLinCGState &state,bool needxrep);
static void LinCGRestart(CLinCGState &state);
private:
static void ClearrFields(CLinCGState &state);
static void UpdateItersData(CLinCGState &state);
};
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
const double CLinCG::m_defaultprecision=1.0E-6;
//+------------------------------------------------------------------+
//| This function initializes linear CG Solver. This solver is used |
//| to solve symmetric positive definite problems. If you want to |
//| solve nonsymmetric (or non-positive definite) problem you may use|
//| LinLSQR solver provided by ALGLIB. |
//| USAGE: |
//| 1. User initializes algorithm state with LinCGCreate() call |
//| 2. User tunes solver parameters with LinCGSetCond() and other |
//| functions |
//| 3. Optionally, user sets starting point with |
//| LinCGSetStartingPoint() |
//| 4. User calls LinCGSolveSparse() function which takes algorithm|
//| state and SparseMatrix object. |
//| 5. User calls LinCGResults() to get solution |
//| 6. Optionally, user may call LinCGSolveSparse() again to solve |
//| another problem with different matrix and/or right part |
//| without reinitializing LinCGState structure. |
//| INPUT PARAMETERS: |
//| N - problem dimension, N > 0 |
//| OUTPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//+------------------------------------------------------------------+
void CLinCG::LinCGCreate(int n,CLinCGState &state)
{
//--- check
if(!CAp::Assert(n>0,__FUNCTION__+": N<=0"))
return;
state.m_n=n;
state.m_prectype=0;
state.m_itsbeforerestart=n;
state.m_itsbeforerupdate=10;
state.m_epsf=m_defaultprecision;
state.m_maxits=0;
state.m_xrep=false;
state.m_running=false;
//--- * allocate arrays
//--- * set RX to NAN (just for the case user calls Results() without
//--- calling SolveSparse()
//--- * set starting point to zero
//--- * we do NOT initialize B here because we assume that user should
//--- initializate it using LinCGSetB() function. In case he forgets
//--- to do so, exception will be thrown in the LinCGIteration().
state.m_rx=vector<double>::Full(state.m_n,AL_NaN);
state.m_startx=vector<double>::Zeros(state.m_n);
state.m_b=vector<double>::Zeros(state.m_n);
state.m_cx=vector<double>::Zeros(state.m_n);
state.m_p=vector<double>::Zeros(state.m_n);
state.m_r=vector<double>::Zeros(state.m_n);
state.m_cr=vector<double>::Zeros(state.m_n);
state.m_z=vector<double>::Zeros(state.m_n);
state.m_cz=vector<double>::Zeros(state.m_n);
state.m_x=vector<double>::Zeros(state.m_n);
state.m_mv=vector<double>::Zeros(state.m_n);
state.m_pv=vector<double>::Zeros(state.m_n);
UpdateItersData(state);
state.m_rstate.ia.Resize(1);
state.m_rstate.ra=vector<double>::Zeros(3);
state.m_rstate.stage=-1;
}
//+------------------------------------------------------------------+
//| This function sets starting point. |
//| By default, zero starting point is used. |
//| INPUT PARAMETERS: |
//| X - starting point, array[N] |
//| OUTPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//+------------------------------------------------------------------+
void CLinCG::LinCGSetStartingPoint(CLinCGState &state,
CRowDouble &x)
{
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": you can not change starting point because LinCGIteration() function is running"))
return;
if(!CAp::Assert(state.m_n<=CAp::Len(x),__FUNCTION__+": Length(X)<N"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(x,state.m_n),__FUNCTION__+": X contains infinite or NaN values!"))
return;
//--- copy
state.m_startx=x;
state.m_startx.Resize(state.m_n);
}
//+------------------------------------------------------------------+
//| This function sets right part. By default, right part is zero. |
//| INPUT PARAMETERS: |
//| B - right part, array[N]. |
//| OUTPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//+------------------------------------------------------------------+
void CLinCG::LinCGSetB(CLinCGState &state,CRowDouble &b)
{
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": you can not set B,because function LinCGIteration is running!"))
return;
if(!CAp::Assert(CAp::Len(b)>=state.m_n,__FUNCTION__+": Length(B)<N"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(b,state.m_n),__FUNCTION__+": B contains infinite or NaN values!"))
return;
//--- copy
state.m_b=b;
state.m_b.Resize(state.m_n);
}
//+------------------------------------------------------------------+
//| This function changes preconditioning settings of |
//| LinCGSolveSparse() function. By default, SolveSparse() uses |
//| diagonal preconditioner, but if you want to use solver without |
//| preconditioning, you can call this function which forces solver |
//| to use unit matrix for preconditioning. |
//| INPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//+------------------------------------------------------------------+
void CLinCG::LinCGSetPrecUnit(CLinCGState &state)
{
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": you can not change preconditioner,because function LinCGIteration is running!"))
return;
state.m_prectype=-1;
}
//+------------------------------------------------------------------+
//| This function changes preconditioning settings of |
//| LinCGSolveSparse() function. LinCGSolveSparse() will use diagonal|
//| of the system matrix as preconditioner. This preconditioning mode|
//| is active by default. |
//| INPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//+------------------------------------------------------------------+
void CLinCG::LinCGSetPrecDiag(CLinCGState &state)
{
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": you can not change preconditioner,because function LinCGIteration is running!"))
return;
state.m_prectype=0;
}
//+------------------------------------------------------------------+
//| This function sets stopping criteria. |
//| INPUT PARAMETERS: |
//| EpsF - algorithm will be stopped if norm of residual is |
//| less than EpsF* || b ||. |
//| MaxIts - algorithm will be stopped if number of iterations |
//| is more than MaxIts. |
//| OUTPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//| NOTES: If both EpsF and MaxIts are zero then small EpsF will be |
//| set to small value. |
//+------------------------------------------------------------------+
void CLinCG::LinCGSetCond(CLinCGState &state,double epsf,int maxits)
{
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": you can not change stopping criteria when LinCGIteration() is running"))
return;
if(!CAp::Assert(MathIsValidNumber(epsf) && epsf>=0.0,__FUNCTION__+": EpsF is negative or contains infinite or NaN values"))
return;
if(!CAp::Assert(maxits>=0,__FUNCTION__+": MaxIts is negative"))
return;
if(epsf==0.0 && maxits==0)
{
state.m_epsf=m_defaultprecision;
state.m_maxits=maxits;
}
else
{
state.m_epsf=epsf;
state.m_maxits=maxits;
}
}
//+------------------------------------------------------------------+
//| Reverse communication version of linear CG. |
//+------------------------------------------------------------------+
bool CLinCG::LinCGIteration(CLinCGState &state)
{
//--- create variables
int i=0;
double uvar=0;
double bnorm=0;
double v=0;
int label=-1;
//--- Reverse communication preparations
//--- This code initializes locals by:
//--- * random values determined during code
//--- generation - on first subroutine call
//--- * values from previous call - on subsequent calls
//
if(state.m_rstate.stage>=0)
{
i=state.m_rstate.ia[0];
uvar=state.m_rstate.ra[0];
bnorm=state.m_rstate.ra[1];
v=state.m_rstate.ra[2];
}
else
{
i=359;
uvar=-58;
bnorm=-919;
v=-909;
}
switch(state.m_rstate.stage)
{
case 0:
label=0;
break;
case 1:
label=1;
break;
case 2:
label=2;
break;
case 3:
label=3;
break;
case 4:
label=4;
break;
case 5:
label=5;
break;
case 6:
label=6;
break;
case 7:
label=7;
break;
default:
//--- Routine body
//--- check
if(!CAp::Assert(CAp::Len(state.m_b)>0,__FUNCTION__+": B is not initialized (you must initialize B by LinCGSetB() call"))
return(false);
state.m_running=true;
state.m_repnmv=0;
ClearrFields(state);
UpdateItersData(state);
//--- Start 0-th iteration
state.m_rx=state.m_startx;
state.m_x=state.m_rx;
state.m_repnmv++;
ClearrFields(state);
state.m_needvmv=true;
state.m_rstate.stage=0;
label=-1;
break;
}
//--- main loop
while(label>=0)
{
switch(label)
{
case 0:
state.m_needvmv=false;
bnorm=0;
state.m_r2=0;
state.m_meritfunction=0;
state.m_r=state.m_b-state.m_mv+0;
state.m_r2=state.m_r.Dot(state.m_r);
state.m_meritfunction= state.m_mv.Dot(state.m_rx)-2*state.m_b.Dot(state.m_rx);
bnorm=state.m_b.Dot(state.m_b);
bnorm=MathSqrt(bnorm);
//--- Output first report
if(!state.m_xrep)
{
label=8;
break;
}
state.m_x=state.m_rx;
ClearrFields(state);
state.m_xupdated=true;
state.m_rstate.stage=1;
label=-1;
break;
case 1:
state.m_xupdated=false;
case 8:
//--- Is x0 a solution?
if(!MathIsValidNumber(state.m_r2) || MathSqrt(state.m_r2)<=(state.m_epsf*bnorm))
{
state.m_running=false;
if(MathIsValidNumber(state.m_r2))
state.m_repterminationtype=1;
else
state.m_repterminationtype=-4;
return(false);
}
//--- Calculate Z and P
state.m_x=state.m_r;
state.m_repnmv++;
ClearrFields(state);
state.m_needprec=true;
state.m_rstate.stage=2;
label=-1;
break;
case 2:
state.m_needprec=false;
state.m_z=state.m_pv;
state.m_p=state.m_z;
//--- Other iterations(1..N)
state.m_repiterationscount=0;
case 10:
state.m_repiterationscount++;
//--- Calculate Alpha
state.m_x=state.m_p;
state.m_repnmv++;
ClearrFields(state);
state.m_needvmv=true;
state.m_rstate.stage=3;
label=-1;
break;
case 3:
state.m_needvmv=false;
if(!MathIsValidNumber(state.m_vmv) || state.m_vmv<=0.0)
{
//--- a) Overflow when calculating VMV
//--- b) non-positive VMV (non-SPD matrix)
state.m_running=false;
if(MathIsValidNumber(state.m_vmv))
state.m_repterminationtype=-5;
else
state.m_repterminationtype=-4;
return(false);
}
state.m_alpha=state.m_r.Dot(state.m_z);
state.m_alpha=state.m_alpha/state.m_vmv;
if(!MathIsValidNumber(state.m_alpha))
{
//--- Overflow when calculating Alpha
state.m_running=false;
state.m_repterminationtype=-4;
return(false);
}
//--- Next step toward solution
state.m_cx=state.m_rx.ToVector()+state.m_p*state.m_alpha;
//--- Calculate R:
//--- * use recurrent relation to update R
//--- * at every ItsBeforeRUpdate-th iteration recalculate it from scratch, using matrix-vector product
//--- in case R grows instead of decreasing, algorithm is terminated with positive completion code
if(!(state.m_itsbeforerupdate==0 || state.m_repiterationscount%state.m_itsbeforerupdate!=0))
{
label=12;
break;
}
//--- Calculate R using recurrent formula
state.m_cr=state.m_r.ToVector()-state.m_mv*state.m_alpha;
state.m_x=state.m_cr;
label=13;
break;
case 12:
//--- Calculate R using matrix-vector multiplication
state.m_x=state.m_cx;
state.m_repnmv++;
ClearrFields(state);
state.m_needmv=true;
state.m_rstate.stage=4;
label=-1;
break;
case 4:
state.m_needmv=false;
state.m_cr=state.m_b-state.m_mv+0;
state.m_x=state.m_cr;
//--- Calculating merit function
//--- Check emergency stopping criterion
v=state.m_mv.Dot(state.m_cx)-2*state.m_b.Dot(state.m_cx);
if(v<state.m_meritfunction)
{
label=14;
break;
}
if(!CApServ::IsFiniteVector(state.m_rx,state.m_n))
{
state.m_running=false;
state.m_repterminationtype=-4;
return(false);
}
//output last report
if(!state.m_xrep)
{
label=16;
break;
}
state.m_x=state.m_rx;
ClearrFields(state);
state.m_xupdated=true;
state.m_rstate.stage=5;
label=-1;
break;
case 5:
state.m_xupdated=false;
case 16:
state.m_running=false;
state.m_repterminationtype=7;
return(false);
case 14:
state.m_meritfunction=v;
case 13:
state.m_rx=state.m_cx;
//--- calculating RNorm
//--- NOTE: monotonic decrease of R2 is not guaranteed by algorithm.
state.m_r2=state.m_cr.Dot(state.m_cr);
//output report
if(!state.m_xrep)
{
label=18;
break;
}
state.m_x=state.m_rx;
ClearrFields(state);
state.m_xupdated=true;
state.m_rstate.stage=6;
label=-1;
break;
case 6:
state.m_xupdated=false;
case 18:
//stopping criterion
//achieved the required precision
if(!MathIsValidNumber(state.m_r2) || (MathSqrt(state.m_r2)<=(state.m_epsf*bnorm)))
{
state.m_running=false;
if(MathIsValidNumber(state.m_r2))
state.m_repterminationtype=1;
else
state.m_repterminationtype=-4;
return(false);
}
if(state.m_repiterationscount>=state.m_maxits && state.m_maxits>0)
{
if(!CApServ::IsFiniteVector(state.m_rx,state.m_n))
{
state.m_running=false;
state.m_repterminationtype=-4;
return(false);
}
//if X is finite number
state.m_running=false;
state.m_repterminationtype=5;
return(false);
}
state.m_x=state.m_cr;
//prepere of parameters for next iteration
state.m_repnmv++;
ClearrFields(state);
state.m_needprec=true;
state.m_rstate.stage=7;
label=-1;
break;
case 7:
state.m_needprec=false;
state.m_cz=state.m_pv;
if(state.m_repiterationscount%state.m_itsbeforerestart!=0)
{
state.m_beta=state.m_cz.Dot(state.m_cr);
uvar=state.m_z.Dot(state.m_r);
//check that UVar is't INF or is't zero
if(!MathIsValidNumber(uvar) || uvar==0.0)
{
state.m_running=false;
state.m_repterminationtype=-4;
return(false);
}
//calculate .BETA
state.m_beta/=uvar;
//check that .BETA neither INF nor NaN
if(!MathIsValidNumber(state.m_beta))
{
state.m_running=false;
state.m_repterminationtype=-1;
return(false);
}
state.m_p=state.m_cz.ToVector()+state.m_p*state.m_beta;
}
else
state.m_p=state.m_cz;
//prepere data for next iteration
//write (k+1)th iteration to (k )th iteration
state.m_r=state.m_cr;
state.m_z=state.m_cz;
label=10;
break;
case 11:
return(false);
}
}
//--- Saving state
state.m_rstate.ia.Set(0,i);
state.m_rstate.ra.Set(0,uvar);
state.m_rstate.ra.Set(1,bnorm);
state.m_rstate.ra.Set(2,v);
//--- return result
return(true);
}
//+------------------------------------------------------------------+
//| Procedure for solution of A*x = b with sparse A. |
//| INPUT PARAMETERS: |
//| State - algorithm state |
//| A - sparse matrix in the CRS format (you MUST |
//| contvert it to CRS format by calling |
//| SparseConvertToCRS() function). |
//| IsUpper - whether upper or lower triangle of A is used: |
//| * IsUpper = True => only upper triangle is used |
//| and lower triangle is not |
//| referenced at all |
//| * IsUpper = False => only lower triangle is used|
//| and upper triangle is not |
//| referenced at all |
//| B - right part, array[N] |
//| RESULT: |
//| This function returns no result. |
//| You can get solution by calling LinCGResults() |
//| NOTE: this function uses lightweight preconditioning - |
//| multiplication by inverse of diag(A). If you want, you can |
//| turn preconditioning off by calling LinCGSetPrecUnit(). |
//| However, preconditioning cost is low and preconditioner is |
//| very important for solution of badly scaled problems. |
//+------------------------------------------------------------------+
void CLinCG::LinCGSolveSparse(CLinCGState &state,CSparseMatrix &a,
bool IsUpper,CRowDouble &b)
{
//--- create variables
int n=state.m_n;
double v=0;
double vmv=0;
//--- check
if(!CAp::Assert(CAp::Len(b)>=state.m_n,__FUNCTION__+": Length(B)<N"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(b,state.m_n),__FUNCTION__+": B contains infinite or NaN values!"))
return;
//--- Allocate temporaries
state.m_tmpd.Resize(n);
//--- Compute diagonal scaling matrix D
if(state.m_prectype==0)
{
//--- Default preconditioner - inverse of matrix diagonal
for(int i=0; i<n; i++)
{
v=CSparse::SparseGetDiagonal(a,i);
if(v>0.0)
state.m_tmpd.Set(i,1/MathSqrt(v));
else
state.m_tmpd.Set(i,1);
}
}
else
{
//--- No diagonal scaling
state.m_tmpd.Fill(1);
}
//--- Solve
LinCGRestart(state);
LinCGSetB(state,b);
while(LinCGIteration(state))
{
//--- Process different requests from optimizer
if(state.m_needmv)
CSparse::SparseSMV(a,IsUpper,state.m_x,state.m_mv);
if(state.m_needvmv)
{
CSparse::SparseSMV(a,IsUpper,state.m_x,state.m_mv);
vmv=state.m_x.Dot(state.m_mv);
state.m_vmv=vmv;
}
if(state.m_needprec)
state.m_pv=state.m_x.ToVector()*state.m_tmpd.Pow(2);
}
}
//+------------------------------------------------------------------+
//| CG - solver: results. |
//| This function must be called after LinCGSolve |
//| INPUT PARAMETERS: |
//| State - algorithm state |
//| OUTPUT PARAMETERS: |
//| X - array[N], solution |
//| Rep - optimization report: |
//| * Rep.TerminationType completetion code: |
//| * -5 input matrix is either not positive |
//| definite, too large or too small |
//| * -4 overflow / underflow during solution |
//| (ill conditioned problem) |
//| * 1 || residual || <= EpsF* || b || |
//| * 5 MaxIts steps was taken |
//| * 7 rounding errors prevent further |
//| progress, best point found is returned |
//| * Rep.IterationsCount contains iterations count |
//| * NMV countains number of matrix - vector |
//| calculations |
//+------------------------------------------------------------------+
void CLinCG::LinCGResult(CLinCGState &state,CRowDouble &x,
CLinCGReport &rep)
{
x.Resize(0);
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": you can not get result,because function LinCGIteration has been launched!"))
return;
if(CAp::Len(x)<state.m_n)
x.Resize(state.m_n);
x=state.m_rx;
rep.m_iterationscount=state.m_repiterationscount;
rep.m_nmv=state.m_repnmv;
rep.m_terminationtype=state.m_repterminationtype;
rep.m_r2=state.m_r2;
}
//+------------------------------------------------------------------+
//| This function sets restart frequency. By default, algorithm is |
//| restarted after N subsequent iterations. |
//+------------------------------------------------------------------+
void CLinCG::LinCGSetRestartFreq(CLinCGState &state,int srf)
{
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": you can not change restart frequency when LinCGIteration() is running"))
return;
if(!CAp::Assert(srf>0,__FUNCTION__+": non-positive SRF"))
return;
state.m_itsbeforerestart=srf;
}
//+------------------------------------------------------------------+
//| This function sets frequency of residual recalculations. |
//| Algorithm updates residual r_k using iterative formula, but |
//| recalculates it from scratch after each 10 iterations. It is done|
//| to avoid accumulation of numerical errors and to stop algorithm |
//| when r_k starts to grow. |
//| Such low update frequence(1 / 10) gives very little overhead, |
//| but makes algorithm a bit more robust against numerical errors. |
//| However, you may change it |
//| INPUT PARAMETERS: |
//| Freq - desired update frequency, Freq >= 0. |
//| Zero value means that no updates will be done. |
//+------------------------------------------------------------------+
void CLinCG::LinCGSetRUpdateFreq(CLinCGState &state,int freq)
{
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": you can not change update frequency when LinCGIteration() is running"))
return;
if(!CAp::Assert(freq>=0,__FUNCTION__+": non-positive Freq"))
return;
state.m_itsbeforerupdate=freq;
}
//+------------------------------------------------------------------+
//| This function turns on / off reporting. |
//| INPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//| NeedXRep - whether iteration reports are needed or not |
//| If NeedXRep is True, algorithm will call rep() callback function |
//| if it is provided to MinCGOptimize(). |
//+------------------------------------------------------------------+
void CLinCG::LinCGSetXRep(CLinCGState &state,bool needxrep)
{
state.m_xrep=needxrep;
}
//+------------------------------------------------------------------+
//| Procedure for restart function LinCGIteration |
//+------------------------------------------------------------------+
void CLinCG::LinCGRestart(CLinCGState &state)
{
state.m_rstate.ia.Resize(1);
state.m_rstate.ra.Resize(3);
state.m_rstate.stage=-1;
ClearrFields(state);
}
//+------------------------------------------------------------------+
//| Clears request fileds (to be sure that we don't forgot to clear |
//| something) |
//+------------------------------------------------------------------+
void CLinCG::ClearrFields(CLinCGState &state)
{
state.m_xupdated=false;
state.m_needmv=false;
state.m_needmtv=false;
state.m_needmv2=false;
state.m_needvmv=false;
state.m_needprec=false;
}
//+------------------------------------------------------------------+
//| Clears request fileds (to be sure that we don't forgot to clear |
//| something) |
//+------------------------------------------------------------------+
void CLinCG::UpdateItersData(CLinCGState &state)
{
state.m_repiterationscount=0;
state.m_repnmv=0;
state.m_repterminationtype=0;
}
//+------------------------------------------------------------------+
//| This object stores state of the LinLSQR method. |
//| You should use ALGLIB functions to work with this object. |
//+------------------------------------------------------------------+
struct CLinLSQRState
{
int m_m;
int m_maxits;
int m_n;
int m_prectype;
int m_repiterationscount;
int m_repnmv;
int m_repterminationtype;
double m_alphai;
double m_alphaip1;
double m_anorm;
double m_betai;
double m_betaip1;
double m_bnorm2;
double m_ci;
double m_dnorm;
double m_epsa;
double m_epsb;
double m_epsc;
double m_lambdai;
double m_phibari;
double m_phibarip1;
double m_phii;
double m_r2;
double m_rhobari;
double m_rhobarip1;
double m_rhoi;
double m_si;
double m_theta;
bool m_needmtv;
bool m_needmv2;
bool m_needmv;
bool m_needprec;
bool m_needvmv;
bool m_running;
bool m_userterminationneeded;
bool m_xrep;
bool m_xupdated;
RCommState m_rstate;
CRowDouble m_b;
CRowDouble m_d;
CRowDouble m_mtv;
CRowDouble m_mv;
CRowDouble m_omegai;
CRowDouble m_omegaip1;
CRowDouble m_rx;
CRowDouble m_tmpd;
CRowDouble m_tmpx;
CRowDouble m_ui;
CRowDouble m_uip1;
CRowDouble m_vi;
CRowDouble m_vip1;
CRowDouble m_x;
CNormEstimatorState m_nes;
//--- constructor / destructor
CLinLSQRState(void);
~CLinLSQRState(void) {}
//---
void Copy(const CLinLSQRState &obj);
//--- overloading
void operator=(const CLinLSQRState &obj) { Copy(obj); }
};
//+------------------------------------------------------------------+
//| Constructor |
//+------------------------------------------------------------------+
CLinLSQRState::CLinLSQRState(void)
{
m_m=0;
m_maxits=0;
m_n=0;
m_prectype=0;
m_repiterationscount=0;
m_repnmv=0;
m_repterminationtype=0;
m_alphai=0;
m_alphaip1=0;
m_anorm=0;
m_betai=0;
m_betaip1=0;
m_bnorm2=0;
m_ci=0;
m_dnorm=0;
m_epsa=0;
m_epsb=0;
m_epsc=0;
m_lambdai=0;
m_phibari=0;
m_phibarip1=0;
m_phii=0;
m_r2=0;
m_rhobari=0;
m_rhobarip1=0;
m_rhoi=0;
m_si=0;
m_theta=0;
m_needmtv=false;
m_needmv2=false;
m_needmv=false;
m_needprec=false;
m_needvmv=false;
m_running=false;
m_userterminationneeded=false;
m_xrep=false;
m_xupdated=false;
}
//+------------------------------------------------------------------+
//| Copy |
//+------------------------------------------------------------------+
void CLinLSQRState::Copy(const CLinLSQRState &obj)
{
m_m=obj.m_m;
m_maxits=obj.m_maxits;
m_n=obj.m_n;
m_prectype=obj.m_prectype;
m_repiterationscount=obj.m_repiterationscount;
m_repnmv=obj.m_repnmv;
m_repterminationtype=obj.m_repterminationtype;
m_alphai=obj.m_alphai;
m_alphaip1=obj.m_alphaip1;
m_anorm=obj.m_anorm;
m_betai=obj.m_betai;
m_betaip1=obj.m_betaip1;
m_bnorm2=obj.m_bnorm2;
m_ci=obj.m_ci;
m_dnorm=obj.m_dnorm;
m_epsa=obj.m_epsa;
m_epsb=obj.m_epsb;
m_epsc=obj.m_epsc;
m_lambdai=obj.m_lambdai;
m_phibari=obj.m_phibari;
m_phibarip1=obj.m_phibarip1;
m_phii=obj.m_phii;
m_r2=obj.m_r2;
m_rhobari=obj.m_rhobari;
m_rhobarip1=obj.m_rhobarip1;
m_rhoi=obj.m_rhoi;
m_si=obj.m_si;
m_theta=obj.m_theta;
m_needmtv=obj.m_needmtv;
m_needmv2=obj.m_needmv2;
m_needmv=obj.m_needmv;
m_needprec=obj.m_needprec;
m_needvmv=obj.m_needvmv;
m_running=obj.m_running;
m_userterminationneeded=obj.m_userterminationneeded;
m_xrep=obj.m_xrep;
m_xupdated=obj.m_xupdated;
m_rstate=obj.m_rstate;
m_b=obj.m_b;
m_d=obj.m_d;
m_mtv=obj.m_mtv;
m_mv=obj.m_mv;
m_omegai=obj.m_omegai;
m_omegaip1=obj.m_omegaip1;
m_rx=obj.m_rx;
m_tmpd=obj.m_tmpd;
m_tmpx=obj.m_tmpx;
m_ui=obj.m_ui;
m_uip1=obj.m_uip1;
m_vi=obj.m_vi;
m_vip1=obj.m_vip1;
m_x=obj.m_x;
m_nes=obj.m_nes;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
struct CLinLSQRReport
{
int m_iterationscount;
int m_nmv;
int m_terminationtype;
//--- constructor / destructor
CLinLSQRReport(void) { ZeroMemory(this); }
~CLinLSQRReport(void) {}
//---
void Copy(const CLinLSQRReport &obj);
//--- overloading
void operator=(const CLinLSQRReport &obj) { Copy(obj); }
};
//+------------------------------------------------------------------+
//| Copy |
//+------------------------------------------------------------------+
void CLinLSQRReport::Copy(const CLinLSQRReport &obj)
{
m_iterationscount=obj.m_iterationscount;
m_nmv=obj.m_nmv;
m_terminationtype=obj.m_terminationtype;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
class CLinLSQR
{
public:
//--- constants
static const double m_atol;
static const double m_btol;
static void LinLSQRCreate(int m,int n,CLinLSQRState &state);
static void LinLSQRCreateBuf(int m,int n,CLinLSQRState &state);
static void LinLSQRSetB(CLinLSQRState &state,CRowDouble &b);
static void LinLSQRSetPrecUnit(CLinLSQRState &state);
static void LinLSQRSetPrecDiag(CLinLSQRState &state);
static void LinLSQRSetLambdaI(CLinLSQRState &state,double lambdai);
static bool LinLSQRIteration(CLinLSQRState &state);
static void LinLSQRSolveSparse(CLinLSQRState &state,CSparseMatrix &a,CRowDouble &b);
static void LinLSQRSetCond(CLinLSQRState &state,double epsa,double epsb,int maxits);
static void LinLSQRResults(CLinLSQRState &state,CRowDouble &x,CLinLSQRReport &rep);
static void LinLSQRSetXRep(CLinLSQRState &state,bool needxrep);
static void CLinLSQR::LinLSQRRestart(CLinLSQRState &state);
static int LinLSQRPeekIterationsCount(CLinLSQRState &s);
static void LinLSQRRequestTermination(CLinLSQRState &state);
private:
static void ClearrFields(CLinLSQRState &state);
};
//+------------------------------------------------------------------+
//| Constants |
//+------------------------------------------------------------------+
const double CLinLSQR::m_atol=1.0E-6;
const double CLinLSQR::m_btol=1.0E-6;
//+------------------------------------------------------------------+
//| This function initializes linear LSQR Solver. This solver is used|
//| to solve non-symmetric (and, possibly, non-square) problems. |
//| Least squares solution is returned for non - compatible systems. |
//| USAGE: |
//| 1. User initializes algorithm state with LinLSQRCreate() call |
//| 2. User tunes solver parameters with LinLSQRSetCond() and other|
//| functions |
//| 3. User calls LinLSQRSolveSparse() function which takes |
//| algorithm state and SparseMatrix object. |
//| 4. User calls LinLSQRResults() to get solution |
//| 5. Optionally, user may call LinLSQRSolveSparse() again to |
//| solve another problem with different matrix and/or right |
//| part without reinitializing LinLSQRState structure. |
//| INPUT PARAMETERS: |
//| M - number of rows in A |
//| N - number of variables, N > 0 |
//| OUTPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//| NOTE: see also LinLSQRCreateBuf() for version which reuses |
//| previously allocated place as much as possible. |
//+------------------------------------------------------------------+
void CLinLSQR::LinLSQRCreate(int m,int n,CLinLSQRState &state)
{
//--- check
if(!CAp::Assert(m>0,__FUNCTION__+": M<=0"))
return;
if(!CAp::Assert(n>0,__FUNCTION__+": N<=0"))
return;
//--- function call
LinLSQRCreateBuf(m,n,state);
}
//+------------------------------------------------------------------+
//| This function initializes linear LSQR Solver. It provides exactly|
//| same functionality as LinLSQRCreate(), but reuses previously |
//| allocated space as much as possible. |
//| INPUT PARAMETERS: |
//| M - number of rows in A |
//| N - number of variables, N > 0 |
//| OUTPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//+------------------------------------------------------------------+
void CLinLSQR::LinLSQRCreateBuf(int m,int n,CLinLSQRState &state)
{
//--- check
if(!CAp::Assert(m>0,__FUNCTION__+": M<=0"))
return;
if(!CAp::Assert(n>0,__FUNCTION__+": N<=0"))
return;
state.m_m=m;
state.m_n=n;
state.m_prectype=0;
state.m_epsa=m_atol;
state.m_epsb=m_btol;
state.m_epsc=1/MathSqrt(CMath::m_machineepsilon);
state.m_maxits=0;
state.m_lambdai=0;
state.m_xrep=false;
state.m_running=false;
state.m_repiterationscount=0;
//--- * allocate arrays
//--- * set RX to NAN (just for the case user calls Results() without
//--- calling SolveSparse()
//--- * set B to zero
CNormEstimator::NormEstimatorCreate(m,n,2,2,state.m_nes);
state.m_rx=vector<double>::Full(state.m_n,AL_NaN);
state.m_ui=vector<double>::Zeros(state.m_m+state.m_n);
state.m_uip1=vector<double>::Zeros(state.m_m+state.m_n);
state.m_vip1=vector<double>::Zeros(state.m_n);
state.m_vi=vector<double>::Zeros(state.m_n);
state.m_omegai=vector<double>::Zeros(state.m_n);
state.m_omegaip1=vector<double>::Zeros(state.m_n);
state.m_d=vector<double>::Zeros(state.m_n);
state.m_x=vector<double>::Zeros(state.m_m+state.m_n);
state.m_mv=vector<double>::Zeros(state.m_m+state.m_n);
state.m_mtv=vector<double>::Zeros(state.m_n);
state.m_b=vector<double>::Zeros(state.m_m);
state.m_rstate.ia.Resize(2);
state.m_rstate.ra=vector<double>::Zeros(1);
state.m_rstate.stage=-1;
}
//+------------------------------------------------------------------+
//| This function sets right part. By default, right part is zero. |
//| INPUT PARAMETERS: |
//| B - right part, array[N]. |
//| OUTPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//+------------------------------------------------------------------+
void CLinLSQR::LinLSQRSetB(CLinLSQRState &state,CRowDouble &b)
{
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": you can not change B when LinLSQRIteration is running"))
return;
if(!CAp::Assert(state.m_m<=CAp::Len(b),__FUNCTION__+": Length(B)<M"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(b,state.m_m),__FUNCTION__+": B contains infinite or NaN values"))
return;
state.m_b=b;
state.m_b.Resize(state.m_m);
state.m_bnorm2=state.m_b.Dot(state.m_b);
}
//+------------------------------------------------------------------+
//| This function changes preconditioning settings of |
//| LinLSQQSolveSparse() function. By default, SolveSparse() uses |
//| diagonal preconditioner, but if you want to use solver without |
//| preconditioning, you can call this function which forces solver |
//| to use unit matrix for preconditioning. |
//| INPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//+------------------------------------------------------------------+
void CLinLSQR::LinLSQRSetPrecUnit(CLinLSQRState &state)
{
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": you can not change preconditioner,because function LinLSQRIteration is running!"))
return;
state.m_prectype=-1;
}
//+------------------------------------------------------------------+
//| This function changes preconditioning settings of |
//| LinCGSolveSparse() function. LinCGSolveSparse() will use diagonal|
//| of the system matrix as preconditioner. This preconditioning mode|
//| is active by default. |
//| INPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//+------------------------------------------------------------------+
void CLinLSQR::LinLSQRSetPrecDiag(CLinLSQRState &state)
{
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": you can not change preconditioner,because function LinCGIteration is running!"))
return;
state.m_prectype=0;
}
//+------------------------------------------------------------------+
//| This function sets optional Tikhonov regularization coefficient. |
//| It is zero by default. |
//| INPUT PARAMETERS: |
//| LambdaI - regularization factor, LambdaI >= 0 |
//| OUTPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//+------------------------------------------------------------------+
void CLinLSQR::LinLSQRSetLambdaI(CLinLSQRState &state,double lambdai)
{
//--- check
if(!CAp::Assert(!state.m_running,"LinLSQRSetLambdaI: you can not set LambdaI,because function LinLSQRIteration is running"))
return;
if(!CAp::Assert(MathIsValidNumber(lambdai) && lambdai>=0.0,"LinLSQRSetLambdaI: LambdaI is infinite or NaN"))
return;
state.m_lambdai=lambdai;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
bool CLinLSQR::LinLSQRIteration(CLinLSQRState &state)
{
//--- create variables
int summn=0;
double bnorm=0;
int i=0;
int i_=0;
int label=-1;
//--- Reverse communication preparations
//--- This code initializes locals by:
//--- * random values determined during code
//--- generation - on first subroutine call
//--- * values from previous call - on subsequent calls
if(state.m_rstate.stage>=0)
{
summn=state.m_rstate.ia[0];
i=state.m_rstate.ia[1];
bnorm=state.m_rstate.ra[0];
}
else
{
summn=359;
i=-58;
bnorm=-919;
}
switch(state.m_rstate.stage)
{
case 0:
label=0;
break;
case 1:
label=1;
break;
case 2:
label=2;
break;
case 3:
label=3;
break;
case 4:
label=4;
break;
case 5:
label=5;
break;
case 6:
label=6;
break;
default:
//--- Routine body
//--- check
if(!CAp::Assert(CAp::Len(state.m_b)>0,__FUNCTION__+": using non-allocated array B"))
return(false);
summn=state.m_m+state.m_n;
bnorm=MathSqrt(state.m_bnorm2);
state.m_userterminationneeded=false;
state.m_running=true;
state.m_repnmv=0;
state.m_repiterationscount=0;
state.m_r2=state.m_bnorm2;
ClearrFields(state);
//estimate for ANorm
CNormEstimator::NormEstimatorRestart(state.m_nes);
label=7;
break;
}
//--- main loop
while(label>=0)
switch(label)
{
case 7:
if(!CNormEstimator::NormEstimatorIteration(state.m_nes))
{
label=8;
break;
}
if(!state.m_nes.m_NeedMv)
{
label=9;
break;
}
for(i_=0; i_<state.m_n; i_++)
state.m_x.Set(i_,state.m_nes.m_X[i_]);
state.m_repnmv++;
ClearrFields(state);
state.m_needmv=true;
state.m_rstate.stage=0;
label=-1;
break;
case 0:
state.m_needmv=false;
state.m_nes.m_Mv=state.m_mv;
label=7;
break;
case 9:
if(!state.m_nes.m_NeedMtv)
{
label=11;
break;
}
for(i_=0; i_<state.m_n; i_++)
state.m_x.Set(i_,state.m_nes.m_X[i_]);
//matrix-vector multiplication
state.m_repnmv++;
ClearrFields(state);
state.m_needmtv=true;
state.m_rstate.stage=1;
label=-1;
break;
case 1:
state.m_needmtv=false;
state.m_nes.m_Mtv=state.m_mtv;
label=7;
break;
case 11:
label=7;
break;
case 8:
CNormEstimator::NormEstimatorResults(state.m_nes,state.m_anorm);
//--- initialize .RX by zeros
state.m_rx.Fill(0);
//output first report
if(!state.m_xrep)
{
label=13;
break;
}
for(i_=0; i_<state.m_n; i_++)
state.m_x.Set(i_,state.m_rx[i_]);
ClearrFields(state);
state.m_xupdated=true;
state.m_rstate.stage=2;
label=-1;
break;
case 2:
state.m_xupdated=false;
case 13:
//--- LSQR, Step 0.
//--- Algorithm outline corresponds to one which was described at p.50 of
//--- "LSQR - an algorithm for sparse linear equations and sparse least
//--- squares" by C.Paige and M.Saunders with one small addition - we
//--- explicitly extend system matrix by additional N lines in order
//--- to handle non-zero lambda, i.e. original A is replaced by
//--- [ A ]
//--- A_mod = [ ]
//--- [ lambda*I ].
//--- Step 0:
//--- x[0] = 0
//--- beta[1]*u[1] = b
//--- alpha[1]*v[1] = A_mod'*u[1]
//--- w[1] = v[1]
//--- phiBar[1] = beta[1]
//--- rhoBar[1] = alpha[1]
//--- d[0] = 0
//--- NOTE:
//--- There are three criteria for stopping:
//--- (S0) maximum number of iterations
//--- (S1) ||Rk||<=EpsB*||B||;
//--- (S2) ||A^T*Rk||/(||A||*||Rk||)<=EpsA.
//--- It is very important that S2 always checked AFTER S1. It is necessary
//--- to avoid division by zero when Rk=0.
state.m_betai=bnorm;
if(state.m_betai==0.0)
{
//--- Zero right part
state.m_running=false;
state.m_repterminationtype=1;
return(false);
}
for(i=0; i<summn; i++)
{
if(i<state.m_m)
state.m_ui.Set(i,state.m_b[i]/state.m_betai);
else
state.m_ui.Set(i,0);
state.m_x.Set(i,state.m_ui[i]);
}
state.m_repnmv++;
ClearrFields(state);
state.m_needmtv=true;
state.m_rstate.stage=3;
label=-1;
break;
case 3:
state.m_needmtv=false;
for(i=0; i<state.m_n; i++)
state.m_mtv.Add(i,state.m_lambdai*state.m_ui[state.m_m+i]);
state.m_alphai=state.m_mtv.Dot(state.m_mtv);
state.m_alphai=MathSqrt(state.m_alphai);
if(state.m_alphai==0.0)
{
//--- Orthogonality stopping criterion is met
state.m_running=false;
state.m_repterminationtype=4;
return(false);
}
state.m_vi=state.m_mtv/state.m_alphai+0;
state.m_omegai=state.m_vi;
state.m_phibari=state.m_betai;
state.m_rhobari=state.m_alphai;
state.m_d.Fill(0);
state.m_dnorm=0;
//--- Steps I=1, 2, ...
case 15:
//--- At I-th step State.RepIterationsCount=I.
state.m_repiterationscount++;
//--- Bidiagonalization part:
//--- beta[i+1]*u[i+1] = A_mod*v[i]-alpha[i]*u[i]
//--- alpha[i+1]*v[i+1] = A_mod'*u[i+1] - beta[i+1]*v[i]
//--- NOTE: beta[i+1]=0 or alpha[i+1]=0 will lead to successful termination
//--- in the end of the current iteration. In this case u/v are zero.
//--- NOTE2: algorithm won't fail on zero alpha or beta (there will be no
//--- division by zero because it will be stopped BEFORE division
//--- occurs). However, near-zero alpha and beta won't stop algorithm
//--- and, although no division by zero will happen, orthogonality
//--- in U and V will be lost.
for(i_=0; i_<state.m_n; i_++)
state.m_x.Set(i_,state.m_vi[i_]);
state.m_repnmv=state.m_repnmv+1;
ClearrFields(state);
state.m_needmv=true;
state.m_rstate.stage=4;
label=-1;
break;
case 4:
state.m_needmv=false;
for(i=0; i<state.m_n; i++)
state.m_mv.Set(state.m_m+i,state.m_lambdai*state.m_vi[i]);
state.m_betaip1=0;
state.m_uip1=state.m_mv.ToVector()-state.m_ui*state.m_alphai;
state.m_betaip1=CAblasF::RDotV2(summn,state.m_uip1);
if(state.m_betaip1!=0.0)
{
state.m_betaip1=MathSqrt(state.m_betaip1);
state.m_uip1/=state.m_betaip1;
}
state.m_x=state.m_uip1;
state.m_repnmv++;
ClearrFields(state);
state.m_needmtv=true;
state.m_rstate.stage=5;
label=-1;
break;
case 5:
state.m_needmtv=false;
for(i=0; i<state.m_n; i++)
state.m_mtv.Add(i,state.m_lambdai*state.m_uip1[state.m_m+i]);
state.m_alphaip1=0;
state.m_vip1=state.m_mtv.ToVector()-state.m_vi*state.m_betaip1 ;
state.m_alphaip1=state.m_vip1.Dot(state.m_vip1);
if(state.m_alphaip1!=0.0)
{
state.m_alphaip1=MathSqrt(state.m_alphaip1);
state.m_vip1/=state.m_alphaip1;
}
//--- Build next orthogonal transformation
state.m_rhoi=CApServ::SafePythag2(state.m_rhobari,state.m_betaip1);
state.m_ci=state.m_rhobari/state.m_rhoi;
state.m_si=state.m_betaip1/state.m_rhoi;
state.m_theta=state.m_si*state.m_alphaip1;
state.m_rhobarip1=-(state.m_ci*state.m_alphaip1);
state.m_phii=state.m_ci*state.m_phibari;
state.m_phibarip1=state.m_si*state.m_phibari;
//--- Update .RNorm
//--- This tricky formula is necessary because simply writing
//--- State.R2:=State.PhiBarIP1*State.PhiBarIP1 does NOT guarantees
//--- monotonic decrease of R2. Roundoff error combined with 80-bit
//--- precision used internally by Intel chips allows R2 to increase
//--- slightly in some rare, but possible cases. This property is
//--- undesirable, so we prefer to guard against R increase.
state.m_r2=MathMin(state.m_r2,state.m_phibarip1*state.m_phibarip1);
//--- Update d and DNorm, check condition-related stopping criteria
state.m_d=(state.m_vi.ToVector()-state.m_d*state.m_theta)/state.m_rhoi;
state.m_dnorm=CAblasF::RDotV2(state.m_n,state.m_d);
if((MathSqrt(state.m_dnorm)*state.m_anorm)>=state.m_epsc)
{
state.m_running=false;
state.m_repterminationtype=7;
return(false);
}
//--- Update x, output report
state.m_rx+=state.m_omegai.ToVector()*state.m_phii/state.m_rhoi;
if(!state.m_xrep)
{
label=17;
break;
}
for(i_=0; i_<state.m_n; i_++)
state.m_x.Set(i_,state.m_rx[i_]);
ClearrFields(state);
state.m_xupdated=true;
state.m_rstate.stage=6;
label=-1;
break;
case 6:
state.m_xupdated=false;
case 17:
//--- Check stopping criteria
//--- 1. achieved required number of iterations;
//--- 2. ||Rk||<=EpsB*||B||;
//--- 3. ||A^T*Rk||/(||A||*||Rk||)<=EpsA;
if(state.m_maxits>0 && state.m_repiterationscount>=state.m_maxits)
{
//--- Achieved required number of iterations
state.m_running=false;
state.m_repterminationtype=5;
return(false);
}
if(state.m_phibarip1<=state.m_epsb*bnorm)
{
//--- ||Rk||<=EpsB*||B||, here ||Rk||=PhiBar
state.m_running=false;
state.m_repterminationtype=1;
return(false);
}
if((state.m_alphaip1*MathAbs(state.m_ci)/state.m_anorm)<=state.m_epsa)
{
//--- ||A^T*Rk||/(||A||*||Rk||)<=EpsA, here ||A^T*Rk||=PhiBar*Alpha[i+1]*|.C|
state.m_running=false;
state.m_repterminationtype=4;
return(false);
}
if(state.m_userterminationneeded)
{
//--- User requested termination
state.m_running=false;
state.m_repterminationtype=8;
return(false);
}
//--- Update omega
state.m_omegaip1=state.m_vip1.ToVector()-state.m_omegai*state.m_theta/state.m_rhoi;
//--- Prepare for the next iteration - rename variables:
//--- u[i] := u[i+1]
//--- v[i] := v[i+1]
//--- rho[i] := rho[i+1]
//--- ...
state.m_ui=state.m_uip1;
state.m_vi=state.m_vip1;
state.m_omegai=state.m_omegaip1;
state.m_alphai=state.m_alphaip1;
state.m_betai=state.m_betaip1;
state.m_phibari=state.m_phibarip1;
state.m_rhobari=state.m_rhobarip1;
label=15;
break;
case 16:
return(false);
}
//--- Saving state
state.m_rstate.ia.Set(0,summn);
state.m_rstate.ia.Set(1,i);
state.m_rstate.ra.Set(0,bnorm);
//--- return result
return(true);
}
//+------------------------------------------------------------------+
//| Procedure for solution of A*x = b with sparse A. |
//| INPUT PARAMETERS: |
//| State - algorithm state |
//| A - sparse M*N matrix in the CRS format (you MUST |
//| contvert it to CRS format by calling |
//| SparseConvertToCRS() function BEFORE you pass it|
//| to this function). |
//| B - right part, array[M] |
//| RESULT: |
//| This function returns no result. |
//| You can get solution by calling LinCGResults() |
//| NOTE: this function uses lightweight preconditioning - |
//| multiplication by inverse of diag(A). If you want, you can |
//| turn preconditioning off by calling LinLSQRSetPrecUnit(). |
//| However, preconditioning cost is low and preconditioner is |
//| very important for solution of badly scaled problems. |
//+------------------------------------------------------------------+
void CLinLSQR::LinLSQRSolveSparse(CLinLSQRState &state,
CSparseMatrix &a,CRowDouble &b)
{
//--- create variables
int n=state.m_n;
int i=0;
int j=0;
int t0=0;
int t1=0;
double v=0;
//--- check
if(!CAp::Assert(!state.m_running,"LinLSQRSolveSparse: you can not call this function when LinLSQRIteration is running"))
return;
if(!CAp::Assert(CAp::Len(b)>=state.m_m,"LinLSQRSolveSparse: Length(B)<M"))
return;
if(!CAp::Assert(CApServ::IsFiniteVector(b,state.m_m),"LinLSQRSolveSparse: B contains infinite or NaN values"))
return;
//--- Allocate temporaries
state.m_tmpd.Resize(n);
state.m_tmpx.Resize(n);
//--- Compute diagonal scaling matrix D
if(state.m_prectype==0)
{
//--- Default preconditioner - inverse of column norms
state.m_tmpd.Fill(0);
t0=0;
t1=0;
while(CSparse::SparseEnumerate(a,t0,t1,i,j,v))
state.m_tmpd.Add(j,CMath::Sqr(v));
for(i=0; i<n; i++)
{
if(state.m_tmpd[i]>0.0)
state.m_tmpd.Set(i,1/MathSqrt(state.m_tmpd[i]));
else
state.m_tmpd.Set(i,1);
}
}
else
{
//--- No diagonal scaling
state.m_tmpd.Fill(1);
}
//--- Solve.
//--- Instead of solving A*x=b we solve preconditioned system (A*D)*(inv(D)*x)=b.
//--- Transformed A is not calculated explicitly, we just modify multiplication
//--- by A or A'. After solution we modify State.RX so it will store untransformed
//--- variables
LinLSQRSetB(state,b);
LinLSQRRestart(state);
while(LinLSQRIteration(state))
{
if(state.m_needmv)
{
for(i=0; i<n; i++)
state.m_tmpx.Set(i,state.m_tmpd[i]*state.m_x[i]);
CSparse::SparseMV(a,state.m_tmpx,state.m_mv);
}
if(state.m_needmtv)
{
CSparse::SparseMTV(a,state.m_x,state.m_mtv);
state.m_mtv*=state.m_tmpd;
}
}
state.m_rx*=state.m_tmpd;
}
//+------------------------------------------------------------------+
//| This function sets stopping criteria. |
//| INPUT PARAMETERS: |
//| EpsA - algorithm will be stopped if |
//| || A^T * Rk || / ( || A || * || Rk ||) <= EpsA. |
//| EpsB - algorithm will be stopped if |
//| || Rk || <= EpsB * || B || |
//| MaxIts - algorithm will be stopped if number of |
//| iterations more than MaxIts. |
//| OUTPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//| NOTE: if EpsA, EpsB, EpsC and MaxIts are zero then these |
//| variables will be setted as default values. |
//+------------------------------------------------------------------+
void CLinLSQR::LinLSQRSetCond(CLinLSQRState &state,double epsa,
double epsb,int maxits)
{
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": you can not call this function when LinLSQRIteration is running"))
return;
if(!CAp::Assert(MathIsValidNumber(epsa) && epsa>=0.0,__FUNCTION__+": EpsA is negative,INF or NAN"))
return;
if(!CAp::Assert(MathIsValidNumber(epsb) && epsb>=0.0,__FUNCTION__+": EpsB is negative,INF or NAN"))
return;
if(!CAp::Assert(maxits>=0,__FUNCTION__+": MaxIts is negative"))
return;
if((epsa==0.0 && epsb==0.0) && maxits==0)
{
state.m_epsa=m_atol;
state.m_epsb=m_btol;
state.m_maxits=state.m_n;
}
else
{
state.m_epsa=epsa;
state.m_epsb=epsb;
state.m_maxits=maxits;
}
}
//+------------------------------------------------------------------+
//| LSQR solver: results. |
//| This function must be called after LinLSQRSolve |
//| INPUT PARAMETERS: |
//| State - algorithm state |
//| OUTPUT PARAMETERS: |
//| X - array[N], solution |
//| Rep - optimization report: |
//| * Rep.TerminationType completetion code: |
//| * 1 || Rk || <= EpsB* || B || |
//| * 4 ||A^T * Rk|| / (||A|| * ||Rk||) <= EpsA |
//| * 5 MaxIts steps was taken |
//| * 7 rounding errors prevent further progress, |
//| X contains best point found so far. |
//| (sometimes returned on singular systems) |
//| * 8 user requested termination via calling |
//| LinLSQRRequestTermination() |
//| * Rep.IterationsCount contains iterations count |
//| * NMV countains number of matrix - vector |
//| calculations |
//+------------------------------------------------------------------+
void CLinLSQR::LinLSQRResults(CLinLSQRState &state,CRowDouble &x,
CLinLSQRReport &rep)
{
x.Resize(0);
//--- check
if(!CAp::Assert(!state.m_running,__FUNCTION__+": you can not call this function when LinLSQRIteration is running"))
return;
x=state.m_rx;
x.Resize(state.m_n);
rep.m_iterationscount=state.m_repiterationscount;
rep.m_nmv=state.m_repnmv;
rep.m_terminationtype=state.m_repterminationtype;
}
//+------------------------------------------------------------------+
//| This function turns on / off reporting. |
//| INPUT PARAMETERS: |
//| State - structure which stores algorithm state |
//| NeedXRep - whether iteration reports are needed or not |
//| If NeedXRep is True, algorithm will call rep() |
//| callback function if it is provided to |
//| MinCGOptimize(). |
//+------------------------------------------------------------------+
void CLinLSQR::LinLSQRSetXRep(CLinLSQRState &state,bool needxrep)
{
state.m_xrep=needxrep;
}
//+------------------------------------------------------------------+
//| This function restarts LinLSQRIteration |
//+------------------------------------------------------------------+
void CLinLSQR::LinLSQRRestart(CLinLSQRState &state)
{
state.m_rstate.ia.Resize(2);
state.m_rstate.ra.Resize(1);
state.m_rstate.stage=-1;
ClearrFields(state);
state.m_repiterationscount=0;
}
//+------------------------------------------------------------------+
//| This function is used to peek into LSQR solver and get current |
//| iteration counter. You can safely "peek" into the solver from |
//| another thread. |
//| INPUT PARAMETERS: |
//| S - solver object |
//| RESULT: |
//| iteration counter, in [0, INF) |
//+------------------------------------------------------------------+
int CLinLSQR::LinLSQRPeekIterationsCount(CLinLSQRState &s)
{
int result=s.m_repiterationscount;
return(result);
}
//+------------------------------------------------------------------+
//| This subroutine submits request for termination of the running |
//| solver. It can be called from some other thread which wants LSQR |
//| solver to terminate (obviously, the thread running LSQR solver |
//| can not request termination because it is already busy working on|
//| LSQR). |
//| As result, solver stops at point which was "current accepted" |
//| when termination request was submitted and returns error code 8 |
//| (successful termination). Such termination is a smooth process |
//| which properly deallocates all temporaries. |
//| INPUT PARAMETERS: |
//| State - solver structure |
//| NOTE: calling this function on solver which is NOT running will |
//| have no effect. |
//| NOTE: multiple calls to this function are possible. First call is|
//| counted, subsequent calls are silently ignored. |
//| NOTE: solver clears termination flag on its start, it means that |
//| if some other thread will request termination too soon, its|
//| request will went unnoticed. |
//+------------------------------------------------------------------+
void CLinLSQR::LinLSQRRequestTermination(CLinLSQRState &state)
{
state.m_userterminationneeded=true;
}
//+------------------------------------------------------------------+
//| Clears request fileds (to be sure that we don't forgot to clear |
//| something) |
//+------------------------------------------------------------------+
void CLinLSQR::ClearrFields(CLinLSQRState &state)
{
state.m_xupdated=false;
state.m_needmv=false;
state.m_needmtv=false;
state.m_needmv2=false;
state.m_needvmv=false;
state.m_needprec=false;
}
//+------------------------------------------------------------------+
Reference in a new issue View git blame Copy permalink