symbfact.c

SuperLU 2.2版本。对大型、稀疏、非对称的线性系统的直接求解

 *   marker: A-row --> A-row/col (0/1)
 *   repfnz: SuperA-col --> PA-row
 *   parent: SuperA-col --> SuperA-col
 *   xplore: SuperA-col --> index to L-structure
 *
 * Return value
 * ============
 *     0  success;
 *   > 0  number of bytes allocated when run out of space.
 *
 */
    NCPformat *Astore;
    int_t     *asub, *xa_begin, *xa_end;
    int_t     jcolp1, jcolm1, jsuper, nsuper, nextl;
    int_t     k, krep, krow, kmark, kperm;
    int_t     fsupc; /* first column of a supernode */
    int_t     myfnz; /* first nonzero column of a U-segment */
    int_t     chperm, chmark, chrep, kchild;
    int_t     xdfs, maxdfs, kpar, oldrep;
    int_t     jptr, jm1ptr;
    int_t     ito, ifrom, istop;	/* used to compress row subscripts */
    int_t     *xsup, *supno, *lsub, *xlsub;
    int_t     nzlmax;
    static int_t first = 1, maxsuper;
    int_t     mem_error;
    
    /* Initializations */
    Astore   = A->Store;
    asub     = Astore->rowind;
    xa_begin = Astore->colbeg;
    xa_end   = Astore->colend;
    xsup     = Glu_persist->xsup;
    supno    = Glu_persist->supno;
    lsub     = Glu_freeable->lsub;
    xlsub    = Glu_freeable->xlsub;
    nzlmax   = Glu_freeable->nzlmax;
    jcolp1   = jcol + 1;
    jcolm1   = jcol - 1;
    jsuper   = nsuper = supno[jcol];
    nextl    = xlsub[jcol];
    if ( first ) {
	maxsuper = sp_ienv_dist(3);
	first = 0;
    }
    
    *nseg = 0;

    /* For each nonzero in A[*,jcol] perform depth-first search. */
    for (k = xa_begin[jcol]; k < xa_end[jcol]; ++k) {
	krow = asub[k];
	kmark = marker[krow];

	/* krow was visited before, go to the next nonzero. */
	if ( kmark == jcol ) continue; 
	
	/* 
	 * For each unmarked neighber krow of jcol ...
	 */
	marker[krow] = jcol;
	kperm = perm_r[krow];

	if ( kperm == EMPTY ) {
	    /* krow is in L:
	     * place it in structure of L[*,jcol].
	     */
	    lsub[nextl++] = krow; 	/* krow is indexed into A */
	    if ( nextl >= nzlmax ) {
		if ( mem_error = symbfact_SubXpand(A->ncol, jcol, nextl, LSUB,
						   &nzlmax, Glu_freeable) )
		    return (mem_error);
		lsub = Glu_freeable->lsub;
	    }
	    if ( kmark != jcolm1 ) jsuper = EMPTY; /* Row index subset test */
	} else {
	    /* krow is in U:
	     * If its supernode krep has been explored, update repfnz[*].
	     */
	    krep = xsup[supno[kperm]+1] - 1;
	    myfnz = repfnz[krep];
	    
	    if ( myfnz != EMPTY ) { /* krep was visited before */
		if ( kperm < myfnz ) repfnz[krep] = kperm;
		/* continue; */
	    } else {
		/* Otherwise perform DFS, starting at krep */
		oldrep = EMPTY;
		parent[krep] = oldrep;
		repfnz[krep] = kperm;
		xdfs = xlsub[krep];
		maxdfs = xprune[krep];
		
		do {
		    /* 
		     * For each unmarked kchild of krep 
		     */
		    while ( xdfs < maxdfs ) {
			kchild = lsub[xdfs++];
			chmark = marker[kchild];
			
			if ( chmark != jcol ) { /* Not reached yet */
			    marker[kchild] = jcol;
			    chperm = perm_r[kchild];
			    
			    /* Case kchild is in L: place it in L[*,k] */
			    if ( chperm == EMPTY ) {
				lsub[nextl++] = kchild;
				if ( nextl >= nzlmax ) {
				    if ( mem_error =
					symbfact_SubXpand(A->ncol, jcol, nextl,
							  LSUB, &nzlmax,
							  Glu_freeable) )
					return (mem_error);
				    lsub = Glu_freeable->lsub;
				}
				if ( chmark != jcolm1 ) jsuper = EMPTY;
			    } else {
				/* Case kchild is in U: 
				 * chrep = its supernode-rep. If its rep 
				 * has been explored, update its repfnz[*].
				 */
				chrep = xsup[supno[chperm]+1] - 1;
				myfnz = repfnz[chrep];
				if ( myfnz != EMPTY ) {/* Visited before */
				    if (chperm < myfnz) repfnz[chrep] = chperm;
				} else {
				    /* Continue DFS at sup-rep of kchild */
				    xplore[krep] = xdfs;
				    oldrep = krep;
				    krep = chrep; /* Go deeper down G(L') */
				    parent[krep] = oldrep;
				    repfnz[krep] = chperm;
				    xdfs = xlsub[krep];     
				    maxdfs = xprune[krep];
				} /* else */
			    } /* else */
			} /* if */
			
		    } /* while */
		    
		    /* krow has no more unexplored neighbors:
		     *    place supernode-rep krep in postorder DFS;
		     *    backtrack DFS to its parent.
		     */
		    segrep[*nseg] = krep;
		    ++(*nseg);
		    kpar = parent[krep]; /* Pop from stack; recurse */
		    if ( kpar == EMPTY ) break; /* DFS done */
		    krep = kpar;
		    xdfs = xplore[krep];
		    maxdfs = xprune[krep];
		} while ( kpar != EMPTY ); /* Until empty stack */
	    } /* else */
	} /* else */
    } /* for each nonzero ... */
    
    /* Check to see if jcol belongs in the same supernode as jcol-1 */
    if ( jcol == 0 ) { /* Do nothing for column 0 */
	nsuper = supno[0] = 0;
    } else {
	fsupc = xsup[nsuper];
	jptr = xlsub[jcol];	/* Not compressed yet */
	jm1ptr = xlsub[jcolm1];
	
#ifdef T2_SUPER
	if ( (nextl-jptr != jptr-jm1ptr-1) ) jsuper = EMPTY;
#endif
	/* Make sure the number of columns in a supernode doesn't
	   exceed threshold. */
	if ( jcol - fsupc >= maxsuper ) jsuper = EMPTY;
	
	/* If jcol starts a new supernode, reclaim storage space in
	 * lsub[*] from the previous supernode. Note we only store
	 * the subscript set of the first and last columns of
	 * a supernode. (first for G(L'), last for pruned graph)
	 */
	if ( jsuper ==EMPTY ) { /* Starts a new supernode */
	    if ( (fsupc < jcolm1-1) ) { /* >= 3 columns in nsuper */
#ifdef CHK_COMPRESS
		printf("  Compress lsub[] at super %d-%d\n",fsupc,jcolm1);
#endif
		ito = xlsub[fsupc+1];
		xlsub[jcolm1] = ito;
		istop = ito + jptr - jm1ptr;
		xprune[jcolm1] = istop; /* Initialize xprune[jcol-1] */
		xlsub[jcol] = istop;
		for (ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito)
		    lsub[ito] = lsub[ifrom];
		nextl = ito;            /* = istop + length(jcol) */
	    }
	    ++nsuper;
	    supno[jcol] = nsuper;
	} /* if a new supernode */
	
    } /* else: jcol > 0 */ 
    
    /* Tidy up the pointers before exit */
    xsup[nsuper+1] = jcolp1;
    supno[jcolp1]  = nsuper;
    xprune[jcol]   = nextl; /* Initialize an upper bound for pruning. */
    xlsub[jcolp1]  = nextl;
    return 0;
} /* COLUMN_DFS */

