kusano 7d535a
kusano 7d535a
/*! @file sgstrf.c
kusano 7d535a
 * \brief Computes an LU factorization of a general sparse 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
 * 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 "slu_sdefs.h"
kusano 7d535a
kusano 7d535a
/*! \brief
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 * Purpose
kusano 7d535a
 * =======
kusano 7d535a
 *
kusano 7d535a
 * SGSTRF computes an LU factorization of a general sparse m-by-n
kusano 7d535a
 * matrix A using partial pivoting with row interchanges.
kusano 7d535a
 * The factorization has the form
kusano 7d535a
 *     Pr * A = L * U
kusano 7d535a
 * where Pr is a row permutation matrix, L is lower triangular with unit
kusano 7d535a
 * diagonal elements (lower trapezoidal if A->nrow > A->ncol), and U is upper 
kusano 7d535a
 * triangular (upper trapezoidal if A->nrow < A->ncol).
kusano 7d535a
 *
kusano 7d535a
 * See supermatrix.h for the definition of 'SuperMatrix' structure.
kusano 7d535a
 *
kusano 7d535a
 * Arguments
kusano 7d535a
 * =========
kusano 7d535a
 *
kusano 7d535a
 * options (input) superlu_options_t*
kusano 7d535a
 *         The structure defines the input parameters to control
kusano 7d535a
 *         how the LU decomposition will be performed.
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 = SLU_NCP; Dtype = SLU_S; Mtype = SLU_GE.
kusano 7d535a
 *
kusano 7d535a
 * relax    (input) int
kusano 7d535a
 *          To control degree of relaxing supernodes. If the number
kusano 7d535a
 *          of nodes (columns) in a subtree of the elimination tree is less
kusano 7d535a
 *          than relax, this subtree is considered as one supernode,
kusano 7d535a
 *          regardless of the row structures of those columns.
kusano 7d535a
 *
kusano 7d535a
 * panel_size (input) int
kusano 7d535a
 *          A panel consists of at most panel_size consecutive columns.
kusano 7d535a
 *
kusano 7d535a
 * etree    (input) int*, dimension (A->ncol)
kusano 7d535a
 *          Elimination tree of A'*A.
kusano 7d535a
 *          Note: etree is a vector of parent pointers for a forest whose
kusano 7d535a
 *          vertices are the integers 0 to A->ncol-1; etree[root]==A->ncol.
kusano 7d535a
 *          On input, the columns of A should be permuted so that the
kusano 7d535a
 *          etree is in a certain postorder.
kusano 7d535a
 *
kusano 7d535a
 * work     (input/output) void*, size (lwork) (in bytes)
kusano 7d535a
 *          User-supplied work space and space for the output data structures.
kusano 7d535a
 *          Not referenced if lwork = 0;
kusano 7d535a
 *
kusano 7d535a
 * lwork   (input) int
kusano 7d535a
 *         Specifies the size of work array in bytes.
kusano 7d535a
 *         = 0:  allocate space internally by system malloc;
kusano 7d535a
 *         > 0:  use user-supplied work array of length lwork in bytes,
kusano 7d535a
 *               returns error if space runs out.
kusano 7d535a
 *         = -1: the routine guesses the amount of space needed without
kusano 7d535a
 *               performing the factorization, and returns it in
kusano 7d535a
 *               *info; no other side effects.
kusano 7d535a
 *
kusano 7d535a
 * perm_c   (input) int*, dimension (A->ncol)
kusano 7d535a
 *	    Column permutation vector, which defines the 
kusano 7d535a
 *          permutation matrix Pc; perm_c[i] = j means column i of A is 
kusano 7d535a
 *          in position j in A*Pc.
kusano 7d535a
 *          When searching for diagonal, perm_c[*] is applied to the
kusano 7d535a
 *          row subscripts of A, so that diagonal threshold pivoting
kusano 7d535a
 *          can find the diagonal of A, rather than that of A*Pc.
kusano 7d535a
 *
kusano 7d535a
 * perm_r   (input/output) int*, dimension (A->nrow)
kusano 7d535a
 *          Row permutation vector which defines the permutation matrix Pr,
kusano 7d535a
 *          perm_r[i] = j means row i of A is in position j in Pr*A.
kusano 7d535a
 *          If options->Fact = SamePattern_SameRowPerm, the pivoting routine
kusano 7d535a
 *             will try to use the input perm_r, unless a certain threshold
kusano 7d535a
 *             criterion is violated. In that case, perm_r is overwritten by
kusano 7d535a
 *             a new permutation determined by partial pivoting or diagonal
kusano 7d535a
 *             threshold pivoting.
kusano 7d535a
 *          Otherwise, perm_r is output argument;
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 = SLU_SC, Dtype = SLU_S, Mtype = SLU_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 = SLU_NC, 
kusano 7d535a
 *          Dtype = SLU_S, Mtype = SLU_TRU.
kusano 7d535a
 *
kusano 7d535a
 * stat     (output) SuperLUStat_t*
kusano 7d535a
 *          Record the statistics on runtime and floating-point operation count.
kusano 7d535a
 *          See slu_util.h for the definition of 'SuperLUStat_t'.
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
 *          > 0: if info = i, and i is
kusano 7d535a
 *             <= A->ncol: U(i,i) is exactly zero. The factorization has
kusano 7d535a
 *                been completed, but the factor U is exactly singular,
kusano 7d535a
 *                and division by zero will occur if it is used to solve a
kusano 7d535a
 *                system of equations.
kusano 7d535a
 *             > A->ncol: number of bytes allocated when memory allocation
kusano 7d535a
 *                failure occurred, plus A->ncol. If lwork = -1, it is
kusano 7d535a
 *                the estimated amount of space needed, plus A->ncol.
kusano 7d535a
 *
kusano 7d535a
 * ======================================================================
kusano 7d535a
 *
kusano 7d535a
 * Local Working Arrays: 
kusano 7d535a
 * ======================
kusano 7d535a
 *   m = number of rows in the matrix
kusano 7d535a
 *   n = number of columns in the matrix
kusano 7d535a
 *
kusano 7d535a
 *   xprune[0:n-1]: xprune[*] points to locations in subscript 
kusano 7d535a
 *	vector lsub[*]. For column i, xprune[i] denotes the point where 
kusano 7d535a
 *	structural pruning begins. I.e. only xlsub[i],..,xprune[i]-1 need 
kusano 7d535a
 *	to be traversed for symbolic factorization.
kusano 7d535a
 *
kusano 7d535a
 *   marker[0:3*m-1]: marker[i] = j means that node i has been 
kusano 7d535a
 *	reached when working on column j.
kusano 7d535a
 *	Storage: relative to original row subscripts
kusano 7d535a
 *	NOTE: There are 3 of them: marker/marker1 are used for panel dfs, 
kusano 7d535a
 *	      see spanel_dfs.c; marker2 is used for inner-factorization,
kusano 7d535a
 *            see scolumn_dfs.c.
kusano 7d535a
 *
kusano 7d535a
 *   parent[0:m-1]: parent vector used during dfs
kusano 7d535a
 *      Storage: relative to new row subscripts
kusano 7d535a
 *
kusano 7d535a
 *   xplore[0:m-1]: xplore[i] gives the location of the next (dfs) 
kusano 7d535a
 *	unexplored neighbor of i in lsub[*]
kusano 7d535a
 *
kusano 7d535a
 *   segrep[0:nseg-1]: contains the list of supernodal representatives
kusano 7d535a
 *	in topological order of the dfs. A supernode representative is the 
kusano 7d535a
 *	last column of a supernode.
kusano 7d535a
 *      The maximum size of segrep[] is n.
kusano 7d535a
 *
kusano 7d535a
 *   repfnz[0:W*m-1]: for a nonzero segment U[*,j] that ends at a 
kusano 7d535a
 *	supernodal representative r, repfnz[r] is the location of the first 
kusano 7d535a
 *	nonzero in this segment.  It is also used during the dfs: repfnz[r]>0
kusano 7d535a
 *	indicates the supernode r has been explored.
kusano 7d535a
 *	NOTE: There are W of them, each used for one column of a panel. 
kusano 7d535a
 *
kusano 7d535a
 *   panel_lsub[0:W*m-1]: temporary for the nonzeros row indices below 
kusano 7d535a
 *      the panel diagonal. These are filled in during spanel_dfs(), and are
kusano 7d535a
 *      used later in the inner LU factorization within the panel.
kusano 7d535a
 *	panel_lsub[]/dense[] pair forms the SPA data structure.
kusano 7d535a
 *	NOTE: There are W of them.
kusano 7d535a
 *
kusano 7d535a
 *   dense[0:W*m-1]: sparse accumulating (SPA) vector for intermediate values;
kusano 7d535a
 *	    	   NOTE: there are W of them.
kusano 7d535a
 *
kusano 7d535a
 *   tempv[0:*]: real temporary used for dense numeric kernels;
kusano 7d535a
 *	The size of this array is defined by NUM_TEMPV() in slu_sdefs.h.
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
kusano 7d535a
void
kusano 7d535a
sgstrf (superlu_options_t *options, SuperMatrix *A,
kusano 7d535a
        int relax, int panel_size, int *etree, void *work, int lwork,
kusano 7d535a
        int *perm_c, int *perm_r, SuperMatrix *L, SuperMatrix *U,
kusano 7d535a
        SuperLUStat_t *stat, int *info)
kusano 7d535a
{
kusano 7d535a
    /* Local working arrays */
