source: trunk/sources/ALGLIB/hessenberg.cs @ 2563

/*************************************************************************
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 hessenberg
    {
        /*************************************************************************
        Reduction of a square matrix to  upper Hessenberg form: Q'*A*Q = H,
        where Q is an orthogonal matrix, H - Hessenberg matrix.

        Input parameters:
            A       -   matrix A with elements [0..N-1, 0..N-1]
            N       -   size of matrix A.

        Output parameters:
            A       -   matrices Q and P in  compact form (see below).
            Tau     -   array of scalar factors which are used to form matrix Q.
                        Array whose index ranges within [0..N-2]

        Matrix H is located on the main diagonal, on the lower secondary  diagonal
        and above the main diagonal of matrix A. The elements which are used to
        form matrix Q are situated in array Tau and below the lower secondary
        diagonal of matrix A as follows:

        Matrix Q is represented as a product of elementary reflections

        Q = H(0)*H(2)*...*H(n-2),

        where each H(i) is given by

        H(i) = 1 - tau * v * (v^T)

        where tau is a scalar stored in Tau[I]; v - is a real vector,
        so that v(0:i) = 0, v(i+1) = 1, v(i+2:n-1) stored in A(i+2:n-1,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 rmatrixhessenberg(ref double[,] a,
            int n,
            ref double[] tau)
        {
            int i = 0;
            double v = 0;
            double[] t = new double[0];
            double[] work = new double[0];
            int i_ = 0;
            int i1_ = 0;

            System.Diagnostics.Debug.Assert(n>=0, "RMatrixHessenberg: incorrect N!");
            
            //
            // Quick return if possible
            //
            if( n<=1 )
            {
                return;
            }
            tau = new double[n-2+1];
            t = new double[n+1];
            work = new double[n-1+1];
            for(i=0; i<=n-2; i++)
            {
                
                //
                // Compute elementary reflector H(i) to annihilate A(i+2:ihi,i)
                //
                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 v);
                i1_ = (1) - (i+1);
                for(i_=i+1; i_<=n-1;i_++)
                {
                    a[i_,i] = t[i_+i1_];
                }
                tau[i] = v;
                t[1] = 1;
                
                //
                // Apply H(i) to A(1:ihi,i+1:ihi) from the right
                //
                reflections.applyreflectionfromtheright(ref a, v, ref t, 0, n-1, i+1, n-1, ref work);
                
                //
                // Apply H(i) to A(i+1:ihi,i+1:n) from the left
                //
                reflections.applyreflectionfromtheleft(ref a, v, ref t, i+1, n-1, i+1, n-1, ref work);
            }
        }


        /*************************************************************************
        Unpacking matrix Q which reduces matrix A to upper Hessenberg form

        Input parameters:
            A   -   output of RMatrixHessenberg subroutine.
            N   -   size of matrix A.
            Tau -   scalar factors which are used to form Q.
                    Output of RMatrixHessenberg subroutine.

        Output parameters:
            Q   -   matrix Q.
                    Array whose indexes range within [0..N-1, 0..N-1].

          -- ALGLIB --
             Copyright 2005 by Bochkanov Sergey
        *************************************************************************/
        public static void rmatrixhessenbergunpackq(ref double[,] a,
            int n,
            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+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
            //
            for(i=0; i<=n-2; 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.applyreflectionfromtheright(ref q, tau[i], ref v, 0, n-1, i+1, n-1, ref work);
            }
        }


        /*************************************************************************
        Unpacking matrix H (the result of matrix A reduction to upper Hessenberg form)

        Input parameters:
            A   -   output of RMatrixHessenberg subroutine.
            N   -   size of matrix A.

        Output parameters:
            H   -   matrix H. Array whose indexes range within [0..N-1, 0..N-1].

          -- ALGLIB --
             Copyright 2005 by Bochkanov Sergey
        *************************************************************************/
        public static void rmatrixhessenbergunpackh(ref double[,] a,
            int n,
            ref double[,] h)
        {
            int i = 0;
            int j = 0;
            double[] v = new double[0];
            double[] work = new double[0];
            int i_ = 0;

            if( n==0 )
            {
                return;
            }
            h = new double[n-1+1, n-1+1];
            for(i=0; i<=n-1; i++)
            {
                for(j=0; j<=i-2; j++)
                {
                    h[i,j] = 0;
                }
                j = Math.Max(0, i-1);
                for(i_=j; i_<=n-1;i_++)
                {
                    h[i,i_] = a[i,i_];
                }
            }
        }


        public static void toupperhessenberg(ref double[,] a,
            int n,
            ref double[] tau)
        {
            int i = 0;
            int ip1 = 0;
            int nmi = 0;
            double v = 0;
            double[] t = new double[0];
            double[] work = new double[0];
            int i_ = 0;
            int i1_ = 0;

            System.Diagnostics.Debug.Assert(n>=0, "ToUpperHessenberg: incorrect N!");
            
            //
            // Quick return if possible
            //
            if( n<=1 )
            {
                return;
            }
            tau = new double[n-1+1];
            t = new double[n+1];
            work = new double[n+1];
            for(i=1; i<=n-1; i++)
            {
                
                //
                // Compute elementary reflector H(i) to annihilate A(i+2:ihi,i)
                //
                ip1 = i+1;
                nmi = n-i;
                i1_ = (ip1) - (1);
                for(i_=1; i_<=nmi;i_++)
                {
                    t[i_] = a[i_+i1_,i];
                }
                reflections.generatereflection(ref t, nmi, ref v);
                i1_ = (1) - (ip1);
                for(i_=ip1; i_<=n;i_++)
                {
                    a[i_,i] = t[i_+i1_];
                }
                tau[i] = v;
                t[1] = 1;
                
                //
                // Apply H(i) to A(1:ihi,i+1:ihi) from the right
                //
                reflections.applyreflectionfromtheright(ref a, v, ref t, 1, n, i+1, n, ref work);
                
                //
                // Apply H(i) to A(i+1:ihi,i+1:n) from the left
                //
                reflections.applyreflectionfromtheleft(ref a, v, ref t, i+1, n, i+1, n, ref work);
            }
        }


        public static void unpackqfromupperhessenberg(ref double[,] a,
            int n,
            ref double[] tau,
            ref double[,] q)
        {
            int i = 0;
            int j = 0;
            double[] v = new double[0];
            double[] work = new double[0];
            int ip1 = 0;
            int nmi = 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
            //
            for(i=1; i<=n-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.applyreflectionfromtheright(ref q, tau[i], ref v, 1, n, i+1, n, ref work);
            }
        }


        public static void unpackhfromupperhessenberg(ref double[,] a,
            int n,
            ref double[] tau,
            ref double[,] h)
        {
            int i = 0;
            int j = 0;
            double[] v = new double[0];
            double[] work = new double[0];
            int i_ = 0;

            if( n==0 )
            {
                return;
            }
            h = new double[n+1, n+1];
            for(i=1; i<=n; i++)
            {
                for(j=1; j<=i-2; j++)
                {
                    h[i,j] = 0;
                }
                j = Math.Max(1, i-1);
                for(i_=j; i_<=n;i_++)
                {
                    h[i,i_] = a[i,i_];
                }
            }
        }
    }
}