#include "stdafx.h"
#include <math.h>
#include <stdlib.h>
#include "hnum_pzsp_defs.h"
namespace harlinn
{
namespace numerics
{
namespace SuperLU
{
namespace DoubleComplex
{
int
pzgstrf_pivotL(
const int pnum, /* process number */
const int jcol, /* current column */
const double u, /* diagonal pivoting threshold */
yes_no_t *usepr, /* re-use the pivot sequence given by
perm_r[]/inv_perm_r[] */
int *perm_r, /* modified - row pivotings */
int *inv_perm_r,/* modified - inverse of perm_r */
int *inv_perm_c,/* in - used to find diagonal of Pc*A*Pc' */
int *pivrow, /* the pivot row for this column */
GlobalLU_t *Glu, /* modified - global LU data structures */
Gstat_t *Gstat /* modified */
)
{
/*
* -- SuperLU MT routine (version 2.0) --
* Lawrence Berkeley National Lab, Univ. of California Berkeley,
* and Xerox Palo Alto Research Center.
* September 10, 2007
*
* Purpose
* =======
* Performs the numerical pivoting on the current column of L,
* and the CDIV operation.
*
* Pivot policy
* ============
* (1) Compute thresh = u * max_(i>=j) abs(A_ij);
* (2) IF user specifies pivot row k and abs(A_kj) >= thresh THEN
* pivot row = k;
* ELSE IF abs(A_jj) >= thresh THEN
* pivot row = j;
* ELSE
* pivot row = m;
*
* Note: If you absolutely want to use a given pivot order, then set u=0.0.
*
* Return value
* ============
* 0 success;
* i > 0 U(i,i) is exactly zero.
*
*/
doublecomplex one = {1.0, 0.0};
register int fsupc; /* first column in the supernode */
register int nsupc; /* no of columns in the supernode */
register int nsupr; /* no of rows in the supernode */
register int lptr; /* the starting subscript of the supernode */
register int pivptr, old_pivptr, diag, diagind;
register int isub, icol, k, itemp;
register double pivmax, rtemp, thresh;
doublecomplex *lu_sup_ptr;
doublecomplex *lu_col_ptr;
doublecomplex temp;
int *lsub_ptr;
int *lsub;
doublecomplex *lusup;
int *xlusup;
/* Initialize pointers */
lsub = Glu->lsub;
lusup = Glu->lusup;
xlusup = Glu->xlusup;
fsupc = Glu->xsup[Glu->supno[jcol]];
nsupc = jcol - fsupc; /* excluding jcol; nsupc >= 0 */
lptr = Glu->xlsub[fsupc];
nsupr = Glu->xlsub_end[fsupc] - lptr;
lu_sup_ptr = &lusup[xlusup[fsupc]]; /* start of the current supernode */
lu_col_ptr = &lusup[xlusup[jcol]]; /* start of jcol in the supernode */
lsub_ptr = &lsub[lptr]; /* start of row indices of the supernode */
#ifdef CHK_PIVOT
printf("Before cdiv: col %d\n", jcol);
for (k = nsupc; k < nsupr; k++)
printf(" lu[%d] %f\n", lsub_ptr[k], lu_col_ptr[k]);
#endif
/* Determine the largest abs numerical value for partial pivoting;
Also search for user-specified pivot, and diagonal element. */
if ( *usepr == YES ) *pivrow = inv_perm_r[jcol];
diagind = inv_perm_c[jcol];
pivmax = 0.0;
pivptr = nsupc;
diag = EMPTY;
old_pivptr = nsupc;
for (isub = nsupc; isub < nsupr; ++isub) {
rtemp = z_abs1 (&lu_col_ptr[isub]);
if ( rtemp > pivmax ) {
pivmax = rtemp;
pivptr = isub;
}
if ( *usepr == YES && lsub_ptr[isub] == *pivrow ) old_pivptr = isub;
if ( lsub_ptr[isub] == diagind ) diag = isub;
}
/* Test for singularity */
if ( pivmax == 0.0 ) {
*pivrow = lsub_ptr[pivptr];
perm_r[*pivrow] = jcol;
inv_perm_r[jcol] = *pivrow;
*usepr = NO;
return (jcol+1);
}
thresh = u * pivmax;
/* Choose appropriate pivotal element by our policy. */
if ( *usepr == YES ) {
rtemp = z_abs1 (&lu_col_ptr[old_pivptr]);
if ( rtemp != 0.0 && rtemp >= thresh )
pivptr = old_pivptr;
else
*usepr = NO;
}
if ( *usepr == NO ) {
/* Can we use diagonal as pivot? */
if ( diag >= 0 ) { /* diagonal exists */
rtemp = z_abs1 (&lu_col_ptr[diag]);
if ( rtemp != 0.0 && rtemp >= thresh ) pivptr = diag;
}
*pivrow = lsub_ptr[pivptr];
}
/* Record pivot row */
perm_r[*pivrow] = jcol;
inv_perm_r[jcol] = *pivrow;
/* Interchange row subscripts */
if ( pivptr != nsupc ) {
itemp = lsub_ptr[pivptr];
lsub_ptr[pivptr] = lsub_ptr[nsupc];
lsub_ptr[nsupc] = itemp;
/* Interchange numerical values as well, for the whole supernode,
* such that L is indexed the same way as A.
*/
k = 0;
for (icol = 0; icol <= nsupc; ++icol, k += nsupr) {
itemp = pivptr + k;
temp = lu_sup_ptr[itemp];
lu_sup_ptr[itemp] = lu_sup_ptr[nsupc + k];
lu_sup_ptr[nsupc + k] = temp;
}
} /* if */
/* CDIV operation */
/* ops[FACT] += nsupr - nsupc;*/
Gstat->procstat[pnum].fcops += 10 * (nsupr - nsupc);
z_div(&temp, &one, &lu_col_ptr[nsupc]);
for (k = nsupc+1; k < nsupr; k++)
zz_mult(&lu_col_ptr[k], &lu_col_ptr[k], &temp);
#ifdef CHK_PIVOT
printf("After cdiv: col %d\n", jcol);
for (k = nsupc; k < nsupr; k++)
printf(" lu[%d] %f\n", lsub_ptr[k], lu_col_ptr[k]);
#endif
return 0;
}
};
};
};
};
Submit an article or tip