kusano 7d535a
    NCPformat *Astore;
kusano 7d535a
    int       *iperm_r = NULL; /* inverse of perm_r; used when 
kusano 7d535a
                                  options->Fact == SamePattern_SameRowPerm */
kusano 7d535a
    int       *iperm_c; /* inverse of perm_c */
kusano 7d535a
    int       *iwork;
kusano 7d535a
    float    *swork;
kusano 7d535a
    int	      *segrep, *repfnz, *parent, *xplore;
kusano 7d535a
    int	      *panel_lsub; /* dense[]/panel_lsub[] pair forms a w-wide SPA */
kusano 7d535a
    int	      *xprune;
kusano 7d535a
    int	      *marker;
kusano 7d535a
    float    *dense, *tempv;
kusano 7d535a
    int       *relax_end;
kusano 7d535a
    float    *a;
kusano 7d535a
    int       *asub;
kusano 7d535a
    int       *xa_begin, *xa_end;
kusano 7d535a
    int       *xsup, *supno;
kusano 7d535a
    int       *xlsub, *xlusup, *xusub;
kusano 7d535a
    int       nzlumax;
kusano 7d535a
    float fill_ratio = sp_ienv(6);  /* estimated fill ratio */
kusano 7d535a
    static    GlobalLU_t Glu; /* persistent to facilitate multiple factors. */