/************************************************************************/
static int_t pivotL
/************************************************************************/
(
 const int_t jcol,     /* current column number     (input)    */
 int_t       *perm_r,  /* row permutation vector    (output)   */
 int_t       *pivrow,  /* the pivot row index       (output)   */
 Glu_persist_t *Glu_persist,   /* global LU data structures (modified) */
 Glu_freeable_t *Glu_freeable
 )
{
/* Purpose
 * =======
 *   pivotL() interchanges row subscripts so that each diagonal block of a
 *   supernode in L has the row subscripts sorted in order of pivots.
 *   The row subscripts in the off-diagonal block are not sorted.
 *
 */
    int_t  fsupc;	/* first column in the supernode */
    int_t  nsupc;	/* number of columns in the supernode */
    int_t  nsupr;       /* number of rows in the supernode */
    int_t  lptr;	/* point_ts to the first subscript of the supernode */
    int_t  diag, diagind;
    int_t  *lsub_ptr;
    int_t  isub, itemp;
    int_t  *lsub, *xlsub;

    /* Initialization. */
    lsub     = Glu_freeable->lsub;
    xlsub    = Glu_freeable->xlsub;
    fsupc    = (Glu_persist->xsup)[(Glu_persist->supno)[jcol]];
    nsupc    = jcol - fsupc; /* excluding jcol; nsupc >= 0 */
    lptr     = xlsub[fsupc];
    nsupr    = xlsub[fsupc+1] - lptr;
    lsub_ptr = &lsub[lptr]; /* start of row indices of the supernode */

    /* Search for diagonal element. */
    /* diagind = iperm_c[jcol];*/
    diagind = jcol;
    diag = EMPTY;
    for (isub = nsupc; isub < nsupr; ++isub)
	if ( lsub_ptr[isub] == diagind ) {
	    diag = isub;
	    break;
	}

    /* Diagonal pivot exists? */
    if ( diag == EMPTY ) {
	printf("At column %d, ", jcol);
	ABORT("pivotL() encounters zero diagonal");
    }

    /* Record pivot row. */
    *pivrow = lsub_ptr[diag];
    perm_r[*pivrow] = jcol;  /* perm_r[] should be Identity. */
    /*assert(*pivrow==jcol);*/
    
    /* Interchange row subscripts. */
    if ( diag != nsupc ) {
	itemp = lsub_ptr[diag];
	lsub_ptr[diag] = lsub_ptr[nsupc];
	lsub_ptr[nsupc] = itemp;
    }

    return 0;
} /* PIVOTL */


