Windows Development in C++, COM API Clients

#include "stdafx.h"



/*
 * -- SuperLU MT routine (version 2.0) --
 * Lawrence Berkeley National Lab, Univ. of California Berkeley,
 * and Xerox Palo Alto Research Center.
 * September 10, 2007
 *
 * History:     Modified from lapack routines DGECON.
 */
#include <math.h>
#include "hnum_pcsp_defs.h"



namespace harlinn
{
    namespace numerics
    {
        namespace SuperLU
        {
            namespace Complex
            {

                /*
                    Purpose   
                    =======   

                    CGSCON estimates the reciprocal of the condition number of a general 
                    real matrix A, in either the 1-norm or the infinity-norm, using   
                    the LU factorization computed by CGETRF.   

                    An estimate is obtained for norm(inv(A)), and the reciprocal of the   
                    condition number is computed as   
                       RCOND = 1 / ( norm(A) * norm(inv(A)) ).   

                    See supermatrix.h for the definition of 'SuperMatrix' structure.
 
                    Arguments   
                    =========   

                    NORM    (input) char*
                            Specifies whether the 1-norm condition number or the   
                            infinity-norm condition number is required:   
                            = '1' or 'O':  1-norm;   
                            = 'I':         Infinity-norm.
	    
                    L       (input) SuperMatrix*
                            The factor L from the factorization Pr*A*Pc=L*U as computed by
                            cgstrf(). Use compressed row subscripts storage for supernodes,
                            i.e., L has types: Stype = SLU_SCP, Dtype = SLU_C, Mtype = SLU_TRLU.
 
                    U       (input) SuperMatrix*
                            The factor U from the factorization Pr*A*Pc=L*U as computed by
                            cgstrf(). Use column-wise storage scheme, i.e., U has types:
                            Stype = SLU_NCP, Dtype = SLU_C, Mtype = SLU_TRU.
	    
                    ANORM   (input) float
                            If NORM = '1' or 'O', the 1-norm of the original matrix A.   
                            If NORM = 'I', the infinity-norm of the original matrix A.
	    
                    RCOND   (output) float*
                            The reciprocal of the condition number of the matrix A,   
                            computed as RCOND = 1/(norm(A) * norm(inv(A))).
	    
                    INFO    (output) int*
                            = 0:  successful exit   
                            < 0:  if INFO = -i, the i-th argument had an illegal value   

                    ===================================================================== 
                */
                void cgscon(char *norm, SuperMatrix *L, SuperMatrix *U, float anorm, float *rcond, int *info)
                {
                    // Local variables 
                    int    kase, kase1, onenrm, i;
                    float ainvnm;
                    complex *work;
    
                    // Test the input parameters. 
                    *info = 0;
                    onenrm = *(unsigned char *)norm == '1' || lsame_(norm, "O");
                    if (! onenrm && ! lsame_(norm, "I")) 
                    {
                        *info = -1;
                    }
                    else if (L->nrow < 0 || L->nrow != L->ncol ||
                             L->Stype != SLU_SCP || L->Dtype != SLU_C || L->Mtype != SLU_TRLU)
                    {
	                 *info = -2;
                    }
                    else if (U->nrow < 0 || U->nrow != U->ncol || U->Stype != SLU_NCP || U->Dtype != SLU_C || U->Mtype != SLU_TRU) 
                    {
	                    *info = -3;
                    }
                    if (*info != 0) 
                    {
	                    i = -(*info);
	                    xerbla_("cgscon", &i);
	                    return;
                    }

                    // Quick return if possible 
                    *rcond = 0.;
                    if ( L->nrow == 0 || U->nrow == 0) 
                    {
	                    *rcond = 1.;
	                    return;
                    }

                    work = complexCalloc( 3*L->nrow );


                    if ( !work )
                    {
	                    SUPERLU_ABORT("Malloc fails for work arrays in cgscon.");
                    }
    
                    // Estimate the norm of inv(A). 
                    ainvnm = 0.;
                    if ( onenrm ) 
                    {
                        kase1 = 1;
                    }
                    else 
                    {
                        kase1 = 2;
                    }
                    kase = 0;

                    do 
                    {
	                    clacon_(&L->nrow, &work[L->nrow], &work[0], &ainvnm, &kase);

	                    if (kase == 0) 
                        {
                            break;
                        }

	                    if (kase == kase1) 
                        {
	                        // Multiply by inv(L). 
	                        sp_ctrsv("Lower", "No transpose", "Unit", L, U, &work[0], info);

	                        // Multiply by inv(U). 
	                        sp_ctrsv("Upper", "No transpose", "Non-unit", L, U, &work[0],info);
	    
	                    } 
                        else 
                        {

	                        // Multiply by inv(U'). 
	                        sp_ctrsv("Upper", "Transpose", "Non-unit", L, U, &work[0], info);

	                        // Multiply by inv(L'). 
	                        sp_ctrsv("Lower", "Transpose", "Unit", L, U, &work[0], info);
	    
	                    }

                    } while ( kase != 0 );

                    // Compute the estimate of the reciprocal condition number. 
                    if (ainvnm != 0.f) 
                    {
                        *rcond = (1.f / ainvnm) / anorm;
                    }

                    SUPERLU_FREE (work);
                    return;
                }
            };
        };
    };
};