kusano 7d535a
kusano 7d535a
    /* Local scalars */
kusano 7d535a
    fact_t    fact = options->Fact;
kusano 7d535a
    double    diag_pivot_thresh = options->DiagPivotThresh;
kusano 7d535a
    int       pivrow;   /* pivotal row number in the original matrix A */
kusano 7d535a
    int       nseg1;	/* no of segments in U-column above panel row jcol */
kusano 7d535a
    int       nseg;	/* no of segments in each U-column */
kusano 7d535a
    register int jcol;	
kusano 7d535a
    register int kcol;	/* end column of a relaxed snode */
kusano 7d535a
    register int icol;
kusano 7d535a
    register int i, k, jj, new_next, iinfo;
kusano 7d535a
    int       m, n, min_mn, jsupno, fsupc, nextlu, nextu;
kusano 7d535a
    int       w_def;	/* upper bound on panel width */
kusano 7d535a
    int       usepr, iperm_r_allocated = 0;
kusano 7d535a
    int       nnzL, nnzU;
kusano 7d535a
    int       *panel_histo = stat->panel_histo;
kusano 7d535a
    flops_t   *ops = stat->ops;
kusano 7d535a
kusano 7d535a
    iinfo    = 0;
kusano 7d535a
    m        = A->nrow;
kusano 7d535a
    n        = A->ncol;
kusano 7d535a
    min_mn   = SUPERLU_MIN(m, n);
kusano 7d535a
    Astore   = A->Store;
kusano 7d535a
    a        = Astore->nzval;
kusano 7d535a
    asub     = Astore->rowind;
kusano 7d535a
    xa_begin = Astore->colbeg;
kusano 7d535a
    xa_end   = Astore->colend;
kusano 7d535a
kusano 7d535a
    /* Allocate storage common to the factor routines */
kusano 7d535a
    *info = sLUMemInit(fact, work, lwork, m, n, Astore->nnz,
kusano 7d535a
                       panel_size, fill_ratio, L, U, &Glu, &iwork, &swork);
kusano 7d535a
    if ( *info ) return;
