kusano 7d535a
kusano 7d535a
/*! @file dgscon.c
kusano 7d535a
 * \brief Estimates reciprocal of the condition number of a general matrix
kusano 7d535a
 * 
kusano 7d535a
 * 
kusano 7d535a
 * -- SuperLU routine (version 3.0) --
kusano 7d535a
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
kusano 7d535a
 * and Lawrence Berkeley National Lab.
kusano 7d535a
 * October 15, 2003
kusano 7d535a
 *
kusano 7d535a
 * Modified from lapack routines DGECON.
kusano 7d535a
 *  
kusano 7d535a
 */
kusano 7d535a
kusano 7d535a
/*
kusano 7d535a
 * File name:	dgscon.c
kusano 7d535a
 * History:     Modified from lapack routines DGECON.
kusano 7d535a
 */
kusano 7d535a
#include <math.h></math.h>
kusano 7d535a
#include "slu_ddefs.h"
kusano 7d535a
kusano 7d535a
/*! \brief
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 *   Purpose   
kusano 7d535a
 *   =======   
kusano 7d535a
 *
kusano 7d535a
 *   DGSCON estimates the reciprocal of the condition number of a general 
kusano 7d535a
 *   real matrix A, in either the 1-norm or the infinity-norm, using   
kusano 7d535a
 *   the LU factorization computed by DGETRF.   *
kusano 7d535a
 *
kusano 7d535a
 *   An estimate is obtained for norm(inv(A)), and the reciprocal of the   
kusano 7d535a
 *   condition number is computed as   
kusano 7d535a
 *      RCOND = 1 / ( norm(A) * norm(inv(A)) ).   
kusano 7d535a
 *
kusano 7d535a
 *   See supermatrix.h for the definition of 'SuperMatrix' structure.
kusano 7d535a
 * 
kusano 7d535a
 *   Arguments   
kusano 7d535a
 *   =========   
kusano 7d535a
 *
kusano 7d535a
 *    NORM    (input) char*
kusano 7d535a
 *            Specifies whether the 1-norm condition number or the   
kusano 7d535a
 *            infinity-norm condition number is required:   
kusano 7d535a
 *            = '1' or 'O':  1-norm;   
kusano 7d535a
 *            = 'I':         Infinity-norm.
kusano 7d535a
 *	    
kusano 7d535a
 *    L       (input) SuperMatrix*
kusano 7d535a
 *            The factor L from the factorization Pr*A*Pc=L*U as computed by
kusano 7d535a
 *            dgstrf(). Use compressed row subscripts storage for supernodes,
kusano 7d535a
 *            i.e., L has types: Stype = SLU_SC, Dtype = SLU_D, Mtype = SLU_TRLU.
kusano 7d535a
 * 
kusano 7d535a
 *    U       (input) SuperMatrix*
kusano 7d535a
 *            The factor U from the factorization Pr*A*Pc=L*U as computed by
kusano 7d535a
 *            dgstrf(). Use column-wise storage scheme, i.e., U has types:
kusano 7d535a
 *            Stype = SLU_NC, Dtype = SLU_D, Mtype = SLU_TRU.
kusano 7d535a
 *	    
kusano 7d535a
 *    ANORM   (input) double
kusano 7d535a
 *            If NORM = '1' or 'O', the 1-norm of the original matrix A.   
kusano 7d535a
 *            If NORM = 'I', the infinity-norm of the original matrix A.
kusano 7d535a
 *	    
kusano 7d535a
 *    RCOND   (output) double*
kusano 7d535a
 *           The reciprocal of the condition number of the matrix A,   
kusano 7d535a
 *           computed as RCOND = 1/(norm(A) * norm(inv(A))).
kusano 7d535a
 *	    
kusano 7d535a
 *    INFO    (output) int*
kusano 7d535a
 *           = 0:  successful exit   
kusano 7d535a
 *           < 0:  if INFO = -i, the i-th argument had an illegal value   
kusano 7d535a
 *
kusano 7d535a
 *    ===================================================================== 
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
kusano 7d535a
void
kusano 7d535a
dgscon(char *norm, SuperMatrix *L, SuperMatrix *U,
kusano 7d535a
       double anorm, double *rcond, SuperLUStat_t *stat, int *info)
kusano 7d535a
{
kusano 7d535a
kusano 7d535a
kusano 7d535a
    /* Local variables */
