kusano 7d535a
kusano 7d535a
/*! @file spanel_dfs.c
kusano 7d535a
 * \brief Peforms a symbolic factorization on a panel of symbols
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 * -- SuperLU routine (version 2.0) --
kusano 7d535a
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
kusano 7d535a
 * and Lawrence Berkeley National Lab.
kusano 7d535a
 * November 15, 1997
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
 *   Performs a symbolic factorization on a panel of columns [jcol, jcol+w).
kusano 7d535a
 *
kusano 7d535a
 *   A supernode representative is the last column of a supernode.
kusano 7d535a
 *   The nonzeros in U[*,j] are segments that end at supernodal
kusano 7d535a
 *   representatives.
kusano 7d535a
 *
kusano 7d535a
 *   The routine returns one list of the supernodal representatives
kusano 7d535a
 *   in topological order of the dfs that generates them. This list is
kusano 7d535a
 *   a superset of the topological order of each individual column within
kusano 7d535a
 *   the panel. 
kusano 7d535a
 *   The location of the first nonzero in each supernodal segment
kusano 7d535a
 *   (supernodal entry location) is also returned. Each column has a 
kusano 7d535a
 *   separate list for this purpose.
kusano 7d535a
 *
kusano 7d535a
 *   Two marker arrays are used for dfs:
kusano 7d535a
 *     marker[i] == jj, if i was visited during dfs of current column jj;
kusano 7d535a
 *     marker1[i] >= jcol, if i was visited by earlier columns in this panel;
kusano 7d535a
 *
kusano 7d535a
 *   marker: A-row --> A-row/col (0/1)
kusano 7d535a
 *   repfnz: SuperA-col --> PA-row
kusano 7d535a
 *   parent: SuperA-col --> SuperA-col
kusano 7d535a
 *   xplore: SuperA-col --> index to L-structure
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
kusano 7d535a
void
kusano 7d535a
spanel_dfs (
kusano 7d535a
	   const int  m,           /* in - number of rows in the matrix */
kusano 7d535a
	   const int  w,           /* in */
kusano 7d535a
	   const int  jcol,        /* in */
kusano 7d535a
	   SuperMatrix *A,       /* in - original matrix */
kusano 7d535a
	   int        *perm_r,     /* in */
kusano 7d535a
	   int        *nseg,	   /* out */
kusano 7d535a
	   float     *dense,      /* out */
kusano 7d535a
	   int        *panel_lsub, /* out */
kusano 7d535a
	   int        *segrep,     /* out */
kusano 7d535a
	   int        *repfnz,     /* out */
kusano 7d535a
	   int        *xprune,     /* out */
kusano 7d535a
	   int        *marker,     /* out */     
kusano 7d535a
	   int        *parent,     /* working array */
kusano 7d535a
	   int        *xplore,     /* working array */
kusano 7d535a
	   GlobalLU_t *Glu         /* modified */
kusano 7d535a
	   )