kusano 7d535a
    
kusano 7d535a
    xsup    = Glu.xsup;
kusano 7d535a
    supno   = Glu.supno;
kusano 7d535a
    xlsub   = Glu.xlsub;
kusano 7d535a
    xlusup  = Glu.xlusup;
kusano 7d535a
    xusub   = Glu.xusub;
kusano 7d535a
    
kusano 7d535a
    SetIWork(m, n, panel_size, iwork, &segrep, &parent, &xplore,
kusano 7d535a
	     &repfnz, &panel_lsub, &xprune, &marker);
kusano 7d535a
    sSetRWork(m, panel_size, swork, &dense, &tempv);
kusano 7d535a
    
kusano 7d535a
    usepr = (fact == SamePattern_SameRowPerm);
kusano 7d535a
    if ( usepr ) {
kusano 7d535a
	/* Compute the inverse of perm_r */
kusano 7d535a
	iperm_r = (int *) intMalloc(m);
kusano 7d535a
	for (k = 0; k < m; ++k) iperm_r[perm_r[k]] = k;
kusano 7d535a
	iperm_r_allocated = 1;
kusano 7d535a
    }
kusano 7d535a
    iperm_c = (int *) intMalloc(n);
kusano 7d535a
    for (k = 0; k < n; ++k) iperm_c[perm_c[k]] = k;
kusano 7d535a
kusano 7d535a
    /* Identify relaxed snodes */
kusano 7d535a
    relax_end = (int *) intMalloc(n);
kusano 7d535a
    if ( options->SymmetricMode == YES ) {
kusano 7d535a
        heap_relax_snode(n, etree, relax, marker, relax_end); 
kusano 7d535a
    } else {
kusano 7d535a
        relax_snode(n, etree, relax, marker, relax_end); 
kusano 7d535a
    }
kusano 7d535a
    
kusano 7d535a
    ifill (perm_r, m, EMPTY);
kusano 7d535a
    ifill (marker, m * NO_MARKER, EMPTY);
kusano 7d535a
    supno[0] = -1;
kusano 7d535a
    xsup[0]  = xlsub[0] = xusub[0] = xlusup[0] = 0;
kusano 7d535a
    w_def    = panel_size;
kusano 7d535a
kusano 7d535a
    /* 
kusano 7d535a
     * Work on one "panel" at a time. A panel is one of the following: 
kusano 7d535a
     *	   (a) a relaxed supernode at the bottom of the etree, or
kusano 7d535a
     *	   (b) panel_size contiguous columns, defined by the user
kusano 7d535a
     */
