source: trunk/sources/ALGLIB/htridiagonal.cs @ 2575

/*************************************************************************
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 htridiagonal
    {
        /*************************************************************************
        Reduction of a Hermitian matrix which is given  by  its  higher  or  lower
        triangular part to a real  tridiagonal  matrix  using  unitary  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 real symmetric matrix T.
                        array with elements [0..N-1].
            E       -   secondary diagonal of real 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 complex scalar, and v is a complex 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 complex scalar, and v is a complex 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 hmatrixtd(ref AP.Complex[,] a,
            int n,
            bool isupper,
            ref AP.Complex[] tau,
            ref double[] d,
            ref double[] e)
        {
            int i = 0;
            AP.Complex alpha = 0;
            AP.Complex taui = 0;
            AP.Complex v = 0;
            AP.Complex[] t = new AP.Complex[0];
            AP.Complex[] t2 = new AP.Complex[0];
            AP.Complex[] t3 = new AP.Complex[0];
            int i_ = 0;
            int i1_ = 0;

            if( n<=0 )
            {
                return;
            }
            for(i=0; i<=n-1; i++)
            {
                System.Diagnostics.Debug.Assert((double)(a[i,i].y)==(double)(0));
            }
            if( n>1 )
            {
                tau = new AP.Complex[n-2+1];
                e = new double[n-2+1];
            }
            d = new double[n-1+1];
            t = new AP.Complex[n-1+1];
            t2 = new AP.Complex[n-1+1];
            t3 = new AP.Complex[n-1+1];
            if( isupper )
            {
                
                //
                // Reduce the upper triangle of A
                //
                a[n-1,n-1] = a[n-1,n-1].x;
                for(i=n-2; i>=0; i--)
                {
                    
                    //
                    // Generate elementary reflector H = I+1 - tau * v * v'
                    //
                    alpha = a[i,i+1];
                    t[1] = alpha;
                    if( i>=1 )
                    {
                        i1_ = (0) - (2);
                        for(i_=2; i_<=i+1;i_++)
                        {
                            t[i_] = a[i_+i1_,i+1];
                        }
                    }
                    creflections.complexgeneratereflection(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_];
                        }
                    }
                    alpha = t[1];
                    e[i] = alpha.x;
                    if( taui!=0 )
                    {
                        
                        //
                        // Apply H(I+1) 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];
                        }
                        hblas.hermitianmatrixvectormultiply(ref a, isupper, 0, i, ref t, taui, ref t2);
                        i1_ = (1) - (0);
                        for(i_=0; i_<=i;i_++)
                        {
                            tau[i_] = t2[i_+i1_];
                        }
                        
                        //
                        // Compute  w := x - 1/2 * tau * (x'*v) * v
                        //
                        v = 0.0;
                        for(i_=0; i_<=i;i_++)
                        {
                            v += AP.Math.Conj(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_];
                        }
                        hblas.hermitianrank2update(ref a, isupper, 0, i, ref t, ref t3, ref t2, -1);
                    }
                    else
                    {
                        a[i,i] = a[i,i].x;
                    }
                    a[i,i+1] = e[i];
                    d[i+1] = a[i+1,i+1].x;
                    tau[i] = taui;
                }
                d[0] = a[0,0].x;
            }
            else
            {
                
                //
                // Reduce the lower triangle of A
                //
                a[0,0] = a[0,0].x;
                for(i=0; i<=n-2; i++)
                {
                    
                    //
                    // Generate elementary reflector H = I - tau * v * v'
                    //
                    i1_ = (i+1) - (1);
                    for(i_=1; i_<=n-i-1;i_++)
                    {
                        t[i_] = a[i_+i1_,i];
                    }
                    creflections.complexgeneratereflection(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].x;
                    if( taui!=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
                        //
                        i1_ = (i+1) - (1);
                        for(i_=1; i_<=n-i-1;i_++)
                        {
                            t[i_] = a[i_+i1_,i];
                        }
                        hblas.hermitianmatrixvectormultiply(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 += AP.Math.Conj(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_];
                        }
                        hblas.hermitianrank2update(ref a, isupper, i+1, n-1, ref t, ref t2, ref t3, -1);
                    }
                    else
                    {
                        a[i+1,i+1] = a[i+1,i+1].x;
                    }
                    a[i+1,i] = e[i];
                    d[i] = a[i,i].x;
                    tau[i] = taui;
                }
                d[n-1] = a[n-1,n-1].x;
            }
        }


        /*************************************************************************
        Unpacking matrix Q which reduces a Hermitian matrix to a real  tridiagonal
        form.

        Input parameters:
            A       -   the result of a HMatrixTD subroutine
            N       -   size of matrix A.
            IsUpper -   storage format (a parameter of HMatrixTD subroutine)
            Tau     -   the result of a HMatrixTD subroutine

        Output parameters:
            Q       -   transformation matrix.
                        array with elements [0..N-1, 0..N-1].

          -- ALGLIB --
             Copyright 2005, 2007, 2008 by Bochkanov Sergey
        *************************************************************************/
        public static void hmatrixtdunpackq(ref AP.Complex[,] a,
            int n,
            bool isupper,
            ref AP.Complex[] tau,
            ref AP.Complex[,] q)
        {
            int i = 0;
            int j = 0;
            AP.Complex[] v = new AP.Complex[0];
            AP.Complex[] work = new AP.Complex[0];
            int i_ = 0;
            int i1_ = 0;

            if( n==0 )
            {
                return;
            }
            
            //
            // init
            //
            q = new AP.Complex[n-1+1, n-1+1];
            v = new AP.Complex[n+1];
            work = new AP.Complex[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;
                    creflections.complexapplyreflectionfromtheleft(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;
                    creflections.complexapplyreflectionfromtheleft(ref q, tau[i], ref v, i+1, n-1, 0, n-1, ref work);
                }
            }
        }


        /*************************************************************************

          -- 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 hermitiantotridiagonal(ref AP.Complex[,] a,
            int n,
            bool isupper,
            ref AP.Complex[] tau,
            ref double[] d,
            ref double[] e)
        {
            int i = 0;
            AP.Complex alpha = 0;
            AP.Complex taui = 0;
            AP.Complex v = 0;
            AP.Complex[] t = new AP.Complex[0];
            AP.Complex[] t2 = new AP.Complex[0];
            AP.Complex[] t3 = new AP.Complex[0];
            int i_ = 0;
            int i1_ = 0;

            if( n<=0 )
            {
                return;
            }
            for(i=1; i<=n; i++)
            {
                System.Diagnostics.Debug.Assert((double)(a[i,i].y)==(double)(0));
            }
            tau = new AP.Complex[Math.Max(1, n-1)+1];
            d = new double[n+1];
            e = new double[Math.Max(1, n-1)+1];
            t = new AP.Complex[n+1];
            t2 = new AP.Complex[n+1];
            t3 = new AP.Complex[n+1];
            if( isupper )
            {
                
                //
                // Reduce the upper triangle of A
                //
                a[n,n] = a[n,n].x;
                for(i=n-1; i>=1; i--)
                {
                    
                    //
                    // Generate elementary reflector H(i) = I - tau * v * v'
                    // to annihilate A(1:i-1,i+1)
                    //
                    alpha = a[i,i+1];
                    t[1] = alpha;
                    if( i>=2 )
                    {
                        i1_ = (1) - (2);
                        for(i_=2; i_<=i;i_++)
                        {
                            t[i_] = a[i_+i1_,i+1];
                        }
                    }
                    creflections.complexgeneratereflection(ref t, i, ref taui);
                    if( i>=2 )
                    {
                        i1_ = (2) - (1);
                        for(i_=1; i_<=i-1;i_++)
                        {
                            a[i_,i+1] = t[i_+i1_];
                        }
                    }
                    alpha = t[1];
                    e[i] = alpha.x;
                    if( taui!=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)
                        //
                        for(i_=1; i_<=i;i_++)
                        {
                            t[i_] = a[i_,i+1];
                        }
                        hblas.hermitianmatrixvectormultiply(ref a, isupper, 1, i, ref t, taui, ref tau);
                        
                        //
                        // Compute  w := x - 1/2 * tau * (x'*v) * v
                        //
                        v = 0.0;
                        for(i_=1; i_<=i;i_++)
                        {
                            v += AP.Math.Conj(tau[i_])*a[i_,i+1];
                        }
                        alpha = -(0.5*taui*v);
                        for(i_=1; 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'
                        //
                        for(i_=1; i_<=i;i_++)
                        {
                            t[i_] = a[i_,i+1];
                        }
                        hblas.hermitianrank2update(ref a, isupper, 1, i, ref t, ref tau, ref t2, -1);
                    }
                    else
                    {
                        a[i,i] = a[i,i].x;
                    }
                    a[i,i+1] = e[i];
                    d[i+1] = a[i+1,i+1].x;
                    tau[i] = taui;
                }
                d[1] = a[1,1].x;
            }
            else
            {
                
                //
                // Reduce the lower triangle of A
                //
                a[1,1] = a[1,1].x;
                for(i=1; i<=n-1; i++)
                {
                    
                    //
                    // Generate elementary reflector H(i) = I - tau * v * v'
                    // to annihilate A(i+2:n,i)
                    //
                    i1_ = (i+1) - (1);
                    for(i_=1; i_<=n-i;i_++)
                    {
                        t[i_] = a[i_+i1_,i];
                    }
                    creflections.complexgeneratereflection(ref t, n-i, ref taui);
                    i1_ = (1) - (i+1);
                    for(i_=i+1; i_<=n;i_++)
                    {
                        a[i_,i] = t[i_+i1_];
                    }
                    e[i] = a[i+1,i].x;
                    if( taui!=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)
                        //
                        i1_ = (i+1) - (1);
                        for(i_=1; i_<=n-i;i_++)
                        {
                            t[i_] = a[i_+i1_,i];
                        }
                        hblas.hermitianmatrixvectormultiply(ref a, isupper, i+1, n, ref t, taui, ref t2);
                        i1_ = (1) - (i);
                        for(i_=i; i_<=n-1;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-1;i_++)
                        {
                            v += AP.Math.Conj(tau[i_])*a[i_+i1_,i];
                        }
                        alpha = -(0.5*taui*v);
                        i1_ = (i+1) - (i);
                        for(i_=i; i_<=n-1;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;i_++)
                        {
                            t[i_] = a[i_+i1_,i];
                        }
                        i1_ = (i) - (1);
                        for(i_=1; i_<=n-i;i_++)
                        {
                            t2[i_] = tau[i_+i1_];
                        }
                        hblas.hermitianrank2update(ref a, isupper, i+1, n, ref t, ref t2, ref t3, -1);
                    }
                    else
                    {
                        a[i+1,i+1] = a[i+1,i+1].x;
                    }
                    a[i+1,i] = e[i];
                    d[i] = a[i,i].x;
                    tau[i] = taui;
                }
                d[n] = a[n,n].x;
            }
        }


        /*************************************************************************

          -- ALGLIB --
             Copyright 2005, 2007 by Bochkanov Sergey
        *************************************************************************/
        public static void unpackqfromhermitiantridiagonal(ref AP.Complex[,] a,
            int n,
            bool isupper,
            ref AP.Complex[] tau,
            ref AP.Complex[,] q)
        {
            int i = 0;
            int j = 0;
            AP.Complex[] v = new AP.Complex[0];
            AP.Complex[] work = new AP.Complex[0];
            int i_ = 0;
            int i1_ = 0;

            if( n==0 )
            {
                return;
            }
            
            //
            // init
            //
            q = new AP.Complex[n+1, n+1];
            v = new AP.Complex[n+1];
            work = new AP.Complex[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)
                    //
                    for(i_=1; i_<=i;i_++)
                    {
                        v[i_] = a[i_,i+1];
                    }
                    v[i] = 1;
                    creflections.complexapplyreflectionfromtheleft(ref q, tau[i], ref v, 1, i, 1, n, ref work);
                }
            }
            else
            {
                for(i=n-1; i>=1; i--)
                {
                    
                    //
                    // Apply H(i)
                    //
                    i1_ = (i+1) - (1);
                    for(i_=1; i_<=n-i;i_++)
                    {
                        v[i_] = a[i_+i1_,i];
                    }
                    v[1] = 1;
                    creflections.complexapplyreflectionfromtheleft(ref q, tau[i], ref v, i+1, n, 1, n, ref work);
                }
            }
        }
    }
}