/************************************************************************* Copyright (c) 1992-2007 The University of Tennessee. All rights reserved. Contributors: * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to pseudocode. See subroutines comments for additional copyrights. >>> SOURCE LICENSE >>> This program 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. A copy of the GNU General Public License is available at http://www.fsf.org/licensing/licenses >>> END OF LICENSE >>> *************************************************************************/ using System; namespace alglib { public class trfac { /************************************************************************* LU decomposition of a general real 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)]. -- ALGLIB routine -- 10.01.2010 Bochkanov Sergey *************************************************************************/ public static void rmatrixlu(ref double[,] a, int m, int n, ref int[] pivots) { System.Diagnostics.Debug.Assert(m>0, "RMatrixLU: incorrect M!"); System.Diagnostics.Debug.Assert(n>0, "RMatrixLU: incorrect N!"); rmatrixplu(ref a, m, n, ref 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)]. -- ALGLIB routine -- 10.01.2010 Bochkanov Sergey *************************************************************************/ public static void cmatrixlu(ref AP.Complex[,] a, int m, int n, ref int[] pivots) { System.Diagnostics.Debug.Assert(m>0, "CMatrixLU: incorrect M!"); System.Diagnostics.Debug.Assert(n>0, "CMatrixLU: incorrect N!"); cmatrixplu(ref a, m, n, ref 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. -- ALGLIB routine -- 15.12.2009 Bochkanov Sergey *************************************************************************/ public static bool hpdmatrixcholesky(ref AP.Complex[,] a, int n, bool isupper) { bool result = new bool(); AP.Complex[] tmp = new AP.Complex[0]; if( n<1 ) { result = false; return result; } tmp = new AP.Complex[2*n]; result = hpdmatrixcholeskyrec(ref a, 0, n, isupper, ref tmp); return result; } /************************************************************************* 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. -- ALGLIB routine -- 15.12.2009 Bochkanov Sergey *************************************************************************/ public static bool spdmatrixcholesky(ref double[,] a, int n, bool isupper) { bool result = new bool(); double[] tmp = new double[0]; if( n<1 ) { result = false; return result; } tmp = new double[2*n]; result = spdmatrixcholeskyrec(ref a, 0, n, isupper, ref tmp); return result; } public static void rmatrixlup(ref double[,] a, int m, int n, ref int[] pivots) { double[] tmp = new double[0]; int i = 0; int j = 0; double mx = 0; double v = 0; int i_ = 0; // // Internal LU decomposition subroutine. // Never call it directly. // System.Diagnostics.Debug.Assert(m>0, "RMatrixLUP: incorrect M!"); System.Diagnostics.Debug.Assert(n>0, "RMatrixLUP: incorrect N!"); // // Scale matrix to avoid overflows, // decompose it, then scale back. // mx = 0; for(i=0; i<=m-1; i++) { for(j=0; j<=n-1; j++) { mx = Math.Max(mx, Math.Abs(a[i,j])); } } if( (double)(mx)!=(double)(0) ) { v = 1/mx; for(i=0; i<=m-1; i++) { for(i_=0; i_<=n-1;i_++) { a[i,i_] = v*a[i,i_]; } } } pivots = new int[Math.Min(m, n)]; tmp = new double[2*Math.Max(m, n)]; rmatrixluprec(ref a, 0, m, n, ref pivots, ref tmp); if( (double)(mx)!=(double)(0) ) { v = mx; for(i=0; i<=m-1; i++) { for(i_=0; i_<=Math.Min(i, n-1);i_++) { a[i,i_] = v*a[i,i_]; } } } } public static void cmatrixlup(ref AP.Complex[,] a, int m, int n, ref int[] pivots) { AP.Complex[] tmp = new AP.Complex[0]; int i = 0; int j = 0; double mx = 0; double v = 0; int i_ = 0; // // Internal LU decomposition subroutine. // Never call it directly. // System.Diagnostics.Debug.Assert(m>0, "CMatrixLUP: incorrect M!"); System.Diagnostics.Debug.Assert(n>0, "CMatrixLUP: incorrect N!"); // // Scale matrix to avoid overflows, // decompose it, then scale back. // mx = 0; for(i=0; i<=m-1; i++) { for(j=0; j<=n-1; j++) { mx = Math.Max(mx, AP.Math.AbsComplex(a[i,j])); } } if( (double)(mx)!=(double)(0) ) { v = 1/mx; for(i=0; i<=m-1; i++) { for(i_=0; i_<=n-1;i_++) { a[i,i_] = v*a[i,i_]; } } } pivots = new int[Math.Min(m, n)]; tmp = new AP.Complex[2*Math.Max(m, n)]; cmatrixluprec(ref a, 0, m, n, ref pivots, ref tmp); if( (double)(mx)!=(double)(0) ) { v = mx; for(i=0; i<=m-1; i++) { for(i_=0; i_<=Math.Min(i, n-1);i_++) { a[i,i_] = v*a[i,i_]; } } } } public static void rmatrixplu(ref double[,] a, int m, int n, ref int[] pivots) { double[] tmp = new double[0]; int i = 0; int j = 0; double mx = 0; double v = 0; int i_ = 0; // // Internal LU decomposition subroutine. // Never call it directly. // System.Diagnostics.Debug.Assert(m>0, "RMatrixPLU: incorrect M!"); System.Diagnostics.Debug.Assert(n>0, "RMatrixPLU: incorrect N!"); tmp = new double[2*Math.Max(m, n)]; pivots = new int[Math.Min(m, n)]; // // Scale matrix to avoid overflows, // decompose it, then scale back. // mx = 0; for(i=0; i<=m-1; i++) { for(j=0; j<=n-1; j++) { mx = Math.Max(mx, Math.Abs(a[i,j])); } } if( (double)(mx)!=(double)(0) ) { v = 1/mx; for(i=0; i<=m-1; i++) { for(i_=0; i_<=n-1;i_++) { a[i,i_] = v*a[i,i_]; } } } rmatrixplurec(ref a, 0, m, n, ref pivots, ref tmp); if( (double)(mx)!=(double)(0) ) { v = mx; for(i=0; i<=Math.Min(m, n)-1; i++) { for(i_=i; i_<=n-1;i_++) { a[i,i_] = v*a[i,i_]; } } } } public static void cmatrixplu(ref AP.Complex[,] a, int m, int n, ref int[] pivots) { AP.Complex[] tmp = new AP.Complex[0]; int i = 0; int j = 0; double mx = 0; AP.Complex v = 0; int i_ = 0; // // Internal LU decomposition subroutine. // Never call it directly. // System.Diagnostics.Debug.Assert(m>0, "CMatrixPLU: incorrect M!"); System.Diagnostics.Debug.Assert(n>0, "CMatrixPLU: incorrect N!"); tmp = new AP.Complex[2*Math.Max(m, n)]; pivots = new int[Math.Min(m, n)]; // // Scale matrix to avoid overflows, // decompose it, then scale back. // mx = 0; for(i=0; i<=m-1; i++) { for(j=0; j<=n-1; j++) { mx = Math.Max(mx, AP.Math.AbsComplex(a[i,j])); } } if( (double)(mx)!=(double)(0) ) { v = 1/mx; for(i=0; i<=m-1; i++) { for(i_=0; i_<=n-1;i_++) { a[i,i_] = v*a[i,i_]; } } } cmatrixplurec(ref a, 0, m, n, ref pivots, ref tmp); if( (double)(mx)!=(double)(0) ) { v = mx; for(i=0; i<=Math.Min(m, n)-1; i++) { for(i_=i; i_<=n-1;i_++) { a[i,i_] = v*a[i,i_]; } } } } /************************************************************************* Recurrent complex LU subroutine. Never call it directly. -- ALGLIB routine -- 04.01.2010 Bochkanov Sergey *************************************************************************/ private static void cmatrixluprec(ref AP.Complex[,] a, int offs, int m, int n, ref int[] pivots, ref AP.Complex[] tmp) { int i = 0; int m1 = 0; int m2 = 0; int i_ = 0; int i1_ = 0; // // Kernel case // if( Math.Min(m, n)<=ablas.ablascomplexblocksize(ref a) ) { cmatrixlup2(ref a, offs, m, n, ref pivots, ref tmp); return; } // // Preliminary step, make N>=M // // ( A1 ) // A = ( ), where A1 is square // ( A2 ) // // Factorize A1, update A2 // if( m>n ) { cmatrixluprec(ref a, offs, n, n, ref pivots, ref tmp); for(i=0; i<=n-1; i++) { i1_ = (offs+n) - (0); for(i_=0; i_<=m-n-1;i_++) { tmp[i_] = a[i_+i1_,offs+i]; } for(i_=offs+n; i_<=offs+m-1;i_++) { a[i_,offs+i] = a[i_,pivots[offs+i]]; } i1_ = (0) - (offs+n); for(i_=offs+n; i_<=offs+m-1;i_++) { a[i_,pivots[offs+i]] = tmp[i_+i1_]; } } ablas.cmatrixrighttrsm(m-n, n, ref a, offs, offs, true, true, 0, ref a, offs+n, offs); return; } // // Non-kernel case // ablas.ablascomplexsplitlength(ref a, m, ref m1, ref m2); cmatrixluprec(ref a, offs, m1, n, ref pivots, ref tmp); if( m2>0 ) { for(i=0; i<=m1-1; i++) { if( offs+i!=pivots[offs+i] ) { i1_ = (offs+m1) - (0); for(i_=0; i_<=m2-1;i_++) { tmp[i_] = a[i_+i1_,offs+i]; } for(i_=offs+m1; i_<=offs+m-1;i_++) { a[i_,offs+i] = a[i_,pivots[offs+i]]; } i1_ = (0) - (offs+m1); for(i_=offs+m1; i_<=offs+m-1;i_++) { a[i_,pivots[offs+i]] = tmp[i_+i1_]; } } } ablas.cmatrixrighttrsm(m2, m1, ref a, offs, offs, true, true, 0, ref a, offs+m1, offs); ablas.cmatrixgemm(m-m1, n-m1, m1, -1.0, ref a, offs+m1, offs, 0, ref a, offs, offs+m1, 0, +1.0, ref a, offs+m1, offs+m1); cmatrixluprec(ref a, offs+m1, m-m1, n-m1, ref pivots, ref tmp); for(i=0; i<=m2-1; i++) { if( offs+m1+i!=pivots[offs+m1+i] ) { i1_ = (offs) - (0); for(i_=0; i_<=m1-1;i_++) { tmp[i_] = a[i_+i1_,offs+m1+i]; } for(i_=offs; i_<=offs+m1-1;i_++) { a[i_,offs+m1+i] = a[i_,pivots[offs+m1+i]]; } i1_ = (0) - (offs); for(i_=offs; i_<=offs+m1-1;i_++) { a[i_,pivots[offs+m1+i]] = tmp[i_+i1_]; } } } } } /************************************************************************* Recurrent real LU subroutine. Never call it directly. -- ALGLIB routine -- 04.01.2010 Bochkanov Sergey *************************************************************************/ private static void rmatrixluprec(ref double[,] a, int offs, int m, int n, ref int[] pivots, ref double[] tmp) { int i = 0; int m1 = 0; int m2 = 0; int i_ = 0; int i1_ = 0; // // Kernel case // if( Math.Min(m, n)<=ablas.ablasblocksize(ref a) ) { rmatrixlup2(ref a, offs, m, n, ref pivots, ref tmp); return; } // // Preliminary step, make N>=M // // ( A1 ) // A = ( ), where A1 is square // ( A2 ) // // Factorize A1, update A2 // if( m>n ) { rmatrixluprec(ref a, offs, n, n, ref pivots, ref tmp); for(i=0; i<=n-1; i++) { if( offs+i!=pivots[offs+i] ) { i1_ = (offs+n) - (0); for(i_=0; i_<=m-n-1;i_++) { tmp[i_] = a[i_+i1_,offs+i]; } for(i_=offs+n; i_<=offs+m-1;i_++) { a[i_,offs+i] = a[i_,pivots[offs+i]]; } i1_ = (0) - (offs+n); for(i_=offs+n; i_<=offs+m-1;i_++) { a[i_,pivots[offs+i]] = tmp[i_+i1_]; } } } ablas.rmatrixrighttrsm(m-n, n, ref a, offs, offs, true, true, 0, ref a, offs+n, offs); return; } // // Non-kernel case // ablas.ablassplitlength(ref a, m, ref m1, ref m2); rmatrixluprec(ref a, offs, m1, n, ref pivots, ref tmp); if( m2>0 ) { for(i=0; i<=m1-1; i++) { if( offs+i!=pivots[offs+i] ) { i1_ = (offs+m1) - (0); for(i_=0; i_<=m2-1;i_++) { tmp[i_] = a[i_+i1_,offs+i]; } for(i_=offs+m1; i_<=offs+m-1;i_++) { a[i_,offs+i] = a[i_,pivots[offs+i]]; } i1_ = (0) - (offs+m1); for(i_=offs+m1; i_<=offs+m-1;i_++) { a[i_,pivots[offs+i]] = tmp[i_+i1_]; } } } ablas.rmatrixrighttrsm(m2, m1, ref a, offs, offs, true, true, 0, ref a, offs+m1, offs); ablas.rmatrixgemm(m-m1, n-m1, m1, -1.0, ref a, offs+m1, offs, 0, ref a, offs, offs+m1, 0, +1.0, ref a, offs+m1, offs+m1); rmatrixluprec(ref a, offs+m1, m-m1, n-m1, ref pivots, ref tmp); for(i=0; i<=m2-1; i++) { if( offs+m1+i!=pivots[offs+m1+i] ) { i1_ = (offs) - (0); for(i_=0; i_<=m1-1;i_++) { tmp[i_] = a[i_+i1_,offs+m1+i]; } for(i_=offs; i_<=offs+m1-1;i_++) { a[i_,offs+m1+i] = a[i_,pivots[offs+m1+i]]; } i1_ = (0) - (offs); for(i_=offs; i_<=offs+m1-1;i_++) { a[i_,pivots[offs+m1+i]] = tmp[i_+i1_]; } } } } } /************************************************************************* Recurrent complex LU subroutine. Never call it directly. -- ALGLIB routine -- 04.01.2010 Bochkanov Sergey *************************************************************************/ private static void cmatrixplurec(ref AP.Complex[,] a, int offs, int m, int n, ref int[] pivots, ref AP.Complex[] tmp) { int i = 0; int n1 = 0; int n2 = 0; int i_ = 0; int i1_ = 0; // // Kernel case // if( Math.Min(m, n)<=ablas.ablascomplexblocksize(ref a) ) { cmatrixplu2(ref a, offs, m, n, ref pivots, ref tmp); return; } // // Preliminary step, make M>=N. // // A = (A1 A2), where A1 is square // Factorize A1, update A2 // if( n>m ) { cmatrixplurec(ref a, offs, m, m, ref pivots, ref tmp); for(i=0; i<=m-1; i++) { i1_ = (offs+m) - (0); for(i_=0; i_<=n-m-1;i_++) { tmp[i_] = a[offs+i,i_+i1_]; } for(i_=offs+m; i_<=offs+n-1;i_++) { a[offs+i,i_] = a[pivots[offs+i],i_]; } i1_ = (0) - (offs+m); for(i_=offs+m; i_<=offs+n-1;i_++) { a[pivots[offs+i],i_] = tmp[i_+i1_]; } } ablas.cmatrixlefttrsm(m, n-m, ref a, offs, offs, false, true, 0, ref a, offs, offs+m); return; } // // Non-kernel case // ablas.ablascomplexsplitlength(ref a, n, ref n1, ref n2); cmatrixplurec(ref a, offs, m, n1, ref pivots, ref tmp); if( n2>0 ) { for(i=0; i<=n1-1; i++) { if( offs+i!=pivots[offs+i] ) { i1_ = (offs+n1) - (0); for(i_=0; i_<=n2-1;i_++) { tmp[i_] = a[offs+i,i_+i1_]; } for(i_=offs+n1; i_<=offs+n-1;i_++) { a[offs+i,i_] = a[pivots[offs+i],i_]; } i1_ = (0) - (offs+n1); for(i_=offs+n1; i_<=offs+n-1;i_++) { a[pivots[offs+i],i_] = tmp[i_+i1_]; } } } ablas.cmatrixlefttrsm(n1, n2, ref a, offs, offs, false, true, 0, ref a, offs, offs+n1); ablas.cmatrixgemm(m-n1, n-n1, n1, -1.0, ref a, offs+n1, offs, 0, ref a, offs, offs+n1, 0, +1.0, ref a, offs+n1, offs+n1); cmatrixplurec(ref a, offs+n1, m-n1, n-n1, ref pivots, ref tmp); for(i=0; i<=n2-1; i++) { if( offs+n1+i!=pivots[offs+n1+i] ) { i1_ = (offs) - (0); for(i_=0; i_<=n1-1;i_++) { tmp[i_] = a[offs+n1+i,i_+i1_]; } for(i_=offs; i_<=offs+n1-1;i_++) { a[offs+n1+i,i_] = a[pivots[offs+n1+i],i_]; } i1_ = (0) - (offs); for(i_=offs; i_<=offs+n1-1;i_++) { a[pivots[offs+n1+i],i_] = tmp[i_+i1_]; } } } } } /************************************************************************* Recurrent real LU subroutine. Never call it directly. -- ALGLIB routine -- 04.01.2010 Bochkanov Sergey *************************************************************************/ private static void rmatrixplurec(ref double[,] a, int offs, int m, int n, ref int[] pivots, ref double[] tmp) { int i = 0; int n1 = 0; int n2 = 0; int i_ = 0; int i1_ = 0; // // Kernel case // if( Math.Min(m, n)<=ablas.ablasblocksize(ref a) ) { rmatrixplu2(ref a, offs, m, n, ref pivots, ref tmp); return; } // // Preliminary step, make M>=N. // // A = (A1 A2), where A1 is square // Factorize A1, update A2 // if( n>m ) { rmatrixplurec(ref a, offs, m, m, ref pivots, ref tmp); for(i=0; i<=m-1; i++) { i1_ = (offs+m) - (0); for(i_=0; i_<=n-m-1;i_++) { tmp[i_] = a[offs+i,i_+i1_]; } for(i_=offs+m; i_<=offs+n-1;i_++) { a[offs+i,i_] = a[pivots[offs+i],i_]; } i1_ = (0) - (offs+m); for(i_=offs+m; i_<=offs+n-1;i_++) { a[pivots[offs+i],i_] = tmp[i_+i1_]; } } ablas.rmatrixlefttrsm(m, n-m, ref a, offs, offs, false, true, 0, ref a, offs, offs+m); return; } // // Non-kernel case // ablas.ablassplitlength(ref a, n, ref n1, ref n2); rmatrixplurec(ref a, offs, m, n1, ref pivots, ref tmp); if( n2>0 ) { for(i=0; i<=n1-1; i++) { if( offs+i!=pivots[offs+i] ) { i1_ = (offs+n1) - (0); for(i_=0; i_<=n2-1;i_++) { tmp[i_] = a[offs+i,i_+i1_]; } for(i_=offs+n1; i_<=offs+n-1;i_++) { a[offs+i,i_] = a[pivots[offs+i],i_]; } i1_ = (0) - (offs+n1); for(i_=offs+n1; i_<=offs+n-1;i_++) { a[pivots[offs+i],i_] = tmp[i_+i1_]; } } } ablas.rmatrixlefttrsm(n1, n2, ref a, offs, offs, false, true, 0, ref a, offs, offs+n1); ablas.rmatrixgemm(m-n1, n-n1, n1, -1.0, ref a, offs+n1, offs, 0, ref a, offs, offs+n1, 0, +1.0, ref a, offs+n1, offs+n1); rmatrixplurec(ref a, offs+n1, m-n1, n-n1, ref pivots, ref tmp); for(i=0; i<=n2-1; i++) { if( offs+n1+i!=pivots[offs+n1+i] ) { i1_ = (offs) - (0); for(i_=0; i_<=n1-1;i_++) { tmp[i_] = a[offs+n1+i,i_+i1_]; } for(i_=offs; i_<=offs+n1-1;i_++) { a[offs+n1+i,i_] = a[pivots[offs+n1+i],i_]; } i1_ = (0) - (offs); for(i_=offs; i_<=offs+n1-1;i_++) { a[pivots[offs+n1+i],i_] = tmp[i_+i1_]; } } } } } /************************************************************************* Complex LUP kernel -- ALGLIB routine -- 10.01.2010 Bochkanov Sergey *************************************************************************/ private static void cmatrixlup2(ref AP.Complex[,] a, int offs, int m, int n, ref int[] pivots, ref AP.Complex[] tmp) { int i = 0; int j = 0; int jp = 0; AP.Complex s = 0; int i_ = 0; int i1_ = 0; // // Quick return if possible // if( m==0 | n==0 ) { return; } // // main cycle // for(j=0; j<=Math.Min(m-1, n-1); j++) { // // Find pivot, swap columns // jp = j; for(i=j+1; i<=n-1; i++) { if( (double)(AP.Math.AbsComplex(a[offs+j,offs+i]))>(double)(AP.Math.AbsComplex(a[offs+j,offs+jp])) ) { jp = i; } } pivots[offs+j] = offs+jp; if( jp!=j ) { i1_ = (offs) - (0); for(i_=0; i_<=m-1;i_++) { tmp[i_] = a[i_+i1_,offs+j]; } for(i_=offs; i_<=offs+m-1;i_++) { a[i_,offs+j] = a[i_,offs+jp]; } i1_ = (0) - (offs); for(i_=offs; i_<=offs+m-1;i_++) { a[i_,offs+jp] = tmp[i_+i1_]; } } // // LU decomposition of 1x(N-J) matrix // if( a[offs+j,offs+j]!=0 & j+1<=n-1 ) { s = 1/a[offs+j,offs+j]; for(i_=offs+j+1; i_<=offs+n-1;i_++) { a[offs+j,i_] = s*a[offs+j,i_]; } } // // Update trailing (M-J-1)x(N-J-1) matrix // if( j(double)(Math.Abs(a[offs+j,offs+jp])) ) { jp = i; } } pivots[offs+j] = offs+jp; if( jp!=j ) { i1_ = (offs) - (0); for(i_=0; i_<=m-1;i_++) { tmp[i_] = a[i_+i1_,offs+j]; } for(i_=offs; i_<=offs+m-1;i_++) { a[i_,offs+j] = a[i_,offs+jp]; } i1_ = (0) - (offs); for(i_=offs; i_<=offs+m-1;i_++) { a[i_,offs+jp] = tmp[i_+i1_]; } } // // LU decomposition of 1x(N-J) matrix // if( (double)(a[offs+j,offs+j])!=(double)(0) & j+1<=n-1 ) { s = 1/a[offs+j,offs+j]; for(i_=offs+j+1; i_<=offs+n-1;i_++) { a[offs+j,i_] = s*a[offs+j,i_]; } } // // Update trailing (M-J-1)x(N-J-1) matrix // if( j(double)(AP.Math.AbsComplex(a[offs+jp,offs+j])) ) { jp = i; } } pivots[offs+j] = offs+jp; if( a[offs+jp,offs+j]!=0 ) { // //Apply the interchange to rows // if( jp!=j ) { for(i=0; i<=n-1; i++) { s = a[offs+j,offs+i]; a[offs+j,offs+i] = a[offs+jp,offs+i]; a[offs+jp,offs+i] = s; } } // //Compute elements J+1:M of J-th column. // if( j+1<=m-1 ) { s = 1/a[offs+j,offs+j]; for(i_=offs+j+1; i_<=offs+m-1;i_++) { a[i_,offs+j] = s*a[i_,offs+j]; } } } if( j(double)(Math.Abs(a[offs+jp,offs+j])) ) { jp = i; } } pivots[offs+j] = offs+jp; if( (double)(a[offs+jp,offs+j])!=(double)(0) ) { // //Apply the interchange to rows // if( jp!=j ) { for(i=0; i<=n-1; i++) { s = a[offs+j,offs+i]; a[offs+j,offs+i] = a[offs+jp,offs+i]; a[offs+jp,offs+i] = s; } } // //Compute elements J+1:M of J-th column. // if( j+1<=m-1 ) { s = 1/a[offs+j,offs+j]; for(i_=offs+j+1; i_<=offs+m-1;i_++) { a[i_,offs+j] = s*a[i_,offs+j]; } } } if( j(double)(0) ) { a[offs,offs] = Math.Sqrt(a[offs,offs].x); result = true; } else { result = false; } return result; } if( n<=ablas.ablascomplexblocksize(ref a) ) { result = hpdmatrixcholesky2(ref a, offs, n, isupper, ref tmp); return result; } // // general case: split task in cache-oblivious manner // result = true; ablas.ablascomplexsplitlength(ref a, n, ref n1, ref n2); result = hpdmatrixcholeskyrec(ref a, offs, n1, isupper, ref tmp); if( !result ) { return result; } if( n2>0 ) { if( isupper ) { ablas.cmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, false, 2, ref a, offs, offs+n1); ablas.cmatrixsyrk(n2, n1, -1.0, ref a, offs, offs+n1, 2, +1.0, ref a, offs+n1, offs+n1, isupper); } else { ablas.cmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, false, 2, ref a, offs+n1, offs); ablas.cmatrixsyrk(n2, n1, -1.0, ref a, offs+n1, offs, 0, +1.0, ref a, offs+n1, offs+n1, isupper); } result = hpdmatrixcholeskyrec(ref a, offs+n1, n2, isupper, ref tmp); if( !result ) { return result; } } return result; } /************************************************************************* Recursive computational subroutine for SPDMatrixCholesky -- ALGLIB routine -- 15.12.2009 Bochkanov Sergey *************************************************************************/ private static bool spdmatrixcholeskyrec(ref double[,] a, int offs, int n, bool isupper, ref double[] tmp) { bool result = new bool(); int n1 = 0; int n2 = 0; // // check N // if( n<1 ) { result = false; return result; } // // special cases // if( n==1 ) { if( (double)(a[offs,offs])>(double)(0) ) { a[offs,offs] = Math.Sqrt(a[offs,offs]); result = true; } else { result = false; } return result; } if( n<=ablas.ablasblocksize(ref a) ) { result = spdmatrixcholesky2(ref a, offs, n, isupper, ref tmp); return result; } // // general case: split task in cache-oblivious manner // result = true; ablas.ablassplitlength(ref a, n, ref n1, ref n2); result = spdmatrixcholeskyrec(ref a, offs, n1, isupper, ref tmp); if( !result ) { return result; } if( n2>0 ) { if( isupper ) { ablas.rmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, false, 1, ref a, offs, offs+n1); ablas.rmatrixsyrk(n2, n1, -1.0, ref a, offs, offs+n1, 1, +1.0, ref a, offs+n1, offs+n1, isupper); } else { ablas.rmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, false, 1, ref a, offs+n1, offs); ablas.rmatrixsyrk(n2, n1, -1.0, ref a, offs+n1, offs, 0, +1.0, ref a, offs+n1, offs+n1, isupper); } result = spdmatrixcholeskyrec(ref a, offs+n1, n2, isupper, ref tmp); if( !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 *************************************************************************/ private static bool hpdmatrixcholesky2(ref AP.Complex[,] aaa, int offs, int n, bool isupper, ref AP.Complex[] tmp) { bool result = new bool(); int i = 0; int j = 0; int k = 0; int j1 = 0; int j2 = 0; double ajj = 0; AP.Complex v = 0; double r = 0; int i_ = 0; int i1_ = 0; result = true; if( n<0 ) { result = false; return result; } // // Quick return if possible // if( n==0 ) { return result; } if( isupper ) { // // Compute the Cholesky factorization A = U'*U. // for(j=0; j<=n-1; j++) { // // Compute U(J,J) and test for non-positive-definiteness. // v = 0.0; for(i_=offs; i_<=offs+j-1;i_++) { v += AP.Math.Conj(aaa[i_,offs+j])*aaa[i_,offs+j]; } ajj = (aaa[offs+j,offs+j]-v).x; if( (double)(ajj)<=(double)(0) ) { aaa[offs+j,offs+j] = ajj; result = false; return result; } ajj = Math.Sqrt(ajj); aaa[offs+j,offs+j] = ajj; // // Compute elements J+1:N-1 of row J. // if( j0 ) { i1_ = (offs) - (0); for(i_=0; i_<=j-1;i_++) { tmp[i_] = -AP.Math.Conj(aaa[i_+i1_,offs+j]); } ablas.cmatrixmv(n-j-1, j, ref aaa, offs, offs+j+1, 1, ref tmp, 0, ref tmp, n); i1_ = (n) - (offs+j+1); for(i_=offs+j+1; i_<=offs+n-1;i_++) { aaa[offs+j,i_] = aaa[offs+j,i_] + tmp[i_+i1_]; } } r = 1/ajj; for(i_=offs+j+1; i_<=offs+n-1;i_++) { aaa[offs+j,i_] = r*aaa[offs+j,i_]; } } } } else { // // Compute the Cholesky factorization A = L*L'. // for(j=0; j<=n-1; j++) { // // Compute L(J+1,J+1) and test for non-positive-definiteness. // v = 0.0; for(i_=offs; i_<=offs+j-1;i_++) { v += AP.Math.Conj(aaa[offs+j,i_])*aaa[offs+j,i_]; } ajj = (aaa[offs+j,offs+j]-v).x; if( (double)(ajj)<=(double)(0) ) { aaa[offs+j,offs+j] = ajj; result = false; return result; } ajj = Math.Sqrt(ajj); aaa[offs+j,offs+j] = ajj; // // Compute elements J+1:N of column J. // if( j0 ) { i1_ = (offs) - (0); for(i_=0; i_<=j-1;i_++) { tmp[i_] = AP.Math.Conj(aaa[offs+j,i_+i1_]); } ablas.cmatrixmv(n-j-1, j, ref aaa, offs+j+1, offs, 0, ref tmp, 0, ref tmp, n); for(i=0; i<=n-j-2; i++) { aaa[offs+j+1+i,offs+j] = (aaa[offs+j+1+i,offs+j]-tmp[n+i])/ajj; } } else { for(i=0; i<=n-j-2; i++) { aaa[offs+j+1+i,offs+j] = aaa[offs+j+1+i,offs+j]/ajj; } } } } } return result; } /************************************************************************* Level-2 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 *************************************************************************/ private static bool spdmatrixcholesky2(ref double[,] aaa, int offs, int n, bool isupper, ref double[] tmp) { bool result = new bool(); int i = 0; int j = 0; int k = 0; int j1 = 0; int j2 = 0; double ajj = 0; double v = 0; double r = 0; int i_ = 0; int i1_ = 0; result = true; if( n<0 ) { result = false; return result; } // // Quick return if possible // if( n==0 ) { return result; } if( isupper ) { // // Compute the Cholesky factorization A = U'*U. // for(j=0; j<=n-1; j++) { // // Compute U(J,J) and test for non-positive-definiteness. // v = 0.0; for(i_=offs; i_<=offs+j-1;i_++) { v += aaa[i_,offs+j]*aaa[i_,offs+j]; } ajj = aaa[offs+j,offs+j]-v; if( (double)(ajj)<=(double)(0) ) { aaa[offs+j,offs+j] = ajj; result = false; return result; } ajj = Math.Sqrt(ajj); aaa[offs+j,offs+j] = ajj; // // Compute elements J+1:N-1 of row J. // if( j0 ) { i1_ = (offs) - (0); for(i_=0; i_<=j-1;i_++) { tmp[i_] = -aaa[i_+i1_,offs+j]; } ablas.rmatrixmv(n-j-1, j, ref aaa, offs, offs+j+1, 1, ref tmp, 0, ref tmp, n); i1_ = (n) - (offs+j+1); for(i_=offs+j+1; i_<=offs+n-1;i_++) { aaa[offs+j,i_] = aaa[offs+j,i_] + tmp[i_+i1_]; } } r = 1/ajj; for(i_=offs+j+1; i_<=offs+n-1;i_++) { aaa[offs+j,i_] = r*aaa[offs+j,i_]; } } } } else { // // Compute the Cholesky factorization A = L*L'. // for(j=0; j<=n-1; j++) { // // Compute L(J+1,J+1) and test for non-positive-definiteness. // v = 0.0; for(i_=offs; i_<=offs+j-1;i_++) { v += aaa[offs+j,i_]*aaa[offs+j,i_]; } ajj = aaa[offs+j,offs+j]-v; if( (double)(ajj)<=(double)(0) ) { aaa[offs+j,offs+j] = ajj; result = false; return result; } ajj = Math.Sqrt(ajj); aaa[offs+j,offs+j] = ajj; // // Compute elements J+1:N of column J. // if( j0 ) { i1_ = (offs) - (0); for(i_=0; i_<=j-1;i_++) { tmp[i_] = aaa[offs+j,i_+i1_]; } ablas.rmatrixmv(n-j-1, j, ref aaa, offs+j+1, offs, 0, ref tmp, 0, ref tmp, n); for(i=0; i<=n-j-2; i++) { aaa[offs+j+1+i,offs+j] = (aaa[offs+j+1+i,offs+j]-tmp[n+i])/ajj; } } else { for(i=0; i<=n-j-2; i++) { aaa[offs+j+1+i,offs+j] = aaa[offs+j+1+i,offs+j]/ajj; } } } } } return result; } } }