kusano 7d535a
    for (jcol = 0; jcol < min_mn; ) {
kusano 7d535a
kusano 7d535a
	if ( relax_end[jcol] != EMPTY ) { /* start of a relaxed snode */
kusano 7d535a
   	    kcol = relax_end[jcol];	  /* end of the relaxed snode */
kusano 7d535a
	    panel_histo[kcol-jcol+1]++;
kusano 7d535a
kusano 7d535a
	    /* --------------------------------------
kusano 7d535a
	     * Factorize the relaxed supernode(jcol:kcol) 
kusano 7d535a
	     * -------------------------------------- */
kusano 7d535a
	    /* Determine the union of the row structure of the snode */
kusano 7d535a
	    if ( (*info = ssnode_dfs(jcol, kcol, asub, xa_begin, xa_end,
kusano 7d535a
				    xprune, marker, &Glu)) != 0 )
kusano 7d535a
		return;
kusano 7d535a
kusano 7d535a
            nextu    = xusub[jcol];
kusano 7d535a
	    nextlu   = xlusup[jcol];
kusano 7d535a
	    jsupno   = supno[jcol];
kusano 7d535a
	    fsupc    = xsup[jsupno];
kusano 7d535a
	    new_next = nextlu + (xlsub[fsupc+1]-xlsub[fsupc])*(kcol-jcol+1);
kusano 7d535a
	    nzlumax = Glu.nzlumax;
kusano 7d535a
	    while ( new_next > nzlumax ) {
kusano 7d535a
		if ( (*info = sLUMemXpand(jcol, nextlu, LUSUP, &nzlumax, &Glu)) )
kusano 7d535a
		    return;
kusano 7d535a
	    }
kusano 7d535a
    
kusano 7d535a
	    for (icol = jcol; icol<= kcol; icol++) {
kusano 7d535a
		xusub[icol+1] = nextu;
kusano 7d535a
		
kusano 7d535a
    		/* Scatter into SPA dense[*] */
kusano 7d535a
    		for (k = xa_begin[icol]; k < xa_end[icol]; k++)
kusano 7d535a
        	    dense[asub[k]] = a[k];
kusano 7d535a
kusano 7d535a
	       	/* Numeric update within the snode */
kusano 7d535a
	        ssnode_bmod(icol, jsupno, fsupc, dense, tempv, &Glu, stat);
kusano 7d535a
kusano 7d535a
		if ( (*info = spivotL(icol, diag_pivot_thresh, &usepr, perm_r,
kusano 7d535a
				      iperm_r, iperm_c, &pivrow, &Glu, stat)) )
kusano 7d535a
		    if ( iinfo == 0 ) iinfo = *info;
kusano 7d535a
		
kusano 7d535a
#ifdef DEBUG
kusano 7d535a
		sprint_lu_col("[1]: ", icol, pivrow, xprune, &Glu);
kusano 7d535a
#endif
kusano 7d535a
kusano 7d535a
	    }
kusano 7d535a
kusano 7d535a
	    jcol = icol;
kusano 7d535a
kusano 7d535a
	} else { /* Work on one panel of panel_size columns */
kusano 7d535a
	    
kusano 7d535a
	    /* Adjust panel_size so that a panel won't overlap with the next 
kusano 7d535a
	     * relaxed snode.
kusano 7d535a
	     */
kusano 7d535a
	    panel_size = w_def;
kusano 7d535a
	    for (k = jcol + 1; k < SUPERLU_MIN(jcol+panel_size, min_mn); k++) 
kusano 7d535a
		if ( relax_end[k] != EMPTY ) {
kusano 7d535a
		    panel_size = k - jcol;
kusano 7d535a
		    break;
kusano 7d535a
		}
kusano 7d535a
	    if ( k == min_mn ) panel_size = min_mn - jcol;
kusano 7d535a
	    panel_histo[panel_size]++;
kusano 7d535a
kusano 7d535a
	    /* symbolic factor on a panel of columns */
kusano 7d535a
	    spanel_dfs(m, panel_size, jcol, A, perm_r, &nseg1,
kusano 7d535a
		      dense, panel_lsub, segrep, repfnz, xprune,
kusano 7d535a
		      marker, parent, xplore, &Glu);
kusano 7d535a
	    
kusano 7d535a
	    /* numeric sup-panel updates in topological order */
kusano 7d535a
	    spanel_bmod(m, panel_size, jcol, nseg1, dense,
kusano 7d535a
		        tempv, segrep, repfnz, &Glu, stat);
kusano 7d535a
	    
kusano 7d535a
	    /* Sparse LU within the panel, and below panel diagonal */
kusano 7d535a
    	    for ( jj = jcol; jj < jcol + panel_size; jj++) {
kusano 7d535a
 		k = (jj - jcol) * m; /* column index for w-wide arrays */
kusano 7d535a
kusano 7d535a
		nseg = nseg1;	/* Begin after all the panel segments */
kusano 7d535a
kusano 7d535a
	    	if ((*info = scolumn_dfs(m, jj, perm_r, &nseg, &panel_lsub[k],
kusano 7d535a
					segrep, &repfnz[k], xprune, marker,
kusano 7d535a
					parent, xplore, &Glu)) != 0) return;
kusano 7d535a
kusano 7d535a
	      	/* Numeric updates */
kusano 7d535a
	    	if ((*info = scolumn_bmod(jj, (nseg - nseg1), &dense[k],
kusano 7d535a
					 tempv, &segrep[nseg1], &repfnz[k],
kusano 7d535a
					 jcol, &Glu, stat)) != 0) return;
kusano 7d535a
		
kusano 7d535a
	        /* Copy the U-segments to ucol[*] */
kusano 7d535a
		if ((*info = scopy_to_ucol(jj, nseg, segrep, &repfnz[k],
kusano 7d535a
					  perm_r, &dense[k], &Glu)) != 0)
kusano 7d535a
		    return;
kusano 7d535a
kusano 7d535a
	    	if ( (*info = spivotL(jj, diag_pivot_thresh, &usepr, perm_r,
kusano 7d535a
				      iperm_r, iperm_c, &pivrow, &Glu, stat)) )
kusano 7d535a
		    if ( iinfo == 0 ) iinfo = *info;
kusano 7d535a
kusano 7d535a
		/* Prune columns (0:jj-1) using column jj */
kusano 7d535a
	    	spruneL(jj, perm_r, pivrow, nseg, segrep,
kusano 7d535a
                        &repfnz[k], xprune, &Glu);
kusano 7d535a
kusano 7d535a
		/* Reset repfnz[] for this column */
kusano 7d535a
	    	resetrep_col (nseg, segrep, &repfnz[k]);
kusano 7d535a
		
kusano 7d535a
#ifdef DEBUG
kusano 7d535a
		sprint_lu_col("[2]: ", jj, pivrow, xprune, &Glu);
kusano 7d535a
#endif
kusano 7d535a
kusano 7d535a
	    }
kusano 7d535a
kusano 7d535a
   	    jcol += panel_size;	/* Move to the next panel */
kusano 7d535a
kusano 7d535a
	} /* else */
kusano 7d535a
kusano 7d535a
    } /* for */
