//+------------------------------------------------------------------+ //| linalg.mqh | //| Copyright 2003-2012 Sergey Bochkanov (ALGLIB project) | //| Copyright 2012-2017, MetaQuotes Software Corp. | //| 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 "alglibinternal.mqh" #include "alglibmisc.mqh" //+------------------------------------------------------------------+ //| Work with matrix forms | //+------------------------------------------------------------------+ class CAblas { private: //--- split lenght static void AblasInternalSplitLength(const int n,const int nb,int &n1,int &n2); //--- real numbers static void RMatrixSyrk2(const int n,const int k,const double alpha,const CMatrixDouble &a,const int ia,const int ja,const int optypea,const double beta,CMatrixDouble &c,const int ic,const int jc,const bool isUpper); static void RMatrixGemmK(const int m,const int n,const int k,const double alpha,const CMatrixDouble &a,const int ia,const int ja,const int optypea,const CMatrixDouble &b,const int ib,const int jb,const int optypeb,const double beta,CMatrixDouble &c,const int ic,const int jc); static void RMatrixRightTrsM2(const int m,const int n,CMatrixDouble &a,const int i1,const int j1,const bool isUpper,const bool isUnit,const int optype,CMatrixDouble &x,const int i2,const int j2); static void RMatrixLeftTrsM2(const int m,const int n,CMatrixDouble &a,const int i1,const int j1,const bool isUpper,const bool isUnit,const int optype,CMatrixDouble &x,const int i2,const int j2); //--- complex numbers static void CMatrixSyrk2(const int n,const int k,const double alpha,const CMatrixComplex &a,const int ia,const int ja,const int optypea,const double beta,CMatrixComplex &c,const int ic,const int jc,const bool isUpper); static void CMatrixGemmk(const int m,const int n,const int k,complex &alpha,const CMatrixComplex &a,const int ia,const int ja,const int optypea,const CMatrixComplex &b,const int ib,const int jb,const int optypeb,complex &beta,CMatrixComplex &c,const int ic,const int jc); static void CMatrixRightTrsM2(const int m,const int n,CMatrixComplex &a,const int i1,const int j1,const bool isUpper,const bool isUnit,const int optype,CMatrixComplex &x,const int i2,const int j2); static void CMatrixLeftTrsM2(const int m,const int n,CMatrixComplex &a,const int i1,const int j1,const bool isUpper,const bool isUnit,const int optype,CMatrixComplex &x,const int i2,const int j2); public: CAblas(void); ~CAblas(void); //--- size static int AblasBlockSize(void) { return(32);} static int AblasMicroBlockSize(void) { return(8); } static int AblasComplexBlockSize(void) { return(24);} //--- split lenght static void AblasSplitLength(const CMatrixDouble &a,const int n,int &n1,int &n2); static void AblasComplexSplitLength(const CMatrixComplex &a,const int n,int &n1,int &n2); //--- real numbers static void RMatrixSyrk(const int n,const int k,const double alpha,const CMatrixDouble &a,const int ia,const int ja,const int optypea,const double beta,CMatrixDouble &c,const int ic,const int jc,const bool isUpper); static void RMatrixGemm(const int m,const int n,const int k,const double alpha,const CMatrixDouble &a,const int ia,const int ja,const int optypea,const CMatrixDouble &b,const int ib,const int jb,const int optypeb,const double beta,CMatrixDouble &c,const int ic,const int jc); static void RMatrixTranspose(const int m,const int n,const CMatrixDouble &a,const int ia,const int ja,CMatrixDouble &b,const int ib,const int jb); static void RMatrixCopy(const int m,const int n,const CMatrixDouble &a,const int ia,const int ja,CMatrixDouble &b,const int ib,const int jb); static void RMatrixRank1(const int m,const int n,CMatrixDouble &a,const int ia,const int ja,const double &u[],const int iu,const double &v[],const int iv); static void RMatrixMVect(const int m,const int n,const CMatrixDouble &a,const int ia,const int ja,const int opa,const double &x[],const int ix,double &y[],const int iy); static void RMatrixRightTrsM(const int m,const int n,CMatrixDouble &a,const int i1,const int j1,const bool isUpper,const bool isUnit,const int optype,CMatrixDouble &x,const int i2,const int j2); static void RMatrixLeftTrsM(const int m,const int n,CMatrixDouble &a,const int i1,const int j1,const bool isUpper,const bool isUnit,const int optype,CMatrixDouble &x,const int i2,const int j2); //--- complex numbers static void CMatrixSyrk(const int n,const int k,const double alpha,CMatrixComplex &a,const int ia,const int ja,const int optypea,const double beta,CMatrixComplex &c,const int ic,const int jc,const bool isUpper); static void CMatrixGemm(const int m,const int n,const int k,complex &alpha,CMatrixComplex &a,const int ia,const int ja,const int optypea,CMatrixComplex &b,const int ib,const int jb,const int optypeb,complex &beta,CMatrixComplex &c,const int ic,const int jc); static void CMatrixTranspose(const int m,const int n,const CMatrixComplex &a,const int ia,const int ja,CMatrixComplex &b,const int ib,const int jb); static void CMatrixCopy(const int m,const int n,const CMatrixComplex &a,const int ia,const int ja,CMatrixComplex &b,const int ib,const int jb); static void CMatrixRank1(const int m,const int n,CMatrixComplex &a,const int ia,const int ja,const complex &u[],const int iu,const complex &v[],const int iv); static void CMatrixMVect(const int m,const int n,const CMatrixComplex &a,const int ia,const int ja,const int opa,const complex &x[],const int ix,complex &y[],const int iy); static void CMatrixRightTrsM(const int m,const int n,CMatrixComplex &a,const int i1,const int j1,const bool isUpper,const bool isUnit,const int optype,CMatrixComplex &x,const int i2,const int j2); static void CMatrixLeftTrsM(const int m,const int n,CMatrixComplex &a,const int i1,const int j1,const bool isUpper,const bool isUnit,const int optype,CMatrixComplex &x,const int i2,const int j2); }; //+------------------------------------------------------------------+ //| Constructor without parameters | //+------------------------------------------------------------------+ CAblas::CAblas(void) { } //+------------------------------------------------------------------+ //| Destructor | //+------------------------------------------------------------------+ CAblas::~CAblas(void) { } //+------------------------------------------------------------------+ //| Same as CMatrixSYRK, but for real matrices | //| OpType may be only 0 or 1. | //+------------------------------------------------------------------+ static void CAblas::RMatrixSyrk(const int n,const int k,const double alpha, const CMatrixDouble &a,const int ia,const int ja, const int optypea,const double beta,CMatrixDouble &c, const int ic,const int jc,const bool isUpper) { //--- create variables int s1=0; int s2=0; int bs=0; //--- calculation size bs=AblasBlockSize(); //--- check if(n<=bs && k<=bs) { RMatrixSyrk2(n,k,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); return; } //--- check if(k>=n) { //--- Split K AblasSplitLength(a,k,s1,s2); //--- check if(optypea==0) { RMatrixSyrk(n,s1,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); RMatrixSyrk(n,s2,alpha,a,ia,ja+s1,optypea,1.0,c,ic,jc,isUpper); } else { RMatrixSyrk(n,s1,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); RMatrixSyrk(n,s2,alpha,a,ia+s1,ja,optypea,1.0,c,ic,jc,isUpper); } } else { //--- Split N AblasSplitLength(a,n,s1,s2); //--- check if(optypea==0 && isUpper) { RMatrixSyrk(s1,k,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); RMatrixGemm(s1,s2,k,alpha,a,ia,ja,0,a,ia+s1,ja,1,beta,c,ic,jc+s1); RMatrixSyrk(s2,k,alpha,a,ia+s1,ja,optypea,beta,c,ic+s1,jc+s1,isUpper); //--- exit the function return; } //--- check if(optypea==0 && !isUpper) { RMatrixSyrk(s1,k,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); RMatrixGemm(s2,s1,k,alpha,a,ia+s1,ja,0,a,ia,ja,1,beta,c,ic+s1,jc); RMatrixSyrk(s2,k,alpha,a,ia+s1,ja,optypea,beta,c,ic+s1,jc+s1,isUpper); //--- exit the function return; } //--- check if(optypea!=0 && isUpper) { RMatrixSyrk(s1,k,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); RMatrixGemm(s1,s2,k,alpha,a,ia,ja,1,a,ia,ja+s1,0,beta,c,ic,jc+s1); RMatrixSyrk(s2,k,alpha,a,ia,ja+s1,optypea,beta,c,ic+s1,jc+s1,isUpper); //--- exit the function return; } //--- check if(optypea!=0 && !isUpper) { RMatrixSyrk(s1,k,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); RMatrixGemm(s2,s1,k,alpha,a,ia,ja+s1,1,a,ia,ja,0,beta,c,ic+s1,jc); RMatrixSyrk(s2,k,alpha,a,ia,ja+s1,optypea,beta,c,ic+s1,jc+s1,isUpper); //--- exit the function return; } } //--- exit the function return; } //+------------------------------------------------------------------+ //| Splits matrix length in two parts, left part should match ABLAS | //| block size | //| INPUT PARAMETERS | //| A - real matrix, is passed to ensure that we didn't split| //| complex matrix using real splitting subroutine. | //| matrix itself is not changed. | //| N - length, N>0 | //| OUTPUT PARAMETERS | //| N1 - length | //| N2 - length | //| N1+N2=N, N1>=N2, N2 may be zero | //+------------------------------------------------------------------+ static void CAblas::AblasSplitLength(const CMatrixDouble &a,const int n, int &n1,int &n2) { //--- initialization n1=0; n2=0; //--- check if(n>AblasBlockSize()) AblasInternalSplitLength(n,AblasBlockSize(),n1,n2); else AblasInternalSplitLength(n,AblasMicroBlockSize(),n1,n2); //--- exit the function return; } //+------------------------------------------------------------------+ //| Complex ABLASSplitLength | //+------------------------------------------------------------------+ static void CAblas::AblasComplexSplitLength(const CMatrixComplex &a,const int n, int &n1,int &n2) { //--- check if(n>AblasComplexBlockSize()) AblasInternalSplitLength(n,AblasComplexBlockSize(),n1,n2); else AblasInternalSplitLength(n,AblasMicroBlockSize(),n1,n2); } //+------------------------------------------------------------------+ //| Same as CMatrixGEMM, but for real numbers. | //| OpType may be only 0 or 1. | //+------------------------------------------------------------------+ static void CAblas::RMatrixGemm(const int m,const int n,const int k,const double alpha, const CMatrixDouble &a,const int ia,const int ja, const int optypea,const CMatrixDouble &b, const int ib,const int jb,const int optypeb, const double beta,CMatrixDouble &c,const int ic, const int jc) { //--- create variables int s1=0; int s2=0; int bs; //--- calculation size bs=AblasBlockSize(); //--- check if(m<=bs && n<=bs && k<=bs) { RMatrixGemmK(m,n,k,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); return; } //--- check if(m>=n && m>=k) { //--- A*B = (A1 A2)^T*B AblasSplitLength(a,m,s1,s2); //--- check if(optypea==0) { RMatrixGemm(s1,n,k,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); RMatrixGemm(s2,n,k,alpha,a,ia+s1,ja,optypea,b,ib,jb,optypeb,beta,c,ic+s1,jc); } else { RMatrixGemm(s1,n,k,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); RMatrixGemm(s2,n,k,alpha,a,ia,ja+s1,optypea,b,ib,jb,optypeb,beta,c,ic+s1,jc); } return; } //--- check if(n>=m && n>=k) { //--- A*B = A*(B1 B2) AblasSplitLength(a,n,s1,s2); //--- check if(optypeb==0) { RMatrixGemm(m,s1,k,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); RMatrixGemm(m,s2,k,alpha,a,ia,ja,optypea,b,ib,jb+s1,optypeb,beta,c,ic,jc+s1); } else { RMatrixGemm(m,s1,k,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); RMatrixGemm(m,s2,k,alpha,a,ia,ja,optypea,b,ib+s1,jb,optypeb,beta,c,ic,jc+s1); } return; } //--- check if(k>=m && k>=n) { //--- A*B = (A1 A2)*(B1 B2)^T AblasSplitLength(a,k,s1,s2); //--- check if(optypea==0 && optypeb==0) { RMatrixGemm(m,n,s1,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); RMatrixGemm(m,n,s2,alpha,a,ia,ja+s1,optypea,b,ib+s1,jb,optypeb,1.0,c,ic,jc); } //--- check if(optypea==0 && optypeb!=0) { RMatrixGemm(m,n,s1,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); RMatrixGemm(m,n,s2,alpha,a,ia,ja+s1,optypea,b,ib,jb+s1,optypeb,1.0,c,ic,jc); } //--- check if(optypea!=0 && optypeb==0) { RMatrixGemm(m,n,s1,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); RMatrixGemm(m,n,s2,alpha,a,ia+s1,ja,optypea,b,ib+s1,jb,optypeb,1.0,c,ic,jc); } //--- check if(optypea!=0 && optypeb!=0) { RMatrixGemm(m,n,s1,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); RMatrixGemm(m,n,s2,alpha,a,ia+s1,ja,optypea,b,ib,jb+s1,optypeb,1.0,c,ic,jc); } } //--- exit the function return; } //+------------------------------------------------------------------+ //| Level 2 subrotuine | //+------------------------------------------------------------------+ static void CAblas::RMatrixSyrk2(const int n,const int k,const double alpha, const CMatrixDouble &a,const int ia,const int ja, const int optypea,const double beta,CMatrixDouble &c, const int ic,const int jc,const bool isUpper) { //--- check if((alpha==0.0 || k==0.0) && beta==1.0) return; //--- create variables int i=0; int j=0; int j1=0; int j2=0; double v=0; int i_=0; int i1_=0; //--- check if(optypea==0) { //--- C=alpha*A*A^H+beta*C for(i=0;i0) { v=0.0; for(i_=ja;i_n) { //--- split AblasComplexSplitLength(a,m,s1,s2); //--- function call CMatrixTranspose(s1,n,a,ia,ja,b,ib,jb); CMatrixTranspose(s2,n,a,ia+s1,ja,b,ib,jb+s1); } else { //--- split AblasComplexSplitLength(a,n,s1,s2); //--- function call CMatrixTranspose(m,s1,a,ia,ja,b,ib,jb); CMatrixTranspose(m,s2,a,ia,ja+s1,b,ib+s1,jb); } } } //+------------------------------------------------------------------+ //| Cache-oblivous real "copy-and-transpose" | //| Input parameters: | //| M - number of rows | //| N - number of columns | //| A - source matrix, MxN submatrix is copied and transposed| //| IA - submatrix offset (row index) | //| JA - submatrix offset (column index) | //| A - destination matrix | //| IB - submatrix offset (row index) | //| JB - submatrix offset (column index) | //+------------------------------------------------------------------+ static void CAblas::RMatrixTranspose(const int m,const int n,const CMatrixDouble &a, const int ia,const int ja,CMatrixDouble &b, const int ib,const int jb) { //--- create variables int i=0; int s1=0; int s2=0; int i_=0; int i1_=0; //--- check if(m<=2*AblasBlockSize() && n<=2*AblasBlockSize()) { //--- base case for(i=0;in) { //--- split AblasSplitLength(a,m,s1,s2); //--- function call RMatrixTranspose(s1,n,a,ia,ja,b,ib,jb); RMatrixTranspose(s2,n,a,ia+s1,ja,b,ib,jb+s1); } else { //--- split AblasSplitLength(a,n,s1,s2); //--- function call RMatrixTranspose(m,s1,a,ia,ja,b,ib,jb); RMatrixTranspose(m,s2,a,ia,ja+s1,b,ib+s1,jb); } } } //+------------------------------------------------------------------+ //| Copy | //| Input parameters: | //| M - number of rows | //| N - number of columns | //| A - source matrix, MxN submatrix is copied and transposed| //| IA - submatrix offset (row index) | //| JA - submatrix offset (column index) | //| B - destination matrix | //| IB - submatrix offset (row index) | //| JB - submatrix offset (column index) | //+------------------------------------------------------------------+ static void CAblas::CMatrixCopy(const int m,const int n,const CMatrixComplex &a, const int ia,const int ja,CMatrixComplex &b, const int ib,const int jb) { //--- create variables int i=0; int i_=0; int i1_=0; //--- copy for(i=0;i=0 | //| N - number of columns of op(A) | //| N>=0 | //| A - target matrix | //| IA - submatrix offset (row index) | //| JA - submatrix offset (column index) | //| OpA - operation type: | //| * OpA=0 => op(A) = A | //| * OpA=1 => op(A) = A^T | //| * OpA=2 => op(A) = A^H | //| X - input vector | //| IX - subvector offset | //| IY - subvector offset | //| OUTPUT PARAMETERS: | //| Y - vector which stores result | //| if M=0, then subroutine does nothing. | //| if N=0, Y is filled by zeros. | //+------------------------------------------------------------------+ static void CAblas::CMatrixMVect(const int m,const int n,const CMatrixComplex &a, const int ia,const int ja,const int opa, const complex &x[],const int ix,complex &y[], const int iy) { //--- create variables int i=0; complex v=0; int i_=0; int i1_=0; //--- check if(m==0) return; //--- check if(n==0) { for(i=0;i op(A) = A | //| * OpA=1 => op(A) = A^T | //| X - input vector | //| IX - subvector offset | //| IY - subvector offset | //| OUTPUT PARAMETERS: | //| Y - vector which stores result | //| if M=0, then subroutine does nothing. | //| if N=0, Y is filled by zeros. | //+------------------------------------------------------------------+ static void CAblas::RMatrixMVect(const int m,const int n,const CMatrixDouble &a, const int ia,const int ja,const int opa, const double &x[],const int ix,double &y[], const int iy) { //--- create variables int i=0; double v=0; int i_=0; int i1_=0; //--- check if(m==0) return; //--- check if(n==0) { for(i=0;i=0 | //| M - matrix size, N>=0 | //| A - matrix, actial matrix is stored in | //| A[I1:I1+N-1,J1:J1+N-1] | //| I1 - submatrix offset | //| J1 - submatrix offset | //| IsUpper - whether matrix is upper triangular | //| IsUnit - whether matrix is unitriangular | //| OpType - transformation type: | //| * 0 - no transformation | //| * 1 - transposition | //| * 2 - conjugate transposition | //| C - matrix, actial matrix is stored in | //| C[I2:I2+M-1,J2:J2+N-1] | //| I2 - submatrix offset | //| J2 - submatrix offset | //+------------------------------------------------------------------+ static void CAblas::CMatrixRightTrsM(const int m,const int n,CMatrixComplex &a, const int i1,const int j1,const bool isUpper, const bool isUnit,const int optype, CMatrixComplex &x,const int i2,const int j2) { //--- create variables complex Alpha(-1,0); complex Beta(1,0); int s1=0; int s2=0; int bs=AblasComplexBlockSize(); //--- check if(m<=bs && n<=bs) { //--- basic algorithm CMatrixRightTrsM2(m,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); //--- exit the function return; } //--- check if(m>=n) { //--- Split X: X*A = (X1 X2)^T*A AblasComplexSplitLength(a,m,s1,s2); CMatrixRightTrsM(s1,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); CMatrixRightTrsM(s2,n,a,i1,j1,isUpper,isUnit,optype,x,i2+s1,j2); } else { //--- Split A: //--- (A1 A12) //--- X*op(A) = X*op( ) //--- ( A2) //--- //--- Different variants depending on //--- IsUpper/OpType combinations AblasComplexSplitLength(a,n,s1,s2); //--- check if(isUpper && optype==0) { //--- (A1 A12)-1 //--- X*A^-1 = (X1 X2)*( ) //--- ( A2) CMatrixRightTrsM(m,s1,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); CMatrixGemm(m,s2,s1,Alpha,x,i2,j2,0,a,i1,j1+s1,0,Beta,x,i2,j2+s1); CMatrixRightTrsM(m,s2,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2,j2+s1); //--- exit the function return; } //--- check if(isUpper && optype!=0) { //--- (A1' )-1 //--- X*A^-1 = (X1 X2)*( ) //--- (A12' A2') CMatrixRightTrsM(m,s2,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2,j2+s1); CMatrixGemm(m,s1,s2,Alpha,x,i2,j2+s1,0,a,i1,j1+s1,optype,Beta,x,i2,j2); CMatrixRightTrsM(m,s1,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); //--- exit the function return; } //--- check if(!isUpper && optype==0) { //--- (A1 )-1 //--- X*A^-1 = (X1 X2)*( ) //--- (A21 A2) CMatrixRightTrsM(m,s2,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2,j2+s1); CMatrixGemm(m,s1,s2,Alpha,x,i2,j2+s1,0,a,i1+s1,j1,0,Beta,x,i2,j2); CMatrixRightTrsM(m,s1,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); //--- exit the function return; } //--- check if(!isUpper && optype!=0) { //--- (A1' A21')-1 //--- X*A^-1 = (X1 X2)*( ) //--- ( A2') CMatrixRightTrsM(m,s1,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); CMatrixGemm(m,s2,s1,Alpha,x,i2,j2,0,a,i1+s1,j1,optype,Beta,x,i2,j2+s1); CMatrixRightTrsM(m,s2,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2,j2+s1); } } //--- exit the function return; } //+------------------------------------------------------------------+ //| This subroutine calculates op(A^-1)*X where: | //| * X is MxN general matrix | //| * A is MxM upper/lower triangular/unitriangular matrix | //| * "op" may be identity transformation, transposition, conjugate | //| transposition | //| Multiplication result replaces X. | //| Cache-oblivious algorithm is used. | //| INPUT PARAMETERS | //| N - matrix size, N>=0 | //| M - matrix size, N>=0 | //| A - matrix, actial matrix is stored in | //| A[I1:I1+M-1,J1:J1+M-1] | //| I1 - submatrix offset | //| J1 - submatrix offset | //| IsUpper - whether matrix is upper triangular | //| IsUnit - whether matrix is unitriangular | //| OpType - transformation type: | //| * 0 - no transformation | //| * 1 - transposition | //| * 2 - conjugate transposition | //| C - matrix, actial matrix is stored in | //| C[I2:I2+M-1,J2:J2+N-1] | //| I2 - submatrix offset | //| J2 - submatrix offset | //+------------------------------------------------------------------+ static void CAblas::CMatrixLeftTrsM(const int m,const int n,CMatrixComplex &a, const int i1,const int j1,const bool isUpper, const bool isUnit,const int optype, CMatrixComplex &x,const int i2,const int j2) { //--- create variables complex Alpha(-1,0); complex Beta(1,0); int s1=0; int s2=0; int bs=AblasComplexBlockSize(); //--- check if(m<=bs && n<=bs) { //--- basic algorithm CMatrixLeftTrsM2(m,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); //--- exit the function return; } //--- check if(n>=m) { //--- Split X: op(A)^-1*X = op(A)^-1*(X1 X2) AblasComplexSplitLength(x,n,s1,s2); CMatrixLeftTrsM(m,s1,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); CMatrixLeftTrsM(m,s2,a,i1,j1,isUpper,isUnit,optype,x,i2,j2+s1); } else { //--- Split A AblasComplexSplitLength(a,m,s1,s2); //--- check if(isUpper && optype==0) { //--- (A1 A12)-1 ( X1 ) //--- A^-1*X* = ( ) *( ) //--- ( A2) ( X2 ) CMatrixLeftTrsM(s2,n,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2+s1,j2); CMatrixGemm(s1,n,s2,Alpha,a,i1,j1+s1,0,x,i2+s1,j2,0,Beta,x,i2,j2); CMatrixLeftTrsM(s1,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); //--- exit the function return; } //--- check if(isUpper && optype!=0) { //--- (A1' )-1 ( X1 ) //--- A^-1*X = ( ) *( ) //--- (A12' A2') ( X2 ) CMatrixLeftTrsM(s1,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); CMatrixGemm(s2,n,s1,Alpha,a,i1,j1+s1,optype,x,i2,j2,0,Beta,x,i2+s1,j2); CMatrixLeftTrsM(s2,n,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2+s1,j2); //--- exit the function return; } //--- check if(!isUpper && optype==0) { //--- (A1 )-1 ( X1 ) //--- A^-1*X = ( ) *( ) //--- (A21 A2) ( X2 ) CMatrixLeftTrsM(s1,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); CMatrixGemm(s2,n,s1,Alpha,a,i1+s1,j1,0,x,i2,j2,0,Beta,x,i2+s1,j2); CMatrixLeftTrsM(s2,n,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2+s1,j2); //--- exit the function return; } //--- check if(!isUpper && optype!=0) { //--- (A1' A21')-1 ( X1 ) //--- A^-1*X = ( ) *( ) //--- ( A2') ( X2 ) CMatrixLeftTrsM(s2,n,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2+s1,j2); CMatrixGemm(s1,n,s2,Alpha,a,i1+s1,j1,optype,x,i2+s1,j2,0,Beta,x,i2,j2); CMatrixLeftTrsM(s1,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); } } //--- exit the function return; } //+------------------------------------------------------------------+ //| Same as CMatrixRightTRSM, but for real matrices | //| OpType may be only 0 or 1. | //+------------------------------------------------------------------+ static void CAblas::RMatrixRightTrsM(const int m,const int n,CMatrixDouble &a, const int i1,const int j1,const bool isUpper, const bool isUnit,const int optype, CMatrixDouble &x,const int i2,const int j2) { //--- create variables int s1=0; int s2=0; int bs=AblasBlockSize(); //--- check if(m<=bs && n<=bs) { //--- basic algorithm RMatrixRightTrsM2(m,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); //--- exit the function return; } //--- check if(m>=n) { //--- Split X: X*A = (X1 X2)^T*A AblasSplitLength(a,m,s1,s2); RMatrixRightTrsM(s1,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); RMatrixRightTrsM(s2,n,a,i1,j1,isUpper,isUnit,optype,x,i2+s1,j2); } else { //--- Split A: //--- (A1 A12) //--- X*op(A) = X*op( ) //--- ( A2) //--- Different variants depending on //--- IsUpper/OpType combinations AblasSplitLength(a,n,s1,s2); //--- check if(isUpper && optype==0) { //--- (A1 A12)-1 //--- X*A^-1 = (X1 X2)*( ) //--- ( A2) RMatrixRightTrsM(m,s1,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); RMatrixGemm(m,s2,s1,-1.0,x,i2,j2,0,a,i1,j1+s1,0,1.0,x,i2,j2+s1); RMatrixRightTrsM(m,s2,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2,j2+s1); //--- exit the function return; } //--- check if(isUpper && optype!=0) { //--- (A1' )-1 //--- X*A^-1 = (X1 X2)*( ) //--- (A12' A2') RMatrixRightTrsM(m,s2,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2,j2+s1); RMatrixGemm(m,s1,s2,-1.0,x,i2,j2+s1,0,a,i1,j1+s1,optype,1.0,x,i2,j2); RMatrixRightTrsM(m,s1,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); //--- exit the function return; } //--- check if(!isUpper && optype==0) { //--- (A1 )-1 //--- X*A^-1 = (X1 X2)*( ) //--- (A21 A2) RMatrixRightTrsM(m,s2,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2,j2+s1); RMatrixGemm(m,s1,s2,-1.0,x,i2,j2+s1,0,a,i1+s1,j1,0,1.0,x,i2,j2); RMatrixRightTrsM(m,s1,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); //--- exit the function return; } //--- check if(!isUpper && optype!=0) { //--- (A1' A21')-1 //--- X*A^-1 = (X1 X2)*( ) //--- ( A2') RMatrixRightTrsM(m,s1,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); RMatrixGemm(m,s2,s1,-1.0,x,i2,j2,0,a,i1+s1,j1,optype,1.0,x,i2,j2+s1); RMatrixRightTrsM(m,s2,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2,j2+s1); } } //--- exit the function return; } //+------------------------------------------------------------------+ //| Same as CMatrixLeftTRSM, but for real matrices | //| OpType may be only 0 or 1. | //+------------------------------------------------------------------+ static void CAblas::RMatrixLeftTrsM(const int m,const int n,CMatrixDouble &a, const int i1,const int j1,const bool isUpper, const bool isUnit,const int optype, CMatrixDouble &x,const int i2,const int j2) { //--- create variables int s1=0; int s2=0; int bs=AblasBlockSize(); //--- check if(m<=bs && n<=bs) { //--- basic algorithm RMatrixLeftTrsM2(m,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); //--- exit the function return; } //--- check if(n>=m) { //--- Split X: op(A)^-1*X = op(A)^-1*(X1 X2) AblasSplitLength(x,n,s1,s2); RMatrixLeftTrsM(m,s1,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); RMatrixLeftTrsM(m,s2,a,i1,j1,isUpper,isUnit,optype,x,i2,j2+s1); } else { //--- Split A AblasSplitLength(a,m,s1,s2); //--- check if(isUpper && optype==0) { //--- (A1 A12)-1 ( X1 ) //--- A^-1*X* = ( ) *( ) //--- ( A2) ( X2 ) RMatrixLeftTrsM(s2,n,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2+s1,j2); RMatrixGemm(s1,n,s2,-1.0,a,i1,j1+s1,0,x,i2+s1,j2,0,1.0,x,i2,j2); RMatrixLeftTrsM(s1,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); //--- exit the function return; } //--- check if(isUpper && optype!=0) { //--- (A1' )-1 ( X1 ) //--- A^-1*X = ( ) *( ) //--- (A12' A2') ( X2 ) RMatrixLeftTrsM(s1,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); RMatrixGemm(s2,n,s1,-1.0,a,i1,j1+s1,optype,x,i2,j2,0,1.0,x,i2+s1,j2); RMatrixLeftTrsM(s2,n,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2+s1,j2); //--- exit the function return; } //--- check if(!isUpper && optype==0) { //--- (A1 )-1 ( X1 ) //--- A^-1*X = ( ) *( ) //--- (A21 A2) ( X2 ) RMatrixLeftTrsM(s1,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); RMatrixGemm(s2,n,s1,-1.0,a,i1+s1,j1,0,x,i2,j2,0,1.0,x,i2+s1,j2); RMatrixLeftTrsM(s2,n,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2+s1,j2); //--- exit the function return; } //--- check if(!isUpper && optype!=0) { //--- (A1' A21')-1 ( X1 ) //--- A^-1*X = ( ) *( ) //--- ( A2') ( X2 ) RMatrixLeftTrsM(s2,n,a,i1+s1,j1+s1,isUpper,isUnit,optype,x,i2+s1,j2); RMatrixGemm(s1,n,s2,-1.0,a,i1+s1,j1,optype,x,i2+s1,j2,0,1.0,x,i2,j2); RMatrixLeftTrsM(s1,n,a,i1,j1,isUpper,isUnit,optype,x,i2,j2); } } //--- exit the function return; } //+------------------------------------------------------------------+ //| This subroutine calculates C=alpha*A*A^H+beta*C or | //| C=alpha*A^H*A+beta*C where: | //| * C is NxN Hermitian matrix given by its upper/lower triangle | //| * A is NxK matrix when A*A^H is calculated, KxN matrix otherwise | //| Additional info: | //| * cache-oblivious algorithm is used. | //| * multiplication result replaces C. If Beta=0, C elements are not| //| used in calculations (not multiplied by zero - just not | //| referenced) | //| * if Alpha=0, A is not used (not multiplied by zero - just not | //| referenced) | //| * if both Beta and Alpha are zero, C is filled by zeros. | //| INPUT PARAMETERS | //| N - matrix size, N>=0 | //| K - matrix size, K>=0 | //| Alpha - coefficient | //| A - matrix | //| IA - submatrix offset | //| JA - submatrix offset | //| OpTypeA - multiplication type: | //| * 0 - A*A^H is calculated | //| * 2 - A^H*A is calculated | //| Beta - coefficient | //| C - matrix | //| IC - submatrix offset | //| JC - submatrix offset | //| IsUpper - whether C is upper triangular or lower triangular| //+------------------------------------------------------------------+ static void CAblas::CMatrixSyrk(const int n,const int k,const double alpha, CMatrixComplex &a,const int ia,const int ja, const int optypea,const double beta,CMatrixComplex &c, const int ic,const int jc,const bool isUpper) { //--- create variables complex Alpha(alpha,0); complex Beta(beta,0); int s1=0; int s2=0; int bs=AblasComplexBlockSize(); //--- check if(n<=bs && k<=bs) { //--- basic algorithm CMatrixSyrk2(n,k,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); //--- exit the function return; } //--- check if(k>=n) { //--- Split K AblasComplexSplitLength(a,k,s1,s2); //--- check if(optypea==0) { CMatrixSyrk(n,s1,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); CMatrixSyrk(n,s2,alpha,a,ia,ja+s1,optypea,1.0,c,ic,jc,isUpper); } else { CMatrixSyrk(n,s1,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); CMatrixSyrk(n,s2,alpha,a,ia+s1,ja,optypea,1.0,c,ic,jc,isUpper); } } else { //--- Split N AblasComplexSplitLength(a,n,s1,s2); //--- check if(optypea==0 && isUpper) { CMatrixSyrk(s1,k,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); CMatrixGemm(s1,s2,k,Alpha,a,ia,ja,0,a,ia+s1,ja,2,Beta,c,ic,jc+s1); CMatrixSyrk(s2,k,alpha,a,ia+s1,ja,optypea,beta,c,ic+s1,jc+s1,isUpper); //--- exit the function return; } //--- check if(optypea==0 && !isUpper) { CMatrixSyrk(s1,k,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); CMatrixGemm(s2,s1,k,Alpha,a,ia+s1,ja,0,a,ia,ja,2,Beta,c,ic+s1,jc); CMatrixSyrk(s2,k,alpha,a,ia+s1,ja,optypea,beta,c,ic+s1,jc+s1,isUpper); //--- exit the function return; } //--- check if(optypea!=0 && isUpper) { CMatrixSyrk(s1,k,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); CMatrixGemm(s1,s2,k,Alpha,a,ia,ja,2,a,ia,ja+s1,0,Beta,c,ic,jc+s1); CMatrixSyrk(s2,k,alpha,a,ia,ja+s1,optypea,beta,c,ic+s1,jc+s1,isUpper); //--- exit the function return; } //--- check if(optypea!=0 && !isUpper) { CMatrixSyrk(s1,k,alpha,a,ia,ja,optypea,beta,c,ic,jc,isUpper); CMatrixGemm(s2,s1,k,Alpha,a,ia,ja+s1,2,a,ia,ja,0,Beta,c,ic+s1,jc); CMatrixSyrk(s2,k,alpha,a,ia,ja+s1,optypea,beta,c,ic+s1,jc+s1,isUpper); } } //--- exit the function return; } //+------------------------------------------------------------------+ //| This subroutine calculates C = alpha*op1(A)*op2(B) +beta*C where:| //| * C is MxN general matrix | //| * op1(A) is MxK matrix | //| * op2(B) is KxN matrix | //| * "op" may be identity transformation, transposition, conjugate | //| transposition | //| Additional info: | //| * cache-oblivious algorithm is used. | //| * multiplication result replaces C. If Beta=0, C elements are not| //| used in calculations (not multiplied by zero - just not | //| referenced) | //| * if Alpha=0, A is not used (not multiplied by zero - just not | //| referenced) | //| * if both Beta and Alpha are zero, C is filled by zeros. | //| INPUT PARAMETERS | //| N - matrix size, N>0 | //| M - matrix size, N>0 | //| K - matrix size, K>0 | //| Alpha - coefficient | //| A - matrix | //| IA - submatrix offset | //| JA - submatrix offset | //| OpTypeA - transformation type: | //| * 0 - no transformation | //| * 1 - transposition | //| * 2 - conjugate transposition | //| B - matrix | //| IB - submatrix offset | //| JB - submatrix offset | //| OpTypeB - transformation type: | //| * 0 - no transformation | //| * 1 - transposition | //| * 2 - conjugate transposition | //| Beta - coefficient | //| C - matrix | //| IC - submatrix offset | //| JC - submatrix offset | //+------------------------------------------------------------------+ static void CAblas::CMatrixGemm(const int m,const int n,const int k,complex &alpha, CMatrixComplex &a,const int ia,const int ja, const int optypea,CMatrixComplex &b,const int ib, const int jb,const int optypeb,complex &beta, CMatrixComplex &c,const int ic,const int jc) { //--- create variables complex Beta(1,0); int s1=0; int s2=0; int bs=AblasComplexBlockSize(); //--- check if(m<=bs && n<=bs && k<=bs) { //--- basic algorithm CMatrixGemmk(m,n,k,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); //--- exit the function return; } //--- check if(m>=n && m>=k) { //--- A*B = (A1 A2)^T*B AblasComplexSplitLength(a,m,s1,s2); CMatrixGemm(s1,n,k,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); //--- check if(optypea==0) CMatrixGemm(s2,n,k,alpha,a,ia+s1,ja,optypea,b,ib,jb,optypeb,beta,c,ic+s1,jc); else CMatrixGemm(s2,n,k,alpha,a,ia,ja+s1,optypea,b,ib,jb,optypeb,beta,c,ic+s1,jc); //--- exit the function return; } //--- check if(n>=m && n>=k) { //---A*B = A*(B1 B2) AblasComplexSplitLength(a,n,s1,s2); //--- check if(optypeb==0) { CMatrixGemm(m,s1,k,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); CMatrixGemm(m,s2,k,alpha,a,ia,ja,optypea,b,ib,jb+s1,optypeb,beta,c,ic,jc+s1); } else { CMatrixGemm(m,s1,k,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); CMatrixGemm(m,s2,k,alpha,a,ia,ja,optypea,b,ib+s1,jb,optypeb,beta,c,ic,jc+s1); } //--- exit the function return; } //--- check if(k>=m && k>=n) { //--- A*B = (A1 A2)*(B1 B2)^T AblasComplexSplitLength(a,k,s1,s2); //--- check if(optypea==0 && optypeb==0) { CMatrixGemm(m,n,s1,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); CMatrixGemm(m,n,s2,alpha,a,ia,ja+s1,optypea,b,ib+s1,jb,optypeb,Beta,c,ic,jc); } //--- check if(optypea==0 && optypeb!=0) { CMatrixGemm(m,n,s1,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); CMatrixGemm(m,n,s2,alpha,a,ia,ja+s1,optypea,b,ib,jb+s1,optypeb,Beta,c,ic,jc); } //--- check if(optypea!=0 && optypeb==0) { CMatrixGemm(m,n,s1,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); CMatrixGemm(m,n,s2,alpha,a,ia+s1,ja,optypea,b,ib+s1,jb,optypeb,Beta,c,ic,jc); } //--- check if(optypea!=0 && optypeb!=0) { CMatrixGemm(m,n,s1,alpha,a,ia,ja,optypea,b,ib,jb,optypeb,beta,c,ic,jc); CMatrixGemm(m,n,s2,alpha,a,ia+s1,ja,optypea,b,ib,jb+s1,optypeb,Beta,c,ic,jc); } } //--- exit the function return; } //+------------------------------------------------------------------+ //| Level 2 variant of CMatrixRightTRSM | //+------------------------------------------------------------------+ static void CAblas::CMatrixRightTrsM2(const int m,const int n,CMatrixComplex &a, const int i1,const int j1,const bool isUpper, const bool isUnit,const int optype, CMatrixComplex &x,const int i2,const int j2) { //--- check if(n*m==0) return; //--- create variables int i=0; int j=0; complex vc=0; complex vd=0; int i_=0; int i1_=0; //--- General case if(isUpper) { //--- Upper triangular matrix if(optype==0) { //--- X*A^(-1) for(i=0;i=0;j--) { vc=0; vd=1; //--- check if(j=0;j--) { vc=0; vd=1; //--- check if(j=0;j--) { //--- check if(isUnit) vd=1; else vd=a[i1+j][j1+j]; //--- change x x[i2+i].Set(j2+j,x[i2+i][j2+j]/vd); //--- check if(j>0) { vc=x[i2+i][j2+j]; i1_=j1-j2; for(i_=j2;i_0) { i1_=j1-j2; vc=0.0; for(i_=j2;i_0) { i1_=j1-j2; vc=0.0; for(i_=j2;i_=0;i--) { for(j=i+1;j<=m-1;j++) { vc=a[i1+i][j1+j]; //--- change x for(i_=j2;i_=0;i--) { //--- check if(isUnit) vd=1; else vd=Beta/a[i1+i][j1+i]; //--- change x for(i_=j2;i_=0;j--) { vc=a[i1+i][j1+j]; for(i_=j2;i_=0;i--) { //--- check if(isUnit) vd=1; else vd=Beta/CMath::Conj(a[i1+i][j1+i]); //--- change x for(i_=j2;i_=0;j--) { vc=CMath::Conj(a[i1+i][j1+j]); for(i_=j2;i_=0;j--) { vr=0; vd=1; //--- check if(j=0;j--) { //--- check if(isUnit) vd=1; else vd=a[i1+j][j1+j]; //--- change x x[i2+i].Set(j2+j,x[i2+i][j2+j]/vd); //--- check if(j>0) { vr=x[i2+i][j2+j]; i1_=j1-j2; //--- change x for(i_=j2;i_0) { i1_=j1-j2; vr=0.0; for(i_=j2;i_=0;i--) { for(j=i+1;j<=m-1;j++) { vr=a[i1+i][j1+j]; for(i_=j2;i_=0;i--) { //--- check if(isUnit) vd=1; else vd=1/a[i1+i][j1+i]; //--- change x for(i_=j2;i_=0;j--) { vr=a[i1+i][j1+j]; for(i_=j2;i_0) { v=0.0; for(i_=ja;i_<=ja+k-1;i_++) v+=a[ia+i][i_]*CMath::Conj(a[ia+j][i_]); } else v=0; //--- check if(beta==0) c[ic+i].Set(jc+j,Alpha*v); else c[ic+i].Set(jc+j,Beta*c[ic+i][jc+j]+Alpha*v); } } //--- exit the function return; } else { //--- C = alpha*A^H*A+beta*C for(i=0;iCAblas::AblasBlockSize()) blocksize=CAblas::AblasBlockSize(); //--- change rowscount=m-blockstart; //--- QR decomposition of submatrix. //--- Matrix is copied to temporary storage to solve //--- some TLB issues arising from non-contiguous memory //--- access pattern. CAblas::RMatrixCopy(rowscount,blocksize,a,blockstart,blockstart,tmpa,0,0); RMatrixQRBaseCase(tmpa,rowscount,blocksize,work,t,taubuf); CAblas::RMatrixCopy(rowscount,blocksize,tmpa,0,0,a,blockstart,blockstart); i1_=-blockstart; for(i_=blockstart;i_<=blockstart+blocksize-1;i_++) tau[i_]=taubuf[i_+i1_]; //--- check if(blockstart+blocksize<=n-1) { //--- Update the rest, choose between: //--- a) Level 2 algorithm (when the rest of the matrix is small enough) //--- b) blocked algorithm, see algorithm 5 from 'A storage efficient WY //--- representation for products of Householder transformations', //--- by R. Schreiber and C. Van Loan. if(n-blockstart-blocksize>=2*CAblas::AblasBlockSize() || rowscount>=4*CAblas::AblasBlockSize()) { //--- Prepare block reflector RMatrixBlockReflector(tmpa,taubuf,true,rowscount,blocksize,tmpt,work); //--- Multiply the rest of A by Q'. //--- Q = E + Y*T*Y' = E + TmpA*TmpT*TmpA' //--- Q' = E + Y*T'*Y' = E + TmpA*TmpT'*TmpA' CAblas::RMatrixGemm(blocksize,n-blockstart-blocksize,rowscount,1.0,tmpa,0,0,1,a,blockstart,blockstart+blocksize,0,0.0,tmpr,0,0); CAblas::RMatrixGemm(blocksize,n-blockstart-blocksize,blocksize,1.0,tmpt,0,0,1,tmpr,0,0,0,0.0,tmpr,blocksize,0); CAblas::RMatrixGemm(rowscount,n-blockstart-blocksize,blocksize,1.0,tmpa,0,0,0,tmpr,blocksize,0,0,1.0,a,blockstart,blockstart+blocksize); } else { //--- Level 2 algorithm for(i=0;iCAblas::AblasBlockSize()) blocksize=CAblas::AblasBlockSize(); //--- change columnscount=n-blockstart; //--- LQ decomposition of submatrix. //--- Matrix is copied to temporary storage to solve //--- some TLB issues arising from non-contiguous memory //--- access pattern. CAblas::RMatrixCopy(blocksize,columnscount,a,blockstart,blockstart,tmpa,0,0); RMatrixLQBaseCase(tmpa,blocksize,columnscount,work,t,taubuf); CAblas::RMatrixCopy(blocksize,columnscount,tmpa,0,0,a,blockstart,blockstart); i1_=-blockstart; for(i_=blockstart;i_<=blockstart+blocksize-1;i_++) { tau[i_]=taubuf[i_+i1_]; } //--- Update the rest, choose between: //--- a) Level 2 algorithm (when the rest of the matrix is small enough) //--- b) blocked algorithm, see algorithm 5 from 'A storage efficient WY //--- representation for products of Householder transformations', //--- by R. Schreiber and C. Van Loan. if(blockstart+blocksize<=m-1) { //--- check if(m-blockstart-blocksize>=2*CAblas::AblasBlockSize()) { //--- prepare RMatrixBlockReflector(tmpa,taubuf,false,columnscount,blocksize,tmpt,work); //--- Multiply the rest of A by Q. //--- Q = E + Y*T*Y' = E + TmpA'*TmpT*TmpA CAblas::RMatrixGemm(m-blockstart-blocksize,blocksize,columnscount,1.0,a,blockstart+blocksize,blockstart,0,tmpa,0,0,1,0.0,tmpr,0,0); CAblas::RMatrixGemm(m-blockstart-blocksize,blocksize,blocksize,1.0,tmpr,0,0,0,tmpt,0,0,0,0.0,tmpr,0,blocksize); CAblas::RMatrixGemm(m-blockstart-blocksize,columnscount,blocksize,1.0,tmpr,0,blocksize,0,tmpa,0,0,0,1.0,a,blockstart+blocksize,blockstart); } else { //--- Level 2 algorithm for(i=0;iCAblas::AblasComplexBlockSize()) blocksize=CAblas::AblasComplexBlockSize(); rowscount=m-blockstart; //--- QR decomposition of submatrix. //--- Matrix is copied to temporary storage to solve //--- some TLB issues arising from non-contiguous memory //--- access pattern. CAblas::CMatrixCopy(rowscount,blocksize,a,blockstart,blockstart,tmpa,0,0); CMatrixQRBaseCase(tmpa,rowscount,blocksize,work,t,taubuf); CAblas::CMatrixCopy(rowscount,blocksize,tmpa,0,0,a,blockstart,blockstart); i1_=-blockstart; for(i_=blockstart;i_<=blockstart+blocksize-1;i_++) tau[i_]=taubuf[i_+i1_]; //--- Update the rest, choose between: //--- a) Level 2 algorithm (when the rest of the matrix is small enough) //--- b) blocked algorithm, see algorithm 5 from 'A storage efficient WY //--- representation for products of Householder transformations', //--- by R. Schreiber and C. Van Loan. if(blockstart+blocksize<=n-1) { //--- check if(n-blockstart-blocksize>=2*CAblas::AblasComplexBlockSize()) { //--- Prepare block reflector CMatrixBlockReflector(tmpa,taubuf,true,rowscount,blocksize,tmpt,work); //--- Multiply the rest of A by Q'. //--- Q = E + Y*T*Y' = E + TmpA*TmpT*TmpA' //--- Q' = E + Y*T'*Y' = E + TmpA*TmpT'*TmpA' CAblas::CMatrixGemm(blocksize,n-blockstart-blocksize,rowscount,Alpha,tmpa,0,0,2,a,blockstart,blockstart+blocksize,0,Beta,tmpr,0,0); CAblas::CMatrixGemm(blocksize,n-blockstart-blocksize,blocksize,Alpha,tmpt,0,0,2,tmpr,0,0,0,Beta,tmpr,blocksize,0); CAblas::CMatrixGemm(rowscount,n-blockstart-blocksize,blocksize,Alpha,tmpa,0,0,0,tmpr,blocksize,0,0,Alpha,a,blockstart,blockstart+blocksize); } else { //--- Level 2 algorithm for(i=0;iCAblas::AblasComplexBlockSize()) blocksize=CAblas::AblasComplexBlockSize(); columnscount=n-blockstart; //--- LQ decomposition of submatrix. //--- Matrix is copied to temporary storage to solve //--- some TLB issues arising from non-contiguous memory //--- access pattern. CAblas::CMatrixCopy(blocksize,columnscount,a,blockstart,blockstart,tmpa,0,0); CMatrixLQBaseCase(tmpa,blocksize,columnscount,work,t,taubuf); CAblas::CMatrixCopy(blocksize,columnscount,tmpa,0,0,a,blockstart,blockstart); i1_=-blockstart; for(i_=blockstart;i_<=blockstart+blocksize-1;i_++) tau[i_]=taubuf[i_+i1_]; //--- Update the rest, choose between: //--- a) Level 2 algorithm (when the rest of the matrix is small enough) //--- b) blocked algorithm, see algorithm 5 from 'A storage efficient WY //--- representation for products of Householder transformations', //--- by R. Schreiber and C. Van Loan. if(blockstart+blocksize<=m-1) { //--- check if(m-blockstart-blocksize>=2*CAblas::AblasComplexBlockSize()) { //--- Prepare block reflector CMatrixBlockReflector(tmpa,taubuf,false,columnscount,blocksize,tmpt,work); //--- Multiply the rest of A by Q. //--- Q = E + Y*T*Y' = E + TmpA'*TmpT*TmpA CAblas::CMatrixGemm(m-blockstart-blocksize,blocksize,columnscount,Alpha,a,blockstart+blocksize,blockstart,0,tmpa,0,0,2,Beta,tmpr,0,0); CAblas::CMatrixGemm(m-blockstart-blocksize,blocksize,blocksize,Alpha,tmpr,0,0,0,tmpt,0,0,0,Beta,tmpr,0,blocksize); CAblas::CMatrixGemm(m-blockstart-blocksize,columnscount,blocksize,Alpha,tmpr,0,blocksize,0,tmpa,0,0,0,Alpha,a,blockstart+blocksize,blockstart); } else { //--- Level 2 algorithm for(i=0;i=0. | //| N - number of columns in given matrix A. N>=0. | //| Tau - scalar factors which are used to form Q. | //| Output of the RMatrixQR subroutine. | //| QColumns - required number of columns of matrix Q. | //| M>=QColumns>=0. | //| Output parameters: | //| Q - first QColumns columns of matrix Q. | //| Array whose indexes range within | //| [0..M-1, 0..QColumns-1]. | //| If QColumns=0, the array remains unchanged. | //+------------------------------------------------------------------+ static void COrtFac::RMatrixQRUnpackQ(CMatrixDouble &a,const int m,const int n, double &tau[],const int qcolumns,CMatrixDouble &q) { //--- check if(!CAp::Assert(qcolumns<=m,__FUNCTION__+": QColumns>M!")) return; //--- check if(m<=0 || n<=0 || qcolumns<=0) return; //--- create arrays double work[]; double t[]; double taubuf[]; //--- create matrix CMatrixDouble tmpa; CMatrixDouble tmpt; CMatrixDouble tmpr; //--- create variables int minmn=MathMin(m,n); int refcnt=MathMin(minmn,qcolumns); int blockstart=CAblas::AblasBlockSize()*(refcnt/CAblas::AblasBlockSize()); int blocksize=refcnt-blockstart; int rowscount=0; int i=0; int j=0; int i_=0; int i1_=0; //--- allocation q.Resize(m,qcolumns); //--- identity matrix for(i=0;i=0) { rowscount=m-blockstart; //--- Copy current block CAblas::RMatrixCopy(rowscount,blocksize,a,blockstart,blockstart,tmpa,0,0); i1_=blockstart; for(i_=0;i_=2*CAblas::AblasBlockSize()) { //--- Prepare block reflector RMatrixBlockReflector(tmpa,taubuf,true,rowscount,blocksize,tmpt,work); //--- Multiply matrix by Q. //--- Q = E + Y*T*Y' = E + TmpA*TmpT*TmpA' CAblas::RMatrixGemm(blocksize,qcolumns,rowscount,1.0,tmpa,0,0,1,q,blockstart,0,0,0.0,tmpr,0,0); CAblas::RMatrixGemm(blocksize,qcolumns,blocksize,1.0,tmpt,0,0,0,tmpr,0,0,0,0.0,tmpr,blocksize,0); CAblas::RMatrixGemm(rowscount,qcolumns,blocksize,1.0,tmpa,0,0,0,tmpr,blocksize,0,0,1.0,q,blockstart,0); } else { //--- Level 2 algorithm for(i=blocksize-1;i>=0;i--) { i1_=i-1; for(i_=1;i_<=rowscount-i;i_++) { t[i_]=tmpa[i_+i1_][i]; } t[1]=1; //--- function call CReflections::ApplyReflectionFromTheLeft(q,taubuf[i],t,blockstart+i,m-1,0,qcolumns-1,work); } } //--- change value blockstart=blockstart-CAblas::AblasBlockSize(); blocksize=CAblas::AblasBlockSize(); } } //+------------------------------------------------------------------+ //| Unpacking of matrix R from the QR decomposition of a matrix A | //| Input parameters: | //| A - matrices Q and R in compact form. | //| Output of RMatrixQR subroutine. | //| M - number of rows in given matrix A. M>=0. | //| N - number of columns in given matrix A. N>=0. | //| Output parameters: | //| R - matrix R, array[0..M-1, 0..N-1]. | //+------------------------------------------------------------------+ static void COrtFac::RMatrixQRUnpackR(CMatrixDouble &a,const int m,const int n,CMatrixDouble &r) { //--- check if(m<=0 || n<=0) return; //--- create variables int i=0; int k=MathMin(m,n); int i_=0; //--- allocation r.Resize(m,n); //--- Prepare matrix for(i=0;i=0. | //| N - number of columns in given matrix A. N>=0. | //| Tau - scalar factors which are used to form Q. | //| Output of the RMatrixLQ subroutine. | //| QRows - required number of rows in matrix Q. N>=QRows>=0.| //| Output parameters: | //| Q - first QRows rows of matrix Q. Array whose indexes| //| range within [0..QRows-1, 0..N-1]. If QRows=0, | //| the array remains unchanged. | //+------------------------------------------------------------------+ static void COrtFac::RMatrixLQUnpackQ(CMatrixDouble &a,const int m,const int n, double &tau[],const int qrows,CMatrixDouble &q) { //--- check if(!CAp::Assert(qrows<=n,__FUNCTION__+": QRows>N!")) return; //--- check if(m<=0 || n<=0 || qrows<=0) return; //--- create arrays double work[]; double t[]; double taubuf[]; //--- create matrix CMatrixDouble tmpa; CMatrixDouble tmpt; CMatrixDouble tmpr; //--- create variables int minmn=MathMin(m,n); int refcnt=MathMin(minmn,qrows); int blockstart=CAblas::AblasBlockSize()*(refcnt/CAblas::AblasBlockSize()); int blocksize=refcnt-blockstart; int columnscount=0; int i=0; int j=0; int i_=0; int i1_=0; //--- allocation ArrayResizeAL(work,MathMax(m,n)+1); ArrayResizeAL(t,MathMax(m,n)+1); ArrayResizeAL(taubuf,minmn); //--- allocation tmpa.Resize(CAblas::AblasBlockSize(),n); tmpt.Resize(CAblas::AblasBlockSize(),2*CAblas::AblasBlockSize()); tmpr.Resize(qrows,2*CAblas::AblasBlockSize()); q.Resize(qrows,n); //--- identity matrix for(i=0;i<=qrows-1;i++) { for(j=0;j=0) { columnscount=n-blockstart; //--- Copy submatrix CAblas::RMatrixCopy(blocksize,columnscount,a,blockstart,blockstart,tmpa,0,0); i1_=blockstart; for(i_=0;i_=2*CAblas::AblasBlockSize()) { //--- Prepare block reflector RMatrixBlockReflector(tmpa,taubuf,false,columnscount,blocksize,tmpt,work); //--- Multiply the rest of A by Q'. //--- Q' = E + Y*T'*Y' = E + TmpA'*TmpT'*TmpA CAblas::RMatrixGemm(qrows,blocksize,columnscount,1.0,q,0,blockstart,0,tmpa,0,0,1,0.0,tmpr,0,0); CAblas::RMatrixGemm(qrows,blocksize,blocksize,1.0,tmpr,0,0,0,tmpt,0,0,1,0.0,tmpr,0,blocksize); CAblas::RMatrixGemm(qrows,columnscount,blocksize,1.0,tmpr,0,blocksize,0,tmpa,0,0,0,1.0,q,0,blockstart); } else { //--- Level 2 algorithm for(i=blocksize-1;i>=0;i--) { i1_=i-1; for(i_=1;i_<=columnscount-i;i_++) t[i_]=tmpa[i][i_+i1_]; t[1]=1; //--- function call CReflections::ApplyReflectionFromTheRight(q,taubuf[i],t,0,qrows-1,blockstart+i,n-1,work); } } //--- change value blockstart=blockstart-CAblas::AblasBlockSize(); blocksize=CAblas::AblasBlockSize(); } } //+------------------------------------------------------------------+ //| Unpacking of matrix L from the LQ decomposition of a matrix A | //| Input parameters: | //| A -matrices Q and L in compact form. | //| Output of RMatrixLQ subroutine. | //| M -number of rows in given matrix A. M>=0. | //| N -number of columns in given matrix A. N>=0. | //| Output parameters: | //| L -matrix L, array[0..M-1,0..N-1]. | //+------------------------------------------------------------------+ static void COrtFac::RMatrixLQUnpackL(CMatrixDouble &a,const int m,const int n,CMatrixDouble &l) { //--- check if(m<=0 || n<=0) return; //--- create variables int i=0; int k=0; int i_=0; //--- allocation l.Resize(m,n); //--- Prepare matrix for(i=0;i=0. | //| N - number of columns in matrix A. N>=0. | //| Tau - scalar factors which are used to form Q. | //| Output of CMatrixQR subroutine . | //| QColumns - required number of columns in matrix Q. | //| M>=QColumns>=0. | //| Output parameters: | //| Q - first QColumns columns of matrix Q. | //| Array whose index ranges within [0..M-1, | //| 0..QColumns-1]. | //| If QColumns=0, array isn't changed. | //+------------------------------------------------------------------+ static void COrtFac::CMatrixQRUnpackQ(CMatrixComplex &a,const int m,const int n, complex &tau[],const int qcolumns,CMatrixComplex &q) { //--- check if(!CAp::Assert(qcolumns<=m,__FUNCTION__+": QColumns>M!")) return; //--- check if(m<=0 || n<=0) return; //--- create arrays complex work[]; complex t[]; complex taubuf[]; //--- create matrix CMatrixComplex tmpa; CMatrixComplex tmpt; CMatrixComplex tmpr; //--- create variables int minmn=MathMin(m,n); int refcnt=MathMin(minmn,qcolumns); int blockstart=CAblas::AblasComplexBlockSize()*(refcnt/CAblas::AblasComplexBlockSize()); int blocksize=refcnt-blockstart; int rowscount=0; int i=0; int j=0; int i_=0; int i1_=0; complex One(1,0); complex Zero(0,0); //--- allocation ArrayResizeAL(work,MathMax(m,n)+1); ArrayResizeAL(t,MathMax(m,n)+1); ArrayResizeAL(taubuf,minmn); //--- allocation tmpa.Resize(m,CAblas::AblasComplexBlockSize()); tmpt.Resize(CAblas::AblasComplexBlockSize(),CAblas::AblasComplexBlockSize()); tmpr.Resize(2*CAblas::AblasComplexBlockSize(),qcolumns); q.Resize(m,qcolumns); for(i=0;i=0) { rowscount=m-blockstart; //--- QR decomposition of submatrix. //--- Matrix is copied to temporary storage to solve //--- some TLB issues arising from non-contiguous memory //--- access pattern. CAblas::CMatrixCopy(rowscount,blocksize,a,blockstart,blockstart,tmpa,0,0); i1_=blockstart; for(i_=0;i_=2*CAblas::AblasComplexBlockSize()) { //--- Prepare block reflector CMatrixBlockReflector(tmpa,taubuf,true,rowscount,blocksize,tmpt,work); //--- Multiply the rest of A by Q. //--- Q = E + Y*T*Y' = E + TmpA*TmpT*TmpA' CAblas::CMatrixGemm(blocksize,qcolumns,rowscount,One,tmpa,0,0,2,q,blockstart,0,0,Zero,tmpr,0,0); CAblas::CMatrixGemm(blocksize,qcolumns,blocksize,One,tmpt,0,0,0,tmpr,0,0,0,Zero,tmpr,blocksize,0); CAblas::CMatrixGemm(rowscount,qcolumns,blocksize,One,tmpa,0,0,0,tmpr,blocksize,0,0,One,q,blockstart,0); } else { //--- Level 2 algorithm for(i=blocksize-1;i>=0;i--) { i1_=i-1; for(i_=1;i_<=rowscount-i;i_++) t[i_]=tmpa[i_+i1_][i]; t[1]=1; //--- function call CComplexReflections::ComplexApplyReflectionFromTheLeft(q,taubuf[i],t,blockstart+i,m-1,0,qcolumns-1,work); } } //--- change value blockstart=blockstart-CAblas::AblasComplexBlockSize(); blocksize=CAblas::AblasComplexBlockSize(); } } //+------------------------------------------------------------------+ //| Unpacking of matrix R from the QR decomposition of a matrix A | //| Input parameters: | //| A - matrices Q and R in compact form. | //| Output of CMatrixQR subroutine. | //| M - number of rows in given matrix A. M>=0. | //| N - number of columns in given matrix A. N>=0. | //| Output parameters: | //| R - matrix R, array[0..M-1, 0..N-1]. | //+------------------------------------------------------------------+ static void COrtFac::CMatrixQRUnpackR(CMatrixComplex &a,const int m,const int n,CMatrixComplex &r) { //--- check if(m<=0 || n<=0) return; //--- create variables complex Zero(0,0); int i=0; int k=MathMin(m,n); int i_=0; //--- allocation r.Resize(m,n); //--- Prepare matrix for(i=0;i=0. | //| N - number of columns in matrix A. N>=0. | //| Tau - scalar factors which are used to form Q. | //| Output of CMatrixLQ subroutine . | //| QRows - required number of rows in matrix Q. | //| N>=QColumns>=0. | //| Output parameters: | //| Q - first QRows rows of matrix Q. | //| Array whose index ranges within [0..QRows-1, | //| 0..N-1]. | //| If QRows=0, array isn't changed. | //+------------------------------------------------------------------+ static void COrtFac::CMatrixLQUnpackQ(CMatrixComplex &a,const int m,const int n, complex &tau[],const int qrows,CMatrixComplex &q) { //--- check if(m<=0 || n<=0) return; //--- create arrays complex work[]; complex t[]; complex taubuf[]; //--- create matrix CMatrixComplex tmpa; CMatrixComplex tmpt; CMatrixComplex tmpr; //--- create variables int minmn=MathMin(m,n); int refcnt=MathMin(minmn,qrows); int blockstart=CAblas::AblasComplexBlockSize()*(refcnt/CAblas::AblasComplexBlockSize()); int blocksize=refcnt-blockstart; int columnscount=0; int i=0; int j=0; int i_=0; int i1_=0; complex One(1,0); complex Zero(0,0); //--- allocation ArrayResizeAL(work,MathMax(m,n)+1); ArrayResizeAL(t,MathMax(m,n)+1); ArrayResizeAL(taubuf,minmn); //--- allocation tmpa.Resize(CAblas::AblasComplexBlockSize(),n); tmpt.Resize(CAblas::AblasComplexBlockSize(),CAblas::AblasComplexBlockSize()); tmpr.Resize(qrows,2*CAblas::AblasComplexBlockSize()); q.Resize(qrows,n); for(i=0;i<=qrows-1;i++) { for(j=0;j=0) { columnscount=n-blockstart; //--- LQ decomposition of submatrix. //--- Matrix is copied to temporary storage to solve //--- some TLB issues arising from non-contiguous memory //--- access pattern. CAblas::CMatrixCopy(blocksize,columnscount,a,blockstart,blockstart,tmpa,0,0); i1_=blockstart; for(i_=0;i_=2*CAblas::AblasComplexBlockSize()) { //--- Prepare block reflector CMatrixBlockReflector(tmpa,taubuf,false,columnscount,blocksize,tmpt,work); //--- Multiply the rest of A by Q'. //--- Q' = E + Y*T'*Y' = E + TmpA'*TmpT'*TmpA CAblas::CMatrixGemm(qrows,blocksize,columnscount,One,q,0,blockstart,0,tmpa,0,0,2,Zero,tmpr,0,0); CAblas::CMatrixGemm(qrows,blocksize,blocksize,One,tmpr,0,0,0,tmpt,0,0,2,Zero,tmpr,0,blocksize); CAblas::CMatrixGemm(qrows,columnscount,blocksize,One,tmpr,0,blocksize,0,tmpa,0,0,0,One,q,0,blockstart); } else { //--- Level 2 algorithm for(i=blocksize-1;i>=0;i--) { i1_=i-1; for(i_=1;i_<=columnscount-i;i_++) t[i_]=CMath::Conj(tmpa[i][i_+i1_]); t[1]=1; //--- function call CComplexReflections::ComplexApplyReflectionFromTheRight(q,CMath::Conj(taubuf[i]),t,0,qrows-1,blockstart+i,n-1,work); } } //--- change value blockstart=blockstart-CAblas::AblasComplexBlockSize(); blocksize=CAblas::AblasComplexBlockSize(); } } //+------------------------------------------------------------------+ //| Unpacking of matrix L from the LQ decomposition of a matrix A | //| Input parameters: | //| A - matrices Q and L in compact form. | //| Output of CMatrixLQ subroutine. | //| M - number of rows in given matrix A. M>=0. | //| N - number of columns in given matrix A. N>=0. | //| Output parameters: | //| L - matrix L, array[0..M-1, 0..N-1]. | //+------------------------------------------------------------------+ static void COrtFac::CMatrixLQUnpackL(CMatrixComplex &a,const int m,const int n,CMatrixComplex &l) { //--- check if(m<=0 || n<=0) return; //--- create variables complex Zero(0,0); int i=0; int k=0; int i_=0; //--- allocation l.Resize(m,n); //--- Prepare matrix for(i=0;i=N, B is the upper bidiagonal MxN matrix and is stored in | //| the corresponding elements of matrix A. Matrix Q is represented | //| as a product of elementary reflections Q = H(0)*H(1)*...*H(n-1), | //| where H(i) = 1-tau*v*v'. Here tau is a scalar which is stored in | //| TauQ[i], and vector v has the following structure: v(0:i-1)=0, | //| v(i)=1, v(i+1:m-1) is stored in elements A(i+1:m-1,i).Matrix P is| //| as follows: P = G(0)*G(1)*...*G(n-2), where G(i) = 1 - tau*u*u'. | //| Tau is stored in TauP[i], u(0:i)=0, u(i+1)=1, u(i+2:n-1) is | //| stored in elements A(i,i+2:n-1). | //| If M n): m=5, n=6 (m < n): | //| ( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 ) | //| ( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 ) | //| ( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 ) | //| ( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 ) | //| ( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 ) | //| ( v1 v2 v3 v4 v5 ) | //| Here vi and ui are vectors which form H(i) and G(i), and d and | //| e - are the diagonal and off-diagonal elements of matrix B. | //+------------------------------------------------------------------+ static void COrtFac::RMatrixBD(CMatrixDouble &a,const int m,const int n,double &tauq[],double &taup[]) { //--- check if(n<=0 || m<=0) return; //--- create arrays double work[]; double t[]; //--- create variables int minmn=0; int maxmn=MathMax(m,n); int i=0; double ltau=0; int i_=0; int i1_=0; //--- allocation ArrayResizeAL(work,maxmn+1); ArrayResizeAL(t,maxmn+1); //--- check if(m>=n) { ArrayResizeAL(tauq,n); ArrayResizeAL(taup,n); } else { ArrayResizeAL(tauq,m); ArrayResizeAL(taup,m); } //--- check if(m>=n) { //--- Reduce to upper bidiagonal form for(i=0;i=QColumns>=0. | //| Output parameters: | //| Q - first QColumns columns of matrix Q. | //| Array[0..M-1, 0..QColumns-1] | //| If QColumns=0, the array is not modified. | //+------------------------------------------------------------------+ static void COrtFac::RMatrixBDUnpackQ(CMatrixDouble &qp,const int m,const int n, double &tauq[],const int qcolumns,CMatrixDouble &q) { //--- check if(!CAp::Assert(qcolumns<=m,__FUNCTION__+": QColumns>M!")) return; //--- check if(!CAp::Assert(qcolumns>=0,__FUNCTION__+": QColumns<0!")) return; //--- check if(m==0 || n==0 || qcolumns==0) return; //--- create variables int i=0; int j=0; //--- allocation q.Resize(m,qcolumns); //--- identity matrix for(i=0;i=n) { //--- setup if(fromtheright) { i1=0; i2=n-1; istep=1; } else { i1=n-1; i2=0; istep=-1; } //--- check if(dotranspose) { i=i1; i1=i2; i2=i; istep=-istep; } //--- Process i=i1; do { i1_=i-1; for(i_=1;i_<=m-i;i_++) v[i_]=qp[i_+i1_][i]; v[1]=1; //--- check if(fromtheright) CReflections::ApplyReflectionFromTheRight(z,tauq[i],v,0,zrows-1,i,m-1,work); else CReflections::ApplyReflectionFromTheLeft(z,tauq[i],v,i,m-1,0,zcolumns-1,work); i=i+istep; } while(i!=i2+istep); } else { //--- setup if(fromtheright) { i1=0; i2=m-2; istep=1; } else { i1=m-2; i2=0; istep=-1; } //--- check if(dotranspose) { i=i1; i1=i2; i2=i; istep=-istep; } //--- Process if(m-1>0) { i=i1; do { i1_=i; for(i_=1;i_<=m-i-1;i_++) v[i_]=qp[i_+i1_][i]; v[1]=1; //--- check if(fromtheright) CReflections::ApplyReflectionFromTheRight(z,tauq[i],v,0,zrows-1,i+1,m-1,work); else CReflections::ApplyReflectionFromTheLeft(z,tauq[i],v,i+1,m-1,0,zcolumns-1,work); i=i+istep; } while(i!=i2+istep); } } } //+------------------------------------------------------------------+ //| Unpacking matrix P which reduces matrix A to bidiagonal form. | //| The subroutine returns transposed matrix P. | //| Input parameters: | //| QP - matrices Q and P in compact form. | //| Output of ToBidiagonal subroutine. | //| M - number of rows in matrix A. | //| N - number of columns in matrix A. | //| TAUP - scalar factors which are used to form P. | //| Output of ToBidiagonal subroutine. | //| PTRows - required number of rows of matrix P^T. | //| N >= PTRows >= 0. | //| Output parameters: | //| PT - first PTRows columns of matrix P^T | //| Array[0..PTRows-1, 0..N-1] | //| If PTRows=0, the array is not modified. | //+------------------------------------------------------------------+ static void COrtFac::RMatrixBDUnpackPT(CMatrixDouble &qp,const int m,const int n, double &taup[],const int ptrows,CMatrixDouble &pt) { //--- check if(!CAp::Assert(ptrows<=n,__FUNCTION__+": PTRows>N!")) return; //--- check if(!CAp::Assert(ptrows>=0,__FUNCTION__+": PTRows<0!")) return; //--- check if(m==0 || n==0 || ptrows==0) return; //--- create variables int i=0; int j=0; //--- allocation pt.Resize(ptrows,n); //--- prepare for(i=0;i<=ptrows-1;i++) { for(j=0;j=n) { //--- setup if(fromtheright) { i1=n-2; i2=0; istep=-1; } else { i1=0; i2=n-2; istep=1; } //--- check if(!dotranspose) { i=i1; i1=i2; i2=i; istep=-istep; } //--- Process if(n-1>0) { i=i1; do { i1_=i; for(i_=1;i_=n) isupper=true; else isupper=false; //--- check if(isupper) { //--- allocation ArrayResizeAL(d,n); ArrayResizeAL(e,n); //--- get result for(i=0;i=0,__FUNCTION__+": incorrect N!")) return; //--- create arrays double t[]; double work[]; //--- create variables int i=0; double v=0; int i_=0; int i1_=0; //--- allocation ArrayResizeAL(tau,n-1); ArrayResizeAL(t,n+1); ArrayResizeAL(work,n); for(i=0;i1) ArrayResizeAL(tau,n-1); ArrayResizeAL(d,n); //--- check if(n>1) ArrayResizeAL(e,n-1); //--- check if(isupper) { //--- Reduce the upper triangle of A for(i=n-2;i>=0;i--) { //--- Generate elementary reflector H() = E - tau * v * v' if(i>=1) { i1_=-2; for(i_=2;i_<=i+1;i_++) t[i_]=a[i_+i1_][i+1]; } t[1]=a[i][i+1]; CReflections::GenerateReflection(t,i+1,taui); //--- check if(i>=1) { i1_=2; for(i_=0;i_=0;i--) { //--- Apply H(i) i1_=i; for(i_=1;i_1) { ArrayResizeAL(tau,n-1); ArrayResizeAL(e,n-1); } ArrayResizeAL(d,n); ArrayResizeAL(t,n); ArrayResizeAL(t2,n); ArrayResizeAL(t3,n); //--- check if(isupper) { //--- Reduce the upper triangle of A a[n-1].Set(n-1,a[n-1][n-1].re); for(i=n-2;i>=0;i--) { //--- Generate elementary reflector H = I+1 - tau * v * v' alpha=a[i][i+1]; t[1]=alpha; //--- check if(i>=1) { i1_=-2; for(i_=2;i_<=i+1;i_++) t[i_]=a[i_+i1_][i+1]; } //--- function call CComplexReflections::ComplexGenerateReflection(t,i+1,taui); //--- check if(i>=1) { i1_=2; for(i_=0;i_=0;i--) { //--- Apply H(i) i1_=i; for(i_=1;i_=0). | //| W - array of the eigenvalues found. | //| Array whose index ranges within [0..M-1]. | //| Z - if ZNeeded is equal to: | //| * 0, Z hasn?t changed; | //| * 1, Z contains eigenvectors. | //| Array whose indexes range within | //| [0..N-1, 0..M-1]. | //| The eigenvectors are stored in the matrix | //| columns. | //| Result: | //| True, if successful. M contains the number of eigenvalues in | //| the given half-interval (could be equal to 0), W contains the| //| eigenvalues, Z contains the eigenvectors (if needed). | //| False, if the bisection method subroutine wasn't able to find| //| the eigenvalues in the given interval or if the inverse | //| iteration subroutine wasn't able to find all the | //| corresponding eigenvectors. In that case, the eigenvalues | //| and eigenvectors are not returned, M is equal to 0. | //+------------------------------------------------------------------+ static bool CEigenVDetect::SMatrixEVDR(CMatrixDouble &ca,const int n,const int zneeded, const bool isupper,const double b1,const double b2, int &m,double &w[],CMatrixDouble &z) { //--- create arrays double tau[]; double e[]; //--- create copy CMatrixDouble a; a=ca; //--- initialization m=0; //--- check if(!CAp::Assert(zneeded==0 || zneeded==1,__FUNCTION__+": incorrect ZNeeded")) return(false); //--- function call COrtFac::SMatrixTD(a,n,isupper,tau,w,e); //--- check if(zneeded==1) { //--- function call COrtFac::SMatrixTDUnpackQ(a,n,isupper,tau,z); } //--- return result return(SMatrixTdEVDR(w,e,n,zneeded,b1,b2,m,z)); } //+------------------------------------------------------------------+ //| Subroutine for finding the eigenvalues and eigenvectors of a | //| symmetric matrix with given indexes by using bisection and | //| inverse iteration methods. | //| Input parameters: | //| A - symmetric matrix which is given by its upper or | //| lower triangular part. Array whose indexes range | //| within [0..N-1, 0..N-1]. | //| N - size of matrix A. | //| ZNeeded - flag controlling whether the eigenvectors are | //| needed or not. | //| If ZNeeded is equal to: | //| * 0, the eigenvectors are not returned; | //| * 1, the eigenvectors are returned. | //| IsUpperA - storage format of matrix A. | //| I1, I2 - index interval for searching (from I1 to I2). | //| 0 <= I1 <= I2 <= N-1. | //| Output parameters: | //| W - array of the eigenvalues found. | //| Array whose index ranges within [0..I2-I1]. | //| Z - if ZNeeded is equal to: | //| * 0, Z hasn?t changed; | //| * 1, Z contains eigenvectors. | //| Array whose indexes range within | //| [0..N-1, 0..I2-I1]. | //| In that case, the eigenvectors are stored in the | //| matrix columns. | //| Result: | //| True, if successful. W contains the eigenvalues, Z contains | //| the eigenvectors (if needed). | //| False, if the bisection method subroutine wasn't able to find| //| the eigenvalues in the given interval or if the inverse | //| iteration subroutine wasn't able to find all the | //| corresponding eigenvectors. In that case, the eigenvalues | //| and eigenvectors are not returned. | //+------------------------------------------------------------------+ static bool CEigenVDetect::SMatrixEVDI(CMatrixDouble &ca,const int n,const int zneeded, const bool isupper,const int i1,const int i2, double &w[],CMatrixDouble &z) { //--- create arrays double tau[]; double e[]; //--- create copy CMatrixDouble a; a=ca; //--- check if(!CAp::Assert(zneeded==0 || zneeded==1,__FUNCTION__+": incorrect ZNeeded")) return(false); //--- function call COrtFac::SMatrixTD(a,n,isupper,tau,w,e); //--- check if(zneeded==1) { //--- function call COrtFac::SMatrixTDUnpackQ(a,n,isupper,tau,z); } //--- return result return(SMatrixTdEVDI(w,e,n,zneeded,i1,i2,z)); } //+------------------------------------------------------------------+ //| Finding the eigenvalues and eigenvectors of a Hermitian matrix | //| The algorithm finds eigen pairs of a Hermitian matrix by reducing| //| it to real tridiagonal form and using the QL/QR algorithm. | //| Input parameters: | //| A - Hermitian matrix which is given by its upper or | //| lower triangular part. | //| Array whose indexes range within | //| [0..N-1, 0..N-1]. | //| N - size of matrix A. | //| IsUpper - storage format. | //| ZNeeded - flag controlling whether the eigenvectors are | //| needed or not. If ZNeeded is equal to: | //| * 0, the eigenvectors are not returned; | //| * 1, the eigenvectors are returned. | //| Output parameters: | //| D - eigenvalues in ascending order. | //| Array whose index ranges within [0..N-1]. | //| Z - if ZNeeded is equal to: | //| * 0, Z hasn?t changed; | //| * 1, Z contains the eigenvectors. | //| Array whose indexes range within | //| [0..N-1, 0..N-1]. | //| The eigenvectors are stored in the matrix | //| columns. | //| Result: | //| True, if the algorithm has converged. | //| False, if the algorithm hasn't converged (rare case). | //| Note: | //| eigenvectors of Hermitian matrix are defined up to | //| multiplication by a complex number L, such that |L|=1. | //+------------------------------------------------------------------+ static bool CEigenVDetect::HMatrixEVD(CMatrixComplex &ca,const int n,int zneeded, const bool isupper,double &d[],CMatrixComplex &z) { //--- create variables int i=0; int k=0; double v=0; int i_=0; bool result; //--- create arrays complex tau[]; double e[]; double work[]; //--- create matrix CMatrixDouble t; CMatrixComplex q; //--- create copy CMatrixComplex a; a=ca; //--- check if(!CAp::Assert(zneeded==0 || zneeded==1,__FUNCTION__+": incorrect ZNeeded")) return(false); //--- Reduce to tridiagonal form COrtFac::HMatrixTD(a,n,isupper,tau,d,e); //--- check if(zneeded==1) { //--- function call COrtFac::HMatrixTDUnpackQ(a,n,isupper,tau,q); zneeded=2; } //--- get result result=SMatrixTdEVD(d,e,n,zneeded,t); //--- Eigenvectors are needed //--- Calculate Z = Q*T = Re(Q)*T + i*Im(Q)*T if(result && zneeded!=0) { ArrayResizeAL(work,n); z.Resize(n,n); for(i=0;i=0 | //| W - array of the eigenvalues found. | //| Array whose index ranges within [0..M-1]. | //| Z - if ZNeeded is equal to: | //| * 0, Z hasn?t changed; | //| * 1, Z contains eigenvectors. | //| Array whose indexes range within | //| [0..N-1, 0..M-1]. | //| The eigenvectors are stored in the matrix | //| columns. | //| Result: | //| True, if successful. M contains the number of eigenvalues | //| in the given half-interval (could be equal to 0), W contains | //| the eigenvalues, Z contains the eigenvectors (if needed). | //| False, if the bisection method subroutine wasn't able to find| //| the eigenvalues in the given interval or if the inverse | //| iteration subroutine wasn't able to find all the | //| corresponding eigenvectors. In that case, the eigenvalues and| //| eigenvectors are not returned, M is equal to 0. | //| Note: | //| eigen vectors of Hermitian matrix are defined up to | //| multiplication by a complex number L, such as |L|=1. | //+------------------------------------------------------------------+ static bool CEigenVDetect::HMatrixEVDR(CMatrixComplex &ca,const int n,int zneeded, bool isupper,const double b1,const double b2, int &m,double &w[],CMatrixComplex &z) { //--- create variables int i=0; int k=0; double v=0; int i_=0; bool result; //--- create arrays complex tau[]; double e[]; double work[]; //--- create matrix CMatrixComplex q; CMatrixDouble t; //--- create copy CMatrixComplex a; a=ca; //--- initialization m=0; //--- check if(!CAp::Assert(zneeded==0 || zneeded==1,__FUNCTION__+": incorrect ZNeeded")) return(false); //--- Reduce to tridiagonal form COrtFac::HMatrixTD(a,n,isupper,tau,w,e); //--- check if(zneeded==1) { //--- function call COrtFac::HMatrixTDUnpackQ(a,n,isupper,tau,q); zneeded=2; } //--- Bisection and inverse iteration result=SMatrixTdEVDR(w,e,n,zneeded,b1,b2,m,t); //--- Eigenvectors are needed //--- Calculate Z = Q*T = Re(Q)*T + i*Im(Q)*T if((result && zneeded!=0) && m!=0) { ArrayResizeAL(work,m); z.Resize(n,m); for(i=0;i1) { i1_=-1; for(i_=1;i_=0. | //| ZNeeded - flag controlling whether the eigenvectors are | //| needed or not. If ZNeeded is equal to: | //| * 0, the eigenvectors are not needed; | //| * 1, the eigenvectors of a tridiagonal matrix | //| are multiplied by the square matrix Z. It is | //| used if the tridiagonal matrix is obtained by | //| the similarity transformation of a symmetric | //| matrix. | //| * 2, the eigenvectors of a tridiagonal matrix | //| replace matrix Z. | //| A, B - half-interval (A, B] to search eigenvalues in. | //| Z - if ZNeeded is equal to: | //| * 0, Z isn't used and remains unchanged; | //| * 1, Z contains the square matrix (array whose | //| indexes range within [0..N-1, 0..N-1]) which | //| reduces the given symmetric matrix to | //| tridiagonal form; | //| * 2, Z isn't used (but changed on the exit). | //| Output parameters: | //| D - array of the eigenvalues found. | //| Array whose index ranges within [0..M-1]. | //| M - number of eigenvalues found in the given | //| half-interval (M>=0). | //| Z - if ZNeeded is equal to: | //| * 0, doesn't contain any information; | //| * 1, contains the product of a given NxN matrix | //| Z (from the left) and NxM matrix of the | //| eigenvectors found (from the right). Array | //| whose indexes range within [0..N-1, 0..M-1]. | //| * 2, contains the matrix of the eigenvectors | //| found. Array whose indexes range within | //| [0..N-1, 0..M-1]. | //| Result: | //| True, if successful. In that case, M contains the number of | //| eigenvalues in the given half-interval (could be equal to 0),| //| D contains the eigenvalues, Z contains the eigenvectors (if | //| needed). It should be noted that the subroutine changes the | //| size of arrays D and Z. | //| False, if the bisection method subroutine wasn't able to find| //| the eigenvalues in the given interval or if the inverse | //| iteration subroutine wasn't able to find all the | //| corresponding eigenvectors. In that case, the eigenvalues and| //| eigenvectors are not returned, M is equal to 0. | //+------------------------------------------------------------------+ static bool CEigenVDetect::SMatrixTdEVDR(double &d[],double &e[],const int n, const int zneeded,const double a, const double b,int &m,CMatrixDouble &z) { //--- create variables bool result; int errorcode=0; int nsplit=0; int i=0; int j=0; int k=0; int cr=0; double v=0; int i_=0; int i1_=0; //--- create arrays int iblock[]; int isplit[]; int ifail[]; double d1[]; double e1[]; double w[]; //--- create matrix CMatrixDouble z2; CMatrixDouble z3; //--- initialization m=0; //--- check if(!CAp::Assert(zneeded>=0 && zneeded<=2,__FUNCTION__+": incorrect ZNeeded!")) return(false); //--- check if(b<=a) { m=0; //--- return result return(true); } //--- check if(n<=0) { m=0; //--- return result return(true); } //--- Copy D,E to D1, E1 ArrayResizeAL(d1,n+1); i1_=-1; for(i_=1;i_<=n;i_++) d1[i_]=d[i_+i1_]; //--- check if(n>1) { ArrayResizeAL(e1,n); i1_=-1; for(i_=1;i_=0. | //| ZNeeded - flag controlling whether the eigenvectors are | //| needed or not. If ZNeeded is equal to: | //| * 0, the eigenvectors are not needed; | //| * 1, the eigenvectors of a tridiagonal matrix | //| are multiplied by the square matrix Z. It is | //| used if the tridiagonal matrix is obtained by | //| the similarity transformation of a symmetric | //| matrix. | //| * 2, the eigenvectors of a tridiagonal matrix | //| replace matrix Z. | //| I1, I2 - index interval for searching (from I1 to I2). | //| 0 <= I1 <= I2 <= N-1. | //| Z - if ZNeeded is equal to: | //| * 0, Z isn't used and remains unchanged; | //| * 1, Z contains the square matrix (array whose | //| indexes range within [0..N-1, 0..N-1]) which | //| reduces the given symmetric matrix to | //| tridiagonal form; | //| * 2, Z isn't used (but changed on the exit). | //| Output parameters: | //| D - array of the eigenvalues found. | //| Array whose index ranges within [0..I2-I1]. | //| Z - if ZNeeded is equal to: | //| * 0, doesn't contain any information; | //| * 1, contains the product of a given NxN matrix | //| Z (from the left) and Nx(I2-I1) matrix of the | //| eigenvectors found (from the right). Array | //| whose indexes range within [0..N-1, 0..I2-I1].| //| * 2, contains the matrix of the eigenvalues | //| found. Array whose indexes range within | //| [0..N-1, 0..I2-I1]. | //| Result: | //| True, if successful. In that case, D contains the | //| eigenvalues, Z contains the eigenvectors (if needed). | //| It should be noted that the subroutine changes the size of | //| arrays D and Z. | //| False, if the bisection method subroutine wasn't able to find| //| the eigenvalues in the given interval or if the inverse | //| iteration subroutine wasn't able to find all the | //| corresponding eigenvectors. In that case, the eigenvalues and| //| eigenvectors are not returned. | //+------------------------------------------------------------------+ static bool CEigenVDetect::SMatrixTdEVDI(double &d[],double &e[],const int n, const int zneeded,const int i1, const int i2,CMatrixDouble &z) { //--- create variables bool result; int errorcode=0; int nsplit=0; int i=0; int j=0; int k=0; int m=0; int cr=0; double v=0; int i_=0; int i1_=0; //--- create arrays int iblock[]; int isplit[]; int ifail[]; double w[]; double d1[]; double e1[]; //--- create matrix CMatrixDouble z2; CMatrixDouble z3; //--- check if(!CAp::Assert((0<=i1 && i1<=i2) && i21) { ArrayResizeAL(e1,n); i1_=-1; for(i_=1;i_0, we have a pair of complex conjugate | //| numbers with positive and negative imaginary | //| parts: the first eigenvalue WR[i] + | //| + sqrt(-1)*WI[i]; the second eigenvalue | //| WR[i+1] + sqrt(-1)*WI[i+1]; | //| WI[i]>0 | //| WI[i+1] = -WI[i] < 0 | //| In that case, the eigenvector corresponding to | //| the first eigenvalue is located in i and i+1 | //| columns of matrices VL/VR (the column number i | //| contains the real part, and the column number | //| i+1 contains the imaginary part), and the vector | //| corresponding to the second eigenvalue is a | //| complex conjugate to the first vector. | //| Arrays whose indexes range within | //| [0..N-1, 0..N-1]. | //| Result: | //| True, if the algorithm has converged. | //| False, if the algorithm has not converged. | //| Note 1: | //| Some users may ask the following question: what if WI[N-1]>0?| //| WI[N] must contain an eigenvalue which is complex conjugate | //| to the N-th eigenvalue, but the array has only size N? | //| The answer is as follows: such a situation cannot occur | //| because the algorithm finds a pairs of eigenvalues, | //| therefore, if WI[i]>0, I is strictly less than N-1. | //| Note 2: | //| The algorithm performance depends on the value of the | //| internal parameter NS of the InternalSchurDecomposition | //| subroutine which defines the number of shifts in the QR | //| algorithm (similarly to the block width in block-matrix | //| algorithms of linear algebra). If you require maximum | //| performance on your machine, it is recommended to adjust | //| this parameter manually. | //| See also the InternalTREVC subroutine. | //| The algorithm is based on the LAPACK 3.0 library. | //+------------------------------------------------------------------+ static bool CEigenVDetect::RMatrixEVD(CMatrixDouble &ca,const int n,const int vneeded, double &wr[],double &wi[], CMatrixDouble &vl,CMatrixDouble &vr) { //--- create variables int i=0; int i_=0; int i1_=0; bool result; //--- create arrays double wr1[]; double wi1[]; //--- create matrix CMatrixDouble a1; CMatrixDouble vl1; CMatrixDouble vr1; //--- create copy CMatrixDouble a; a=ca; //--- check if(!CAp::Assert(vneeded>=0 && vneeded<=3,__FUNCTION__+": incorrect VNeeded!")) return(false); a1.Resize(n+1,n+1); for(i=1;i<=n;i++) { i1_=-1; for(i_=1;i_<=n;i_++) a1[i].Set(i_,a[i-1][i_+i1_]); } //--- get result result=NonSymmetricEVD(a1,n,vneeded,wr1,wi1,vl1,vr1); //--- check if(result) { //--- allocation ArrayResizeAL(wr,n); ArrayResizeAL(wi,n); i1_=1; for(i_=0;i_=0 && zneeded<=3,__FUNCTION__+": Incorrent ZNeeded")) return(false); //--- check if(zneeded<0 || zneeded>3) return(false); //--- initialization result=true; //--- check if(n==0) return(result); //--- check if(n==1) { //--- check if(zneeded==2 || zneeded==3) { z.Resize(2,2); z[1].Set(1,1); } //--- return result return(result); } //--- initialization maxit=30; //--- allocation ArrayResizeAL(wtemp,n+1); ArrayResizeAL(work1,n); ArrayResizeAL(work2,n); ArrayResizeAL(workc,n+1); ArrayResizeAL(works,n+1); //--- Determine the unit roundoff and over/underflow thresholds. eps=CMath::m_machineepsilon; eps2=CMath::Sqr(eps); safmin=CMath::m_minrealnumber; safmax=CMath::m_maxrealnumber; ssfmax=MathSqrt(safmax)/3; ssfmin=MathSqrt(safmin)/eps2; //--- Here we are using transposition to get rid of column operations wastranspose=false; zrows=0; //--- check if(zneeded==1) zrows=n; //--- check if(zneeded==2) zrows=n; //--- check if(zneeded==3) zrows=1; //--- check if(zneeded==1) { wastranspose=true; //--- function call CBlas::InplaceTranspose(z,1,n,1,n,wtemp); } //--- check if(zneeded==2) { wastranspose=true; z.Resize(n+1,n+1); for(i=1;i<=n;i++) { for(j=1;j<=n;j++) { //--- check if(i==j) z[i].Set(j,1); else z[i].Set(j,0); } } } //--- check if(zneeded==3) { wastranspose=false; z.Resize(2,n+1); for(j=1;j<=n;j++) { //--- check if(j==1) z[1].Set(j,1); else z[1].Set(j,0); } } //--- initialization nmaxit=n*maxit; jtot=0; //--- Determine where the matrix splits and choose QL or QR iteration //--- for each block, according to whether top or bottom diagonal //--- element is smaller. l1=1; nm1=n-1; while(true) { //--- check if(l1>n) break; //--- check if(l1>1) e[l1-1]=0; gotoflag=false; m=l1; //--- check if(l1<=nm1) { for(m=l1;m<=nm1;m++) { tst=MathAbs(e[m]); //--- check if(tst==0.0) { gotoflag=true; //--- break the cycle break; } //--- check if(tst<=MathSqrt(MathAbs(d[m]))*MathSqrt(MathAbs(d[m+1]))*eps) { e[m]=0; gotoflag=true; //--- break the cycle break; } } } //--- check if(!gotoflag) m=n; //--- change values l=l1; lsv=l; lend=m; lendsv=lend; l1=m+1; //--- check if(lend==l) continue; //--- Scale submatrix in rows and columns L to LEND if(l==lend) anorm=MathAbs(d[l]); else { anorm=MathMax(MathAbs(d[l])+MathAbs(e[l]),MathAbs(e[lend-1])+MathAbs(d[lend])); for(i=l+1;i<=lend-1;i++) anorm=MathMax(anorm,MathAbs(d[i])+MathAbs(e[i])+MathAbs(e[i-1])); } iscale=0; //--- check if(anorm==0.0) continue; //--- check if(anorm>(double)(ssfmax)) { iscale=1; tmp=ssfmax/anorm; tmpint=lend-1; for(i_=l;i_<=lend;i_++) d[i_]=tmp*d[i_]; for(i_=l;i_<=tmpint;i_++) e[i_]=tmp*e[i_]; } //--- check if(anorml) { //--- QL Iteration //--- Look for small subdiagonal element. while(true) { gotoflag=false; //--- check if(l!=lend) { lendm1=lend-1; for(m=l;m<=lendm1;m++) { tst=CMath::Sqr(MathAbs(e[m])); //--- check if(tst<=eps2*MathAbs(d[m])*MathAbs(d[m+1])+safmin) { gotoflag=true; //--- break the cycle break; } } } //--- check if(!gotoflag) m=lend; //--- check if(m0) { //--- function call TdEVDEv2(d[l],e[l],d[l+1],rt1,rt2,c,s); //--- change values work1[l]=c; work2[l]=s; workc[1]=work1[l]; works[1]=work2[l]; //--- check if(!wastranspose) CRotations::ApplyRotationsFromTheRight(false,1,zrows,l,l+1,workc,works,z,wtemp); else CRotations::ApplyRotationsFromTheLeft(false,l,l+1,1,zrows,workc,works,z,wtemp); } else //--- function call TdEVDE2(d[l],e[l],d[l+1],rt1,rt2); //--- change values d[l]=rt1; d[l+1]=rt2; e[l]=0; l=l+2; //--- check if(l<=lend) continue; //--- break the cycle break; } //--- check if(jtot==nmaxit) break; jtot=jtot+1; //--- Form shift. g=(d[l+1]-p)/(2*e[l]); //--- function call r=TdEVDPythag(g,1); g=d[m]-p+e[l]/(g+TdEVDExtSign(r,g)); s=1; c=1; p=0; //--- Inner loop mm1=m-1; for(i=mm1;i>=l;i--) { f=s*e[i]; b=c*e[i]; //--- function call CRotations::GenerateRotation(g,f,c,s,r); //--- check if(i!=m-1) e[i+1]=r; g=d[i+1]-p; r=(d[i]-g)*s+2*c*b; p=s*r; d[i+1]=g+p; g=c*r-b; //--- If eigenvectors are desired, then save rotations. if(zneeded>0) { work1[i]=c; work2[i]=-s; } } //--- If eigenvectors are desired, then apply saved rotations. if(zneeded>0) { for(i=l;i=lendp1;m--) { tst=CMath::Sqr(MathAbs(e[m-1])); //--- check if(tst<=(double)(eps2*MathAbs(d[m])*MathAbs(d[m-1])+safmin)) { gotoflag=true; //--- break the cycle break; } } } //--- check if(!gotoflag) m=lend; //--- check if(m>lend) e[m-1]=0; p=d[l]; //--- check if(m!=l) { //--- If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 //--- to compute its eigensystem. if(m==l-1) { //--- check if(zneeded>0) { //--- function call TdEVDEv2(d[l-1],e[l-1],d[l],rt1,rt2,c,s); work1[m]=c; work2[m]=s; workc[1]=c; works[1]=s; //--- check if(!wastranspose) { //--- function call CRotations::ApplyRotationsFromTheRight(true,1,zrows,l-1,l,workc,works,z,wtemp); } else { //--- function call CRotations::ApplyRotationsFromTheLeft(true,l-1,l,1,zrows,workc,works,z,wtemp); } } else { //--- function call TdEVDE2(d[l-1],e[l-1],d[l],rt1,rt2); } d[l-1]=rt1; d[l]=rt2; e[l-1]=0; l=l-2; //--- check if(l>=lend) continue; //--- break the cycle break; } //--- check if(jtot==nmaxit) break; jtot=jtot+1; //--- Form shift. g=(d[l-1]-p)/(2*e[l-1]); //--- function call r=TdEVDPythag(g,1); g=d[m]-p+e[l-1]/(g+TdEVDExtSign(r,g)); s=1; c=1; p=0; //--- Inner loop lm1=l-1; for(i=m;i<=lm1;i++) { f=s*e[i]; b=c*e[i]; //--- function call CRotations::GenerateRotation(g,f,c,s,r); //--- check if(i!=m) e[i-1]=r; //--- change values g=d[i]-p; r=(d[i+1]-g)*s+2*c*b; p=s*r; d[i]=g+p; g=c*r-b; //--- If eigenvectors are desired, then save rotations. if(zneeded>0) { work1[i]=c; work2[i]=s; } } //--- If eigenvectors are desired, then apply saved rotations. if(zneeded>0) { mm=l-m+1; for(i=m;i<=l-1;i++) { workc[i-m+1]=work1[i]; works[i-m+1]=work2[i]; } //--- check if(!wastranspose) { //--- function call CRotations::ApplyRotationsFromTheRight(true,1,zrows,m,l,workc,works,z,wtemp); } else { //--- function call CRotations::ApplyRotationsFromTheLeft(true,m,l,1,zrows,workc,works,z,wtemp); } } d[l]=d[l]-p; e[lm1]=g; continue; } //--- Eigenvalue found. d[l]=p; l=l-1; //--- check if(l>=lend) continue; //--- break the cycle break; } } //--- Undo scaling if necessary if(iscale==1) { tmp=anorm/ssfmax; tmpint=lendsv-1; for(i_=lsv;i_<=lendsv;i_++) d[i_]=tmp*d[i_]; for(i_=lsv;i_<=tmpint;i_++) e[i_]=tmp*e[i_]; } //--- check if(iscale==2) { tmp=anorm/ssfmin; tmpint=lendsv-1; for(i_=lsv;i_<=lendsv;i_++) d[i_]=tmp*d[i_]; for(i_=lsv;i_<=tmpint;i_++) e[i_]=tmp*e[i_]; } //--- Check for no convergence to an eigenvalue after a total //--- of N*MAXIT iterations. if(jtot>=nmaxit) { result=false; //--- check if(wastranspose) { //--- function call CBlas::InplaceTranspose(z,1,n,1,n,wtemp); } //--- return result return(result); } } //--- Order eigenvalues and eigenvectors. if(zneeded==0) { //--- Sort if(n==1) return(result); //--- check if(n==2) { //--- check if(d[1]>d[2]) { tmp=d[1]; d[1]=d[2]; d[2]=tmp; } //--- return result return(result); } i=2; do { t=i; while(t!=1) { k=t/2; //--- check if(d[k]>=d[t]) t=1; else { //--- change values tmp=d[k]; d[k]=d[t]; d[t]=tmp; t=k; } } i=i+1; } //--- cycle while(i<=n); i=n-1; do { tmp=d[i+1]; d[i+1]=d[1]; d[1]=tmp; t=1; while(t!=0) { k=2*t; //--- check if(k>i) t=0; else { //--- check if(kd[k]) k=k+1; } //--- check if(d[t]>=d[k]) t=0; else { //--- change values tmp=d[k]; d[k]=d[t]; d[t]=tmp; t=k; } } } i=i-1; } while(i>=1); } else { //--- Use Selection Sort to minimize swaps of eigenvectors for(ii=2;ii<=n;ii++) { i=ii-1; k=i; p=d[i]; for(j=ii;j<=n;j++) { //--- check if(d[j]MathAbs(c)) { acmx=a; acmn=c; } else { acmx=c; acmn=a; } //--- check if(adf>ab) { rt=adf*MathSqrt(1+CMath::Sqr(ab/adf)); } else { //--- check if(adf0.0) { rt1=0.5*(sm+rt); //--- Order of execution important. //--- To get fully accurate smaller eigenvalue, //--- next line needs to be executed in higher precision. rt2=acmx/rt1*acmn-b/rt1*b; } else { //--- Includes case RT1 = RT2 = 0 rt1=0.5*rt; rt2=-(0.5*rt); } } } //+------------------------------------------------------------------+ //| DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric | //| matrix | //| [ A B ] | //| [ B C ]. | //| On return, RT1 is the eigenvalue of larger absolute value, RT2 is| //| the eigenvalue of smaller absolute value, and (CS1,SN1) is the | //| unit right eigenvector for RT1, giving the decomposition | //| [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ] | //| [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]. | //| -- LAPACK auxiliary routine (version 3.0) -- | //| Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., | //| Courant Institute, Argonne National Lab, and Rice University| //| October 31, 1992 | //+------------------------------------------------------------------+ static void CEigenVDetect::TdEVDEv2(const double a,const double b,const double c, double &rt1,double &rt2,double &cs1,double &sn1) { //--- create variables int sgn1=0; int sgn2=0; double ab=0; double acmn=0; double acmx=0; double acs=0; double adf=0; double cs=0; double ct=0; double df=0; double rt=0; double sm=0; double tb=0; double tn=0; //--- initialization rt1=0; rt2=0; cs1=0; sn1=0; //--- Compute the eigenvalues sm=a+c; df=a-c; adf=MathAbs(df); tb=b+b; ab=MathAbs(tb); //--- check if(MathAbs(a)>MathAbs(c)) { acmx=a; acmn=c; } else { acmx=c; acmn=a; } //--- check if(adf>ab) rt=adf*MathSqrt(1+CMath::Sqr(ab/adf)); else { //--- check if(adf0.0) { rt1=0.5*(sm+rt); sgn1=1; //--- Order of execution important. //--- To get fully accurate smaller eigenvalue, //--- next line needs to be executed in higher precision. rt2=acmx/rt1*acmn-b/rt1*b; } else { //--- Includes case RT1 = RT2 = 0 rt1=0.5*rt; rt2=-(0.5*rt); sgn1=1; } } //--- Compute the eigenvector if(df>=0.0) { cs=df+rt; sgn2=1; } else { cs=df-rt; sgn2=-1; } acs=MathAbs(cs); //--- check if(acs>ab) { ct=-(tb/cs); sn1=1/MathSqrt(1+ct*ct); cs1=ct*sn1; } else { //--- check if(ab==0.0) { cs1=1; sn1=0; } else { tn=-(cs/tb); cs1=1/MathSqrt(1+tn*tn); sn1=tn*cs1; } } //--- check if(sgn1==sgn2) { tn=cs1; cs1=-sn1; sn1=tn; } } //+------------------------------------------------------------------+ //| Internal routine | //+------------------------------------------------------------------+ static double CEigenVDetect::TdEVDPythag(const double a,const double b) { //--- create variables double result=0; //--- check if(MathAbs(a)=0.0) result=MathAbs(a); else result=-MathAbs(a); //--- return result return(result); } //+------------------------------------------------------------------+ //| Internal subroutine | //+------------------------------------------------------------------+ static bool CEigenVDetect::InternalBisectionEigenValues(double &cd[],double &ce[], const int n,int irange, const int iorder,const double vl, const double vu,const int il, const int iu,const double abstol, double &w[],int &m, int &nsplit,int &iblock[], int &isplit[],int &errorcode) { //--- create variables bool result; double fudge=0; double relfac=0; bool ncnvrg; bool toofew; int ib=0; int ibegin=0; int idiscl=0; int idiscu=0; int ie=0; int iend=0; int iinfo=0; int im=0; int iin=0; int ioff=0; int iout=0; int itmax=0; int iw=0; int iwoff=0; int j=0; int itmp1=0; int jb=0; int jdisc=0; int je=0; int nwl=0; int nwu=0; int tmpi=0; double atoli=0; double bnorm=0; double gl=0; double gu=0; double pivmin=0; double rtoli=0; double safemn=0; double tmp1=0; double tmp2=0; double tnorm=0; double ulp=0; double wkill=0; double wl=0; double wlu=0; double wu=0; double wul=0; double scalefactor=0; double t=0; //--- create arrays int idumma[]; double work[]; int iwork[]; int ia1s2[]; double ra1s2[]; double ra1siin[]; double ra2siin[]; double ra3siin[]; double ra4siin[]; int iworkspace[]; double rworkspace[]; //--- create matrix CMatrixDouble ra1s2x2; CMatrixInt ia1s2x2; CMatrixDouble ra1siinx2; CMatrixInt ia1siinx2; //--- create copy double d[]; ArrayResizeAL(d,ArraySize(cd)); ArrayCopy(d,cd); //--- create copy double e[]; ArrayResizeAL(e,ArraySize(ce)); ArrayCopy(e,ce); //--- initialization m=0; nsplit=0; errorcode=0; //--- check if(n==0) return(true); //--- Get machine constants //--- NB is the minimum vector length for vector bisection, or 0 //--- if only scalar is to be done. fudge=2; relfac=2; safemn=CMath::m_minrealnumber; ulp=2*CMath::m_machineepsilon; rtoli=ulp*relfac; //--- allocation ArrayResizeAL(idumma,2); ArrayResizeAL(work,4*n+1); ArrayResizeAL(iwork,3*n+1); ArrayResizeAL(w,n+1); ArrayResizeAL(iblock,n+1); ArrayResizeAL(isplit,n+1); ArrayResizeAL(ia1s2,3); ArrayResizeAL(ra1s2,3); ArrayResizeAL(ra1siin,n+1); ArrayResizeAL(ra2siin,n+1); ArrayResizeAL(ra3siin,n+1); ArrayResizeAL(ra4siin,n+1); ArrayResizeAL(iworkspace,n+1); ArrayResizeAL(rworkspace,n+1); ra1siinx2.Resize(n+1,3); ia1siinx2.Resize(n+1,3); ra1s2x2.Resize(3,3); ia1s2x2.Resize(3,3); //--- initialization wlu=0; wul=0; //--- Check for Errors result=false; errorcode=0; //--- check if(irange<=0 || irange>=4) errorcode=-4; //--- check if(iorder<=0 || iorder>=3) errorcode=-5; //--- check if(n<0) errorcode=-3; //--- check if(irange==2 && vl>=vu) errorcode=-6; //--- check if(irange==3 && (il<1 || il>MathMax(1,n))) errorcode=-8; //--- check if(irange==3 && (iun)) errorcode=-9; //--- check if(errorcode!=0) return(result); //--- Initialize error flags ncnvrg=false; toofew=false; //--- Simplifications: if(irange==3 && il==1 && iu==n) irange=1; //--- Special Case when N=1 if(n==1) { nsplit=1; isplit[1]=1; //--- check if((irange==2 && vl>=d[1]) || vuMathSqrt(MathSqrt(CMath::m_minrealnumber))*MathSqrt(CMath::m_maxrealnumber)) scalefactor=t; //--- check if(ttmp1) { isplit[nsplit]=j-1; nsplit=nsplit+1; work[j-1]=0; } else { work[j-1]=tmp1; pivmin=MathMax(pivmin,tmp1); } } isplit[nsplit]=n; pivmin=pivmin*safemn; //--- Compute Interval and ATOLI if(irange==3) { //--- RANGE='I': Compute the interval containing eigenvalues //--- IL through IU. //--- Compute Gershgorin interval for entire (split) matrix //--- and use it as the initial interval gu=d[1]; gl=d[1]; tmp1=0; for(j=1;j=n || nwu<1 || nwu>n) { errorcode=4; return(false); } } else { //--- RANGE='A' or 'V' -- Set ATOLI tnorm=MathMax(MathAbs(d[1])+MathAbs(e[1]),MathAbs(d[n])+MathAbs(e[n-1])); for(j=2;j=d[ibegin]-pivmin) nwl=nwl+1; //--- check if(irange==1 || wu>=d[ibegin]-pivmin) nwu=nwu+1; //--- check if((irange==1 || wl=d[ibegin]-pivmin) { m=m+1; w[m]=d[ibegin]; iblock[m]=jb; } } else { //--- General Case -- IIN > 1 //--- Compute Gershgorin Interval //--- and use it as the initial interval gu=d[ibegin]; gl=d[ibegin]; tmp1=0; for(j=ibegin;j<=iend-1;j++) { //--- change values tmp2=MathAbs(e[j]); gu=MathMax(gu,d[j]+tmp1+tmp2); gl=MathMin(gl,d[j]-tmp1-tmp2); tmp1=tmp2; } //--- change values gu=MathMax(gu,d[iend]+tmp1); gl=MathMin(gl,d[iend]-tmp1); bnorm=MathMax(MathAbs(gl),MathAbs(gu)); gl=gl-fudge*bnorm*ulp*iin-fudge*pivmin; gu=gu+fudge*bnorm*ulp*iin+fudge*pivmin; //--- Compute ATOLI for the current submatrix if(abstol<=0.0) atoli=ulp*MathMax(MathAbs(gl),MathAbs(gu)); else atoli=abstol; //--- check if(irange>1) { //--- check if(gu=gu) continue; } //--- Set Up Initial Interval work[n+1]=gl; work[n+iin+1]=gu; //--- Calling DLAEBZ //--- CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, //--- D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), //--- IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM, //--- IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) for(tmpi=1;tmpi<=iin;tmpi++) { ra1siin[tmpi]=d[ibegin-1+tmpi]; //--- check if(ibegin-1+tmpiiout-iinfo) { ncnvrg=true; ib=-jb; } else ib=jb; for(je=iwork[j]+1+iwoff;je<=iwork[j+iin]+iwoff;je++) { w[je]=tmp1; iblock[je]=ib; } } m=m+im; } } //--- If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU //--- If NWL+1 < IL or NWU > IU, discard extra eigenvalues. if(irange==3) { im=0; idiscl=il-1-nwl; idiscu=nwu-iu; //--- check if(idiscl>0 || idiscu>0) { for(je=1;je<=m;je++) { //--- check if(w[je]<=wlu && idiscl>0) { idiscl=idiscl-1; } else { //--- check if(w[je]>=wul && idiscu>0) idiscu=idiscu-1; else { im=im+1; w[im]=w[je]; iblock[im]=iblock[je]; } } } m=im; } //--- check if(idiscl>0 || idiscu>0) { //--- Code to deal with effects of bad arithmetic: //--- Some low eigenvalues to be discarded are not in (WL,WLU], //--- or high eigenvalues to be discarded are not in (WUL,WU] //--- so just kill off the smallest IDISCL/largest IDISCU //--- eigenvalues, by simply finding the smallest/largest //--- eigenvalue(s). //--- (If N(w) is monotone non-decreasing, this should never //--- happen.) if(idiscl>0) { wkill=wu; for(jdisc=1;jdisc<=idiscl;jdisc++) { iw=0; for(je=1;je<=m;je++) { //--- check if(iblock[je]!=0 && (w[je]<(double)(wkill) || iw==0)) { iw=je; wkill=w[je]; } } iblock[iw]=0; } } //--- check if(idiscu>0) { wkill=wl; for(jdisc=1;jdisc<=idiscu;jdisc++) { iw=0; for(je=1;je<=m;je++) { //--- check if(iblock[je]!=0 && (w[je]>(double)(wkill) || iw==0)) { iw=je; wkill=w[je]; } } iblock[iw]=0; } } im=0; for(je=1;je<=m;je++) { //--- check if(iblock[je]!=0) { im=im+1; w[im]=w[je]; iblock[im]=iblock[je]; } } m=im; } //--- check if(idiscl<0 || idiscu<0) toofew=true; } //--- If ORDER='B', do nothing -- the eigenvalues are already sorted //--- by block. //--- If ORDER='E', sort the eigenvalues from smallest to largest if(iorder==1 && nsplit>1) { for(je=1;je<=m-1;je++) { ie=0; tmp1=w[je]; for(j=je+1;j<=m;j++) { //--- check if(w[j]n) { info=-4; return; } for(j=2;j<=m;j++) { //--- check if(iblock[j]1) { eps1=MathAbs(eps*xj); pertol=10*eps1; sep=xj-xjm; //--- check if(sepmaxits) { //--- If stopping criterion was not satisfied, update info and //--- store eigenvector number in array ifail. info=info+1; ifail[info]=j; break; } //--- Normalize and scale the righthand side vector Pb. v=0; for(ti=1;ti<=blksiz;ti++) v=v+MathAbs(work1[ti]); scl=blksiz*onenrm*MathMax(eps,MathAbs(work4[blksiz]))/v; for(i_=1;i_<=blksiz;i_++) work1[i_]=scl*work1[i_]; //--- Solve the system LU = Pb. TdIninternalDLAGTS(blksiz,work4,work2,work3,work5,iwork,work1,tol,iinfo); //--- Reorthogonalize by modified Gram-Schmidt if eigenvalues are //--- close enough. if(jblk!=1) { //--- check if(MathAbs(xj-xjm)>ortol) gpind=j; //--- check if(gpind!=j) { for(i=gpind;i1) tol=MathMax(tol,MathMax(MathAbs(a[2]),MathAbs(b[1]))); for(k=3;k<=n;k++) tol=MathMax(tol,MathMax(MathAbs(a[k]),MathMax(MathAbs(b[k-1]),MathAbs(d[k-2])))); tol=tol*eps; //--- check if(tol==0.0) tol=eps; } for(k=2;k<=n;k++) { //--- check if(iin[k-1]==0) y[k]=y[k]-c[k-1]*y[k-1]; else { temp=y[k-1]; y[k-1]=y[k]; y[k]=temp-c[k-1]*y[k]; } } for(k=n;k>=1;k--) { //--- check if(k<=n-2) temp=y[k]-b[k]*y[k+1]-d[k]*y[k+2]; else { //--- check if(k==n-1) temp=y[k]-b[k]*y[k+1]; else temp=y[k]; } ak=a[k]; pert=MathAbs(tol); //--- check if(ak<0.0) pert=-pert; while(true) { absak=MathAbs(ak); //--- check if(absak<1.0) { //--- check if(absakabsak) { ak=ak+pert; pert=2*pert; continue; } else { temp=temp*bignum; ak=ak*bignum; } } else { //--- check if(MathAbs(temp)>absak*bignum) { ak=ak+pert; pert=2*pert; continue; } } } //--- break the cycle break; } y[k]=temp/ak; } } //+------------------------------------------------------------------+ //| Internal subroutine | //+------------------------------------------------------------------+ static void CEigenVDetect::InternalDLAEBZ(const int ijob,const int nitmax, const int n,const int mmax,const int minp, const double abstol,const double reltol, const double pivmin,double &d[], double &e[],double &e2[],int &nval[], CMatrixDouble &ab,double &c[],int &mout, CMatrixInt &nab,double &work[], int &iwork[],int &info) { //--- create variables int itmp1=0; int itmp2=0; int j=0; int ji=0; int jit=0; int jp=0; int kf=0; int kfnew=0; int kl=0; int klnew=0; double tmp1=0; double tmp2=0; //--- initialization mout=0; info=0; //--- check if(ijob<1 || ijob>3) { info=-1; //--- exit the function return; } //--- Initialize NAB if(ijob==1) { //--- Compute the number of eigenvalues in the initial intervals. mout=0; //--- DIR$ NOVECTOR for(ji=1;ji<=minp;ji++) { for(jp=1;jp<=2;jp++) { tmp1=d[1]-ab[ji][jp]; //--- check if(MathAbs(tmp1)=nval[ji]) { ab[ji].Set(2,tmp1); nab[ji].Set(2,itmp1); } } } kl=klnew; //--- Check for convergence kfnew=kf; for(ji=kf;ji<=kl;ji++) { tmp1=MathAbs(ab[ji][2]-ab[ji][1]); tmp2=MathMax(MathAbs(ab[ji][2]),MathAbs(ab[ji][1])); //--- check if(tmp1<(double)(MathMax(abstol,MathMax(pivmin,reltol*tmp2))) || nab[ji][1]>=nab[ji][2]) { //--- Converged -- Swap with position KFNEW, //--- then increment KFNEW if(ji>kfnew) { tmp1=ab[ji][1]; tmp2=ab[ji][2]; itmp1=nab[ji][1]; itmp2=nab[ji][2]; //--- change values ab[ji].Set(1,ab[kfnew][1]); ab[ji].Set(2,ab[kfnew][2]); nab[ji].Set(1,nab[kfnew][1]); nab[ji].Set(2,nab[kfnew][2]); ab[kfnew].Set(1,tmp1); ab[kfnew].Set(2,tmp2); nab[kfnew].Set(1,itmp1); nab[kfnew].Set(2,itmp2); //--- check if(ijob==3) { itmp1=nval[ji]; nval[ji]=nval[kfnew]; nval[kfnew]=itmp1; } } kfnew=kfnew+1; } } kf=kfnew; //--- Choose Midpoints for(ji=kf;ji<=kl;ji++) c[ji]=0.5*(ab[ji][1]+ab[ji][2]); //--- If no more intervals to refine, quit. if(kf>kl) break; } //--- Converged info=(int)MathMax(kl+1-kf,0); mout=kl; } //+------------------------------------------------------------------+ //| Internal subroutine | //| -- LAPACK routine (version 3.0) -- | //| Univ. of Tennessee,Univ. of California Berkeley,NAG Ltd., | //| Courant Institute,Argonne National Lab, and Rice University | //| June 30,1999 | //+------------------------------------------------------------------+ static void CEigenVDetect::InternalTREVC(CMatrixDouble &t,const int n,const int side, const int howmny,bool &cvselect[],CMatrixDouble &vl, CMatrixDouble &vr,int &m,int &info) { //--- create variables bool allv; bool bothv; bool leftv; bool over; bool pair; bool rightv; bool somev; int i=0; int ierr=0; int ii=0; int ip=0; int iis=0; int j=0; int j1=0; int j2=0; int jnxt=0; int k=0; int ki=0; int n2=0; double beta=0; double bignum=0; double emax=0; double ovfl=0; double rec=0; double remax=0; double scl=0; double smin=0; double smlnum=0; double ulp=0; double unfl=0; double vcrit=0; double vmax=0; double wi=0; double wr=0; double xnorm=0; bool skipflag; int k1=0; int k2=0; int k3=0; int k4=0; double vt=0; int i_=0; int i1_=0; //--- create arrays double work[]; double temp[]; bool rswap4[]; bool zswap4[]; double civ4[]; double crv4[]; //--- create matrix CMatrixDouble x; CMatrixDouble temp11; CMatrixDouble temp22; CMatrixDouble temp11b; CMatrixDouble temp21b; CMatrixDouble temp12b; CMatrixDouble temp22b; CMatrixInt ipivot44; //--- create copy double vselect[]; ArrayResizeAL(vselect,ArraySize(cvselect)); ArrayCopy(vselect,cvselect); //--- initialization m=0; info=0; //--- allocation x.Resize(3,3); temp11.Resize(2,2); temp11b.Resize(2,2); temp21b.Resize(3,2); temp12b.Resize(2,3); temp22b.Resize(3,3); temp22.Resize(3,3); ArrayResizeAL(work,3*n+1); ArrayResizeAL(temp,n+1); ArrayResizeAL(rswap4,5); ArrayResizeAL(zswap4,5); ArrayResizeAL(civ4,5); ArrayResizeAL(crv4,5); ipivot44.Resize(5,5); //--- check if(howmny!=1) { //--- check if(side==1 || side==3) vr.Resize(n+1,n+1); //--- check if(side==2 || side==3) vl.Resize(n+1,n+1); } //--- Decode and test the input parameters bothv=side==3; rightv=side==1 || bothv; leftv=side==2 || bothv; allv=howmny==2; over=howmny==1; somev=howmny==3; info=0; //--- check if(n<0) { info=-2; return; } //--- check if(!rightv && !leftv) { info=-3; return; } //--- check if((!allv && !over) && !somev) { info=-4; return; } //--- Set M to the number of columns required to store the selected //--- eigenvectors, standardize the array SELECT if necessary, and //--- test MM. if(somev) { m=0; pair=false; for(j=1;j<=n;j++) { //--- check if(pair) { pair=false; vselect[j]=false; } else { //--- check if(j=1;ki--) { skipflag=false; //--- check if(ip==1) skipflag=true; else { //--- check if(ki!=1) { //--- check if(t[ki][ki-1]!=0.0) ip=-1; } //--- check if(somev) { //--- check if(ip==0) { //--- check if(!vselect[ki]) skipflag=true; } else { //--- check if(!vselect[ki-1]) skipflag=true; } } } //--- check if(!skipflag) { //--- Compute the KI-th eigenvalue (WR,WI). wr=t[ki][ki]; wi=0; //--- check if(ip!=0) wi=MathSqrt(MathAbs(t[ki][ki-1]))*MathSqrt(MathAbs(t[ki-1][ki])); smin=MathMax(ulp*(MathAbs(wr)+MathAbs(wi)),smlnum); //--- check if(ip==0) { //--- Real right eigenvector work[ki+n]=1; //--- Form right-hand side for(k=1;k<=ki-1;k++) work[k+n]=-t[k][ki]; //--- Solve the upper quasi-triangular system: //--- (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK. jnxt=ki-1; for(j=ki-1;j>=1;j--) { //--- check if(j>jnxt) continue; j1=j; j2=j; jnxt=j-1; //--- check if(j>1) { //--- check if(t[j][j-1]!=0.0) { j1=j-1; jnxt=j-2; } } //--- check if(j1==j2) { //--- 1-by-1 diagonal block temp11[1].Set(1,t[j][j]); temp11b[1].Set(1,work[j+n]); //--- function call InternalHsEVDLALN2(false,1,1,smin,1,temp11,1.0,1.0,temp11b,wr,0.0,rswap4,zswap4,ipivot44,civ4,crv4,x,scl,xnorm,ierr); //--- Scale X(1,1) to avoid overflow when updating //--- the right-hand side. if(xnorm>1.0) { //--- check if(work[j]>bignum/xnorm) { x[1].Set(1,x[1][1]/xnorm); scl=scl/xnorm; } } //--- Scale if necessary if(scl!=1.0) { k1=n+1; k2=n+ki; for(i_=k1;i_<=k2;i_++) work[i_]=scl*work[i_]; } work[j+n]=x[1][1]; //--- Update right-hand side k1=1+n; k2=j-1+n; k3=j-1; vt=-x[1][1]; i1_=1-k1; for(i_=k1;i_<=k2;i_++) { work[i_]=work[i_]+vt*t[i_+i1_][j]; } } else { //--- 2-by-2 diagonal block temp22[1].Set(1,t[j-1][j-1]); temp22[1].Set(2,t[j-1][j]); temp22[2].Set(1,t[j][j-1]); temp22[2].Set(2,t[j][j]); temp21b[1].Set(1,work[j-1+n]); temp21b[2].Set(1,work[j+n]); //--- function call InternalHsEVDLALN2(false,2,1,smin,1.0,temp22,1.0,1.0,temp21b,wr,0,rswap4,zswap4,ipivot44,civ4,crv4,x,scl,xnorm,ierr); //--- Scale X(1,1) and X(2,1) to avoid overflow when //--- updating the right-hand side. if(xnorm>1.0) { beta=MathMax(work[j-1],work[j]); //--- check if(beta>bignum/xnorm) { x[1].Set(1,x[1][1]/xnorm); x[2].Set(1,x[2][1]/xnorm); scl=scl/xnorm; } } //--- Scale if necessary if(scl!=1.0) { k1=1+n; k2=ki+n; for(i_=k1;i_<=k2;i_++) work[i_]=scl*work[i_]; } work[j-1+n]=x[1][1]; work[j+n]=x[2][1]; //--- Update right-hand side k1=1+n; k2=j-2+n; k3=j-2; k4=j-1; vt=-x[1][1]; i1_=1-k1; for(i_=k1;i_<=k2;i_++) work[i_]=work[i_]+vt*t[i_+i1_][k4]; vt=-x[2][1]; i1_=1-k1; for(i_=k1;i_<=k2;i_++) work[i_]=work[i_]+vt*t[i_+i1_][j]; } } //--- Copy the vector x or Q*x to VR and normalize. if(!over) { k1=1+n; k2=ki+n; i1_=k1-1; for(i_=1;i_<=ki;i_++) vr[i_].Set(iis,work[i_+i1_]); //--- function call ii=CBlas::ColumnIdxAbsMax(vr,1,ki,iis); remax=1/MathAbs(vr[ii][iis]); for(i_=1;i_<=ki;i_++) vr[i_].Set(iis,remax*vr[i_][iis]); for(k=ki+1;k<=n;k++) vr[k].Set(iis,0); } else { //--- check if(ki>1) { for(i_=1;i_<=n;i_++) temp[i_]=vr[i_][ki]; //--- function call CBlas::MatrixVectorMultiply(vr,1,n,1,ki-1,false,work,1+n,ki-1+n,1.0,temp,1,n,work[ki+n]); for(i_=1;i_<=n;i_++) vr[i_].Set(ki,temp[i_]); } //--- function call ii=CBlas::ColumnIdxAbsMax(vr,1,n,ki); remax=1/MathAbs(vr[ii][ki]); for(i_=1;i_<=n;i_++) vr[i_].Set(ki,remax*vr[i_][ki]); } } else { //--- Complex right eigenvector. //--- Initial solve //--- [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0. //--- [ (T(KI,KI-1) T(KI,KI) ) ] if(MathAbs(t[ki-1][ki])>=MathAbs(t[ki][ki-1])) { work[ki-1+n]=1; work[ki+n2]=wi/t[ki-1][ki]; } else { work[ki-1+n]=-(wi/t[ki][ki-1]); work[ki+n2]=1; } work[ki+n]=0; work[ki-1+n2]=0; //--- Form right-hand side for(k=1;k<=ki-2;k++) { work[k+n]=-(work[ki-1+n]*t[k][ki-1]); work[k+n2]=-(work[ki+n2]*t[k][ki]); } //--- Solve upper quasi-triangular system: //--- (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2) jnxt=ki-2; for(j=ki-2;j>=1;j--) { //--- check if(j>jnxt) continue; j1=j; j2=j; jnxt=j-1; //--- check if(j>1) { //--- check if(t[j][j-1]!=0.0) { j1=j-1; jnxt=j-2; } } //--- check if(j1==j2) { //--- 1-by-1 diagonal block temp11[1].Set(1,t[j][j]); temp12b[1].Set(1,work[j+n]); temp12b[1].Set(2,work[j+n+n]); //--- function call InternalHsEVDLALN2(false,1,2,smin,1.0,temp11,1.0,1.0,temp12b,wr,wi,rswap4,zswap4,ipivot44,civ4,crv4,x,scl,xnorm,ierr); //--- Scale X(1,1) and X(1,2) to avoid overflow when //--- updating the right-hand side. if(xnorm>1.0) { //--- check if(work[j]>bignum/xnorm) { x[1].Set(1,x[1][1]/xnorm); x[1].Set(2,x[1][2]/xnorm); scl=scl/xnorm; } } //--- Scale if necessary if(scl!=1.0) { k1=1+n; k2=ki+n; for(i_=k1;i_<=k2;i_++) work[i_]=scl*work[i_]; k1=1+n2; k2=ki+n2; for(i_=k1;i_<=k2;i_++) work[i_]=scl*work[i_]; } work[j+n]=x[1][1]; work[j+n2]=x[1][2]; //--- Update the right-hand side k1=1+n; k2=j-1+n; k3=1; k4=j-1; vt=-x[1][1]; i1_=k3-k1; for(i_=k1;i_<=k2;i_++) work[i_]=work[i_]+vt*t[i_+i1_][j]; //--- change values k1=1+n2; k2=j-1+n2; k3=1; k4=j-1; vt=-x[1][2]; i1_=k3-k1; for(i_=k1;i_<=k2;i_++) work[i_]=work[i_]+vt*t[i_+i1_][j]; } else { //--- 2-by-2 diagonal block temp22[1].Set(1,t[j-1][j-1]); temp22[1].Set(2,t[j-1][j]); temp22[2].Set(1,t[j][j-1]); temp22[2].Set(2,t[j][j]); temp22b[1].Set(1,work[j-1+n]); temp22b[1].Set(2,work[j-1+n+n]); temp22b[2].Set(1,work[j+n]); temp22b[2].Set(2,work[j+n+n]); //--- function call InternalHsEVDLALN2(false,2,2,smin,1.0,temp22,1.0,1.0,temp22b,wr,wi,rswap4,zswap4,ipivot44,civ4,crv4,x,scl,xnorm,ierr); //--- Scale X to avoid overflow when updating //--- the right-hand side. if(xnorm>1.0) { beta=MathMax(work[j-1],work[j]); //--- check if(beta>bignum/xnorm) { rec=1/xnorm; x[1].Set(1,x[1][1]*rec); x[1].Set(2,x[1][2]*rec); x[2].Set(1,x[2][1]*rec); x[2].Set(2,x[2][2]*rec); scl=scl*rec; } } //--- Scale if necessary if(scl!=1.0) { for(i_=1+n;i_<=ki+n;i_++) work[i_]=scl*work[i_]; for(i_=1+n2;i_<=ki+n2;i_++) work[i_]=scl*work[i_]; } //--- change values work[j-1+n]=x[1][1]; work[j+n]=x[2][1]; work[j-1+n2]=x[1][2]; work[j+n2]=x[2][2]; //--- Update the right-hand side vt=-x[1][1]; i1_=-n; for(i_=n+1;i_<=n+j-2;i_++) work[i_]=work[i_]+vt*t[i_+i1_][j-1]; vt=-x[2][1]; i1_=-n; for(i_=n+1;i_<=n+j-2;i_++) work[i_]=work[i_]+vt*t[i_+i1_][j]; vt=-x[1][2]; i1_=-n2; for(i_=n2+1;i_<=n2+j-2;i_++) work[i_]=work[i_]+vt*t[i_+i1_][j-1]; vt=-x[2][2]; i1_=-n2; for(i_=n2+1;i_<=n2+j-2;i_++) work[i_]=work[i_]+vt*t[i_+i1_][j]; } } //--- Copy the vector x or Q*x to VR and normalize. if(!over) { i1_=n; for(i_=1;i_<=ki;i_++) vr[i_].Set(iis-1,work[i_+i1_]); i1_=n2; for(i_=1;i_<=ki;i_++) vr[i_].Set(iis,work[i_+i1_]); emax=0; for(k=1;k<=ki;k++) emax=MathMax(emax,MathAbs(vr[k][iis-1])+MathAbs(vr[k][iis])); remax=1/emax; //--- copy for(i_=1;i_<=ki;i_++) vr[i_].Set(iis-1,remax*vr[i_][iis-1]); for(i_=1;i_<=ki;i_++) vr[i_].Set(iis,remax*vr[i_][iis]); for(k=ki+1;k<=n;k++) vr[k].Set(iis-1,0); vr[k].Set(iis,0); } else { //--- check if(ki>2) { for(i_=1;i_<=n;i_++) temp[i_]=vr[i_][ki-1]; //--- function call CBlas::MatrixVectorMultiply(vr,1,n,1,ki-2,false,work,1+n,ki-2+n,1.0,temp,1,n,work[ki-1+n]); for(i_=1;i_<=n;i_++) vr[i_].Set(ki-1,temp[i_]); for(i_=1;i_<=n;i_++) temp[i_]=vr[i_][ki]; //--- function call CBlas::MatrixVectorMultiply(vr,1,n,1,ki-2,false,work,1+n2,ki-2+n2,1.0,temp,1,n,work[ki+n2]); for(i_=1;i_<=n;i_++) vr[i_].Set(ki,temp[i_]); } else { vt=work[ki-1+n]; //--- copy for(i_=1;i_<=n;i_++) vr[i_].Set(ki-1,vt*vr[i_][ki-1]); vt=work[ki+n2]; for(i_=1;i_<=n;i_++) vr[i_].Set(ki,vt*vr[i_][ki]); } emax=0; for(k=1;k<=n;k++) emax=MathMax(emax,MathAbs(vr[k][ki-1])+MathAbs(vr[k][ki])); remax=1/emax; //--- copy for(i_=1;i_<=n;i_++) vr[i_].Set(ki-1,remax*vr[i_][ki-1]); for(i_=1;i_<=n;i_++) vr[i_].Set(ki,remax*vr[i_][ki]); } } iis=iis-1; //--- check if(ip!=0) iis=iis-1; } //--- check if(ip==1) ip=0; //--- check if(ip==-1) ip=1; } } //--- check if(leftv) { //--- Compute left eigenvectors. ip=0; iis=1; for(ki=1;ki<=n;ki++) { skipflag=false; //--- check if(ip==-1) skipflag=true; else { //--- check if(ki!=n) { //--- check if(t[ki+1][ki]!=0.0) ip=1; } //--- check if(somev) { //--- check if(!vselect[ki]) skipflag=true; } } //--- check if(!skipflag) { //--- Compute the KI-th eigenvalue (WR,WI). wr=t[ki][ki]; wi=0; //--- check if(ip!=0) wi=MathSqrt(MathAbs(t[ki][ki+1]))*MathSqrt(MathAbs(t[ki+1][ki])); smin=MathMax(ulp*(MathAbs(wr)+MathAbs(wi)),smlnum); //--- check if(ip==0) { //--- Real left eigenvector. work[ki+n]=1; //--- Form right-hand side for(k=ki+1;k<=n;k++) work[k+n]=-t[ki][k]; //--- Solve the quasi-triangular system: //--- (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK vmax=1; vcrit=bignum; jnxt=ki+1; for(j=ki+1;j<=n;j++) { //--- check if(jvcrit) { rec=1/vmax; for(i_=ki+n;i_<=n+n;i_++) { work[i_]=rec*work[i_]; } vmax=1; vcrit=bignum; } i1_=n; vt=0.0; for(i_=ki+1;i_vcrit) { rec=1/vmax; for(i_=ki+n;i_<=n+n;i_++) work[i_]=rec*work[i_]; vmax=1; vcrit=bignum; } i1_=n; vt=0.0; for(i_=ki+1;i_=MathAbs(t[ki+1][ki])) { work[ki+n]=wi/t[ki][ki+1]; work[ki+1+n2]=1; } else { work[ki+n]=1; work[ki+1+n2]=-(wi/t[ki+1][ki]); } work[ki+1+n]=0; work[ki+n2]=0; //--- Form right-hand side for(k=ki+2;k<=n;k++) { work[k+n]=-(work[ki+n]*t[ki][k]); work[k+n2]=-(work[ki+1+n2]*t[ki+1][k]); } //--- Solve complex quasi-triangular system: //--- ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2 vmax=1; vcrit=bignum; jnxt=ki+2; for(j=ki+2;j<=n;j++) { //--- check if(jvcrit) { rec=1/vmax; for(i_=ki+n;i_<=n+n;i_++) work[i_]=rec*work[i_]; for(i_=ki+n2;i_<=n+n2;i_++) work[i_]=rec*work[i_]; vmax=1; vcrit=bignum; } i1_=n; vt=0.0; for(i_=ki+2;i_vcrit) { rec=1/vmax; for(i_=ki+n;i_<=n+n;i_++) work[i_]=rec*work[i_]; for(i_=ki+n2;i_<=n+n2;i_++) work[i_]=rec*work[i_]; vmax=1; vcrit=bignum; } i1_=n; vt=0.0; for(i_=ki+2;i_1.0) { //--- check if(bnorm>bignum*cnorm) scl=1/bnorm; } //--- Compute X x[1].Set(1,b[1][1]*scl/csr); xnorm=MathAbs(x[1][1]); } else { //--- Complex 1x1 system (w is complex) //--- C = ca A - w D csr=ca*a[1][1]-wr*d1; csi=-(wi*d1); cnorm=MathAbs(csr)+MathAbs(csi); //--- If | C | < SMINI, use C = SMINI if(cnorm1.0) { //--- check if(bnorm>bignum*cnorm) scl=1/bnorm; } //--- Compute X InternalHsEVDLADIV(scl*b[1][1],scl*b[1][2],csr,csi,tmp1,tmp2); x[1].Set(1,tmp1); x[1].Set(2,tmp2); xnorm=MathAbs(x[1][1])+MathAbs(x[1][2]); } } else { //--- 2x2 System //--- Compute the real part of C = ca A - w D (or ca A' - w D ) crv4[1+0]=ca*a[1][1]-wr*d1; crv4[2+2]=ca*a[2][2]-wr*d2; //--- check if(ltrans) { crv4[1+2]=ca*a[2][1]; crv4[2+0]=ca*a[1][2]; } else { crv4[2+0]=ca*a[2][1]; crv4[1+2]=ca*a[1][2]; } //--- check if(nw==1) { //--- Real 2x2 system (w is real) //--- Find the largest element in C cmax=0; icmax=0; for(j=1;j<=4;j++) { //--- check if(MathAbs(crv4[j])>cmax) { cmax=MathAbs(crv4[j]); icmax=j; } } //--- If norm(C) < SMINI, use SMINI*identity. if(cmax1.0) { //--- check if(bnorm>bignum*smini) scl=1/bnorm; } //--- change values temp=scl/smini; x[1].Set(1,temp*b[1][1]); x[2].Set(1,temp*b[2][1]); xnorm=temp*bnorm; info=1; //--- exit the function return; } //--- Gaussian elimination with complete pivoting. ur11=crv4[icmax]; cr21=crv4[ipivot44[2][icmax]]; ur12=crv4[ipivot44[3][icmax]]; cr22=crv4[ipivot44[4][icmax]]; ur11r=1/ur11; lr21=ur11r*cr21; ur22=cr22-ur12*lr21; //--- If smaller pivot < SMINI, use SMINI if(MathAbs(ur22)1.0 && MathAbs(ur22)<1.0) { //--- check if(bbnd>=bignum*MathAbs(ur22)) scl=1/bbnd; } xr2=br2*scl/ur22; xr1=scl*br1*ur11r-xr2*(ur11r*ur12); //--- check if(zswap4[icmax]) { x[1].Set(1,xr2); x[2].Set(1,xr1); } else { x[1].Set(1,xr1); x[2].Set(1,xr2); } xnorm=MathMax(MathAbs(xr1),MathAbs(xr2)); //--- Further scaling if norm(A) norm(X) > overflow if(xnorm>1.0 && cmax>1.0) { //--- check if(xnorm>bignum/cmax) { temp=cmax/bignum; x[1].Set(1,temp*x[1][1]); x[2].Set(1,temp*x[2][1]); xnorm=temp*xnorm; scl=temp*scl; } } } else { //--- Complex 2x2 system (w is complex) //--- Find the largest element in C civ4[1+0]=-(wi*d1); civ4[2+0]=0; civ4[1+2]=0; civ4[2+2]=-(wi*d2); cmax=0; icmax=0; for(j=1;j<=4;j++) { //--- check if(MathAbs(crv4[j])+MathAbs(civ4[j])>cmax) { cmax=MathAbs(crv4[j])+MathAbs(civ4[j]); icmax=j; } } //--- If norm(C) < SMINI, use SMINI*identity. if(cmax1.0) { //--- check if(bnorm>bignum*smini) scl=1/bnorm; } //--- change values temp=scl/smini; x[1].Set(1,temp*b[1][1]); x[2].Set(1,temp*b[2][1]); x[1].Set(2,temp*b[1][2]); x[2].Set(2,temp*b[2][2]); xnorm=temp*bnorm; info=1; //--- exit the function return; } //--- Gaussian elimination with complete pivoting. ur11=crv4[icmax]; ui11=civ4[icmax]; cr21=crv4[ipivot44[2][icmax]]; ci21=civ4[ipivot44[2][icmax]]; ur12=crv4[ipivot44[3][icmax]]; ui12=civ4[ipivot44[3][icmax]]; cr22=crv4[ipivot44[4][icmax]]; ci22=civ4[ipivot44[4][icmax]]; //--- check if(icmax==1 || icmax==4) { //--- Code when off-diagonals of pivoted C are real if(MathAbs(ur11)>MathAbs(ui11)) { temp=ui11/ur11; ur11r=1/(ur11*(1+CMath::Sqr(temp))); ui11r=-(temp*ur11r); } else { temp=ur11/ui11; ui11r=-(1/(ui11*(1+CMath::Sqr(temp)))); ur11r=-(temp*ui11r); } //--- change values lr21=cr21*ur11r; li21=cr21*ui11r; ur12s=ur12*ur11r; ui12s=ur12*ui11r; ur22=cr22-ur12*lr21; ui22=ci22-ur12*li21; } else { //--- Code when diagonals of pivoted C are real ur11r=1/ur11; ui11r=0; lr21=cr21*ur11r; li21=ci21*ur11r; ur12s=ur12*ur11r; ui12s=ui12*ur11r; ur22=cr22-ur12*lr21+ui12*li21; ui22=-(ur12*li21)-ui12*lr21; } u22abs=MathAbs(ur22)+MathAbs(ui22); //--- If smaller pivot < SMINI, use SMINI if(u22abs1.0 && u22abs<1.0) { //--- check if(bbnd>=bignum*u22abs) { //--- change values scl=1/bbnd; br1=scl*br1; bi1=scl*bi1; br2=scl*br2; bi2=scl*bi2; } } //--- function call InternalHsEVDLADIV(br2,bi2,ur22,ui22,xr2,xi2); xr1=ur11r*br1-ui11r*bi1-ur12s*xr2+ui12s*xi2; xi1=ui11r*br1+ur11r*bi1-ui12s*xr2-ur12s*xi2; //--- check if(zswap4[icmax]) { x[1].Set(1,xr2); x[2].Set(1,xr1); x[1].Set(2,xi2); x[2].Set(2,xi1); } else { x[1].Set(1,xr1); x[2].Set(1,xr2); x[1].Set(2,xi1); x[2].Set(2,xi2); } xnorm=MathMax(MathAbs(xr1)+MathAbs(xi1),MathAbs(xr2)+MathAbs(xi2)); //--- Further scaling if norm(A) norm(X) > overflow if(xnorm>1.0 && cmax>1.0) { //--- check if(xnorm>bignum/cmax) { //--- change values temp=cmax/bignum; x[1].Set(1,temp*x[1][1]); x[2].Set(1,temp*x[2][1]); x[1].Set(2,temp*x[1][2]); x[2].Set(2,temp*x[2][2]); xnorm=temp*xnorm; scl=temp*scl; } } } } } //+------------------------------------------------------------------+ //| performs complex division in real arithmetic | //| a + i*b | //| p + i*q = --------- | //| c + i*d | //| The algorithm is due to Robert L. Smith and can be found | //| in D. Knuth, The art of Computer Programming, Vol.2, p.195 | //| -- LAPACK auxiliary routine (version 3.0) -- | //| Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., | //| Courant Institute, Argonne National Lab, and Rice University| //| October 31, 1992 | //+------------------------------------------------------------------+ static void CEigenVDetect::InternalHsEVDLADIV(const double a,const double b, const double c,const double d, double &p,double &q) { //--- create variables double e=0; double f=0; //--- initialization p=0; q=0; //--- check if(MathAbs(d)=0 && vneeded<=3,__FUNCTION__+": incorrect VNeeded!")) return(false); //--- check if(vneeded==0) { //--- Eigen values only ToUpperHessenberg(a,n,tau); //--- function call CHsSchur::InternalSchurDecomposition(a,n,0,0,wr,wi,s,info); //--- get result result=info==0; //--- return result return(result); } //--- Eigen values and vectors ToUpperHessenberg(a,n,tau); //--- function call UnpackQFromUpperHessenberg(a,n,tau,s); //--- function call CHsSchur::InternalSchurDecomposition(a,n,1,1,wr,wi,s,info); //--- get result result=info==0; //--- check if(!result) return(result); //--- check if(vneeded==1 || vneeded==3) { vr.Resize(n+1,n+1); for(i=1;i<=n;i++) { for(i_=1;i_<=n;i_++) vr[i].Set(i_,s[i][i_]); } } //--- check if(vneeded==2 || vneeded==3) { vl.Resize(n+1,n+1); for(i=1;i<=n;i++) { for(i_=1;i_<=n;i_++) vl[i].Set(i_,s[i][i_]); } } //--- function call InternalTREVC(a,n,vneeded,1,sel,vl,vr,m,info); //--- get result result=info==0; //--- return result return(result); } //+------------------------------------------------------------------+ //| Upper Hessenberg form | //+------------------------------------------------------------------+ static void CEigenVDetect::ToUpperHessenberg(CMatrixDouble &a,const int n,double &tau[]) { //--- create variables int i=0; int ip1=0; int nmi=0; double v=0; int i_=0; int i1_=0; //--- create arrays double t[]; double work[]; //--- check if(!CAp::Assert(n>=0,__FUNCTION__+": incorrect N!")) return; //--- check if(n<=1) return; //--- allocation ArrayResizeAL(tau,n); ArrayResizeAL(t,n+1); ArrayResizeAL(work,n+1); //--- calculations for(i=1;i=1 | //| OUTPUT PARAMETERS: | //| A - orthogonal NxN matrix, array[0..N-1,0..N-1] | //+------------------------------------------------------------------+ static void CMatGen::RMatrixRndOrthogonal(const int n,CMatrixDouble &a) { //--- create variables int i=0; int j=0; //--- check if(!CAp::Assert(n>=1,__FUNCTION__+": N<1!")) return; //--- allocation a.Resize(n,n); for(i=0;i=1&&c>=1.0,__FUNCTION__+": N<1 or C<1!")) return; //--- allocation a.Resize(n,n); //--- check if(n==1) { a[0].Set(0,2*CMath::RandomInteger(2)-1); //--- exit the function return; } //--- initialization l1=0; l2=MathLog(1/c); for(i=0;i=1 | //| OUTPUT PARAMETERS: | //| A - orthogonal NxN matrix, array[0..N-1,0..N-1] | //+------------------------------------------------------------------+ static void CMatGen::CMatrixRndOrthogonal(const int n,CMatrixComplex &a) { //--- create variables int i=0; int j=0; //--- check if(!CAp::Assert(n>=1,__FUNCTION__+": N<1!")) return; //--- allocation a.Resize(n,n); for(i=0;i=1&&c>=1.0,__FUNCTION__+": N<1 or C<1!")) return; //--- allocation a.Resize(n,n); //--- check if(n==1) { //--- function call CHighQualityRand::HQRndRandomize(state); //--- function call CHighQualityRand::HQRndUnit2(state,v.re,v.im); a[0].Set(0,v); return; } //--- initialization l1=0; l2=MathLog(1/c); for(i=0;i=1&&c>=1.0,__FUNCTION__+": N<1 or C<1!")) return; //--- allocation a.Resize(n,n); //--- check if(n==1) { a[0].Set(0,2*CMath::RandomInteger(2)-1); return; } //--- Prepare matrix l1=0; l2=MathLog(1/c); for(i=0;i=1&&c>=1.0,__FUNCTION__+": N<1 or C<1!")) return; //--- allocation a.Resize(n,n); //--- check if(n==1) { a[0].Set(0,2*CMath::RandomInteger(2)-1); return; } //--- Prepare matrix l1=0; l2=MathLog(1/c); for(i=0;i=1 && m>=1,__FUNCTION__+": N<1 or M<1!")) return; //--- check if(n==1) { tau=2*CMath::RandomInteger(2)-1; for(i=0;i=1 && m>=1,__FUNCTION__+": N<1 or M<1!")) return; //--- check if(m==1) { tau=2*CMath::RandomInteger(2)-1; for(j=0;j=1 && m>=1,__FUNCTION__+": N<1 or M<1!")) return; //--- check if(n==1) { //--- function call CHighQualityRand::HQRndRandomize(state); //--- function call CHighQualityRand::HQRndUnit2(state,tau.re,tau.im); for(i=0;i=1 && m>=1,__FUNCTION__+": N<1 or M<1!")) return; //--- check if(m==1) { //--- function call CHighQualityRand::HQRndRandomize(state); //--- function call CHighQualityRand::HQRndUnit2(state,tau.re,tau.im); for(j=0;j>N | //| * worst case - N>>M, small M, large N, matrix does not fit in CPU| //| cache | //| INPUT PARAMETERS: | //| A - array[0..M-1, 0..N-1]. | //| M - number of rows in matrix A. | //| N - number of columns in matrix A. | //| OUTPUT PARAMETERS: | //| A - matrices L and U in compact form: | //| * L is stored under main diagonal | //| * U is stored on and above main diagonal | //| Pivots - permutation matrix in compact form. | //| array[0..Min(M-1,N-1)]. | //+------------------------------------------------------------------+ static void CTrFac::RMatrixLU(CMatrixDouble &a,const int m,const int n,int &pivots[]) { //--- check if(!CAp::Assert(m>0,__FUNCTION__+": incorrect M!")) return; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!")) return; //--- function call RMatrixPLU(a,m,n,pivots); } //+------------------------------------------------------------------+ //| LU decomposition of a general complex matrix with row pivoting | //| A is represented as A = P*L*U, where: | //| * L is lower unitriangular matrix | //| * U is upper triangular matrix | //| * P = P0*P1*...*PK, K=min(M,N)-1, | //| Pi - permutation matrix for I and Pivots[I] | //| This is cache-oblivous implementation of LU decomposition. It is | //| optimized for square matrices. As for rectangular matrices: | //| * best case - M>>N | //| * worst case - N>>M, small M, large N, matrix does not fit in CPU| //| cache | //| INPUT PARAMETERS: | //| A - array[0..M-1, 0..N-1]. | //| M - number of rows in matrix A. | //| N - number of columns in matrix A. | //| OUTPUT PARAMETERS: | //| A - matrices L and U in compact form: | //| * L is stored under main diagonal | //| * U is stored on and above main diagonal | //| Pivots - permutation matrix in compact form. | //| array[0..Min(M-1,N-1)]. | //+------------------------------------------------------------------+ static void CTrFac::CMatrixLU(CMatrixComplex &a,const int m,const int n,int &pivots[]) { //--- check if(!CAp::Assert(m>0,__FUNCTION__+": incorrect M!")) return; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!")) return; //--- function call CMatrixPLU(a,m,n,pivots); } //+------------------------------------------------------------------+ //| Cache-oblivious Cholesky decomposition | //| The algorithm computes Cholesky decomposition of a Hermitian | //| positive - definite matrix. The result of an algorithm is a | //| representation of A as A=U'*U or A=L*L' (here X' detones | //| conj(X^T)). | //| INPUT PARAMETERS: | //| A - upper or lower triangle of a factorized matrix. | //| array with elements [0..N-1, 0..N-1]. | //| N - size of matrix A. | //| IsUpper - if IsUpper=True, then A contains an upper | //| triangle of a symmetric matrix, otherwise A | //| contains a lower one. | //| OUTPUT PARAMETERS: | //| A - the result of factorization. If IsUpper=True, | //| then the upper triangle contains matrix U, so | //| that A = U'*U, and the elements below the main | //| diagonal are not modified. Similarly, if | //| IsUpper = False. | //| RESULT: | //| If the matrix is positive-definite, the function returns | //| True. Otherwise, the function returns False. Contents of A is| //| not determined in such case. | //+------------------------------------------------------------------+ static bool CTrFac::HPDMatrixCholesky(CMatrixComplex &a,const int n,const bool isupper) { //--- create array complex tmp[]; //--- check if(n<1) return(false); //--- return result return(HPDMatrixCholeskyRec(a,0,n,isupper,tmp)); } //+------------------------------------------------------------------+ //| Cache-oblivious Cholesky decomposition | //| The algorithm computes Cholesky decomposition of a symmetric | //| positive - definite matrix. The result of an algorithm is a | //| representation of A as A=U^T*U or A=L*L^T | //| INPUT PARAMETERS: | //| A - upper or lower triangle of a factorized matrix. | //| array with elements [0..N-1, 0..N-1]. | //| N - size of matrix A. | //| IsUpper - if IsUpper=True, then A contains an upper | //| triangle of a symmetric matrix, otherwise A | //| contains a lower one. | //| OUTPUT PARAMETERS: | //| A - the result of factorization. If IsUpper=True, | //| then the upper triangle contains matrix U, so | //| that A = U^T*U, and the elements below the main | //| diagonal are not modified. Similarly, if | //| IsUpper = False. | //| RESULT: | //| If the matrix is positive-definite, the function returns | //| True. Otherwise, the function returns False. Contents of A is| //| not determined in such case. | //+------------------------------------------------------------------+ static bool CTrFac::SPDMatrixCholesky(CMatrixDouble &a,const int n,const bool isupper) { //--- create array double tmp[]; //--- check if(n<1) return(false); //--- return result return(SPDMatrixCholeskyRec(a,0,n,isupper,tmp)); } //+------------------------------------------------------------------+ //| LUP decomposition of general real matrix | //+------------------------------------------------------------------+ static void CTrFac::RMatrixLUP(CMatrixDouble &a,const int m,const int n,int &pivots[]) { //--- create variables int i=0; int j=0; double mx=0; double v=0; int i_=0; //--- create array double tmp[]; //--- Internal LU decomposition subroutine. //--- Never call it directly. if(!CAp::Assert(m>0,__FUNCTION__+": incorrect M!")) return; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!")) return; //--- Scale matrix to avoid overflows, //--- decompose it, then scale back. mx=0; for(i=0;i0,__FUNCTION__+": incorrect M!")) return; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!")) return; //--- Scale matrix to avoid overflows, //--- decompose it, then scale back. mx=0; for(i=0;i0,__FUNCTION__+": incorrect M!")) return; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!")) return; //--- allocation ArrayResizeAL(tmp,2*MathMax(m,n)); ArrayResizeAL(pivots,MathMin(m,n)); //--- Scale matrix to avoid overflows, //--- decompose it, then scale back. mx=0; for(i=0;i0,__FUNCTION__+": incorrect M!")) return; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": incorrect N!")) return; //--- allocation ArrayResizeAL(tmp,2*MathMax(m,n)); ArrayResizeAL(pivots,MathMin(m,n)); //--- Scale matrix to avoid overflows, //--- decompose it, then scale back. mx=0; for(i=0;i0.0) { a[offs].Set(offs,MathSqrt(a[offs][offs])); //--- return result return(true); } else { //--- return result return(false); } } //--- check if(n<=CAblas::AblasBlockSize()) { //--- return result return(SPDMatrixCholesky2(a,offs,n,isupper,tmp)); } //--- general case: split task in cache-oblivious manner result=true; CAblas::AblasSplitLength(a,n,n1,n2); result=SPDMatrixCholeskyRec(a,offs,n1,isupper,tmp); //--- check if(!result) { //--- return result return(result); } //--- check if(n2>0) { if(isupper) { //--- function call CAblas::RMatrixLeftTrsM(n1,n2,a,offs,offs,isupper,false,1,a,offs,offs+n1); //--- function call CAblas::RMatrixSyrk(n2,n1,-1.0,a,offs,offs+n1,1,1.0,a,offs+n1,offs+n1,isupper); } else { //--- function call CAblas::RMatrixRightTrsM(n2,n1,a,offs,offs,isupper,false,1,a,offs+n1,offs); //--- function call CAblas::RMatrixSyrk(n2,n1,-1.0,a,offs+n1,offs,0,1.0,a,offs+n1,offs+n1,isupper); } result=SPDMatrixCholeskyRec(a,offs+n1,n2,isupper,tmp); //--- check if(!result) { //--- return result return(result); } } //--- return result return(result); } //+------------------------------------------------------------------+ //| Recurrent complex LU subroutine. | //| Never call it directly. | //+------------------------------------------------------------------+ static void CTrFac::CMatrixLUPRec(CMatrixComplex &a,const int offs,const int m, const int n,int &pivots[],complex &tmp[]) { //--- create variables int i=0; int m1=0; int m2=0; int i_=0; int i1_=0; complex One(1.0,0.0); complex _One(-1.0,0.0); //--- Kernel case if(MathMin(m,n)<=CAblas::AblasComplexBlockSize()) { CMatrixLUP2(a,offs,m,n,pivots,tmp); //--- exit the function return; } //--- Preliminary step, make N>=M //--- ( A1 ) //--- A = ( ), where A1 is square //--- ( A2 ) //--- Factorize A1, update A2 if(m>n) { //--- function call CMatrixLUPRec(a,offs,n,n,pivots,tmp); for(i=0;i0) { for(i=0;i=M //--- ( A1 ) //--- A = ( ), where A1 is square //--- ( A2 ) //--- Factorize A1, update A2 if(m>n) { //--- function call RMatrixLUPRec(a,offs,n,n,pivots,tmp); for(i=0;i0) { for(i=0;i=N. //--- A = (A1 A2), where A1 is square //--- Factorize A1, update A2 if(n>m) { //--- function call CMatrixPLURec(a,offs,m,m,pivots,tmp); for(i=0;i0) { for(i=0;i=N. //--- A = (A1 A2), where A1 is square //--- Factorize A1, update A2 if(n>m) { //--- function call RMatrixPLURec(a,offs,m,m,pivots,tmp); for(i=0;i0) { for(i=0;iCMath::AbsComplex(a[offs+j][offs+jp])) jp=i; } pivots[offs+j]=offs+jp; //--- check if(jp!=j) { i1_=offs; for(i_=0;i_MathAbs(a[offs+j][offs+jp])) jp=i; } pivots[offs+j]=offs+jp; //--- check if(jp!=j) { i1_=offs; for(i_=0;i_CMath::AbsComplex(a[offs+jp][offs+j])) jp=i; } pivots[offs+j]=offs+jp; if(a[offs+jp][offs+j]!=zero) { //--- Apply the interchange to rows if(jp!=j) { for(i=0;iMathAbs(a[offs+jp][offs+j])) jp=i; } pivots[offs+j]=offs+jp; //--- check if(a[offs+jp][offs+j]!=0.0) { //--- Apply the interchange to rows if(jp!=j) { for(i=0;i0.0) { a[offs].Set(offs,MathSqrt(a[offs][offs].re)); result=true; } else { result=false; } //--- return result return(result); } //--- check if(n<=CAblas::AblasComplexBlockSize()) { result=HPDMatrixCholesky2(a,offs,n,isupper,tmp); //--- return result return(result); } //--- general case: split task in cache-oblivious manner result=true; //--- function call CAblas::AblasComplexSplitLength(a,n,n1,n2); result=HPDMatrixCholeskyRec(a,offs,n1,isupper,tmp); //--- check if(!result) return(result); //--- check if(n2>0) { //--- check if(isupper) { //--- function call CAblas::CMatrixLeftTrsM(n1,n2,a,offs,offs,isupper,false,2,a,offs,offs+n1); //--- function call CAblas::CMatrixSyrk(n2,n1,-1.0,a,offs,offs+n1,2,1.0,a,offs+n1,offs+n1,isupper); } else { //--- function call CAblas::CMatrixRightTrsM(n2,n1,a,offs,offs,isupper,false,2,a,offs+n1,offs); //--- function call CAblas::CMatrixSyrk(n2,n1,-1.0,a,offs+n1,offs,0,1.0,a,offs+n1,offs+n1,isupper); } result=HPDMatrixCholeskyRec(a,offs+n1,n2,isupper,tmp); //--- check if(!result) return(result); } //--- return result return(result); } //+------------------------------------------------------------------+ //| Level-2 Hermitian Cholesky subroutine. | //| -- LAPACK routine (version 3.0) -- | //| Univ. of Tennessee,Univ. of California Berkeley,NAG Ltd., | //| Courant Institute,Argonne National Lab,and Rice University | //| February 29,1992 | //+------------------------------------------------------------------+ static bool CTrFac::HPDMatrixCholesky2(CMatrixComplex &aaa,const int offs, const int n,const bool isupper,complex &tmp[]) { //--- create variables bool result; int i=0; int j=0; double ajj=0; complex v=0; complex cR; double r=0; int i_=0; int i1_=0; //--- initialization result=true; //--- check if(n<0) return(false); //--- check if(n==0) return(result); //--- check if(isupper) { //--- Compute the Cholesky factorization A = U'*U. for(j=0;j0) { i1_=offs; for(i_=0;i_0) { i1_=offs; for(i_=0;i_0) { i1_=offs; for(i_=0;i_0) { i1_=offs; for(i_=0;i_=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- allocation ArrayResizeAL(t,n); //--- fiiling array for(i=0;i=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- change values nrm=0; for(i=0;i=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- allocation ArrayResizeAL(t,n); //--- fiiling array for(i=0;i=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- change values nrm=0; for(i=0;i=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- allocation ArrayResizeAL(t,n); //--- fiiling array for(i=0;i=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- change values nrm=0; for(i=0;i=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- function call CMatrixRCondLUInternal(lua,n,true,false,0.0,v); //--- return result return(v); } //+------------------------------------------------------------------+ //| Estimate of the condition number of a matrix given by its LU | //| decomposition (infinity norm). | //| The algorithm calculates a lower bound of the condition number. | //| In this case, the algorithm does not return a lower bound of the | //| condition number, but an inverse number (to avoid an overflow in | //| case of a singular matrix). | //| Input parameters: | //| LUA - LU decomposition of a matrix in compact form. | //| Output of the CMatrixLU subroutine. | //| N - size of matrix A. | //| Result: 1/LowerBound(cond(A)) | //| NOTE: | //| if k(A) is very large, then matrix is assumed degenerate, | //| k(A)=INF, 0.0 is returned in such cases. | //+------------------------------------------------------------------+ static double CRCond::CMatrixLURCondInf(CMatrixComplex &lua,const int n) { //--- create a variable double v=0; //--- check if(!CAp::Assert(n>=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- function call CMatrixRCondLUInternal(lua,n,false,false,0.0,v); //--- return result return(v); } //+------------------------------------------------------------------+ //| Triangular matrix: estimate of a condition number (1-norm) | //| The algorithm calculates a lower bound of the condition number. | //| In this case, the algorithm does not return a lower bound of the | //| condition number, but an inverse number (to avoid an overflow in | //| case of a singular matrix). | //| Input parameters: | //| A - matrix. Array[0..N-1, 0..N-1]. | //| N - size of A. | //| IsUpper - True, if the matrix is upper triangular. | //| IsUnit - True, if the matrix has a unit diagonal. | //| Result: 1/LowerBound(cond(A)) | //| NOTE: | //| if k(A) is very large, then matrix is assumed degenerate, | //| k(A)=INF, 0.0 is returned in such cases. | //+------------------------------------------------------------------+ static double CRCond::CMatrixTrRCond1(CMatrixComplex &a,const int n, const bool isupper,const bool isunit) { //--- create variables int i=0; int j=0; double v=0; double nrm=0; int j1=0; int j2=0; //--- create arrays int pivots[]; double t[]; //--- check if(!CAp::Assert(n>=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- allocation ArrayResizeAL(t,n); //--- fiiling array for(i=0;i=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- change values nrm=0; for(i=0;i0 after this block if(anorm==0.0) return; //--- check if(n==1) { rc=1; return; } //--- Estimate the norm of inv(A). ainvnm=0; kase=0; while(true) { //--- function call RMatrixEstimateNorm(n,ev,ex,iwork,ainvnm,kase); //--- check if(kase==0) break; //--- from 1-based array to 0-based for(i=0;i=0;i--) ex[i+1]=ex[i]; } //--- Compute the estimate of the reciprocal condition number. if(ainvnm!=0.0) { rc=1/ainvnm; rc=rc/anorm; //--- check if(rc=0;i--) ex[i+1]=ex[i]; } //--- Compute the estimate of the reciprocal condition number. if(ainvnm!=0.0) { rc=1/ainvnm; rc=rc/anorm; //--- check if(rc=1)) return; //--- allocation ArrayResizeAL(tmp,n); //--- initialization rc=0; //--- prepare parameters for triangular solver maxgrowth=1/RCondThreshold(); sa=0; //--- check if(isupper) { for(i=0;i=1;i--) { i1_=1; v=0.0; for(i_=0;i_=0;i--) ex[i+1]=ex[i]; } //--- Compute the estimate of the reciprocal condition number. if(ainvnm!=0.0) { v=1/ainvnm; rc=v/anorm; //--- check if(rc=1)) return; //--- allocation ArrayResizeAL(tmp,n); //--- initialization rc=0; //--- prepare parameters for triangular solver maxgrowth=1/RCondThreshold(); sa=0; //--- check if(isupper) { for(i=0;i=1;i--) { i1_=1; v=0.0; for(i_=0;i_0 if(anorm==0.0) return; //--- check if(n==1) { rc=1; return; } //--- Estimate the norm of inv(A). ainvnm=0; kase=0; while(true) { //--- function call CMatrixEstimateNorm(n,ev,ex,ainvnm,kase,isave,rsave); //--- check if(kase==0) break; for(i=0;i=0;i--) ex[i+1]=ex[i]; } //--- Compute the estimate of the reciprocal condition number. if(ainvnm!=0.0) { rc=1/ainvnm; rc=rc/anorm; //--- check if(rc=1;i--) { //--- check if(i>1) { i1_=1; v=0.0; for(i_=0;i_<=i-2;i_++) v+=lua[i-1][i_]*ex[i_+i1_]; } else v=0; ex[i]=ex[i]+v; } } else { //--- Multiply by L' for(i=0;i=1) { for(i_=0;i_0 after this block if(anorm==0.0) return; //--- check if(n==1) { rc=1; return; } //--- Estimate the norm of inv(A). ainvnm=0; kase=0; while(true) { //--- function call RMatrixEstimateNorm(n,ev,ex,iwork,ainvnm,kase); //--- check if(kase==0) break; //--- from 1-based array to 0-based for(i=0;i=0;i--) ex[i+1]=ex[i]; } //--- Compute the estimate of the reciprocal condition number. if(ainvnm!=0.0) { rc=1/ainvnm; rc=rc/anorm; //--- check if(rc=1;i--) { v=0; //--- check if(i>1) { i1_=1; v=0.0; for(i_=0;i_<=i-2;i_++) v+=lua[i-1][i_]*ex[i_+i1_]; } ex[i]=v+ex[i]; } } else { //--- Multiply by L' for(i=1;i<=n;i++) cwork2[i]=0; for(i=1;i<=n;i++) { v=ex[i]; //--- check if(i>1) { i1_=-1; for(i_=1;i_=0;i--) ex[i+1]=ex[i]; } //--- Compute the estimate of the reciprocal condition number. if(ainvnm!=0.0) { rc=1/ainvnm; rc=rc/anorm; //--- check if(rc=0.0) x[i]=1; else x[i]=-1; //--- check if(x[i]>0) isgn[i]=1; //--- check if(x[i]<0) isgn[i]=-1; //--- check if(x[i]==0) isgn[i]=0; } kase=2; isgn[posjump]=2; //--- exit the function return; } //--- ................ entry (jump = 2) //--- first iteration. x has been overwritten by trandpose(a)*x. if(isgn[posjump]==2) { isgn[posj]=1; for(i=2;i<=n;i++) { //--- check if(MathAbs(x[i])>MathAbs(x[isgn[posj]])) isgn[posj]=i; } isgn[positer]=2; //--- main loop - iterations 2,3,...,itmax. for(i=1;i<=n;i++) x[i]=0; x[isgn[posj]]=1; kase=1; isgn[posjump]=3; //--- exit the function return; } //--- ................ entry (jump = 3) //--- x has been overwritten by a*x. if(isgn[posjump]==3) { for(i_=1;i_<=n;i_++) v[i_]=x[i_]; v[posestold]=est; //--- change value est=0; for(i=1;i<=n;i++) est=est+MathAbs(v[i]); flg=false; for(i=1;i<=n;i++) { if(((x[i]>=0.0)&&(isgn[i]<0))||((x[i]<0.0)&&(isgn[i]>=0))) flg=true; } //--- repeated sign vector detected, hence algorithm has converged. //--- or may be cycling. if(!flg || est<=v[posestold]) { v[posaltsgn]=1; for(i=1;i<=n;i++) { x[i]=v[posaltsgn]*(1+(double)(i-1)/(double)(n-1)); v[posaltsgn]=-v[posaltsgn]; } kase=1; isgn[posjump]=5; //--- exit the function return; } for(i=1;i<=n;i++) { //--- check if(x[i]>=0.0) { x[i]=1; isgn[i]=1; } else { x[i]=-1; isgn[i]=-1; } } kase=2; isgn[posjump]=4; //--- exit the function return; } //--- ................ entry (jump = 4) //--- x has been overwritten by trandpose(a)*x. if(isgn[posjump]==4) { isgn[posjlast]=isgn[posj]; isgn[posj]=1; for(i=2;i<=n;i++) { //--- check if(MathAbs(x[i])>MathAbs(x[isgn[posj]])) isgn[posj]=i; } //--- check if(x[isgn[posjlast]]!=MathAbs(x[isgn[posj]]) && isgn[positer]est) { for(i_=1;i_<=n;i_++) v[i_]=x[i_]; est=v[postemp]; } kase=0; //--- exit the function return; } } //+------------------------------------------------------------------+ //| Internal subroutine | //+------------------------------------------------------------------+ static void CRCond::CMatrixEstimateNorm(const int n,complex &v[],complex &x[], double &est,int &kase,int &isave[], double &rsave[]) { //--- create variables int itmax=0; int i=0; int iter=0; int j=0; int jlast=0; int jump=0; double absxi=0; double altsgn=0; double estold=0; double safmin=0; double temp=0; int i_=0; //--- Executable Statements itmax=5; safmin=CMath::m_minrealnumber; //--- check if(kase==0) { //--- allocation ArrayResizeAL(v,n+1); ArrayResizeAL(x,n+1); ArrayResizeAL(isave,5); ArrayResizeAL(rsave,4); for(i=1;i<=n;i++) x[i]=1.0/(double)n; kase=1; jump=1; //--- function call InternalComplexRCondSaveAll(isave,rsave,i,iter,j,jlast,jump,absxi,altsgn,estold,temp); //--- exit the function return; } //--- function call InternalComplexRCondLoadAll(isave,rsave,i,iter,j,jlast,jump,absxi,altsgn,estold,temp); //--- entry (jump = 1) //--- first iteration. x has been overwritten by a*x. if(jump==1) { //--- check if(n==1) { v[1]=x[1]; //--- function call est=CMath::AbsComplex(v[1]); kase=0; //--- function call InternalComplexRCondSaveAll(isave,rsave,i,iter,j,jlast,jump,absxi,altsgn,estold,temp); //--- exit the function return; } //--- function call est=InternalComplexRCondScSum1(x,n); for(i=1;i<=n;i++) { //--- function call absxi=CMath::AbsComplex(x[i]); //--- check if(absxi>safmin) x[i]=x[i]/absxi; else x[i]=1; } kase=2; jump=2; //--- function call InternalComplexRCondSaveAll(isave,rsave,i,iter,j,jlast,jump,absxi,altsgn,estold,temp); //--- exit the function return; } //--- entry (jump = 2) //--- first iteration. x has been overwritten by ctrans(a)*x. if(jump==2) { j=InternalComplexRCondIcMax1(x,n); iter=2; //--- main loop - iterations 2,3,...,itmax. for(i=1;i<=n;i++) x[i]=0; x[j]=1; kase=1; jump=3; //--- function call InternalComplexRCondSaveAll(isave,rsave,i,iter,j,jlast,jump,absxi,altsgn,estold,temp); //--- exit the function return; } //--- entry (jump = 3) //--- x has been overwritten by a*x. if(jump==3) { for(i_=1;i_<=n;i_++) v[i_]=x[i_]; estold=est; //--- function call est=InternalComplexRCondScSum1(v,n); //--- test for cycling. if(est<=estold) { //--- iteration complete. final stage. altsgn=1; for(i=1;i<=n;i++) { x[i]=altsgn*(1+(double)(i-1)/(double)(n-1)); altsgn=-altsgn; } kase=1; jump=5; //--- function call InternalComplexRCondSaveAll(isave,rsave,i,iter,j,jlast,jump,absxi,altsgn,estold,temp); //--- exit the function return; } for(i=1;i<=n;i++) { absxi=CMath::AbsComplex(x[i]); //--- check if(absxi>safmin) x[i]=x[i]/absxi; else x[i]=1; } kase=2; jump=4; //--- function call InternalComplexRCondSaveAll(isave,rsave,i,iter,j,jlast,jump,absxi,altsgn,estold,temp); //--- exit the function return; } //--- entry (jump = 4) //--- x has been overwritten by ctrans(a)*x. if(jump==4) { jlast=j; j=InternalComplexRCondIcMax1(x,n); //--- check if(CMath::AbsComplex(x[jlast])!=CMath::AbsComplex(x[j]) && iterest) { for(i_=1;i_<=n;i_++) v[i_]=x[i_]; est=temp; } kase=0; //--- function call InternalComplexRCondSaveAll(isave,rsave,i,iter,j,jlast,jump,absxi,altsgn,estold,temp); //--- exit the function return; } } //+------------------------------------------------------------------+ //| Internal subroutine | //+------------------------------------------------------------------+ static double CRCond::InternalComplexRCondScSum1(complex &x[],const int n) { //--- create variables double result=0; int i=0; //--- get result for(i=1;i<=n;i++) result=result+CMath::AbsComplex(x[i]); //--- return result return(result); } //+------------------------------------------------------------------+ //| Internal subroutine | //+------------------------------------------------------------------+ static int CRCond::InternalComplexRCondIcMax1(complex &x[],const int n) { //--- create variables int result=0; int i=0; double m=0; //--- get result result=1; m=CMath::AbsComplex(x[1]); for(i=2;i<=n;i++) { //--- check if(CMath::AbsComplex(x[i])>m) { result=i; m=CMath::AbsComplex(x[i]); } } //--- return result return(result); } //+------------------------------------------------------------------+ //| Internal subroutine | //+------------------------------------------------------------------+ static void CRCond::InternalComplexRCondSaveAll(int &isave[],double &rsave[], int &i,int &iter,int &j, int &jlast,int &jump, double &absxi,double &altsgn, double &estold,double &temp) { //--- copy isave[0]=i; isave[1]=iter; isave[2]=j; isave[3]=jlast; isave[4]=jump; //--- copy rsave[0]=absxi; rsave[1]=altsgn; rsave[2]=estold; rsave[3]=temp; } //+------------------------------------------------------------------+ //| Internal subroutine | //+------------------------------------------------------------------+ static void CRCond::InternalComplexRCondLoadAll(int &isave[],double &rsave[], int &i,int &iter,int &j, int &jlast,int &jump, double &absxi,double &altsgn, double &estold,double &temp) { //--- get i=isave[0]; iter=isave[1]; j=isave[2]; jlast=isave[3]; jump=isave[4]; //--- get absxi=rsave[0]; altsgn=rsave[1]; estold=rsave[2]; temp=rsave[3]; } //+------------------------------------------------------------------+ //| Matrix inverse report: | //| * R1 reciprocal of condition number in 1-norm | //| * RInf reciprocal of condition number in inf-norm | //+------------------------------------------------------------------+ class CMatInvReport { public: double m_r1; double m_rinf; //--- constructor, destructor CMatInvReport(void); ~CMatInvReport(void); }; //+------------------------------------------------------------------+ //| Constructor without parameters | //+------------------------------------------------------------------+ CMatInvReport::CMatInvReport(void) { } //+------------------------------------------------------------------+ //| Destructor | //+------------------------------------------------------------------+ CMatInvReport::~CMatInvReport(void) { } //+------------------------------------------------------------------+ //| Matrix inverse report: | //| * R1 reciprocal of condition number in 1-norm | //| * RInf reciprocal of condition number in inf-norm | //+------------------------------------------------------------------+ class CMatInvReportShell { private: CMatInvReport m_innerobj; public: //--- constructors, destructor CMatInvReportShell(void); CMatInvReportShell(CMatInvReport &obj); ~CMatInvReportShell(void); //--- methods double GetR1(void); void SetR1(double r); double GetRInf(void); void SetRInf(double r); CMatInvReport *GetInnerObj(void); }; //+------------------------------------------------------------------+ //| Constructor without parameters | //+------------------------------------------------------------------+ CMatInvReportShell::CMatInvReportShell(void) { } //+------------------------------------------------------------------+ //| Copy | //+------------------------------------------------------------------+ CMatInvReportShell::CMatInvReportShell(CMatInvReport &obj) { m_innerobj.m_r1=obj.m_r1; m_innerobj.m_rinf=obj.m_rinf; } //+------------------------------------------------------------------+ //| Destructor | //+------------------------------------------------------------------+ CMatInvReportShell::~CMatInvReportShell(void) { } //+------------------------------------------------------------------+ //| Returns the value of the variable r1 | //+------------------------------------------------------------------+ double CMatInvReportShell::GetR1(void) { return(m_innerobj.m_r1); } //+------------------------------------------------------------------+ //| Changing the value of the variable r1 | //+------------------------------------------------------------------+ void CMatInvReportShell::SetR1(double r) { m_innerobj.m_r1=r; } //+------------------------------------------------------------------+ //| Returns the value of the variable rinf | //+------------------------------------------------------------------+ double CMatInvReportShell::GetRInf(void) { return(m_innerobj.m_rinf); } //+------------------------------------------------------------------+ //| Changing the value of the variable rint | //+------------------------------------------------------------------+ void CMatInvReportShell::SetRInf(double r) { m_innerobj.m_rinf=r; } //+------------------------------------------------------------------+ //| Return object of class | //+------------------------------------------------------------------+ CMatInvReport *CMatInvReportShell::GetInnerObj(void) { return(GetPointer(m_innerobj)); } //+------------------------------------------------------------------+ //| Inverse matrix | //+------------------------------------------------------------------+ class CMatInv { private: //--- private methods static void RMatrixTrInverseRec(CMatrixDouble &a,const int offs,const int n,const bool isupper,const bool isunit,double &tmp[],int &info,CMatInvReport &rep); static void CMatrixTrInverseRec(CMatrixComplex &a,const int offs,const int n,const bool isupper,const bool isunit,complex &tmp[],int &info,CMatInvReport &rep); static void RMatrixLUInverseRec(CMatrixDouble &a,const int offs,const int n,double &work[],int &info,CMatInvReport &rep); static void CMatrixLUInverseRec(CMatrixComplex &a,const int offs,const int n,complex &work[],int &info,CMatInvReport &rep); static void SPDMatrixCholeskyInverseRec(CMatrixDouble &a,const int offs,const int n,const bool isupper,double &tmp[]); static void HPDMatrixCholeskyInverseRec(CMatrixComplex &a,const int offs,const int n,const bool isupper,complex &tmp[]); public: CMatInv(void); ~CMatInv(void); //--- public methods static void RMatrixLUInverse(CMatrixDouble &a,int &pivots[],const int n,int &info,CMatInvReport &rep); static void RMatrixInverse(CMatrixDouble &a,const int n,int &info,CMatInvReport &rep); static void SPDMatrixCholeskyInverse(CMatrixDouble &a,const int n,const bool isupper,int &info,CMatInvReport &rep); static void SPDMatrixInverse(CMatrixDouble &a,const int n,const bool isupper,int &info,CMatInvReport &rep); static void RMatrixTrInverse(CMatrixDouble &a,const int n,const bool isupper,const bool isunit,int &info,CMatInvReport &rep); static void CMatrixLUInverse(CMatrixComplex &a,int &pivots[],const int n,int &info,CMatInvReport &rep); static void CMatrixInverse(CMatrixComplex &a,const int n,int &info,CMatInvReport &rep); static void HPDMatrixCholeskyInverse(CMatrixComplex &a,const int n,const bool isupper,int &info,CMatInvReport &rep); static void HPDMatrixInverse(CMatrixComplex &a,const int n,const bool isupper,int &info,CMatInvReport &rep); static void CMatrixTrInverse(CMatrixComplex &a,const int n,const bool isupper,const bool isunit,int &info,CMatInvReport &rep); }; //+------------------------------------------------------------------+ //| Constructor without parameters | //+------------------------------------------------------------------+ CMatInv::CMatInv(void) { } //+------------------------------------------------------------------+ //| Destructor | //+------------------------------------------------------------------+ CMatInv::~CMatInv(void) { } //+------------------------------------------------------------------+ //| Inversion of a matrix given by its LU decomposition. | //| INPUT PARAMETERS: | //| A - LU decomposition of the matrix | //| (output of RMatrixLU subroutine). | //| Pivots - table of permutations | //| (the output of RMatrixLU subroutine). | //| N - size of matrix A (optional) : | //| * if given, only principal NxN submatrix is | //| processed and overwritten. other elements are | //| unchanged. | //| * if not given, size is automatically determined | //| from matrix size (A must be square matrix) | //| OUTPUT PARAMETERS: | //| Info - return code: | //| * -3 A is singular, or VERY close to singular.| //| it is filled by zeros in such cases. | //| * 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 | //| A - inverse of matrix A. | //| Array whose indexes range within | //| [0..N-1, 0..N-1]. | //| 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. | //+------------------------------------------------------------------+ static void CMatInv::RMatrixLUInverse(CMatrixDouble &a,int &pivots[], const int n,int &info,CMatInvReport &rep) { //--- create variables int i=0; int j=0; int k=0; double v=0; //--- create array double work[]; //--- initialization info=0; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": N<=0!")) return; //--- check if(!CAp::Assert(CAp::Cols(a)>=n,__FUNCTION__+": cols(A)=n,__FUNCTION__+": rows(A)=n,__FUNCTION__+": len(Pivots)n-1 || pivots[i]0,__FUNCTION__+": incorrect Pivots array!")) return; //--- calculate condition numbers rep.m_r1=CRCond::RMatrixLURCond1(a,n); rep.m_rinf=CRCond::RMatrixLURCondInf(a,n); //--- check if(rep.m_r1=0;j--) { k=pivots[j]; v=a[i][j]; a[i].Set(j,a[i][k]); a[i].Set(k,v); } } } //+------------------------------------------------------------------+ //| Inversion of a general matrix. | //| Input parameters: | //| A - matrix. | //| N - size of matrix A (optional) : | //| * if given, only principal NxN submatrix is | //| processed and overwritten. other elements are | //| unchanged. | //| * if not given, size is automatically determined | //| from matrix size (A must be square matrix) | //| Output parameters: | //| Info - return code, same as in RMatrixLUInverse | //| Rep - solver report, same as in RMatrixLUInverse | //| A - inverse of matrix A, same as in RMatrixLUInverse | //| Result: | //| True, if the matrix is not singular. | //| False, if the matrix is singular. | //+------------------------------------------------------------------+ static void CMatInv::RMatrixInverse(CMatrixDouble &a,const int n,int &info, CMatInvReport &rep) { //--- create array int pivots[]; //--- initialization info=0; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": N<=0!")) return; //--- check if(!CAp::Assert(CAp::Cols(a)>=n,__FUNCTION__+": cols(A)=n,__FUNCTION__+": rows(A)0,__FUNCTION__+": N<=0!")) return; //--- check if(!CAp::Assert(CAp::Cols(a)>=n,__FUNCTION__+": cols(A)=n,__FUNCTION__+": rows(A)=n,__FUNCTION__+": len(Pivots)n-1 || pivots[i]0,__FUNCTION__+": incorrect Pivots array!")) return; //--- calculate condition numbers rep.m_r1=CRCond::CMatrixLURCond1(a,n); rep.m_rinf=CRCond::CMatrixLURCondInf(a,n); //--- check if(rep.m_r1=0;j--) { k=pivots[j]; v=a[i][j]; a[i].Set(j,a[i][k]); a[i].Set(k,v); } } } //+------------------------------------------------------------------+ //| Inversion of a general matrix. | //| Input parameters: | //| A - matrix | //| N - size of matrix A (optional) : | //| * if given, only principal NxN submatrix is | //| processed and overwritten. other elements are | //| unchanged. | //| * if not given, size is automatically determined | //| from matrix size (A must be square matrix) | //| Output parameters: | //| Info - return code, same as in RMatrixLUInverse | //| Rep - solver report, same as in RMatrixLUInverse | //| A - inverse of matrix A, same as in RMatrixLUInverse | //+------------------------------------------------------------------+ static void CMatInv::CMatrixInverse(CMatrixComplex &a,const int n,int &info, CMatInvReport &rep) { //--- create array int pivots[]; //--- initialization info=0; //--- check if(!CAp::Assert(n>0,__FUNCTION__+": N<=0!")) return; //--- check if(!CAp::Assert(CAp::Cols(a)>=n,__FUNCTION__+": cols(A)=n,__FUNCTION__+": rows(A)0,__FUNCTION__+": N<=0!")) return; //--- check if(!CAp::Assert(CAp::Cols(a)>=n,__FUNCTION__+": cols(A)=n,__FUNCTION__+": rows(A)0,__FUNCTION__+": N<=0!")) return; //--- check if(!CAp::Assert(CAp::Cols(a)>=n,__FUNCTION__+": cols(A)=n,__FUNCTION__+": rows(A)0,__FUNCTION__+": N<=0!")) return; //--- check if(!CAp::Assert(CAp::Cols(a)>=n,__FUNCTION__+": cols(A)=n,__FUNCTION__+": rows(A)0,__FUNCTION__+": N<=0!")) return; //--- check if(!CAp::Assert(CAp::Cols(a)>=n,__FUNCTION__+": cols(A)=n,__FUNCTION__+": rows(A)0,__FUNCTION__+": N<=0!")) return; //--- check if(!CAp::Assert(CAp::Cols(a)>=n,__FUNCTION__+": cols(A)=n,__FUNCTION__+": rows(A)0,__FUNCTION__+": N<=0!")) return; //--- check if(!CAp::Assert(CAp::Cols(a)>=n,__FUNCTION__+": cols(A)=n,__FUNCTION__+": rows(A)0) { i1_=offs; for(i_=0;i_=0;j--) { //--- check if(!isunit) { //--- check if(a[offs+j][offs+j]==0.0) { info=-3; //--- exit the function return; } a[offs+j].Set(offs+j,1/a[offs+j][offs+j]); ajj=-a[offs+j][offs+j]; } else ajj=-1; //--- check if(jj+1) { i1_=-offs; v=0.0; for(i_=offs+j+1;i_0) { //--- check if(isupper) { for(i=0;i0) { i1_=offs; for(i_=0;i_=0;j--) { //--- check if(!isunit) { //--- check if(a[offs+j][offs+j]==0) { info=-3; //--- exit the function return; } a[offs+j].Set(offs+j,One/a[offs+j][offs+j]); ajj=-a[offs+j][offs+j]; } else ajj=-1; //--- check if(jj+1) { i1_=-offs; v=0.0; for(i_=offs+j+1;i_0) { //--- check if(isupper) { for(i=0;i=0;j--) { //--- Copy current column of L to WORK and replace with zeros. for(i=j+1;i0,__FUNCTION__+": internal error!")) return; //--- X :=inv(U1)*U12*inv(U2) CAblas::RMatrixLeftTrsM(n1,n2,a,offs,offs,true,false,0,a,offs,offs+n1); CAblas::RMatrixRightTrsM(n1,n2,a,offs+n1,offs+n1,true,false,0,a,offs,offs+n1); //--- Y :=inv(L2)*L12*inv(L1) CAblas::RMatrixLeftTrsM(n2,n1,a,offs+n1,offs+n1,false,true,0,a,offs+n1,offs); CAblas::RMatrixRightTrsM(n2,n1,a,offs,offs,false,true,0,a,offs+n1,offs); //--- W :=inv(L1*U1)+X*Y RMatrixLUInverseRec(a,offs,n1,work,info,rep); //--- check if(info<=0) return; //--- function call CAblas::RMatrixGemm(n1,n1,n2,1.0,a,offs,offs+n1,0,a,offs+n1,offs,0,1.0,a,offs,offs); //--- X :=-X*inv(L2) //--- Y :=-inv(U2)*Y CAblas::RMatrixRightTrsM(n1,n2,a,offs+n1,offs+n1,false,true,0,a,offs,offs+n1); for(i=0;i=0;j--) { //--- Copy current column of L to WORK and replace with zeros. for(i=j+1;i0,__FUNCTION__+": internal error!")) return; //--- X :=inv(U1)*U12*inv(U2) CAblas::CMatrixLeftTrsM(n1,n2,a,offs,offs,true,false,0,a,offs,offs+n1); CAblas::CMatrixRightTrsM(n1,n2,a,offs+n1,offs+n1,true,false,0,a,offs,offs+n1); //--- Y :=inv(L2)*L12*inv(L1) CAblas::CMatrixLeftTrsM(n2,n1,a,offs+n1,offs+n1,false,true,0,a,offs+n1,offs); CAblas::CMatrixRightTrsM(n2,n1,a,offs,offs,false,true,0,a,offs+n1,offs); //--- W :=inv(L1*U1)+X*Y CMatrixLUInverseRec(a,offs,n1,work,info,rep); //--- check if(info<=0) return; CAblas::CMatrixGemm(n1,n1,n2,One,a,offs,offs+n1,0,a,offs+n1,offs,0,One,a,offs,offs); //--- X :=-X*inv(L2) //--- Y :=-inv(U2)*Y CAblas::CMatrixRightTrsM(n1,n2,a,offs+n1,offs+n1,false,true,0,a,offs,offs+n1); for(i=0;i0, contains matrix U*Q. | //| VT - if NCVT>0, contains matrix (P^T)*VT. | //| C - if NCC>0, contains matrix Q'*C. | //| Result: | //| True, if the algorithm has converged. | //| False, if the algorithm hasn't converged (rare case). | //| Additional information: | //| The type of convergence is controlled by the internal | //| parameter TOL. If the parameter is greater than 0, the | //| singular values will have relative accuracy TOL. If TOL<0, | //| the singular values will have absolute accuracy | //| ABS(TOL)*norm(B). By default, |TOL| falls within the range of| //| 10*Epsilon and 100*Epsilon, where Epsilon is the machine | //| precision. It is not recommended to use TOL less than | //| 10*Epsilon since this will considerably slow down the | //| algorithm and may not lead to error decreasing. | //| History: | //| * 31 March, 2007. | //| changed MAXITR from 6 to 12. | //| -- LAPACK routine (version 3.0) -- | //| Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., | //| Courant Institute, Argonne National Lab, and Rice University| //| October 31, 1999. | //+------------------------------------------------------------------+ static bool CBdSingValueDecompose::RMatrixBdSVD(double &d[],double &ce[], const int n,const bool isupper, const bool isfractionalaccuracyrequired, CMatrixDouble &u,const int nru, CMatrixDouble &c,const int ncc, CMatrixDouble &vt,const int ncvt) { //--- create variables bool result; //--- create arrays double d1[]; double e1[]; int i_=0; int i1_=0; //--- create copy double e[]; ArrayCopy(e,ce); //--- allocation ArrayResizeAL(d1,n+1); //--- change values i1_=-1; for(i_=1;i_<=n;i_++) d1[i_]=d[i_+i1_]; //--- check if(n>1) { //--- allocation ArrayResizeAL(e1,n); //--- change values i1_=-1; for(i_=1;i_0) for(i_=vstart;i_<=vstart+ncvt-1;i_++) vt[vstart].Set(i_,-1*vt[vstart][i_]); } //--- return result return(result); } //--- initialization ll=0; oldsn=0; //--- allocation ArrayResizeAL(work0,n); ArrayResizeAL(work1,n); ArrayResizeAL(work2,n); ArrayResizeAL(work3,n); uend=ustart+(int)MathMax(nru-1,0); vend=vstart+(int)MathMax(ncvt-1,0); cend=cstart+(int)MathMax(ncc-1,0); ArrayResizeAL(utemp,uend+1); ArrayResizeAL(vttemp,vend+1); ArrayResizeAL(ctemp,cend+1); //--- initialization maxitr=12; rightside=true; fwddir=true; //--- resize E from N-1 to N ArrayResizeAL(etemp,n+1); for(i=1;i0) { //--- function call CRotations::ApplyRotationsFromTheRight(fwddir,ustart,uend,1+ustart-1,n+ustart-1,work0,work1,u,utemp); } //--- check if(ncc>0) { //--- function call CRotations::ApplyRotationsFromTheLeft(fwddir,1+cstart-1,n+cstart-1,cstart,cend,work0,work1,c,ctemp); } } //--- Compute singular values to relative accuracy TOL //--- (By setting TOL to be negative, algorithm will compute //--- singular values to absolute accuracy ABS(TOL)*norm(input matrix)) tolmul=MathMax(10,MathMin(100,MathPow(eps,-0.125))); tol=tolmul*eps; //--- check if(!isfractionalaccuracyrequired) tol=-tol; //--- Compute approximate maximum, minimum singular values smax=0; for(i=1;i<=n;i++) smax=MathMax(smax,MathAbs(d[i])); for(i=1;i=0.0) { //--- Relative accuracy desired sminoa=MathAbs(d[1]); //--- check if(sminoa!=0.0) { mu=sminoa; for(i=2;i<=n;i++) { mu=MathAbs(d[i])*(mu/(mu+MathAbs(e[i-1]))); sminoa=MathMin(sminoa,mu); //--- check if(sminoa==0.0) break; } } //--- change values sminoa=sminoa/MathSqrt(n); thresh=MathMax(tol*sminoa,maxitr*n*n*unfl); } else { //--- Absolute accuracy desired thresh=MathMax(MathAbs(tol)*smax,maxitr*n*n*unfl); } //--- Prepare for main iteration loop for the singular values //--- (MAXIT is the maximum number of passes through the inner //--- loop permitted before nonconvergence signalled.) maxit=maxitr*n*n; iter=0; oldll=-1; oldm=-1; //--- M points to last element of unconverged part of matrix m=n; //--- Begin main iteration loop while(true) { //--- Check for convergence or exceeding iteration count if(m<=1) break; //--- check if(iter>maxit) return(false); //--- Find diagonal block of matrix to work on if(tol<0.0 && MathAbs(d[m])<=thresh) d[m]=0; //--- change values smax=MathAbs(d[m]); smin=smax; matrixsplitflag=false; for(lll=1;lll<=m-1;lll++) { ll=m-lll; abss=MathAbs(d[ll]); abse=MathAbs(e[ll]); //--- check if(tol<0.0 && abss<=thresh) d[ll]=0; //--- check if(abse<=thresh) { matrixsplitflag=true; break; } //--- change values smin=MathMin(smin,abss); smax=MathMax(smax,MathMax(abss,abse)); } //--- check if(!matrixsplitflag) ll=0; else { //--- Matrix splits since E(LL) = 0 e[ll]=0; //--- check if(ll==m-1) { //--- Convergence of bottom singular value, return to top of loop m=m-1; continue; } } ll=ll+1; //--- E(LL) through E(M-1) are nonzero, E(LL-1) is zero if(ll==m-1) { //--- 2 by 2 block, handle separately SVDV2x2(d[m-1],e[m-1],d[m],sigmn,sigmx,sinr,cosr,sinl,cosl); d[m-1]=sigmx; e[m-1]=0; d[m]=sigmn; //--- Compute singular vectors, if desired if(ncvt>0) { mm0=m+(vstart-1); mm1=m-1+(vstart-1); //--- swap for(i_=vstart;i_<=vend;i_++) vttemp[i_]=cosr*vt[mm1][i_]; for(i_=vstart;i_<=vend;i_++) vttemp[i_]=vttemp[i_]+sinr*vt[mm0][i_]; for(i_=vstart;i_<=vend;i_++) vt[mm0].Set(i_,cosr*vt[mm0][i_]); for(i_=vstart;i_<=vend;i_++) vt[mm0].Set(i_,vt[mm0][i_]-sinr*vt[mm1][i_]); for(i_=vstart;i_<=vend;i_++) vt[mm1].Set(i_,vttemp[i_]); } //--- check if(nru>0) { mm0=m+ustart-1; mm1=m-1+ustart-1; //--- swap for(i_=ustart;i_<=uend;i_++) utemp[i_]=cosl*u[i_][mm1]; for(i_=ustart;i_<=uend;i_++) utemp[i_]=utemp[i_]+sinl*u[i_][mm0]; for(i_=ustart;i_<=uend;i_++) u[i_].Set(mm0,cosl*u[i_][mm0]); for(i_=ustart;i_<=uend;i_++) u[i_].Set(mm0,u[i_][mm0]-sinl*u[i_][mm1]); for(i_=ustart;i_<=uend;i_++) u[i_].Set(mm1,utemp[i_]); } //--- check if(ncc>0) { mm0=m+cstart-1; mm1=m-1+cstart-1; //--- swap for(i_=cstart;i_<=cend;i_++) ctemp[i_]=cosl*c[mm1][i_]; for(i_=cstart;i_<=cend;i_++) ctemp[i_]=ctemp[i_]+sinl*c[mm0][i_]; for(i_=cstart;i_<=cend;i_++) c[mm0].Set(i_,cosl*c[mm0][i_]); for(i_=cstart;i_<=cend;i_++) c[mm0].Set(i_,c[mm0][i_]-sinl*c[mm1][i_]); for(i_=cstart;i_<=cend;i_++) c[mm1].Set(i_,ctemp[i_]); } m=m-2; continue; } //--- If working on new submatrix, choose shift direction //--- (from larger end diagonal element towards smaller) //--- Previously was //--- "if (LL>OLDM) or (M //--- Very strange that LAPACK still contains it. bchangedir=false; //--- check if(idir==1 && MathAbs(d[ll])<1.0E-3*MathAbs(d[m])) bchangedir=true; //--- check if(idir==2 && MathAbs(d[m])<1.0E-3*MathAbs(d[ll])) bchangedir=true; //--- check if(ll!=oldll || m!=oldm || bchangedir) { //--- check if(MathAbs(d[ll])>=MathAbs(d[m])) { //--- Chase bulge from top(big end) to bottom(small end) idir=1; } else { //--- Chase bulge from bottom (big end) to top (small end) idir=2; } } //--- Apply convergence tests if(idir==1) { //--- Run convergence test in forward direction //--- First apply standard test to bottom of matrix if(MathAbs(e[m-1])<=MathAbs(tol)*MathAbs(d[m]) || (tol<0.0 && MathAbs(e[m-1])<=thresh)) { e[m-1]=0; continue; } //--- check if(tol>=0.0) { //--- If relative accuracy desired, //--- apply convergence criterion forward mu=MathAbs(d[ll]); sminl=mu; iterflag=false; for(lll=ll;lll<=m-1;lll++) { //--- check if(MathAbs(e[lll])<=tol*mu) { e[lll]=0; iterflag=true; //--- break the cycle break; } //--- change values sminlo=sminl; mu=MathAbs(d[lll+1])*(mu/(mu+MathAbs(e[lll]))); sminl=MathMin(sminl,mu); } //--- check if(iterflag) continue; } } else { //--- Run convergence test in backward direction //--- First apply standard test to top of matrix if(MathAbs(e[ll])<=MathAbs(tol)*MathAbs(d[ll]) || (tol<0.0 && MathAbs(e[ll])<=thresh)) { e[ll]=0; continue; } //--- check if(tol>=0.0) { //--- If relative accuracy desired, //--- apply convergence criterion backward mu=MathAbs(d[m]); sminl=mu; iterflag=false; for(lll=m-1;lll>=ll;lll--) { //--- check if(MathAbs(e[lll])<=(double)(tol*mu)) { e[lll]=0; iterflag=true; //--- break the cycle break; } sminlo=sminl; mu=MathAbs(d[lll])*(mu/(mu+MathAbs(e[lll]))); sminl=MathMin(sminl,mu); } //--- check if(iterflag) continue; } } //--- change values oldll=ll; oldm=m; //--- Compute shift. First, test if shifting would ruin relative //--- accuracy, and if so set the shift to zero. if(tol>=0.0 && n*tol*(sminl/smax)<=MathMax(eps,0.01*tol)) { //--- Use a zero shift to avoid loss of relative accuracy shift=0; } else { //--- Compute the shift from 2-by-2 block at end of matrix if(idir==1) { sll=MathAbs(d[ll]); SVD2x2(d[m-1],e[m-1],d[m],shift,r); } else { sll=MathAbs(d[m]); SVD2x2(d[ll],e[ll],d[ll+1],shift,r); } //--- Test if shift negligible, and if so set to zero if(sll>0.0) { //--- check if(CMath::Sqr(shift/sll)ll) e[i-1]=oldsn*r; //--- function call CRotations::GenerateRotation(oldcs*r,d[i+1]*sn,oldcs,oldsn,tmp); //--- change values d[i]=tmp; work0[i-ll+1]=cs; work1[i-ll+1]=sn; work2[i-ll+1]=oldcs; work3[i-ll+1]=oldsn; } //--- change values h=d[m]*cs; d[m]=h*oldcs; e[m-1]=h*oldsn; //--- Update singular vectors if(ncvt>0) { //--- function call CRotations::ApplyRotationsFromTheLeft(fwddir,ll+vstart-1,m+vstart-1,vstart,vend,work0,work1,vt,vttemp); } //--- check if(nru>0) { //--- function call CRotations::ApplyRotationsFromTheRight(fwddir,ustart,uend,ll+ustart-1,m+ustart-1,work2,work3,u,utemp); } //--- check if(ncc>0) { //--- function call CRotations::ApplyRotationsFromTheLeft(fwddir,ll+cstart-1,m+cstart-1,cstart,cend,work2,work3,c,ctemp); } //--- Check for convergence if(MathAbs(e[m-1])<=thresh) e[m-1]=0; } else { //--- Chase bulge from bottom to top //--- Save cosines and sines for later singular vector updates cs=1; oldcs=1; for(i=m;i>=ll+1;i--) { //--- function call CRotations::GenerateRotation(d[i]*cs,e[i-1],cs,sn,r); //--- check if(i0) { //--- function call CRotations::ApplyRotationsFromTheLeft(!fwddir,ll+vstart-1,m+vstart-1,vstart,vend,work2,work3,vt,vttemp); } //--- check if(nru>0) { //--- function call CRotations::ApplyRotationsFromTheRight(!fwddir,ustart,uend,ll+ustart-1,m+ustart-1,work0,work1,u,utemp); } //--- check if(ncc>0) { //--- function call CRotations::ApplyRotationsFromTheLeft(!fwddir,ll+cstart-1,m+cstart-1,cstart,cend,work0,work1,c,ctemp); } //--- Check for convergence if(MathAbs(e[ll])<=thresh) e[ll]=0; } } else { //--- Use nonzero shift if(idir==1) { //--- Chase bulge from top to bottom //--- Save cosines and sines for later singular vector updates f=(MathAbs(d[ll])-shift)*(ExtSignBdSQR(1,d[ll])+shift/d[ll]); g=e[ll]; for(i=ll;ill) e[i-1]=r; //--- change values f=cosr*d[i]+sinr*e[i]; e[i]=cosr*e[i]-sinr*d[i]; g=sinr*d[i+1]; d[i+1]=cosr*d[i+1]; //--- function call CRotations::GenerateRotation(f,g,cosl,sinl,r); //--- change values d[i]=r; f=cosl*e[i]+sinl*d[i+1]; d[i+1]=cosl*d[i+1]-sinl*e[i]; //--- check if(i0) { //--- function call CRotations::ApplyRotationsFromTheLeft(fwddir,ll+vstart-1,m+vstart-1,vstart,vend,work0,work1,vt,vttemp); } //--- check if(nru>0) { //--- function call CRotations::ApplyRotationsFromTheRight(fwddir,ustart,uend,ll+ustart-1,m+ustart-1,work2,work3,u,utemp); } //--- check if(ncc>0) { //--- function call CRotations::ApplyRotationsFromTheLeft(fwddir,ll+cstart-1,m+cstart-1,cstart,cend,work2,work3,c,ctemp); } //--- Check for convergence if(MathAbs(e[m-1])<=thresh) e[m-1]=0; } else { //--- Chase bulge from bottom to top //--- Save cosines and sines for later singular vector updates f=(MathAbs(d[m])-shift)*(ExtSignBdSQR(1,d[m])+shift/d[m]); g=e[m-1]; for(i=m;i>=ll+1;i--) { //--- function call CRotations::GenerateRotation(f,g,cosr,sinr,r); //--- check if(ill+1) { g=sinl*e[i-2]; e[i-2]=cosl*e[i-2]; } //--- change values work0[i-ll]=cosr; work1[i-ll]=-sinr; work2[i-ll]=cosl; work3[i-ll]=-sinl; } e[ll]=f; //--- Check for convergence if(MathAbs(e[ll])<=thresh) e[ll]=0; //--- Update singular vectors if desired if(ncvt>0) { //--- function call CRotations::ApplyRotationsFromTheLeft(!fwddir,ll+vstart-1,m+vstart-1,vstart,vend,work2,work3,vt,vttemp); } //--- check if(nru>0) { //--- function call CRotations::ApplyRotationsFromTheRight(!fwddir,ustart,uend,ll+ustart-1,m+ustart-1,work0,work1,u,utemp); } //--- check if(ncc>0) { //--- function call CRotations::ApplyRotationsFromTheLeft(!fwddir,ll+cstart-1,m+cstart-1,cstart,cend,work0,work1,c,ctemp); } } } //--- QR iteration finished, go back and check convergence continue; } //--- All singular values converged, so make them positive for(i=1;i<=n;i++) { //--- check if(d[i]<0.0) { d[i]=-d[i]; //--- Change sign of singular vectors, if desired if(ncvt>0) { for(i_=vstart;i_<=vend;i_++) vt[i+vstart-1].Set(i_,-1*vt[i+vstart-1][i_]); } } } //--- Sort the singular values into decreasing order (insertion sort on //--- singular values, but only one transposition per singular vector) for(i=1;i0) { j=n+1-i; //--- swap for(i_=vstart;i_<=vend;i_++) vttemp[i_]=vt[isub+vstart-1][i_]; for(i_=vstart;i_<=vend;i_++) vt[isub+vstart-1].Set(i_,vt[j+vstart-1][i_]); for(i_=vstart;i_<=vend;i_++) vt[j+vstart-1].Set(i_,vttemp[i_]); } if(nru>0) { j=n+1-i; //--- swap for(i_=ustart;i_<=uend;i_++) utemp[i_]=u[i_][isub+ustart-1]; for(i_=ustart;i_<=uend;i_++) u[i_].Set(isub+ustart-1,u[i_][j+ustart-1]); for(i_=ustart;i_<=uend;i_++) u[i_].Set(j+ustart-1,utemp[i_]); } //--- check if(ncc>0) { j=n+1-i; //--- swap for(i_=cstart;i_<=cend;i_++) ctemp[i_]=c[isub+cstart-1][i_]; for(i_=cstart;i_<=cend;i_++) c[isub+cstart-1].Set(i_,c[j+cstart-1][i_]); for(i_=cstart;i_<=cend;i_++) c[j+cstart-1].Set(i_,ctemp[i_]); } } } //--- return result return(result); } //+------------------------------------------------------------------+ //| Internal subroutine | //+------------------------------------------------------------------+ static double CBdSingValueDecompose::ExtSignBdSQR(const double a, const double b) { //--- create a variable double result=0; //--- check if(b>=0.0) result=MathAbs(a); else result=-MathAbs(a); //--- return result return(result); } //+------------------------------------------------------------------+ //| Internal subroutine | //+------------------------------------------------------------------+ static void CBdSingValueDecompose::SVD2x2(const double f,const double g, const double h,double &ssmin, double &ssmax) { //--- create variables double aas=0; double at=0; double au=0; double c=0; double fa=0; double fhmn=0; double fhmx=0; double ga=0; double ha=0; //--- initialization ssmin=0; ssmax=0; fa=MathAbs(f); ga=MathAbs(g); ha=MathAbs(h); fhmn=MathMin(fa,ha); fhmx=MathMax(fa,ha); //--- check if(fhmn==0.0) { ssmin=0; //--- check if(fhmx==0.0) ssmax=ga; else ssmax=MathMax(fhmx,ga)*MathSqrt(1+CMath::Sqr(MathMin(fhmx,ga)/MathMax(fhmx,ga))); } else { //--- check if(gafa; //--- check if(swp) { //--- Now FA .ge. HA pmax=3; temp=ft; ft=ht; ht=temp; temp=fa; fa=ha; ha=temp; } gt=g; ga=MathAbs(gt); //--- check if(ga==0.0) { //--- Diagonal matrix ssmin=ha; ssmax=fa; clt=1; crt=1; slt=0; srt=0; } else { gasmal=true; //--- check if(ga>fa) { pmax=2; //--- check if(fa/ga1.0) { v=ga/ha; ssmin=fa/v; } else { v=fa/ga; ssmin=v*ha; } //--- change values clt=1; slt=ht/gt; srt=1; crt=ft/gt; } } //--- check if(gasmal) { //--- Normal case d=fa-ha; //--- check if(d==fa) l=1; else l=d/fa; //--- change values m=gt/ft; t=2-l; mm=m*m; tt=t*t; s=MathSqrt(tt+mm); //--- check if(l==0.0) r=MathAbs(m); else r=MathSqrt(l*l+mm); //--- change values a=0.5*(s+r); ssmin=ha/a; ssmax=fa*a; //--- check if(mm==0.0) { //--- Note that M is very tiny if(l==0.0) t=ExtSignBdSQR(2,ft)*ExtSignBdSQR(1,gt); else t=gt/ExtSignBdSQR(d,ft)+m/t; } else t=(m/(s+t)+m/(r+l))*(1+a); //--- change values l=MathSqrt(t*t+4); crt=2/l; srt=t/l; clt=(crt+srt*m)/a; v=ht/ft; slt=v*srt/a; } } //--- check if(swp) { csl=srt; snl=crt; csr=slt; snr=clt; } else { csl=clt; snl=slt; csr=crt; snr=srt; } //--- Correct signs of SSMAX and SSMIN if(pmax==1) tsign=ExtSignBdSQR(1,csr)*ExtSignBdSQR(1,csl)*ExtSignBdSQR(1,f); //--- check if(pmax==2) tsign=ExtSignBdSQR(1,snr)*ExtSignBdSQR(1,csl)*ExtSignBdSQR(1,g); //--- check if(pmax==3) tsign=ExtSignBdSQR(1,snr)*ExtSignBdSQR(1,snl)*ExtSignBdSQR(1,h); //--- get result ssmax=ExtSignBdSQR(ssmax,tsign); ssmin=ExtSignBdSQR(ssmin,tsign*ExtSignBdSQR(1,f)*ExtSignBdSQR(1,h)); } //+------------------------------------------------------------------+ //| Singular value decomposition | //+------------------------------------------------------------------+ class CSingValueDecompose { public: //--- constructor, destructor CSingValueDecompose(void); ~CSingValueDecompose(void); //--- method static bool RMatrixSVD(CMatrixDouble &ca,const int m,const int n,const int uneeded,const int vtneeded,const int additionalmemory,double &w[],CMatrixDouble &u,CMatrixDouble &vt); }; //+------------------------------------------------------------------+ //| Constructor without parameters | //+------------------------------------------------------------------+ CSingValueDecompose::CSingValueDecompose(void) { } //+------------------------------------------------------------------+ //| Destructor | //+------------------------------------------------------------------+ CSingValueDecompose::~CSingValueDecompose(void) { } //+------------------------------------------------------------------+ //| Singular value decomposition of a rectangular matrix. | //| The algorithm calculates the singular value decomposition of a | //| matrix of size MxN: A = U * S * V^T | //| The algorithm finds the singular values and, optionally, matrices| //| U and V^T. The algorithm can find both first min(M,N) columns of | //| matrix U and rows of matrix V^T (singular vectors), and matrices | //| U and V^T wholly (of sizes MxM and NxN respectively). | //| Take into account that the subroutine does not return matrix V | //| but V^T. | //| Input parameters: | //| A - matrix to be decomposed. | //| Array whose indexes range within | //| [0..M-1, 0..N-1]. | //| M - number of rows in matrix A. | //| N - number of columns in matrix A. | //| UNeeded - 0, 1 or 2. See the description of the | //| parameter U. | //| VTNeeded - 0, 1 or 2. See the description of the | //| parameter VT. | //| AdditionalMemory - | //| If the parameter: | //| * equals 0, the algorithm doesn?t use | //| additional memory (lower requirements, | //| lower performance). | //| * equals 1, the algorithm uses additional | //| memory of size min(M,N)*min(M,N) of real | //| numbers. It often speeds up the algorithm.| //| * equals 2, the algorithm uses additional | //| memory of size M*min(M,N) of real numbers.| //| It allows to get a maximum performance. | //| The recommended value of the parameter is 2. | //| Output parameters: | //| W - contains singular values in descending order.| //| U - if UNeeded=0, U isn't changed, the left | //| singular vectors are not calculated. | //| if Uneeded=1, U contains left singular | //| vectors (first min(M,N) columns of matrix U).| //| Array whose indexes range within | //| [0..M-1, 0..Min(M,N)-1]. if UNeeded=2, U | //| contains matrix U wholly. Array whose indexes| //| range within [0..M-1, 0..M-1]. | //| VT - if VTNeeded=0, VT isn?t changed, the right | //| singular vectors are not calculated. | //| if VTNeeded=1, VT contains right singular | //| vectors (first min(M,N) rows of matrix V^T). | //| Array whose indexes range within | //| [0..min(M,N)-1, 0..N-1]. if VTNeeded=2, VT | //| contains matrix V^T wholly. Array whose | //| indexes range within [0..N-1, 0..N-1]. | //+------------------------------------------------------------------+ static bool CSingValueDecompose::RMatrixSVD(CMatrixDouble &ca,const int m, const int n,const int uneeded, const int vtneeded, const int additionalmemory, double &w[],CMatrixDouble &u, CMatrixDouble &vt) { //--- create variables bool result; bool isupper; int minmn=0; int ncu=0; int nrvt=0; int nru=0; int ncvt=0; int i=0; int j=0; //--- create arrays double tauq[]; double taup[]; double tau[]; double e[]; double work[]; //--- create matrix CMatrixDouble t2; //--- create copy CMatrixDouble a; a=ca; //--- initialization result=true; //--- check if(m==0 || n==0) return(true); //--- check if(!CAp::Assert(uneeded>=0 && uneeded<=2,__FUNCTION__+": wrong parameters!")) return(false); //--- check if(!CAp::Assert(vtneeded>=0 && vtneeded<=2,__FUNCTION__+": wrong parameters!")) return(false); //--- check if(!CAp::Assert(additionalmemory>=0 && additionalmemory<=2,__FUNCTION__+": wrong parameters!")) return(false); //--- initialization minmn=MathMin(m,n); ArrayResizeAL(w,minmn+1); ncu=0; nru=0; //--- check if(uneeded==1) { nru=m; ncu=minmn; u.Resize(nru,ncu); } //--- check if(uneeded==2) { nru=m; ncu=m; u.Resize(nru,ncu); } nrvt=0; ncvt=0; //--- check if(vtneeded==1) { nrvt=minmn; ncvt=n; vt.Resize(nrvt,ncvt); } //--- check if(vtneeded==2) { nrvt=n; ncvt=n; vt.Resize(nrvt,ncvt); } //--- M much larger than N //--- Use bidiagonal reduction with QR-decomposition if((double)(m)>(double)(1.6*n)) { //--- check if(uneeded==0) { //-- No left singular vectors to be computed COrtFac::RMatrixQR(a,m,n,tau); for(i=0;i(double)(1.6*m)) { //--- check if(vtneeded==0) { //--- No right singular vectors to be computed COrtFac::RMatrixLQ(a,m,n,tau); for(i=0;i<=m-1;i++) { for(j=i+1;j<=m-1;j++) a[i].Set(j,0); } //--- function call COrtFac::RMatrixBD(a,m,m,tauq,taup); //--- function call COrtFac::RMatrixBDUnpackQ(a,m,m,tauq,ncu,u); //--- function call COrtFac::RMatrixBDUnpackDiagonals(a,m,m,isupper,w,e); ArrayResizeAL(work,m+1); //--- function call CBlas::InplaceTranspose(u,0,nru-1,0,ncu-1,work); //--- get result result=CBdSingValueDecompose::RMatrixBdSVD(w,e,m,isupper,false,a,0,u,nru,vt,0); //--- function call CBlas::InplaceTranspose(u,0,nru-1,0,ncu-1,work); //--- return result return(result); } else { //--- Right singular vectors (may be full matrix VT) to be computed COrtFac::RMatrixLQ(a,m,n,tau); //--- function call COrtFac::RMatrixLQUnpackQ(a,m,n,tau,nrvt,vt); for(i=0;i<=m-1;i++) { for(j=i+1;j<=m-1;j++) a[i].Set(j,0); } //--- function call COrtFac::RMatrixBD(a,m,m,tauq,taup); //--- function call COrtFac::RMatrixBDUnpackQ(a,m,m,tauq,ncu,u); //--- function call COrtFac::RMatrixBDUnpackDiagonals(a,m,m,isupper,w,e); ArrayResizeAL(work,MathMax(m,n)+1); //--- function call CBlas::InplaceTranspose(u,0,nru-1,0,ncu-1,work); //--- check if(additionalmemory<1) { //--- No additional memory available COrtFac::RMatrixBDMultiplyByP(a,m,m,taup,vt,m,n,false,true); //--- get result result=CBdSingValueDecompose::RMatrixBdSVD(w,e,m,isupper,false,a,0,u,nru,vt,n); } else { //--- Large VT. Transforming intermediate matrix T2 COrtFac::RMatrixBDUnpackPT(a,m,m,taup,m,t2); //--- get result result=CBdSingValueDecompose::RMatrixBdSVD(w,e,m,isupper,false,a,0,u,nru,t2,m); //--- function call CBlas::CopyMatrix(vt,0,m-1,0,n-1,a,0,m-1,0,n-1); //--- function call CBlas::MatrixMatrixMultiply(t2,0,m-1,0,m-1,false,a,0,m-1,0,n-1,false,1.0,vt,0,m-1,0,n-1,0.0,work); } //--- function call CBlas::InplaceTranspose(u,0,nru-1,0,ncu-1,work); //--- return result return(result); } } //--- M<=N //--- We can use inplace transposition of U to get rid of columnwise operations if(m<=n) { //--- function call COrtFac::RMatrixBD(a,m,n,tauq,taup); //--- function call COrtFac::RMatrixBDUnpackQ(a,m,n,tauq,ncu,u); //--- function call COrtFac::RMatrixBDUnpackPT(a,m,n,taup,nrvt,vt); //--- function call COrtFac::RMatrixBDUnpackDiagonals(a,m,n,isupper,w,e); ArrayResizeAL(work,m+1); //--- function call CBlas::InplaceTranspose(u,0,nru-1,0,ncu-1,work); //--- get result result=CBdSingValueDecompose::RMatrixBdSVD(w,e,minmn,isupper,false,a,0,u,nru,vt,ncvt); //--- function call CBlas::InplaceTranspose(u,0,nru-1,0,ncu-1,work); //--- return result return(result); } //--- Simple bidiagonal reduction COrtFac::RMatrixBD(a,m,n,tauq,taup); //--- function call COrtFac::RMatrixBDUnpackQ(a,m,n,tauq,ncu,u); //--- function call COrtFac::RMatrixBDUnpackPT(a,m,n,taup,nrvt,vt); //--- function call COrtFac::RMatrixBDUnpackDiagonals(a,m,n,isupper,w,e); //--- check if(additionalmemory<2 || uneeded==0) { //--- We cant use additional memory or there is no need in such operations result=CBdSingValueDecompose::RMatrixBdSVD(w,e,minmn,isupper,false,u,nru,a,0,vt,ncvt); } else { //--- We can use additional memory t2.Resize(minmn,m); //--- function call CBlas::CopyAndTranspose(u,0,m-1,0,minmn-1,t2,0,minmn-1,0,m-1); //--- get result result=CBdSingValueDecompose::RMatrixBdSVD(w,e,minmn,isupper,false,u,0,t2,m,vt,ncvt); //--- function call CBlas::CopyAndTranspose(t2,0,minmn-1,0,m-1,u,0,m-1,0,minmn-1); } //--- return result return(result); } //+------------------------------------------------------------------+ //| Structure which stores state of linear CG solver between | //| subsequent calls of FBLSCgIteration(). Initialized with | //| FBLSCGCreate(). | //| USAGE: | //| 1. call to FBLSCGCreate() | //| 2. F:=FBLSCgIteration(State) | //| 3. if F is False, iterations are over | //| 4. otherwise, fill State.AX with A*x, State.XAX with x'*A*x | //| 5. goto 2 | //| If you want to rerminate iterations, pass zero or negative value | //| to XAX. | //| FIELDS: | //| E1 - 2-norm of residual at the start | //| E2 - 2-norm of residual at the end | //| X - on return from FBLSCgIteration() it contains | //| vector for matrix-vector product | //| AX - must be filled with A*x if FBLSCgIteration() | //| returned True | //| XAX - must be filled with x'*A*x | //| XK - contains result (if FBLSCgIteration() returned | //| False) | //| Other fields are private and should not be used by outsiders. | //+------------------------------------------------------------------+ class CFblsLinCgState { public: //--- variables double m_e1; double m_e2; double m_x[]; double m_ax[]; double m_xax; double m_xk[]; int m_n; double m_rk[]; double m_rk1[]; double m_xk1[]; double m_pk[]; double m_pk1[]; double m_b[]; RCommState m_rstate; double m_tmp2[]; //--- constructor, destructor CFblsLinCgState(void); ~CFblsLinCgState(void); }; //+------------------------------------------------------------------+ //| Constructor without parameters | //+------------------------------------------------------------------+ CFblsLinCgState::CFblsLinCgState(void) { } //+------------------------------------------------------------------+ //| Destructor | //+------------------------------------------------------------------+ CFblsLinCgState::~CFblsLinCgState(void) { } //+------------------------------------------------------------------+ //| Fast basic linear solutions | //+------------------------------------------------------------------+ class CFbls { private: //--- auxiliary functions for FblsCGiteration static void Func_lbl_rcomm(CFblsLinCgState &state,int n,int k,double rk2,double rk12,double pap,double s,double betak,double v1,double v2); static bool Func_lbl_3(CFblsLinCgState &state,int &n,int &k,double &rk2,double &rk12,double &pap,double &s,double &betak,double &v1,double &v2); static bool Func_lbl_5(CFblsLinCgState &state,int &n,int &k,double &rk2,double &rk12,double &pap,double &s,double &betak,double &v1,double &v2); public: //--- constructor, destructor CFbls(void); ~CFbls(void); //--- methods static void FblsCholeskySolve(CMatrixDouble &cha,const double sqrtscalea,const int n,const bool isupper,double &xb[],double &tmp[]); static void FblsSolveCGx(CMatrixDouble &a,const int m,const int n,const double alpha,const double &b[],double &x[],double &buf[]); static void FblsCGCreate(double &x[],double &b[],const int n,CFblsLinCgState &state); static bool FblsCGIteration(CFblsLinCgState &state); }; //+------------------------------------------------------------------+ //| Constructor without parameters | //+------------------------------------------------------------------+ CFbls::CFbls(void) { } //+------------------------------------------------------------------+ //| Destructor | //+------------------------------------------------------------------+ CFbls::~CFbls(void) { } //+------------------------------------------------------------------+ //| Basic Cholesky solver for ScaleA*Cholesky(A)'*x = y. | //| This subroutine assumes that: | //| * A*ScaleA is well scaled | //| * A is well-conditioned, so no zero divisions or overflow may | //| occur | //| INPUT PARAMETERS: | //| CHA - Cholesky decomposition of A | //| SqrtScaleA- square root of scale factor ScaleA | //| N - matrix size | //| IsUpper - storage type | //| XB - right part | //| Tmp - buffer; function automatically allocates it, if | //| it is too small. It can be reused if function is | //| called several times. | //| OUTPUT PARAMETERS: | //| XB - solution | //| NOTES: no assertion or tests are done during algorithm operation | //+------------------------------------------------------------------+ static void CFbls::FblsCholeskySolve(CMatrixDouble &cha,const double sqrtscalea, const int n,const bool isupper,double &xb[], double &tmp[]) { //--- create variables int i=0; double v=0; int i_=0; //--- check if(CAp::Len(tmp)=0;i--) { //--- check if(i0) { for(i_=0;i_=0;i--) { xb[i]=xb[i]/(sqrtscalea*cha[i][i]); //--- check if(i>0) { v=xb[i]; for(i_=0;i_0 is a scalar | //| * I is NxN identity matrix | //| * b is Nx1 vector | //| * X is Nx1 unknown vector. | //| N iterations of linear conjugate gradient are used to solve | //| problem. | //| INPUT PARAMETERS: | //| A - array[M,N], matrix | //| M - number of rows | //| N - number of unknowns | //| B - array[N], right part | //| X - initial approxumation, array[N] | //| Buf - buffer; function automatically allocates it, if it | //| is too small. It can be reused if function is called | //| several times with same M and N. | //| OUTPUT PARAMETERS: | //| X - improved solution | //| NOTES: | //| * solver checks quality of improved solution. If (because of | //| problem condition number, numerical noise, etc.) new solution| //| is WORSE than original approximation, then original | //| approximation is returned. | //| * solver assumes that both A, B, Alpha are well scaled (i.e. | //| they are less than sqrt(overflow) and greater than | //| sqrt(underflow)). | //+------------------------------------------------------------------+ static void CFbls::FblsSolveCGx(CMatrixDouble &a,const int m,const int n, const double alpha,const double &b[], double &x[],double &buf[]) { //--- create variables int k=0; int offsrk=0; int offsrk1=0; int offsxk=0; int offsxk1=0; int offspk=0; int offspk1=0; int offstmp1=0; int offstmp2=0; int bs=0; double e1=0; double e2=0; double rk2=0; double rk12=0; double pap=0; double s=0; double betak=0; double v1=0; double v2=0; int i_=0; int i1_=0; //--- Test for special case: B=0 v1=0.0; for(i_=0;i_=0) { //--- initialization n=state.m_rstate.ia[0]; k=state.m_rstate.ia[1]; rk2=state.m_rstate.ra[0]; rk12=state.m_rstate.ra[1]; pap=state.m_rstate.ra[2]; s=state.m_rstate.ra[3]; betak=state.m_rstate.ra[4]; v1=state.m_rstate.ra[5]; v2=state.m_rstate.ra[6]; } else { //--- initialization n=-983; k=-989; rk2=-834; rk12=900; pap=-287; s=364; betak=214; v1=-338; v2=-686; } //--- check if(state.m_rstate.stage==0) { //--- copy for(i_=0;i_n-1) return(Func_lbl_5(state,n,k,rk2,rk12,pap,s,betak,v1,v2)); //--- Calculate A*p(k) - store in State.Tmp2 //--- and p(k)'*A*p(k) - store in PAP //--- If PAP=0,break (iteration is over) for(int i_=0;i_=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- check if(!CAp::Assert(CAp::Len(pivots)>=n,__FUNCTION__+": Pivots array is too short!")) return(EMPTY_VALUE); //--- check if(!CAp::Assert(CAp::Rows(a)>=n,__FUNCTION__+": rows(A)=n,__FUNCTION__+": cols(A)=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- check if(!CAp::Assert(CAp::Rows(a)>=n,__FUNCTION__+": rows(A)=n,__FUNCTION__+": cols(A)=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- check if(!CAp::Assert(CAp::Len(pivots)>=n,__FUNCTION__+": Pivots array is too short!")) return(EMPTY_VALUE); //--- check if(!CAp::Assert(CAp::Rows(a)>=n,__FUNCTION__+": rows(A)=n,__FUNCTION__+": cols(A)=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- check if(!CAp::Assert(CAp::Rows(a)>=n,__FUNCTION__+": rows(A)=n,__FUNCTION__+": cols(A)=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- check if(!CAp::Assert(CAp::Rows(a)>=n,__FUNCTION__+": rows(A)=n,__FUNCTION__+": cols(A)=1,__FUNCTION__+": N<1!")) return(EMPTY_VALUE); //--- check if(!CAp::Assert(CAp::Rows(a)>=n,__FUNCTION__+": rows(A)=n,__FUNCTION__+": cols(A)0,__FUNCTION__+": N<=0!")) return(false); //--- check if(!CAp::Assert(problemtype==1 || problemtype==2 || problemtype==3,__FUNCTION__+": incorrect ProblemType!")) return(false); //--- initialization result=true; //--- Problem 1: A*x = lambda*B*x //--- Reducing to: //--- C*y = lambda*y //--- C = L^(-1) * A * L^(-T) //--- x = L^(-T) * y if(problemtype==1) { //--- Factorize B in T: B = LL' t.Resize(n,n); //--- check if(isupperb) { for(i=0;i=0 && updrow=0 && updcolumn=0. | //| Output parameters: | //| A - contains matrix T. | //| Array whose indexes range within [0..N-1, 0..N-1]. | //| S - contains Schur vectors. | //| Array whose indexes range within [0..N-1, 0..N-1]. | //| Note 1: | //| The block structure of matrix T can be easily recognized: | //| since all the elements below the blocks are zeros, the | //| elements a[i+1,i] which are equal to 0 show the block border.| //| Note 2: | //| The algorithm performance depends on the value of the | //| internal parameter NS of the InternalSchurDecomposition | //| subroutine which defines the number of shifts in the QR | //| algorithm (similarly to the block width in block-matrix | //| algorithms in linear algebra). If you require maximum | //| performance on your machine, it is recommended to adjust | //| this parameter manually. | //| Result: | //| True, | //| if the algorithm has converged and parameters A and S | //| contain the result. | //| False, | //| if the algorithm has not converged. | //| Algorithm implemented on the basis of the DHSEQR subroutine | //| (LAPACK 3.0 library). | //+------------------------------------------------------------------+ static bool CSchur::RMatrixSchur(CMatrixDouble &a,const int n,CMatrixDouble &s) { //--- create variables bool result; int info=0; int i=0; int j=0; //--- create arrays double tau[]; double wi[]; double wr[]; //--- create matrix CMatrixDouble a1; CMatrixDouble s1; //--- Upper Hessenberg form of the 0-based matrix COrtFac::RMatrixHessenberg(a,n,tau); COrtFac::RMatrixHessenbergUnpackQ(a,n,tau,s); //--- Convert from 0-based arrays to 1-based, //--- then call InternalSchurDecomposition //--- Awkward, of course, but Schur decompisiton subroutine //--- is too complex to fix it. a1.Resize(n+1,n+1); s1.Resize(n+1,n+1); for(i=1;i<=n;i++) { for(j=1;j<=n;j++) { a1[i].Set(j,a[i-1][j-1]); s1[i].Set(j,s[i-1][j-1]); } } //--- function call CHsSchur::InternalSchurDecomposition(a1,n,1,1,wr,wi,s1,info); result=info==0; //--- convert from 1-based arrays to -based for(i=1;i<=n;i++) { for(j=1;j<=n;j++) { a[i-1].Set(j-1,a1[i][j]); s[i-1].Set(j-1,s1[i][j]); } } //--- return result return(result); } //+------------------------------------------------------------------+