/************************************************************************/
static int_t set_usub
/************************************************************************/
(
 const int_t n,       /* total number of columns (input) */
 const int_t jcol,    /* current column number (input) */
 const int_t nseg,    /* number of supernodal segments in U[*,jcol] (input) */
 int_t       *segrep, /* list of U-segment representatives (output) */
 int_t       *repfnz, /* list of first nonzeros in the U-segments (output) */
 Glu_persist_t *Glu_persist,   /* global LU data structures (modified) */
 Glu_freeable_t *Glu_freeable
 )
{
/* 
 * Purpose
 * =======
 *   set_usub() sets up data structure to store supernodal segments in U.
 *   The supernodal segments in each column are stored in topological order.
 *   
 * NOTE
 * ====
 *   For each supernodal segment, we only store the index of the first
 *   nonzero index, rather than the indices of the whole segment, because
 *   those indices can be generated from first nonzero and supnodal
 *   representative.
 *   Therefore, for G(U), we store the "skeleton" of it.
 *   
 */
    int_t ksub, krep, ksupno;
    int_t k, kfnz;
    int_t jsupno, nextu;
    int_t new_next, mem_error;
    int_t *supno;
    int_t *usub, *xusub;
    int_t nzumax;

    supno   = Glu_persist->supno;
    usub    = Glu_freeable->usub;
    xusub   = Glu_freeable->xusub;
    nzumax  = Glu_freeable->nzumax;
    jsupno  = supno[jcol];
    nextu   = xusub[jcol];

    new_next = nextu + nseg;
    while ( new_next > nzumax ) {
	if (mem_error = symbfact_SubXpand(n, jcol, nextu, USUB, &nzumax,
					  Glu_freeable))
	    return (mem_error);
	usub = Glu_freeable->usub;
    }

    /* We store U-segments in topological order. */
    k = nseg - 1;
    for (ksub = 0; ksub < nseg; ++ksub) {
	krep = segrep[k--];
	ksupno = supno[krep];

	if ( ksupno != jsupno ) { /* Should go into usub[*] */
	    kfnz = repfnz[krep];
	    if ( kfnz != EMPTY ) { /* Nonzero U-segment */
		usub[nextu++] = kfnz;

/*	    	fsupc = xsup[ksupno];
	        isub = xlsub[fsupc] + kfnz - fsupc;
		irow = lsub[isub];
		usub[nextu++] = perm_r[irow];*/
	    } /* if ... */
	} /* if ... */
    } /* for each segment... */

    xusub[jcol + 1] = nextu; /* Close U[*,jcol] */
    return 0;
} /* SET_USUB */


