kusano 7d535a
kusano 7d535a
/*! @file dpivotL.c
kusano 7d535a
 * \brief Performs numerical pivoting
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
 * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
kusano 7d535a
 *
kusano 7d535a
 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
kusano 7d535a
 * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
kusano 7d535a
 * 
kusano 7d535a
 * Permission is hereby granted to use or copy this program for any
kusano 7d535a
 * purpose, provided the above notices are retained on all copies.
kusano 7d535a
 * Permission to modify the code and to distribute modified code is
kusano 7d535a
 * granted, provided the above notices are retained, and a notice that
kusano 7d535a
 * the code was modified is included with the above copyright notice.
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
kusano 7d535a
kusano 7d535a
#include <math.h></math.h>
kusano 7d535a
#include <stdlib.h></stdlib.h>
kusano 7d535a
#include "slu_ddefs.h"
kusano 7d535a
kusano 7d535a
#undef DEBUG
kusano 7d535a
kusano 7d535a
/*! \brief
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 * Purpose
kusano 7d535a
 * =======
kusano 7d535a
 *   Performs the numerical pivoting on the current column of L,
kusano 7d535a
 *   and the CDIV operation.
kusano 7d535a
 *
kusano 7d535a
 *   Pivot policy:
kusano 7d535a
 *   (1) Compute thresh = u * max_(i>=j) abs(A_ij);
kusano 7d535a
 *   (2) IF user specifies pivot row k and abs(A_kj) >= thresh THEN
kusano 7d535a
 *           pivot row = k;
kusano 7d535a
 *       ELSE IF abs(A_jj) >= thresh THEN
kusano 7d535a
 *           pivot row = j;
kusano 7d535a
 *       ELSE
kusano 7d535a
 *           pivot row = m;
kusano 7d535a
 * 
kusano 7d535a
 *   Note: If you absolutely want to use a given pivot order, then set u=0.0.
kusano 7d535a
 *
kusano 7d535a
 *   Return value: 0      success;
kusano 7d535a
 *                 i > 0  U(i,i) is exactly zero.
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
kusano 7d535a
int
kusano 7d535a
dpivotL(
kusano 7d535a
        const int  jcol,     /* in */
kusano 7d535a
        const double u,      /* in - diagonal pivoting threshold */
kusano 7d535a
        int        *usepr,   /* re-use the pivot sequence given by perm_r/iperm_r */
kusano 7d535a
        int        *perm_r,  /* may be modified */
kusano 7d535a
        int        *iperm_r, /* in - inverse of perm_r */
kusano 7d535a
        int        *iperm_c, /* in - used to find diagonal of Pc*A*Pc' */
kusano 7d535a
        int        *pivrow,  /* out */
kusano 7d535a
        GlobalLU_t *Glu,     /* modified - global LU data structures */
kusano 7d535a
	SuperLUStat_t *stat  /* output */
kusano 7d535a
       )
