kusano 7d535a
kusano 7d535a
/*! @file dpivotgrowth.c
kusano 7d535a
 * \brief Computes the reciprocal pivot growth factor
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 * -- SuperLU routine (version 2.0) --
kusano 7d535a
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
kusano 7d535a
 * and Lawrence Berkeley National Lab.
kusano 7d535a
 * November 15, 1997
kusano 7d535a
 * 
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
 * Compute the reciprocal pivot growth factor of the leading ncols columns
kusano 7d535a
 * of the matrix, using the formula:
kusano 7d535a
 *     min_j ( max_i(abs(A_ij)) / max_i(abs(U_ij)) )
kusano 7d535a
 *
kusano 7d535a
 * Arguments
kusano 7d535a
 * =========
kusano 7d535a
 *
kusano 7d535a
 * ncols    (input) int
kusano 7d535a
 *          The number of columns of matrices A, L and U.
kusano 7d535a
 *
kusano 7d535a
 * A        (input) SuperMatrix*
kusano 7d535a
 *	    Original matrix A, permuted by columns, of dimension
kusano 7d535a
 *          (A->nrow, A->ncol). The type of A can be:
kusano 7d535a
 *          Stype = NC; Dtype = SLU_D; Mtype = GE.
kusano 7d535a
 *
kusano 7d535a
 * L        (output) SuperMatrix*
kusano 7d535a
 *          The factor L from the factorization Pr*A=L*U; use compressed row 
kusano 7d535a
 *          subscripts storage for supernodes, i.e., L has type: 
kusano 7d535a
 *          Stype = SC; Dtype = SLU_D; Mtype = TRLU.
kusano 7d535a
 *
kusano 7d535a
 * U        (output) SuperMatrix*
kusano 7d535a
 *	    The factor U from the factorization Pr*A*Pc=L*U. Use column-wise
kusano 7d535a
 *          storage scheme, i.e., U has types: Stype = NC;
kusano 7d535a
 *          Dtype = SLU_D; Mtype = TRU.
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
kusano 7d535a
double
kusano 7d535a
dPivotGrowth(int ncols, SuperMatrix *A, int *perm_c, 
kusano 7d535a
             SuperMatrix *L, SuperMatrix *U)
kusano 7d535a
{
kusano 7d535a
kusano 7d535a
    NCformat *Astore;
kusano 7d535a
    SCformat *Lstore;
kusano 7d535a
    NCformat *Ustore;
kusano 7d535a
    double  *Aval, *Lval, *Uval;
kusano 7d535a
    int      fsupc, nsupr, luptr, nz_in_U;
kusano 7d535a
    int      i, j, k, oldcol;
kusano 7d535a
    int      *inv_perm_c;
kusano 7d535a
    double   rpg, maxaj, maxuj;
kusano 7d535a
    double   smlnum;
kusano 7d535a
    double   *luval;
kusano 7d535a
   
kusano 7d535a
    /* Get machine constants. */
kusano 7d535a
    smlnum = dlamch_("S");
kusano 7d535a
    rpg = 1. / smlnum;
kusano 7d535a
kusano 7d535a
    Astore = A->Store;
kusano 7d535a
    Lstore = L->Store;
kusano 7d535a
    Ustore = U->Store;
kusano 7d535a
    Aval = Astore->nzval;
kusano 7d535a
    Lval = Lstore->nzval;
kusano 7d535a
    Uval = Ustore->nzval;
kusano 7d535a
    
kusano 7d535a
    inv_perm_c = (int *) SUPERLU_MALLOC(A->ncol*sizeof(int));
kusano 7d535a
    for (j = 0; j < A->ncol; ++j) inv_perm_c[perm_c[j]] = j;
kusano 7d535a
kusano 7d535a
    for (k = 0; k <= Lstore->nsuper; ++k) {
kusano 7d535a
	fsupc = L_FST_SUPC(k);
kusano 7d535a
	nsupr = L_SUB_START(fsupc+1) - L_SUB_START(fsupc);
kusano 7d535a
	luptr = L_NZ_START(fsupc);
kusano 7d535a
	luval = &Lval[luptr];
kusano 7d535a
	nz_in_U = 1;
kusano 7d535a
	
kusano 7d535a
	for (j = fsupc; j < L_FST_SUPC(k+1) && j < ncols; ++j) {
kusano 7d535a
	    maxaj = 0.;
kusano 7d535a
            oldcol = inv_perm_c[j];
kusano 7d535a
	    for (i = Astore->colptr[oldcol]; i < Astore->colptr[oldcol+1]; ++i)
kusano 7d535a
		maxaj = SUPERLU_MAX( maxaj, fabs(Aval[i]) );
kusano 7d535a
	
kusano 7d535a
	    maxuj = 0.;
kusano 7d535a
	    for (i = Ustore->colptr[j]; i < Ustore->colptr[j+1]; i++)
kusano 7d535a
		maxuj = SUPERLU_MAX( maxuj, fabs(Uval[i]) );
kusano 7d535a
	    
kusano 7d535a
	    /* Supernode */
kusano 7d535a
	    for (i = 0; i < nz_in_U; ++i)
kusano 7d535a
		maxuj = SUPERLU_MAX( maxuj, fabs(luval[i]) );
kusano 7d535a
kusano 7d535a
	    ++nz_in_U;
kusano 7d535a
	    luval += nsupr;
kusano 7d535a
kusano 7d535a
	    if ( maxuj == 0. )
kusano 7d535a
		rpg = SUPERLU_MIN( rpg, 1.);
kusano 7d535a
	    else
kusano 7d535a
		rpg = SUPERLU_MIN( rpg, maxaj / maxuj );
kusano 7d535a
	}
kusano 7d535a
	
kusano 7d535a
	if ( j >= ncols ) break;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    SUPERLU_FREE(inv_perm_c);
kusano 7d535a
    return (rpg);
kusano 7d535a
}