source: trunk/sources/ALGLIB/ctrlinsolve.cs @ 2636

/*************************************************************************
This file is a part of 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 ctrlinsolve
    {
        /*************************************************************************
        Utility subroutine performing the "safe" solution of a  system  of  linear
        equations with triangular complex coefficient matrices.

        The feature of an algorithm is that it could not cause an  overflow  or  a
        division by zero regardless of the matrix used as the input. If an overflow
        is possible, an error code is returned.

        The algorithm can solve systems of equations with upper/lower triangular
        matrices,  with/without unit diagonal, and systems of types A*x=b, A^T*x=b,
        A^H*x=b.

        Input parameters:
            A       -   system matrix.
                        Array whose indexes range within [1..N, 1..N].
            N       -   size of matrix A.
            X       -   right-hand member of a system.
                        Array whose index ranges within [1..N].
            IsUpper -   matrix type. If it is True, the system matrix is the upper
                        triangular matrix and is located in the corresponding part
                        of matrix A.
            Trans   -   problem type.
                        If Trans is:
                            * 0, A*x=b
                            * 1, A^T*x=b
                            * 2, A^H*x=b
            Isunit  -   matrix type. If it is True, the system matrix has  a  unit
                        diagonal (the elements on the main diagonal are  not  used
                        in the calculation process), otherwise the matrix is
                        considered to be a general triangular matrix.
            CNORM   -   array which is stored in norms of rows and columns of  the
                        matrix. If the array hasn't been filled up during previous
                        executions  of  an  algorithm  with the same matrix as the
                        input,  it  will  be  filled  up by the subroutine. If the
                        array is filled up, the subroutine uses it without filling
                        it up again.
            NORMIN  -   flag defining the state of column norms array. If True, the
                        array is filled up.
            WORKA   -   working array whose index ranges within [1..N].
            WORKX   -   working array whose index ranges within [1..N].

        Output parameters (if the result is True):
            X       -   solution. Array whose index ranges within [1..N].
            CNORM   -   array of column norms whose index ranges within [1..N].

        Result:
            True, if the matrix is not singular  and  the  algorithm  finished its
                work correctly without causing an overflow.
            False, if  the  matrix  is  singular  or  the  algorithm was cancelled
                because of an overflow possibility.

        Note:
            The disadvantage of an algorithm is that  sometimes  it  overestimates
            an overflow probability. This is not a problem when  solving  ordinary
            systems. If the elements of the matrix used as the input are close  to
            MaxRealNumber, a false overflow detection is possible, but in practice
            such matrices can rarely be found.
            You can find more reliable subroutines in the LAPACK library
            (xLATRS subroutine ).

          -- ALGLIB --
             Copyright 31.03.2006 by Bochkanov Sergey
        *************************************************************************/
        public static bool complexsafesolvetriangular(ref AP.Complex[,] a,
            int n,
            ref AP.Complex[] x,
            bool isupper,
            int trans,
            bool isunit,
            ref AP.Complex[] worka,
            ref AP.Complex[] workx)
        {
            bool result = new bool();
            int i = 0;
            int l = 0;
            int j = 0;
            bool dolswp = new bool();
            double ma = 0;
            double mx = 0;
            double v = 0;
            AP.Complex t = 0;
            AP.Complex r = 0;
            int i_ = 0;
            int i1_ = 0;

            System.Diagnostics.Debug.Assert(trans>=0 & trans<=2, "ComplexSafeSolveTriangular: incorrect parameters!");
            result = true;
            
            //
            // Quick return if possible
            //
            if( n<=0 )
            {
                return result;
            }
            
            //
            // Main cycle
            //
            for(l=1; l<=n; l++)
            {
                
                //
                // Prepare subtask L
                //
                dolswp = false;
                if( trans==0 )
                {
                    if( isupper )
                    {
                        i = n+1-l;
                        i1_ = (i) - (1);
                        for(i_=1; i_<=l;i_++)
                        {
                            worka[i_] = a[i,i_+i1_];
                        }
                        i1_ = (i) - (1);
                        for(i_=1; i_<=l;i_++)
                        {
                            workx[i_] = x[i_+i1_];
                        }
                        dolswp = true;
                    }
                    if( !isupper )
                    {
                        i = l;
                        for(i_=1; i_<=l;i_++)
                        {
                            worka[i_] = a[i,i_];
                        }
                        for(i_=1; i_<=l;i_++)
                        {
                            workx[i_] = x[i_];
                        }
                    }
                }
                if( trans==1 )
                {
                    if( isupper )
                    {
                        i = l;
                        for(i_=1; i_<=l;i_++)
                        {
                            worka[i_] = a[i_,i];
                        }
                        for(i_=1; i_<=l;i_++)
                        {
                            workx[i_] = x[i_];
                        }
                    }
                    if( !isupper )
                    {
                        i = n+1-l;
                        i1_ = (i) - (1);
                        for(i_=1; i_<=l;i_++)
                        {
                            worka[i_] = a[i_+i1_,i];
                        }
                        i1_ = (i) - (1);
                        for(i_=1; i_<=l;i_++)
                        {
                            workx[i_] = x[i_+i1_];
                        }
                        dolswp = true;
                    }
                }
                if( trans==2 )
                {
                    if( isupper )
                    {
                        i = l;
                        for(i_=1; i_<=l;i_++)
                        {
                            worka[i_] = AP.Math.Conj(a[i_,i]);
                        }
                        for(i_=1; i_<=l;i_++)
                        {
                            workx[i_] = x[i_];
                        }
                    }
                    if( !isupper )
                    {
                        i = n+1-l;
                        i1_ = (i) - (1);
                        for(i_=1; i_<=l;i_++)
                        {
                            worka[i_] = AP.Math.Conj(a[i_+i1_,i]);
                        }
                        i1_ = (i) - (1);
                        for(i_=1; i_<=l;i_++)
                        {
                            workx[i_] = x[i_+i1_];
                        }
                        dolswp = true;
                    }
                }
                if( dolswp )
                {
                    t = workx[l];
                    workx[l] = workx[1];
                    workx[1] = t;
                    t = worka[l];
                    worka[l] = worka[1];
                    worka[1] = t;
                }
                if( isunit )
                {
                    worka[l] = 1;
                }
                
                //
                // Test if workA[L]=0
                //
                if( worka[l]==0 )
                {
                    result = false;
                    return result;
                }
                
                //
                // Now we have:
                //
                //  workA[1:L]*workX[1:L] = b[I]
                //
                // with known workA[1:L] and workX[1:L-1]
                // and unknown workX[L]
                //
                t = 0;
                if( l>=2 )
                {
                    ma = 0;
                    for(j=1; j<=l-1; j++)
                    {
                        ma = Math.Max(ma, AP.Math.AbsComplex(worka[j]));
                    }
                    mx = 0;
                    for(j=1; j<=l-1; j++)
                    {
                        mx = Math.Max(mx, AP.Math.AbsComplex(workx[j]));
                    }
                    if( (double)(Math.Max(ma, mx))>(double)(1) )
                    {
                        v = AP.Math.MaxRealNumber/Math.Max(ma, mx);
                        v = v/(l-1);
                        if( (double)(v)<(double)(Math.Min(ma, mx)) )
                        {
                            result = false;
                            return result;
                        }
                    }
                    t = 0.0;
                    for(i_=1; i_<=l-1;i_++)
                    {
                        t += worka[i_]*workx[i_];
                    }
                }
                
                //
                // Now we have:
                //
                //  workA[L]*workX[L] + T = b[I]
                //
                if( (double)(Math.Max(AP.Math.AbsComplex(t), AP.Math.AbsComplex(x[i])))>=(double)(0.5*AP.Math.MaxRealNumber) )
                {
                    result = false;
                    return result;
                }
                r = x[i]-t;
                
                //
                // Now we have:
                //
                //  workA[L]*workX[L] = R
                //
                if( r!=0 )
                {
                    if( (double)(Math.Log(AP.Math.AbsComplex(r))-Math.Log(AP.Math.AbsComplex(worka[l])))>=(double)(Math.Log(AP.Math.MaxRealNumber)) )
                    {
                        result = false;
                        return result;
                    }
                }
                
                //
                // X[I]
                //
                x[i] = r/worka[l];
            }
            return result;
        }
    }
}