/************************************************************************/
static void pruneL
/************************************************************************/
(
 const int_t  jcol,    /* in */
 const int_t  *perm_r, /* in */
 const int_t  pivrow,  /* in */
 const int_t  nseg,    /* in */
 const int_t  *segrep, /* in */
 const int_t  *repfnz, /* in */
 int_t  *xprune,       /* out */
 Glu_persist_t *Glu_persist,   /* global LU data structures (modified) */
 Glu_freeable_t *Glu_freeable
 )
{
/*
 * Purpose
 * =======
 *   pruneL() prunes the L-structure of supernodes whose L-structure
 *   contains the current pivot row "pivrow".
 *
 */
    int_t  jsupno, irep, irep1, kmin, kmax, krow;
    int_t  i, ktemp;
    int_t  do_prune; /* logical variable */
    int_t  *supno;
    int_t  *lsub, *xlsub;

    supno  = Glu_persist->supno;
    lsub   = Glu_freeable->lsub;
    xlsub  = Glu_freeable->xlsub;
    
    /*
     * For each supernode-rep irep in U[*,j]
     */
    jsupno = supno[jcol];
    for (i = 0; i < nseg; i++) {
	irep = segrep[i];
	irep1 = irep + 1;

	/* Do not prune with a zero U-segment */
 	if ( repfnz[irep] == EMPTY ) continue;

	/*
	 * If irep has not been pruned & it has a nonzero in row L[pivrow,i]
	 */
	do_prune = FALSE;
	if ( supno[irep] != jsupno ) {
	    if ( xprune[irep] >= xlsub[irep1] ) {
		kmin = xlsub[irep];
		kmax = xlsub[irep1] - 1;
		for (krow = kmin; krow <= kmax; ++krow) 
		    if ( lsub[krow] == pivrow ) {
			do_prune = TRUE;
			break;
		    }
	    }
	    
    	    if ( do_prune ) {
	     	/* Do a quicksort-type partition. */
	        while ( kmin <= kmax ) {
	    	    if ( perm_r[lsub[kmax]] == EMPTY ) 
			kmax--;
		    else if ( perm_r[lsub[kmin]] != EMPTY )
			kmin++;
		    else { /* kmin below pivrow, and kmax above pivrow: 
		            * 	   interchange the two subscripts
			    */
		        ktemp = lsub[kmin];
		        lsub[kmin] = lsub[kmax];
		        lsub[kmax] = ktemp;
		        kmin++;
		        kmax--;
		    }
	        } /* while */
	        xprune[irep] = kmin; /* Pruning */
#if ( DEBUGlevel>=3 )
		printf(".. pruneL(): use col %d: xprune[%d] = %d\n",
		       jcol, irep, kmin);
#endif
	    } /* if do_prune */
	} /* if */
    } /* for each U-segment ... */
} /* PRUNEL */