/************************************************************************* 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 tridiagonal { /************************************************************************* Reduction of a symmetric matrix which is given by its higher or lower triangular part to a tridiagonal matrix using orthogonal similarity transformation: Q'*A*Q=T. Input parameters: A - matrix to be transformed array with elements [0..N-1, 0..N-1]. N - size of matrix A. IsUpper - storage format. If IsUpper = True, then matrix A is given by its upper triangle, and the lower triangle is not used and not modified by the algorithm, and vice versa if IsUpper = False. Output parameters: A - matrices T and Q in compact form (see lower) Tau - array of factors which are forming matrices H(i) array with elements [0..N-2]. D - main diagonal of symmetric matrix T. array with elements [0..N-1]. E - secondary diagonal of symmetric matrix T. array with elements [0..N-2]. If IsUpper=True, the matrix Q is represented as a product of elementary reflectors Q = H(n-2) . . . H(2) H(0). Each H(i) has the form H(i) = I - tau * v * v' where tau is a real scalar, and v is a real vector with v(i+1:n-1) = 0, v(i) = 1, v(0:i-1) is stored on exit in A(0:i-1,i+1), and tau in TAU(i). If IsUpper=False, the matrix Q is represented as a product of elementary reflectors Q = H(0) H(2) . . . H(n-2). Each H(i) has the form H(i) = I - tau * v * v' where tau is a real scalar, and v is a real vector with v(0:i) = 0, v(i+1) = 1, v(i+2:n-1) is stored on exit in A(i+2:n-1,i), and tau in TAU(i). The contents of A on exit are illustrated by the following examples with n = 5: if UPLO = 'U': if UPLO = 'L': ( d e v1 v2 v3 ) ( d ) ( d e v2 v3 ) ( e d ) ( d e v3 ) ( v0 e d ) ( d e ) ( v0 v1 e d ) ( d ) ( v0 v1 v2 e d ) where d and e denote diagonal and off-diagonal elements of T, and vi denotes an element of the vector defining H(i). -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University October 31, 1992 *************************************************************************/ public static void smatrixtd(ref double[,] a, int n, bool isupper, ref double[] tau, ref double[] d, ref double[] e) { int i = 0; double alpha = 0; double taui = 0; double v = 0; double[] t = new double[0]; double[] t2 = new double[0]; double[] t3 = new double[0]; int i_ = 0; int i1_ = 0; if( n<=0 ) { return; } t = new double[n+1]; t2 = new double[n+1]; t3 = new double[n+1]; if( n>1 ) { tau = new double[n-2+1]; } d = new double[n-1+1]; if( n>1 ) { e = new double[n-2+1]; } 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_ = (0) - (2); for(i_=2; i_<=i+1;i_++) { t[i_] = a[i_+i1_,i+1]; } } t[1] = a[i,i+1]; reflections.generatereflection(ref t, i+1, ref taui); if( i>=1 ) { i1_ = (2) - (0); for(i_=0; i_<=i-1;i_++) { a[i_,i+1] = t[i_+i1_]; } } a[i,i+1] = t[1]; e[i] = a[i,i+1]; if( (double)(taui)!=(double)(0) ) { // // Apply H from both sides to A // a[i,i+1] = 1; // // Compute x := tau * A * v storing x in TAU // i1_ = (0) - (1); for(i_=1; i_<=i+1;i_++) { t[i_] = a[i_+i1_,i+1]; } sblas.symmetricmatrixvectormultiply(ref a, isupper, 0, i, ref t, taui, ref t3); i1_ = (1) - (0); for(i_=0; i_<=i;i_++) { tau[i_] = t3[i_+i1_]; } // // Compute w := x - 1/2 * tau * (x'*v) * v // v = 0.0; for(i_=0; i_<=i;i_++) { v += tau[i_]*a[i_,i+1]; } alpha = -(0.5*taui*v); for(i_=0; i_<=i;i_++) { tau[i_] = tau[i_] + alpha*a[i_,i+1]; } // // Apply the transformation as a rank-2 update: // A := A - v * w' - w * v' // i1_ = (0) - (1); for(i_=1; i_<=i+1;i_++) { t[i_] = a[i_+i1_,i+1]; } i1_ = (0) - (1); for(i_=1; i_<=i+1;i_++) { t3[i_] = tau[i_+i1_]; } sblas.symmetricrank2update(ref a, isupper, 0, i, ref t, ref t3, ref t2, -1); a[i,i+1] = e[i]; } d[i+1] = a[i+1,i+1]; tau[i] = taui; } d[0] = a[0,0]; } else { // // Reduce the lower triangle of A // for(i=0; i<=n-2; i++) { // // Generate elementary reflector H = E - tau * v * v' // i1_ = (i+1) - (1); for(i_=1; i_<=n-i-1;i_++) { t[i_] = a[i_+i1_,i]; } reflections.generatereflection(ref t, n-i-1, ref taui); i1_ = (1) - (i+1); for(i_=i+1; i_<=n-1;i_++) { a[i_,i] = t[i_+i1_]; } e[i] = a[i+1,i]; if( (double)(taui)!=(double)(0) ) { // // Apply H from both sides to A // a[i+1,i] = 1; // // Compute x := tau * A * v storing y in TAU // i1_ = (i+1) - (1); for(i_=1; i_<=n-i-1;i_++) { t[i_] = a[i_+i1_,i]; } sblas.symmetricmatrixvectormultiply(ref a, isupper, i+1, n-1, ref t, taui, ref t2); i1_ = (1) - (i); for(i_=i; i_<=n-2;i_++) { tau[i_] = t2[i_+i1_]; } // // Compute w := x - 1/2 * tau * (x'*v) * v // i1_ = (i+1)-(i); v = 0.0; for(i_=i; i_<=n-2;i_++) { v += tau[i_]*a[i_+i1_,i]; } alpha = -(0.5*taui*v); i1_ = (i+1) - (i); for(i_=i; i_<=n-2;i_++) { tau[i_] = tau[i_] + alpha*a[i_+i1_,i]; } // // Apply the transformation as a rank-2 update: // A := A - v * w' - w * v' // // i1_ = (i+1) - (1); for(i_=1; i_<=n-i-1;i_++) { t[i_] = a[i_+i1_,i]; } i1_ = (i) - (1); for(i_=1; i_<=n-i-1;i_++) { t2[i_] = tau[i_+i1_]; } sblas.symmetricrank2update(ref a, isupper, i+1, n-1, ref t, ref t2, ref t3, -1); a[i+1,i] = e[i]; } d[i] = a[i,i]; tau[i] = taui; } d[n-1] = a[n-1,n-1]; } } /************************************************************************* Unpacking matrix Q which reduces symmetric matrix to a tridiagonal form. Input parameters: A - the result of a SMatrixTD subroutine N - size of matrix A. IsUpper - storage format (a parameter of SMatrixTD subroutine) Tau - the result of a SMatrixTD subroutine Output parameters: Q - transformation matrix. array with elements [0..N-1, 0..N-1]. -- ALGLIB -- Copyright 2005-2008 by Bochkanov Sergey *************************************************************************/ public static void smatrixtdunpackq(ref double[,] a, int n, bool isupper, ref double[] tau, ref double[,] q) { int i = 0; int j = 0; double[] v = new double[0]; double[] work = new double[0]; int i_ = 0; int i1_ = 0; if( n==0 ) { return; } // // init // q = new double[n-1+1, n-1+1]; v = new double[n+1]; work = new double[n-1+1]; for(i=0; i<=n-1; i++) { for(j=0; j<=n-1; j++) { if( i==j ) { q[i,j] = 1; } else { q[i,j] = 0; } } } // // unpack Q // if( isupper ) { for(i=0; i<=n-2; i++) { // // Apply H(i) // i1_ = (0) - (1); for(i_=1; i_<=i+1;i_++) { v[i_] = a[i_+i1_,i+1]; } v[i+1] = 1; reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, 0, i, 0, n-1, ref work); } } else { for(i=n-2; i>=0; i--) { // // Apply H(i) // i1_ = (i+1) - (1); for(i_=1; i_<=n-i-1;i_++) { v[i_] = a[i_+i1_,i]; } v[1] = 1; reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i+1, n-1, 0, n-1, ref work); } } } public static void totridiagonal(ref double[,] a, int n, bool isupper, ref double[] tau, ref double[] d, ref double[] e) { int i = 0; int ip1 = 0; int im1 = 0; int nmi = 0; int nm1 = 0; double alpha = 0; double taui = 0; double v = 0; double[] t = new double[0]; double[] t2 = new double[0]; double[] t3 = new double[0]; int i_ = 0; int i1_ = 0; if( n<=0 ) { return; } t = new double[n+1]; t2 = new double[n+1]; t3 = new double[n+1]; tau = new double[Math.Max(1, n-1)+1]; d = new double[n+1]; e = new double[Math.Max(1, n-1)+1]; if( isupper ) { // // Reduce the upper triangle of A // for(i=n-1; i>=1; i--) { // // Generate elementary reflector H(i) = I - tau * v * v' // to annihilate A(1:i-1,i+1) // // DLARFG( I, A( I, I+1 ), A( 1, I+1 ), 1, TAUI ); // ip1 = i+1; im1 = i-1; if( i>=2 ) { i1_ = (1) - (2); for(i_=2; i_<=i;i_++) { t[i_] = a[i_+i1_,ip1]; } } t[1] = a[i,ip1]; reflections.generatereflection(ref t, i, ref taui); if( i>=2 ) { i1_ = (2) - (1); for(i_=1; i_<=im1;i_++) { a[i_,ip1] = t[i_+i1_]; } } a[i,ip1] = t[1]; e[i] = a[i,i+1]; if( (double)(taui)!=(double)(0) ) { // // Apply H(i) from both sides to A(1:i,1:i) // a[i,i+1] = 1; // // Compute x := tau * A * v storing x in TAU(1:i) // // DSYMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO, TAU, 1 ); // ip1 = i+1; for(i_=1; i_<=i;i_++) { t[i_] = a[i_,ip1]; } sblas.symmetricmatrixvectormultiply(ref a, isupper, 1, i, ref t, taui, ref tau); // // Compute w := x - 1/2 * tau * (x'*v) * v // ip1 = i+1; v = 0.0; for(i_=1; i_<=i;i_++) { v += tau[i_]*a[i_,ip1]; } alpha = -(0.5*taui*v); for(i_=1; i_<=i;i_++) { tau[i_] = tau[i_] + alpha*a[i_,ip1]; } // // Apply the transformation as a rank-2 update: // A := A - v * w' - w * v' // // DSYR2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A, LDA ); // for(i_=1; i_<=i;i_++) { t[i_] = a[i_,ip1]; } sblas.symmetricrank2update(ref a, isupper, 1, i, ref t, ref tau, ref t2, -1); a[i,i+1] = e[i]; } d[i+1] = a[i+1,i+1]; tau[i] = taui; } d[1] = a[1,1]; } else { // // Reduce the lower triangle of A // for(i=1; i<=n-1; i++) { // // Generate elementary reflector H(i) = I - tau * v * v' // to annihilate A(i+2:n,i) // //DLARFG( N-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1, TAUI ); // nmi = n-i; ip1 = i+1; i1_ = (ip1) - (1); for(i_=1; i_<=nmi;i_++) { t[i_] = a[i_+i1_,i]; } reflections.generatereflection(ref t, nmi, ref taui); i1_ = (1) - (ip1); for(i_=ip1; i_<=n;i_++) { a[i_,i] = t[i_+i1_]; } e[i] = a[i+1,i]; if( (double)(taui)!=(double)(0) ) { // // Apply H(i) from both sides to A(i+1:n,i+1:n) // a[i+1,i] = 1; // // Compute x := tau * A * v storing y in TAU(i:n-1) // //DSYMV( UPLO, N-I, TAUI, A( I+1, I+1 ), LDA, A( I+1, I ), 1, ZERO, TAU( I ), 1 ); // ip1 = i+1; nmi = n-i; nm1 = n-1; i1_ = (ip1) - (1); for(i_=1; i_<=nmi;i_++) { t[i_] = a[i_+i1_,i]; } sblas.symmetricmatrixvectormultiply(ref a, isupper, i+1, n, ref t, taui, ref t2); i1_ = (1) - (i); for(i_=i; i_<=nm1;i_++) { tau[i_] = t2[i_+i1_]; } // // Compute w := x - 1/2 * tau * (x'*v) * v // nm1 = n-1; ip1 = i+1; i1_ = (ip1)-(i); v = 0.0; for(i_=i; i_<=nm1;i_++) { v += tau[i_]*a[i_+i1_,i]; } alpha = -(0.5*taui*v); i1_ = (ip1) - (i); for(i_=i; i_<=nm1;i_++) { tau[i_] = tau[i_] + alpha*a[i_+i1_,i]; } // // Apply the transformation as a rank-2 update: // A := A - v * w' - w * v' // //DSYR2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1, A( I+1, I+1 ), LDA ); // nm1 = n-1; nmi = n-i; ip1 = i+1; i1_ = (ip1) - (1); for(i_=1; i_<=nmi;i_++) { t[i_] = a[i_+i1_,i]; } i1_ = (i) - (1); for(i_=1; i_<=nmi;i_++) { t2[i_] = tau[i_+i1_]; } sblas.symmetricrank2update(ref a, isupper, i+1, n, ref t, ref t2, ref t3, -1); a[i+1,i] = e[i]; } d[i] = a[i,i]; tau[i] = taui; } d[n] = a[n,n]; } } public static void unpackqfromtridiagonal(ref double[,] a, int n, bool isupper, ref double[] tau, ref double[,] q) { int i = 0; int j = 0; int ip1 = 0; int nmi = 0; double[] v = new double[0]; double[] work = new double[0]; int i_ = 0; int i1_ = 0; if( n==0 ) { return; } // // init // q = new double[n+1, n+1]; v = new double[n+1]; work = new double[n+1]; for(i=1; i<=n; i++) { for(j=1; j<=n; j++) { if( i==j ) { q[i,j] = 1; } else { q[i,j] = 0; } } } // // unpack Q // if( isupper ) { for(i=1; i<=n-1; i++) { // // Apply H(i) // ip1 = i+1; for(i_=1; i_<=i;i_++) { v[i_] = a[i_,ip1]; } v[i] = 1; reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, 1, i, 1, n, ref work); } } else { for(i=n-1; i>=1; i--) { // // Apply H(i) // ip1 = i+1; nmi = n-i; i1_ = (ip1) - (1); for(i_=1; i_<=nmi;i_++) { v[i_] = a[i_+i1_,i]; } v[1] = 1; reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i+1, n, 1, n, ref work); } } } } }