kusano 7d535a
kusano 7d535a
    *info = iinfo;
kusano 7d535a
    
kusano 7d535a
    if ( m > n ) {
kusano 7d535a
	k = 0;
kusano 7d535a
        for (i = 0; i < m; ++i) 
kusano 7d535a
            if ( perm_r[i] == EMPTY ) {
kusano 7d535a
    		perm_r[i] = n + k;
kusano 7d535a
		++k;
kusano 7d535a
	    }
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    countnz(min_mn, xprune, &nnzL, &nnzU, &Glu);
kusano 7d535a
    fixupL(min_mn, perm_r, &Glu);
kusano 7d535a
kusano 7d535a
    sLUWorkFree(iwork, swork, &Glu); /* Free work space and compress storage */
kusano 7d535a
kusano 7d535a
    if ( fact == SamePattern_SameRowPerm ) {
kusano 7d535a
        /* L and U structures may have changed due to possibly different
kusano 7d535a
	   pivoting, even though the storage is available.
kusano 7d535a
	   There could also be memory expansions, so the array locations
kusano 7d535a
           may have changed, */
kusano 7d535a
        ((SCformat *)L->Store)->nnz = nnzL;
kusano 7d535a
	((SCformat *)L->Store)->nsuper = Glu.supno[n];
kusano 7d535a
	((SCformat *)L->Store)->nzval = Glu.lusup;
kusano 7d535a
	((SCformat *)L->Store)->nzval_colptr = Glu.xlusup;
kusano 7d535a
	((SCformat *)L->Store)->rowind = Glu.lsub;
kusano 7d535a
	((SCformat *)L->Store)->rowind_colptr = Glu.xlsub;
kusano 7d535a
	((NCformat *)U->Store)->nnz = nnzU;
kusano 7d535a
	((NCformat *)U->Store)->nzval = Glu.ucol;
kusano 7d535a
	((NCformat *)U->Store)->rowind = Glu.usub;
kusano 7d535a
	((NCformat *)U->Store)->colptr = Glu.xusub;
kusano 7d535a
    } else {
kusano 7d535a
        sCreate_SuperNode_Matrix(L, A->nrow, min_mn, nnzL, Glu.lusup, 
kusano 7d535a
	                         Glu.xlusup, Glu.lsub, Glu.xlsub, Glu.supno,
kusano 7d535a
			         Glu.xsup, SLU_SC, SLU_S, SLU_TRLU);
kusano 7d535a
    	sCreate_CompCol_Matrix(U, min_mn, min_mn, nnzU, Glu.ucol, 
kusano 7d535a
			       Glu.usub, Glu.xusub, SLU_NC, SLU_S, SLU_TRU);
kusano 7d535a
    }
kusano 7d535a
    
kusano 7d535a
    ops[FACT] += ops[TRSV] + ops[GEMV];	
kusano 7d535a
    stat->expansions = --(Glu.num_expansions);
kusano 7d535a
    
kusano 7d535a
    if ( iperm_r_allocated ) SUPERLU_FREE (iperm_r);
kusano 7d535a
    SUPERLU_FREE (iperm_c);
kusano 7d535a
    SUPERLU_FREE (relax_end);
kusano 7d535a
kusano 7d535a
}