kusano 7d535a
    int    kase, kase1, onenrm, i;
kusano 7d535a
    double ainvnm;
kusano 7d535a
    double *work;
kusano 7d535a
    int    *iwork;
kusano 7d535a
    extern int drscl_(int *, double *, double *, int *);
kusano 7d535a
kusano 7d535a
    extern int dlacon_(int *, double *, double *, int *, double *, int *);
kusano 7d535a
kusano 7d535a
    
kusano 7d535a
    /* Test the input parameters. */
kusano 7d535a
    *info = 0;
kusano 7d535a
    onenrm = *(unsigned char *)norm == '1' || lsame_(norm, "O");
kusano 7d535a
    if (! onenrm && ! lsame_(norm, "I")) *info = -1;
kusano 7d535a
    else if (L->nrow < 0 || L->nrow != L->ncol ||
kusano 7d535a
             L->Stype != SLU_SC || L->Dtype != SLU_D || L->Mtype != SLU_TRLU)
kusano 7d535a
	 *info = -2;
kusano 7d535a
    else if (U->nrow < 0 || U->nrow != U->ncol ||
kusano 7d535a
             U->Stype != SLU_NC || U->Dtype != SLU_D || U->Mtype != SLU_TRU) 
kusano 7d535a
	*info = -3;
kusano 7d535a
    if (*info != 0) {
kusano 7d535a
	i = -(*info);
kusano 7d535a
	xerbla_("dgscon", &i);
kusano 7d535a
	return;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    /* Quick return if possible */
kusano 7d535a
    *rcond = 0.;
kusano 7d535a
    if ( L->nrow == 0 || U->nrow == 0) {
kusano 7d535a
	*rcond = 1.;
kusano 7d535a
	return;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    work = doubleCalloc( 3*L->nrow );
kusano 7d535a
    iwork = intMalloc( L->nrow );
kusano 7d535a
kusano 7d535a
kusano 7d535a
    if ( !work || !iwork )
kusano 7d535a
	ABORT("Malloc fails for work arrays in dgscon.");
kusano 7d535a
    
kusano 7d535a
    /* Estimate the norm of inv(A). */
kusano 7d535a
    ainvnm = 0.;
kusano 7d535a
    if ( onenrm ) kase1 = 1;
kusano 7d535a
    else kase1 = 2;
kusano 7d535a
    kase = 0;
kusano 7d535a
kusano 7d535a
    do {
kusano 7d535a
	dlacon_(&L->nrow, &work[L->nrow], &work[0], &iwork[0], &ainvnm, &kase);
kusano 7d535a
kusano 7d535a
	if (kase == 0) break;
kusano 7d535a
kusano 7d535a
	if (kase == kase1) {
kusano 7d535a
	    /* Multiply by inv(L). */
kusano 7d535a
	    sp_dtrsv("L", "No trans", "Unit", L, U, &work[0], stat, info);
kusano 7d535a
kusano 7d535a
	    /* Multiply by inv(U). */
kusano 7d535a
	    sp_dtrsv("U", "No trans", "Non-unit", L, U, &work[0], stat, info);
kusano 7d535a
	    
kusano 7d535a
	} else {
kusano 7d535a
kusano 7d535a
	    /* Multiply by inv(U'). */
kusano 7d535a
	    sp_dtrsv("U", "Transpose", "Non-unit", L, U, &work[0], stat, info);
kusano 7d535a
kusano 7d535a
	    /* Multiply by inv(L'). */
kusano 7d535a
	    sp_dtrsv("L", "Transpose", "Unit", L, U, &work[0], stat, info);
kusano 7d535a
	    
kusano 7d535a
	}
kusano 7d535a
kusano 7d535a
    } while ( kase != 0 );
kusano 7d535a
kusano 7d535a
    /* Compute the estimate of the reciprocal condition number. */
kusano 7d535a
    if (ainvnm != 0.) *rcond = (1. / ainvnm) / anorm;
kusano 7d535a
kusano 7d535a
    SUPERLU_FREE (work);
kusano 7d535a
    SUPERLU_FREE (iwork);
kusano 7d535a
    return;
kusano 7d535a
kusano 7d535a
} /* dgscon */
kusano 7d535a