|
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 |
|