cgstrf.c

LU矩阵分解单机版最新版本

    int       *panel_histo = stat->panel_histo;
    flops_t   *ops = stat->ops;

    iinfo    = 0;
    m        = A->nrow;
    n        = A->ncol;
    min_mn   = SUPERLU_MIN(m, n);
    Astore   = A->Store;
    a        = Astore->nzval;
    asub     = Astore->rowind;
    xa_begin = Astore->colbeg;
    xa_end   = Astore->colend;

    /* Allocate storage common to the factor routines */
    *info = cLUMemInit(fact, work, lwork, m, n, Astore->nnz,
                       panel_size, L, U, &Glu, &iwork, &cwork);
    if ( *info ) return;
    
    xsup    = Glu.xsup;
    supno   = Glu.supno;
    xlsub   = Glu.xlsub;
    xlusup  = Glu.xlusup;
    xusub   = Glu.xusub;
    
    SetIWork(m, n, panel_size, iwork, &segrep, &parent, &xplore,
	     &repfnz, &panel_lsub, &xprune, &marker);
    cSetRWork(m, panel_size, cwork, &dense, &tempv);
    
    usepr = (fact == SamePattern_SameRowPerm);
    if ( usepr ) {
	/* Compute the inverse of perm_r */
	iperm_r = (int *) intMalloc(m);
	for (k = 0; k < m; ++k) iperm_r[perm_r[k]] = k;
	iperm_r_allocated = 1;
    }
    iperm_c = (int *) intMalloc(n);
    for (k = 0; k < n; ++k) iperm_c[perm_c[k]] = k;

    /* Identify relaxed snodes */
    relax_end = (int *) intMalloc(n);
    if ( options->SymmetricMode == YES ) {
        heap_relax_snode(n, etree, relax, marker, relax_end); 
    } else {
        relax_snode(n, etree, relax, marker, relax_end); 
    }
    
    ifill (perm_r, m, EMPTY);
    ifill (marker, m * NO_MARKER, EMPTY);
    supno[0] = -1;
    xsup[0]  = xlsub[0] = xusub[0] = xlusup[0] = 0;
    w_def    = panel_size;

    /* 
     * Work on one "panel" at a time. A panel is one of the following: 
     *	   (a) a relaxed supernode at the bottom of the etree, or
     *	   (b) panel_size contiguous columns, defined by the user
     */
    for (jcol = 0; jcol < min_mn; ) {

	if ( relax_end[jcol] != EMPTY ) { /* start of a relaxed snode */
   	    kcol = relax_end[jcol];	  /* end of the relaxed snode */
	    panel_histo[kcol-jcol+1]++;

	    /* --------------------------------------
	     * Factorize the relaxed supernode(jcol:kcol) 
	     * -------------------------------------- */
	    /* Determine the union of the row structure of the snode */
	    if ( (*info = csnode_dfs(jcol, kcol, asub, xa_begin, xa_end,
				    xprune, marker, &Glu)) != 0 )
		return;

            nextu    = xusub[jcol];
	    nextlu   = xlusup[jcol];
	    jsupno   = supno[jcol];
	    fsupc    = xsup[jsupno];
	    new_next = nextlu + (xlsub[fsupc+1]-xlsub[fsupc])*(kcol-jcol+1);
	    nzlumax = Glu.nzlumax;
	    while ( new_next > nzlumax ) {
		if ( (*info = cLUMemXpand(jcol, nextlu, LUSUP, &nzlumax, &Glu)) )
		    return;
	    }
    
	    for (icol = jcol; icol<= kcol; icol++) {
		xusub[icol+1] = nextu;
		
    		/* Scatter into SPA dense[*] */
    		for (k = xa_begin[icol]; k < xa_end[icol]; k++)
        	    dense[asub[k]] = a[k];

	       	/* Numeric update within the snode */
	        csnode_bmod(icol, jsupno, fsupc, dense, tempv, &Glu, stat);

		if ( (*info = cpivotL(icol, diag_pivot_thresh, &usepr, perm_r,
				      iperm_r, iperm_c, &pivrow, &Glu, stat)) )
		    if ( iinfo == 0 ) iinfo = *info;
		
#ifdef DEBUG
		cprint_lu_col("[1]: ", icol, pivrow, xprune, &Glu);
#endif

	    }

	    jcol = icol;

	} else { /* Work on one panel of panel_size columns */
	    
	    /* Adjust panel_size so that a panel won't overlap with the next 
	     * relaxed snode.
	     */
	    panel_size = w_def;
	    for (k = jcol + 1; k < SUPERLU_MIN(jcol+panel_size, min_mn); k++) 
		if ( relax_end[k] != EMPTY ) {
		    panel_size = k - jcol;
		    break;
		}
	    if ( k == min_mn ) panel_size = min_mn - jcol;
	    panel_histo[panel_size]++;

	    /* symbolic factor on a panel of columns */
	    cpanel_dfs(m, panel_size, jcol, A, perm_r, &nseg1,
		      dense, panel_lsub, segrep, repfnz, xprune,
		      marker, parent, xplore, &Glu);
	    
	    /* numeric sup-panel updates in topological order */
	    cpanel_bmod(m, panel_size, jcol, nseg1, dense,
		        tempv, segrep, repfnz, &Glu, stat);
	    
	    /* Sparse LU within the panel, and below panel diagonal */
    	    for ( jj = jcol; jj < jcol + panel_size; jj++) {
 		k = (jj - jcol) * m; /* column index for w-wide arrays */

		nseg = nseg1;	/* Begin after all the panel segments */

	    	if ((*info = ccolumn_dfs(m, jj, perm_r, &nseg, &panel_lsub[k],
					segrep, &repfnz[k], xprune, marker,
					parent, xplore, &Glu)) != 0) return;

	      	/* Numeric updates */
	    	if ((*info = ccolumn_bmod(jj, (nseg - nseg1), &dense[k],
					 tempv, &segrep[nseg1], &repfnz[k],
					 jcol, &Glu, stat)) != 0) return;
		
	        /* Copy the U-segments to ucol[*] */
		if ((*info = ccopy_to_ucol(jj, nseg, segrep, &repfnz[k],
					  perm_r, &dense[k], &Glu)) != 0)
		    return;

	    	if ( (*info = cpivotL(jj, diag_pivot_thresh, &usepr, perm_r,
				      iperm_r, iperm_c, &pivrow, &Glu, stat)) )
		    if ( iinfo == 0 ) iinfo = *info;

		/* Prune columns (0:jj-1) using column jj */
	    	cpruneL(jj, perm_r, pivrow, nseg, segrep,
                        &repfnz[k], xprune, &Glu);

		/* Reset repfnz[] for this column */
	    	resetrep_col (nseg, segrep, &repfnz[k]);
		
#ifdef DEBUG
		cprint_lu_col("[2]: ", jj, pivrow, xprune, &Glu);
#endif

	    }

   	    jcol += panel_size;	/* Move to the next panel */

	} /* else */

    } /* for */

    *info = iinfo;
    
    if ( m > n ) {
	k = 0;
        for (i = 0; i < m; ++i) 
            if ( perm_r[i] == EMPTY ) {
    		perm_r[i] = n + k;
		++k;
	    }
    }

    countnz(min_mn, xprune, &nnzL, &nnzU, &Glu);
    fixupL(min_mn, perm_r, &Glu);

    cLUWorkFree(iwork, cwork, &Glu); /* Free work space and compress storage */

    if ( fact == SamePattern_SameRowPerm ) {
        /* L and U structures may have changed due to possibly different
	   pivoting, even though the storage is available.
	   There could also be memory expansions, so the array locations
           may have changed, */
        ((SCformat *)L->Store)->nnz = nnzL;
	((SCformat *)L->Store)->nsuper = Glu.supno[n];
	((SCformat *)L->Store)->nzval = Glu.lusup;
	((SCformat *)L->Store)->nzval_colptr = Glu.xlusup;
	((SCformat *)L->Store)->rowind = Glu.lsub;
	((SCformat *)L->Store)->rowind_colptr = Glu.xlsub;
	((NCformat *)U->Store)->nnz = nnzU;
	((NCformat *)U->Store)->nzval = Glu.ucol;
	((NCformat *)U->Store)->rowind = Glu.usub;
	((NCformat *)U->Store)->colptr = Glu.xusub;
    } else {
        cCreate_SuperNode_Matrix(L, A->nrow, min_mn, nnzL, Glu.lusup, 
	                         Glu.xlusup, Glu.lsub, Glu.xlsub, Glu.supno,
			         Glu.xsup, SLU_SC, SLU_C, SLU_TRLU);
    	cCreate_CompCol_Matrix(U, min_mn, min_mn, nnzU, Glu.ucol, 
			       Glu.usub, Glu.xusub, SLU_NC, SLU_C, SLU_TRU);
    }
    
    ops[FACT] += ops[TRSV] + ops[GEMV];	
    
    if ( iperm_r_allocated ) SUPERLU_FREE (iperm_r);
    SUPERLU_FREE (iperm_c);
    SUPERLU_FREE (relax_end);

}