kusano 7d535a
{
kusano 7d535a
kusano 7d535a
    NCPformat *Astore;
kusano 7d535a
    float    *a;
kusano 7d535a
    int       *asub;
kusano 7d535a
    int       *xa_begin, *xa_end;
kusano 7d535a
    int	      krep, chperm, chmark, chrep, oldrep, kchild, myfnz;
kusano 7d535a
    int       k, krow, kmark, kperm;
kusano 7d535a
    int       xdfs, maxdfs, kpar;
kusano 7d535a
    int       jj;	   /* index through each column in the panel */
kusano 7d535a
    int       *marker1;	   /* marker1[jj] >= jcol if vertex jj was visited 
kusano 7d535a
			      by a previous column within this panel.   */
kusano 7d535a
    int       *repfnz_col; /* start of each column in the panel */
kusano 7d535a
    float    *dense_col;  /* start of each column in the panel */
kusano 7d535a
    int       nextl_col;   /* next available position in panel_lsub[*,jj] */
kusano 7d535a
    int       *xsup, *supno;
kusano 7d535a
    int       *lsub, *xlsub;
kusano 7d535a
kusano 7d535a
    /* Initialize pointers */
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
    marker1    = marker + m;
kusano 7d535a
    repfnz_col = repfnz;
kusano 7d535a
    dense_col  = dense;
kusano 7d535a
    *nseg      = 0;
kusano 7d535a
    xsup       = Glu->xsup;
kusano 7d535a
    supno      = Glu->supno;
kusano 7d535a
    lsub       = Glu->lsub;
kusano 7d535a
    xlsub      = Glu->xlsub;
kusano 7d535a
kusano 7d535a
    /* For each column in the panel */
kusano 7d535a
    for (jj = jcol; jj < jcol + w; jj++) {
kusano 7d535a
	nextl_col = (jj - jcol) * m;
kusano 7d535a
kusano 7d535a
#ifdef CHK_DFS
kusano 7d535a
	printf("\npanel col %d: ", jj);
kusano 7d535a
#endif
kusano 7d535a
kusano 7d535a
	/* For each nonz in A[*,jj] do dfs */
kusano 7d535a
	for (k = xa_begin[jj]; k < xa_end[jj]; k++) {
kusano 7d535a
	    krow = asub[k];
kusano 7d535a
            dense_col[krow] = a[k];
kusano 7d535a
	    kmark = marker[krow];    	
kusano 7d535a
	    if ( kmark == jj ) 
kusano 7d535a
		continue;     /* krow visited before, go to the next nonzero */
kusano 7d535a
kusano 7d535a
	    /* For each unmarked nbr krow of jj
kusano 7d535a
	     * krow is in L: place it in structure of L[*,jj]
kusano 7d535a
	     */
kusano 7d535a
	    marker[krow] = jj;
kusano 7d535a
	    kperm = perm_r[krow];
kusano 7d535a
	    
kusano 7d535a
	    if ( kperm == EMPTY ) {
kusano 7d535a
		panel_lsub[nextl_col++] = krow; /* krow is indexed into A */
kusano 7d535a
	    }
kusano 7d535a
	    /* 
kusano 7d535a
	     * krow is in U: if its supernode-rep krep
kusano 7d535a
	     * has been explored, update repfnz[*]
kusano 7d535a
	     */
kusano 7d535a
	    else {
kusano 7d535a
		
kusano 7d535a
		krep = xsup[supno[kperm]+1] - 1;
kusano 7d535a
		myfnz = repfnz_col[krep];
kusano 7d535a
		
kusano 7d535a
#ifdef CHK_DFS
kusano 7d535a
		printf("krep %d, myfnz %d, perm_r[%d] %d\n", krep, myfnz, krow, kperm);
kusano 7d535a
#endif
kusano 7d535a
		if ( myfnz != EMPTY ) {	/* Representative visited before */
kusano 7d535a
		    if ( myfnz > kperm ) repfnz_col[krep] = kperm;
kusano 7d535a
		    /* continue; */
kusano 7d535a
		}
kusano 7d535a
		else {
kusano 7d535a
		    /* Otherwise, perform dfs starting at krep */
kusano 7d535a
		    oldrep = EMPTY;
kusano 7d535a
		    parent[krep] = oldrep;
kusano 7d535a
		    repfnz_col[krep] = kperm;
kusano 7d535a
		    xdfs = xlsub[krep];
kusano 7d535a
		    maxdfs = xprune[krep];
kusano 7d535a
		    
kusano 7d535a
#ifdef CHK_DFS 
kusano 7d535a
		    printf("  xdfs %d, maxdfs %d: ", xdfs, maxdfs);
kusano 7d535a
		    for (i = xdfs; i < maxdfs; i++) printf(" %d", lsub[i]);
kusano 7d535a
		    printf("\n");
kusano 7d535a
#endif
kusano 7d535a
		    do {
kusano 7d535a
			/* 
kusano 7d535a
			 * For each unmarked kchild of krep 
kusano 7d535a
			 */
kusano 7d535a
			while ( xdfs < maxdfs ) {
kusano 7d535a
			    
kusano 7d535a
			    kchild = lsub[xdfs];
kusano 7d535a
			    xdfs++;
kusano 7d535a
			    chmark = marker[kchild];
kusano 7d535a
			    
kusano 7d535a
			    if ( chmark != jj ) { /* Not reached yet */
kusano 7d535a
				marker[kchild] = jj;
kusano 7d535a
				chperm = perm_r[kchild];
kusano 7d535a
			      
kusano 7d535a
				/* Case kchild is in L: place it in L[*,j] */
kusano 7d535a
				if ( chperm == EMPTY ) {
kusano 7d535a
				    panel_lsub[nextl_col++] = kchild;
kusano 7d535a
				} 
kusano 7d535a
				/* Case kchild is in U: 
kusano 7d535a
				 *   chrep = its supernode-rep. If its rep has 
kusano 7d535a
				 *   been explored, update its repfnz[*]
kusano 7d535a
				 */
kusano 7d535a
				else {
kusano 7d535a
				    
kusano 7d535a
				    chrep = xsup[supno[chperm]+1] - 1;
kusano 7d535a
				    myfnz = repfnz_col[chrep];
kusano 7d535a
#ifdef CHK_DFS
kusano 7d535a
				    printf("chrep %d,myfnz %d,perm_r[%d] %d\n",chrep,myfnz,kchild,chperm);
kusano 7d535a
#endif
kusano 7d535a
				    if ( myfnz != EMPTY ) { /* Visited before */
kusano 7d535a
					if ( myfnz > chperm )
kusano 7d535a
					    repfnz_col[chrep] = chperm;
kusano 7d535a
				    }
kusano 7d535a
				    else {
kusano 7d535a
					/* Cont. dfs at snode-rep of kchild */
kusano 7d535a
					xplore[krep] = xdfs;	
kusano 7d535a
					oldrep = krep;
kusano 7d535a
					krep = chrep; /* Go deeper down G(L) */
kusano 7d535a
					parent[krep] = oldrep;
kusano 7d535a
					repfnz_col[krep] = chperm;
kusano 7d535a
					xdfs = xlsub[krep];     
kusano 7d535a
					maxdfs = xprune[krep];
kusano 7d535a
#ifdef CHK_DFS 
kusano 7d535a
					printf("  xdfs %d, maxdfs %d: ", xdfs, maxdfs);
kusano 7d535a
					for (i = xdfs; i < maxdfs; i++) printf(" %d", lsub[i]);	
kusano 7d535a
					printf("\n");
kusano 7d535a
#endif
kusano 7d535a
				    } /* else */
kusano 7d535a
				  
kusano 7d535a
				} /* else */
kusano 7d535a
			      
kusano 7d535a
			    } /* if... */
kusano 7d535a
			    
kusano 7d535a
			} /* while xdfs < maxdfs */
kusano 7d535a
			
kusano 7d535a
			/* krow has no more unexplored nbrs:
kusano 7d535a
			 *    Place snode-rep krep in postorder DFS, if this 
kusano 7d535a
			 *    segment is seen for the first time. (Note that
kusano 7d535a
			 *    "repfnz[krep]" may change later.)
kusano 7d535a
			 *    Backtrack dfs to its parent.
kusano 7d535a
			 */
kusano 7d535a
			if ( marker1[krep] < jcol ) {
kusano 7d535a
			    segrep[*nseg] = krep;
kusano 7d535a
			    ++(*nseg);
kusano 7d535a
			    marker1[krep] = jj;
kusano 7d535a
			}
kusano 7d535a
			
kusano 7d535a
			kpar = parent[krep]; /* Pop stack, mimic recursion */
kusano 7d535a
			if ( kpar == EMPTY ) break; /* dfs done */
kusano 7d535a
			krep = kpar;
kusano 7d535a
			xdfs = xplore[krep];
kusano 7d535a
			maxdfs = xprune[krep];
kusano 7d535a
			
kusano 7d535a
#ifdef CHK_DFS 
kusano 7d535a
			printf("  pop stack: krep %d,xdfs %d,maxdfs %d: ", krep,xdfs,maxdfs);
kusano 7d535a
			for (i = xdfs; i < maxdfs; i++) printf(" %d", lsub[i]);
kusano 7d535a
			printf("\n");
kusano 7d535a
#endif
kusano 7d535a
		    } while ( kpar != EMPTY ); /* do-while - until empty stack */
kusano 7d535a
		    
kusano 7d535a
		} /* else */
kusano 7d535a
		
kusano 7d535a
	    } /* else */
kusano 7d535a
	    
kusano 7d535a
	} /* for each nonz in A[*,jj] */
kusano 7d535a
	
kusano 7d535a
	repfnz_col += m;    /* Move to next column */
kusano 7d535a
        dense_col += m;
kusano 7d535a
	
kusano 7d535a
    } /* for jj ... */
kusano 7d535a
    
kusano 7d535a
}