kusano 7d535a
{
kusano 7d535a
kusano 7d535a
    int          fsupc;	    /* first column in the supernode */
kusano 7d535a
    int          nsupc;	    /* no of columns in the supernode */
kusano 7d535a
    int          nsupr;     /* no of rows in the supernode */
kusano 7d535a
    int          lptr;	    /* points to the starting subscript of the supernode */
kusano 7d535a
    int          pivptr, old_pivptr, diag, diagind;
kusano 7d535a
    double       pivmax, rtemp, thresh;
kusano 7d535a
    double       temp;
kusano 7d535a
    double       *lu_sup_ptr; 
kusano 7d535a
    double       *lu_col_ptr;
kusano 7d535a
    int          *lsub_ptr;
kusano 7d535a
    int          isub, icol, k, itemp;
kusano 7d535a
    int          *lsub, *xlsub;
kusano 7d535a
    double       *lusup;
kusano 7d535a
    int          *xlusup;
kusano 7d535a
    flops_t      *ops = stat->ops;
kusano 7d535a
kusano 7d535a
    /* Initialize pointers */
kusano 7d535a
    lsub       = Glu->lsub;
kusano 7d535a
    xlsub      = Glu->xlsub;
kusano 7d535a
    lusup      = Glu->lusup;
kusano 7d535a
    xlusup     = Glu->xlusup;
kusano 7d535a
    fsupc      = (Glu->xsup)[(Glu->supno)[jcol]];
kusano 7d535a
    nsupc      = jcol - fsupc;	        /* excluding jcol; nsupc >= 0 */
kusano 7d535a
    lptr       = xlsub[fsupc];
kusano 7d535a
    nsupr      = xlsub[fsupc+1] - lptr;
kusano 7d535a
    lu_sup_ptr = &lusup[xlusup[fsupc]];	/* start of the current supernode */
kusano 7d535a
    lu_col_ptr = &lusup[xlusup[jcol]];	/* start of jcol in the supernode */
kusano 7d535a
    lsub_ptr   = &lsub[lptr];	/* start of row indices of the supernode */
kusano 7d535a
kusano 7d535a
#ifdef DEBUG
kusano 7d535a
if ( jcol == MIN_COL ) {
kusano 7d535a
    printf("Before cdiv: col %d\n", jcol);
kusano 7d535a
    for (k = nsupc; k < nsupr; k++) 
kusano 7d535a
	printf("  lu[%d] %f\n", lsub_ptr[k], lu_col_ptr[k]);
kusano 7d535a
}
kusano 7d535a
#endif
kusano 7d535a
    
kusano 7d535a
    /* Determine the largest abs numerical value for partial pivoting;
kusano 7d535a
       Also search for user-specified pivot, and diagonal element. */
kusano 7d535a
    if ( *usepr ) *pivrow = iperm_r[jcol];
kusano 7d535a
    diagind = iperm_c[jcol];
kusano 7d535a
    pivmax = 0.0;
kusano 7d535a
    pivptr = nsupc;
kusano 7d535a
    diag = EMPTY;
kusano 7d535a
    old_pivptr = nsupc;
kusano 7d535a
    for (isub = nsupc; isub < nsupr; ++isub) {
kusano 7d535a
	rtemp = fabs (lu_col_ptr[isub]);
kusano 7d535a
	if ( rtemp > pivmax ) {
kusano 7d535a
	    pivmax = rtemp;
kusano 7d535a
	    pivptr = isub;
kusano 7d535a
	}
kusano 7d535a
	if ( *usepr && lsub_ptr[isub] == *pivrow ) old_pivptr = isub;
kusano 7d535a
	if ( lsub_ptr[isub] == diagind ) diag = isub;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    /* Test for singularity */
kusano 7d535a
    if ( pivmax == 0.0 ) {
kusano 7d535a
#if 1
kusano 7d535a
	*pivrow = lsub_ptr[pivptr];
kusano 7d535a
	perm_r[*pivrow] = jcol;
kusano 7d535a
#else
kusano 7d535a
	perm_r[diagind] = jcol;
kusano 7d535a
#endif
kusano 7d535a
	*usepr = 0;
kusano 7d535a
	return (jcol+1);
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    thresh = u * pivmax;
kusano 7d535a
    
kusano 7d535a
    /* Choose appropriate pivotal element by our policy. */
kusano 7d535a
    if ( *usepr ) {
kusano 7d535a
        rtemp = fabs (lu_col_ptr[old_pivptr]);
kusano 7d535a
	if ( rtemp != 0.0 && rtemp >= thresh )
kusano 7d535a
	    pivptr = old_pivptr;
kusano 7d535a
	else
kusano 7d535a
	    *usepr = 0;
kusano 7d535a
    }
kusano 7d535a
    if ( *usepr == 0 ) {
kusano 7d535a
	/* Use diagonal pivot? */
kusano 7d535a
	if ( diag >= 0 ) { /* diagonal exists */
kusano 7d535a
	    rtemp = fabs (lu_col_ptr[diag]);
kusano 7d535a
	    if ( rtemp != 0.0 && rtemp >= thresh ) pivptr = diag;
kusano 7d535a
        }
kusano 7d535a
	*pivrow = lsub_ptr[pivptr];
kusano 7d535a
    }
kusano 7d535a
    
kusano 7d535a
    /* Record pivot row */
kusano 7d535a
    perm_r[*pivrow] = jcol;
kusano 7d535a
    
kusano 7d535a
    /* Interchange row subscripts */
kusano 7d535a
    if ( pivptr != nsupc ) {
kusano 7d535a
	itemp = lsub_ptr[pivptr];
kusano 7d535a
	lsub_ptr[pivptr] = lsub_ptr[nsupc];
kusano 7d535a
	lsub_ptr[nsupc] = itemp;
kusano 7d535a
kusano 7d535a
	/* Interchange numerical values as well, for the whole snode, such 
kusano 7d535a
	 * that L is indexed the same way as A.
kusano 7d535a
 	 */
kusano 7d535a
	for (icol = 0; icol <= nsupc; icol++) {
kusano 7d535a
	    itemp = pivptr + icol * nsupr;
kusano 7d535a
	    temp = lu_sup_ptr[itemp];
kusano 7d535a
	    lu_sup_ptr[itemp] = lu_sup_ptr[nsupc + icol*nsupr];
kusano 7d535a
	    lu_sup_ptr[nsupc + icol*nsupr] = temp;
kusano 7d535a
	}
kusano 7d535a
    } /* if */
kusano 7d535a
kusano 7d535a
    /* cdiv operation */
kusano 7d535a
    ops[FACT] += nsupr - nsupc;
kusano 7d535a
kusano 7d535a
    temp = 1.0 / lu_col_ptr[nsupc];
kusano 7d535a
    for (k = nsupc+1; k < nsupr; k++) 
kusano 7d535a
	lu_col_ptr[k] *= temp;
kusano 7d535a
kusano 7d535a
    return 0;
kusano 7d535a
}
kusano 7d535a