//+------------------------------------------------------------------+ //| 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_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; in-1 || p[i]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_r10.0)) return; //--- prepare: check inputs,allocate space... if(n<=0 || m<=0) { info=-1; return; } for(i=0; in-1 || p[i]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=0; i--) { //--- calculation for(i_=i+1; i_=0; i--) { //--- check if(i0) { 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=0; i--) { for(i_=i+1; i_=0; i--) { //--- check if(i0) { 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)=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,__FUNCTION__+": Length(X)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)=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_=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_ 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)0) { c=matrix::Zeros(ne,ne); c.Set(0,ne-1,-a[0]); for(i=1; i 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)::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)::Zeros(n); return; } CAblasF::RCopyAllocV(n,b,x); for(int i=0; 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)::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)::Zeros(n); return; } x=b; x.Resize(n); for(int i=0; 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,__FUNCTION__+": length(P)=n,__FUNCTION__+": length(Q)=i && p[i]=i && q[i]::Zeros(n); return; } } x=b; x.Resize(n); for(i=0; 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)=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)=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)=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,__FUNCTION__+": length(B) 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=(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::Full(state.m_n,AL_NaN); state.m_startx=vector::Zeros(state.m_n); state.m_b=vector::Zeros(state.m_n); state.m_cx=vector::Zeros(state.m_n); state.m_p=vector::Zeros(state.m_n); state.m_r=vector::Zeros(state.m_n); state.m_cr=vector::Zeros(state.m_n); state.m_z=vector::Zeros(state.m_n); state.m_cz=vector::Zeros(state.m_n); state.m_x=vector::Zeros(state.m_n); state.m_mv=vector::Zeros(state.m_n); state.m_pv=vector::Zeros(state.m_n); UpdateItersData(state); state.m_rstate.ia.Resize(1); state.m_rstate.ra=vector::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)=state.m_n,__FUNCTION__+": Length(B)=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_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)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)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::Full(state.m_n,AL_NaN); state.m_ui=vector::Zeros(state.m_m+state.m_n); state.m_uip1=vector::Zeros(state.m_m+state.m_n); state.m_vip1=vector::Zeros(state.m_n); state.m_vi=vector::Zeros(state.m_n); state.m_omegai=vector::Zeros(state.m_n); state.m_omegaip1=vector::Zeros(state.m_n); state.m_d=vector::Zeros(state.m_n); state.m_x=vector::Zeros(state.m_m+state.m_n); state.m_mv=vector::Zeros(state.m_m+state.m_n); state.m_mtv=vector::Zeros(state.m_n); state.m_b=vector::Zeros(state.m_m); state.m_rstate.ia.Resize(2); state.m_rstate.ra=vector::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)= 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_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_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)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=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; } //+------------------------------------------------------------------+