source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/spdgevd.cs @ 3021

/*************************************************************************
Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).

>>> 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 spdgevd
    {
        /*************************************************************************
        Algorithm for solving the following generalized symmetric positive-definite
        eigenproblem:
            A*x = lambda*B*x (1) or
            A*B*x = lambda*x (2) or
            B*A*x = lambda*x (3).
        where A is a symmetric matrix, B - symmetric positive-definite matrix.
        The problem is solved by reducing it to an ordinary  symmetric  eigenvalue
        problem.

        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 matrices A and B.
            IsUpperA    -   storage format of matrix A.
            B           -   symmetric positive-definite matrix which is given by
                            its upper or lower triangular part.
                            Array whose indexes range within [0..N-1, 0..N-1].
            IsUpperB    -   storage format of matrix B.
            ZNeeded     -   if ZNeeded is equal to:
                             * 0, the eigenvectors are not returned;
                             * 1, the eigenvectors are returned.
            ProblemType -   if ProblemType is equal to:
                             * 1, the following problem is solved: A*x = lambda*B*x;
                             * 2, the following problem is solved: A*B*x = lambda*x;
                             * 3, the following problem is solved: B*A*x = lambda*x.

        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 eigenvectors.
                            Array whose indexes range within [0..N-1, 0..N-1].
                            The eigenvectors are stored in matrix columns. It should
                            be noted that the eigenvectors in such problems do not
                            form an orthogonal system.

        Result:
            True, if the problem was solved successfully.
            False, if the error occurred during the Cholesky decomposition of matrix
            B (the matrix isn’t positive-definite) or during the work of the iterative
            algorithm for solving the symmetric eigenproblem.

        See also the GeneralizedSymmetricDefiniteEVDReduce subroutine.

          -- ALGLIB --
             Copyright 1.28.2006 by Bochkanov Sergey
        *************************************************************************/
        public static bool smatrixgevd(double[,] a,
            int n,
            bool isuppera,
            ref double[,] b,
            bool isupperb,
            int zneeded,
            int problemtype,
            ref double[] d,
            ref double[,] z)
        {
            bool result = new bool();
            double[,] r = new double[0,0];
            double[,] t = new double[0,0];
            bool isupperr = new bool();
            int j1 = 0;
            int j2 = 0;
            int j1inc = 0;
            int j2inc = 0;
            int i = 0;
            int j = 0;
            double v = 0;
            int i_ = 0;

            a = (double[,])a.Clone();

            
            //
            // Reduce and solve
            //
            result = smatrixgevdreduce(ref a, n, isuppera, ref b, isupperb, problemtype, ref r, ref isupperr);
            if( !result )
            {
                return result;
            }
            result = sevd.smatrixevd(a, n, zneeded, isuppera, ref d, ref t);
            if( !result )
            {
                return result;
            }
            
            //
            // Transform eigenvectors if needed
            //
            if( zneeded!=0 )
            {
                
                //
                // fill Z with zeros
                //
                z = new double[n-1+1, n-1+1];
                for(j=0; j<=n-1; j++)
                {
                    z[0,j] = 0.0;
                }
                for(i=1; i<=n-1; i++)
                {
                    for(i_=0; i_<=n-1;i_++)
                    {
                        z[i,i_] = z[0,i_];
                    }
                }
                
                //
                // Setup R properties
                //
                if( isupperr )
                {
                    j1 = 0;
                    j2 = n-1;
                    j1inc = +1;
                    j2inc = 0;
                }
                else
                {
                    j1 = 0;
                    j2 = 0;
                    j1inc = 0;
                    j2inc = +1;
                }
                
                //
                // Calculate R*Z
                //
                for(i=0; i<=n-1; i++)
                {
                    for(j=j1; j<=j2; j++)
                    {
                        v = r[i,j];
                        for(i_=0; i_<=n-1;i_++)
                        {
                            z[i,i_] = z[i,i_] + v*t[j,i_];
                        }
                    }
                    j1 = j1+j1inc;
                    j2 = j2+j2inc;
                }
            }
            return result;
        }


        /*************************************************************************
        Algorithm for reduction of the following generalized symmetric positive-
        definite eigenvalue problem:
            A*x = lambda*B*x (1) or
            A*B*x = lambda*x (2) or
            B*A*x = lambda*x (3)
        to the symmetric eigenvalues problem C*y = lambda*y (eigenvalues of this and
        the given problems are the same, and the eigenvectors of the given problem
        could be obtained by multiplying the obtained eigenvectors by the
        transformation matrix x = R*y).

        Here A is a symmetric matrix, B - symmetric positive-definite matrix.

        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 matrices A and B.
            IsUpperA    -   storage format of matrix A.
            B           -   symmetric positive-definite matrix which is given by
                            its upper or lower triangular part.
                            Array whose indexes range within [0..N-1, 0..N-1].
            IsUpperB    -   storage format of matrix B.
            ProblemType -   if ProblemType is equal to:
                             * 1, the following problem is solved: A*x = lambda*B*x;
                             * 2, the following problem is solved: A*B*x = lambda*x;
                             * 3, the following problem is solved: B*A*x = lambda*x.

        Output parameters:
            A           -   symmetric matrix which is given by its upper or lower
                            triangle depending on IsUpperA. Contains matrix C.
                            Array whose indexes range within [0..N-1, 0..N-1].
            R           -   upper triangular or low triangular transformation matrix
                            which is used to obtain the eigenvectors of a given problem
                            as the product of eigenvectors of C (from the right) and
                            matrix R (from the left). If the matrix is upper
                            triangular, the elements below the main diagonal
                            are equal to 0 (and vice versa). Thus, we can perform
                            the multiplication without taking into account the
                            internal structure (which is an easier though less
                            effective way).
                            Array whose indexes range within [0..N-1, 0..N-1].
            IsUpperR    -   type of matrix R (upper or lower triangular).

        Result:
            True, if the problem was reduced successfully.
            False, if the error occurred during the Cholesky decomposition of
                matrix B (the matrix is not positive-definite).

          -- ALGLIB --
             Copyright 1.28.2006 by Bochkanov Sergey
        *************************************************************************/
        public static bool smatrixgevdreduce(ref double[,] a,
            int n,
            bool isuppera,
            ref double[,] b,
            bool isupperb,
            int problemtype,
            ref double[,] r,
            ref bool isupperr)
        {
            bool result = new bool();
            double[,] t = new double[0,0];
            double[] w1 = new double[0];
            double[] w2 = new double[0];
            double[] w3 = new double[0];
            int i = 0;
            int j = 0;
            double v = 0;
            int i_ = 0;
            int i1_ = 0;

            System.Diagnostics.Debug.Assert(n>0, "SMatrixGEVDReduce: N<=0!");
            System.Diagnostics.Debug.Assert(problemtype==1 | problemtype==2 | problemtype==3, "SMatrixGEVDReduce: incorrect ProblemType!");
            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 = new double[n-1+1, n-1+1];
                if( isupperb )
                {
                    for(i=0; i<=n-1; i++)
                    {
                        for(i_=i; i_<=n-1;i_++)
                        {
                            t[i_,i] = b[i,i_];
                        }
                    }
                }
                else
                {
                    for(i=0; i<=n-1; i++)
                    {
                        for(i_=0; i_<=i;i_++)
                        {
                            t[i,i_] = b[i,i_];
                        }
                    }
                }
                if( !trfac.spdmatrixcholesky(ref t, n, false) )
                {
                    result = false;
                    return result;
                }
                
                //
                // Invert L in T
                //
                if( !trinverse.rmatrixtrinverse(ref t, n, false, false) )
                {
                    result = false;
                    return result;
                }
                
                //
                // Build L^(-1) * A * L^(-T) in R
                //
                w1 = new double[n+1];
                w2 = new double[n+1];
                r = new double[n-1+1, n-1+1];
                for(j=1; j<=n; j++)
                {
                    
                    //
                    // Form w2 = A * l'(j) (here l'(j) is j-th column of L^(-T))
                    //
                    i1_ = (0) - (1);
                    for(i_=1; i_<=j;i_++)
                    {
                        w1[i_] = t[j-1,i_+i1_];
                    }
                    sblas.symmetricmatrixvectormultiply(ref a, isuppera, 0, j-1, ref w1, 1.0, ref w2);
                    if( isuppera )
                    {
                        blas.matrixvectormultiply(ref a, 0, j-1, j, n-1, true, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0);
                    }
                    else
                    {
                        blas.matrixvectormultiply(ref a, j, n-1, 0, j-1, false, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0);
                    }
                    
                    //
                    // Form l(i)*w2 (here l(i) is i-th row of L^(-1))
                    //
                    for(i=1; i<=n; i++)
                    {
                        i1_ = (1)-(0);
                        v = 0.0;
                        for(i_=0; i_<=i-1;i_++)
                        {
                            v += t[i-1,i_]*w2[i_+i1_];
                        }
                        r[i-1,j-1] = v;
                    }
                }
                
                //
                // Copy R to A
                //
                for(i=0; i<=n-1; i++)
                {
                    for(i_=0; i_<=n-1;i_++)
                    {
                        a[i,i_] = r[i,i_];
                    }
                }
                
                //
                // Copy L^(-1) from T to R and transpose
                //
                isupperr = true;
                for(i=0; i<=n-1; i++)
                {
                    for(j=0; j<=i-1; j++)
                    {
                        r[i,j] = 0;
                    }
                }
                for(i=0; i<=n-1; i++)
                {
                    for(i_=i; i_<=n-1;i_++)
                    {
                        r[i,i_] = t[i_,i];
                    }
                }
                return result;
            }
            
            //
            // Problem 2:  A*B*x = lambda*x
            // or
            // problem 3:  B*A*x = lambda*x
            //
            // Reducing to:
            //     C*y = lambda*y
            //     C = U * A * U'
            //     B = U'* U
            //
            if( problemtype==2 | problemtype==3 )
            {
                
                //
                // Factorize B in T: B = U'*U
                //
                t = new double[n-1+1, n-1+1];
                if( isupperb )
                {
                    for(i=0; i<=n-1; i++)
                    {
                        for(i_=i; i_<=n-1;i_++)
                        {
                            t[i,i_] = b[i,i_];
                        }
                    }
                }
                else
                {
                    for(i=0; i<=n-1; i++)
                    {
                        for(i_=i; i_<=n-1;i_++)
                        {
                            t[i,i_] = b[i_,i];
                        }
                    }
                }
                if( !trfac.spdmatrixcholesky(ref t, n, true) )
                {
                    result = false;
                    return result;
                }
                
                //
                // Build U * A * U' in R
                //
                w1 = new double[n+1];
                w2 = new double[n+1];
                w3 = new double[n+1];
                r = new double[n-1+1, n-1+1];
                for(j=1; j<=n; j++)
                {
                    
                    //
                    // Form w2 = A * u'(j) (here u'(j) is j-th column of U')
                    //
                    i1_ = (j-1) - (1);
                    for(i_=1; i_<=n-j+1;i_++)
                    {
                        w1[i_] = t[j-1,i_+i1_];
                    }
                    sblas.symmetricmatrixvectormultiply(ref a, isuppera, j-1, n-1, ref w1, 1.0, ref w3);
                    i1_ = (1) - (j);
                    for(i_=j; i_<=n;i_++)
                    {
                        w2[i_] = w3[i_+i1_];
                    }
                    i1_ = (j-1) - (j);
                    for(i_=j; i_<=n;i_++)
                    {
                        w1[i_] = t[j-1,i_+i1_];
                    }
                    if( isuppera )
                    {
                        blas.matrixvectormultiply(ref a, 0, j-2, j-1, n-1, false, ref w1, j, n, 1.0, ref w2, 1, j-1, 0.0);
                    }
                    else
                    {
                        blas.matrixvectormultiply(ref a, j-1, n-1, 0, j-2, true, ref w1, j, n, 1.0, ref w2, 1, j-1, 0.0);
                    }
                    
                    //
                    // Form u(i)*w2 (here u(i) is i-th row of U)
                    //
                    for(i=1; i<=n; i++)
                    {
                        i1_ = (i)-(i-1);
                        v = 0.0;
                        for(i_=i-1; i_<=n-1;i_++)
                        {
                            v += t[i-1,i_]*w2[i_+i1_];
                        }
                        r[i-1,j-1] = v;
                    }
                }
                
                //
                // Copy R to A
                //
                for(i=0; i<=n-1; i++)
                {
                    for(i_=0; i_<=n-1;i_++)
                    {
                        a[i,i_] = r[i,i_];
                    }
                }
                if( problemtype==2 )
                {
                    
                    //
                    // Invert U in T
                    //
                    if( !trinverse.rmatrixtrinverse(ref t, n, true, false) )
                    {
                        result = false;
                        return result;
                    }
                    
                    //
                    // Copy U^-1 from T to R
                    //
                    isupperr = true;
                    for(i=0; i<=n-1; i++)
                    {
                        for(j=0; j<=i-1; j++)
                        {
                            r[i,j] = 0;
                        }
                    }
                    for(i=0; i<=n-1; i++)
                    {
                        for(i_=i; i_<=n-1;i_++)
                        {
                            r[i,i_] = t[i,i_];
                        }
                    }
                }
                else
                {
                    
                    //
                    // Copy U from T to R and transpose
                    //
                    isupperr = false;
                    for(i=0; i<=n-1; i++)
                    {
                        for(j=i+1; j<=n-1; j++)
                        {
                            r[i,j] = 0;
                        }
                    }
                    for(i=0; i<=n-1; i++)
                    {
                        for(i_=i; i_<=n-1;i_++)
                        {
                            r[i_,i] = t[i,i_];
                        }
                    }
                }
            }
            return result;
        }
    }
}