|
Blender
V2.59
|
00001 00006 /* 00007 * -- SuperLU routine (version 2.0) -- 00008 * Univ. of California Berkeley, Xerox Palo Alto Research Center, 00009 * and Lawrence Berkeley National Lab. 00010 * November 15, 1997 00011 * 00012 */ 00013 /* 00014 Copyright (c) 1994 by Xerox Corporation. All rights reserved. 00015 00016 THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY 00017 EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. 00018 00019 Permission is hereby granted to use or copy this program for any 00020 purpose, provided the above notices are retained on all copies. 00021 Permission to modify the code and to distribute modified code is 00022 granted, provided the above notices are retained, and a notice that 00023 the code was modified is included with the above copyright notice. 00024 */ 00025 00026 #include "ssp_defs.h" 00027 #include "util.h" 00028 00029 void 00030 spanel_dfs ( 00031 const int m, /* in - number of rows in the matrix */ 00032 const int w, /* in */ 00033 const int jcol, /* in */ 00034 SuperMatrix *A, /* in - original matrix */ 00035 int *perm_r, /* in */ 00036 int *nseg, /* out */ 00037 float *dense, /* out */ 00038 int *panel_lsub, /* out */ 00039 int *segrep, /* out */ 00040 int *repfnz, /* out */ 00041 int *xprune, /* out */ 00042 int *marker, /* out */ 00043 int *parent, /* working array */ 00044 int *xplore, /* working array */ 00045 GlobalLU_t *Glu /* modified */ 00046 ) 00047 { 00048 /* 00049 * Purpose 00050 * ======= 00051 * 00052 * Performs a symbolic factorization on a panel of columns [jcol, jcol+w). 00053 * 00054 * A supernode representative is the last column of a supernode. 00055 * The nonzeros in U[*,j] are segments that end at supernodal 00056 * representatives. 00057 * 00058 * The routine returns one list of the supernodal representatives 00059 * in topological order of the dfs that generates them. This list is 00060 * a superset of the topological order of each individual column within 00061 * the panel. 00062 * The location of the first nonzero in each supernodal segment 00063 * (supernodal entry location) is also returned. Each column has a 00064 * separate list for this purpose. 00065 * 00066 * Two marker arrays are used for dfs: 00067 * marker[i] == jj, if i was visited during dfs of current column jj; 00068 * marker1[i] >= jcol, if i was visited by earlier columns in this panel; 00069 * 00070 * marker: A-row --> A-row/col (0/1) 00071 * repfnz: SuperA-col --> PA-row 00072 * parent: SuperA-col --> SuperA-col 00073 * xplore: SuperA-col --> index to L-structure 00074 * 00075 */ 00076 NCPformat *Astore; 00077 float *a; 00078 int *asub; 00079 int *xa_begin, *xa_end; 00080 int krep, chperm, chmark, chrep, oldrep, kchild, myfnz; 00081 int k, krow, kmark, kperm; 00082 int xdfs, maxdfs, kpar; 00083 int jj; /* index through each column in the panel */ 00084 int *marker1; /* marker1[jj] >= jcol if vertex jj was visited 00085 by a previous column within this panel. */ 00086 int *repfnz_col; /* start of each column in the panel */ 00087 float *dense_col; /* start of each column in the panel */ 00088 int nextl_col; /* next available position in panel_lsub[*,jj] */ 00089 int *xsup, *supno; 00090 int *lsub, *xlsub; 00091 00092 /* Initialize pointers */ 00093 Astore = A->Store; 00094 a = Astore->nzval; 00095 asub = Astore->rowind; 00096 xa_begin = Astore->colbeg; 00097 xa_end = Astore->colend; 00098 marker1 = marker + m; 00099 repfnz_col = repfnz; 00100 dense_col = dense; 00101 *nseg = 0; 00102 xsup = Glu->xsup; 00103 supno = Glu->supno; 00104 lsub = Glu->lsub; 00105 xlsub = Glu->xlsub; 00106 00107 /* For each column in the panel */ 00108 for (jj = jcol; jj < jcol + w; jj++) { 00109 nextl_col = (jj - jcol) * m; 00110 00111 #ifdef CHK_DFS 00112 printf("\npanel col %d: ", jj); 00113 #endif 00114 00115 /* For each nonz in A[*,jj] do dfs */ 00116 for (k = xa_begin[jj]; k < xa_end[jj]; k++) { 00117 krow = asub[k]; 00118 dense_col[krow] = a[k]; 00119 kmark = marker[krow]; 00120 if ( kmark == jj ) 00121 continue; /* krow visited before, go to the next nonzero */ 00122 00123 /* For each unmarked nbr krow of jj 00124 * krow is in L: place it in structure of L[*,jj] 00125 */ 00126 marker[krow] = jj; 00127 kperm = perm_r[krow]; 00128 00129 if ( kperm == EMPTY ) { 00130 panel_lsub[nextl_col++] = krow; /* krow is indexed into A */ 00131 } 00132 /* 00133 * krow is in U: if its supernode-rep krep 00134 * has been explored, update repfnz[*] 00135 */ 00136 else { 00137 00138 krep = xsup[supno[kperm]+1] - 1; 00139 myfnz = repfnz_col[krep]; 00140 00141 #ifdef CHK_DFS 00142 printf("krep %d, myfnz %d, perm_r[%d] %d\n", krep, myfnz, krow, kperm); 00143 #endif 00144 if ( myfnz != EMPTY ) { /* Representative visited before */ 00145 if ( myfnz > kperm ) repfnz_col[krep] = kperm; 00146 /* continue; */ 00147 } 00148 else { 00149 /* Otherwise, perform dfs starting at krep */ 00150 oldrep = EMPTY; 00151 parent[krep] = oldrep; 00152 repfnz_col[krep] = kperm; 00153 xdfs = xlsub[krep]; 00154 maxdfs = xprune[krep]; 00155 00156 #ifdef CHK_DFS 00157 printf(" xdfs %d, maxdfs %d: ", xdfs, maxdfs); 00158 for (i = xdfs; i < maxdfs; i++) printf(" %d", lsub[i]); 00159 printf("\n"); 00160 #endif 00161 do { 00162 /* 00163 * For each unmarked kchild of krep 00164 */ 00165 while ( xdfs < maxdfs ) { 00166 00167 kchild = lsub[xdfs]; 00168 xdfs++; 00169 chmark = marker[kchild]; 00170 00171 if ( chmark != jj ) { /* Not reached yet */ 00172 marker[kchild] = jj; 00173 chperm = perm_r[kchild]; 00174 00175 /* Case kchild is in L: place it in L[*,j] */ 00176 if ( chperm == EMPTY ) { 00177 panel_lsub[nextl_col++] = kchild; 00178 } 00179 /* Case kchild is in U: 00180 * chrep = its supernode-rep. If its rep has 00181 * been explored, update its repfnz[*] 00182 */ 00183 else { 00184 00185 chrep = xsup[supno[chperm]+1] - 1; 00186 myfnz = repfnz_col[chrep]; 00187 #ifdef CHK_DFS 00188 printf("chrep %d,myfnz %d,perm_r[%d] %d\n",chrep,myfnz,kchild,chperm); 00189 #endif 00190 if ( myfnz != EMPTY ) { /* Visited before */ 00191 if ( myfnz > chperm ) 00192 repfnz_col[chrep] = chperm; 00193 } 00194 else { 00195 /* Cont. dfs at snode-rep of kchild */ 00196 xplore[krep] = xdfs; 00197 oldrep = krep; 00198 krep = chrep; /* Go deeper down G(L) */ 00199 parent[krep] = oldrep; 00200 repfnz_col[krep] = chperm; 00201 xdfs = xlsub[krep]; 00202 maxdfs = xprune[krep]; 00203 #ifdef CHK_DFS 00204 printf(" xdfs %d, maxdfs %d: ", xdfs, maxdfs); 00205 for (i = xdfs; i < maxdfs; i++) printf(" %d", lsub[i]); 00206 printf("\n"); 00207 #endif 00208 } /* else */ 00209 00210 } /* else */ 00211 00212 } /* if... */ 00213 00214 } /* while xdfs < maxdfs */ 00215 00216 /* krow has no more unexplored nbrs: 00217 * Place snode-rep krep in postorder DFS, if this 00218 * segment is seen for the first time. (Note that 00219 * "repfnz[krep]" may change later.) 00220 * Backtrack dfs to its parent. 00221 */ 00222 if ( marker1[krep] < jcol ) { 00223 segrep[*nseg] = krep; 00224 ++(*nseg); 00225 marker1[krep] = jj; 00226 } 00227 00228 kpar = parent[krep]; /* Pop stack, mimic recursion */ 00229 if ( kpar == EMPTY ) break; /* dfs done */ 00230 krep = kpar; 00231 xdfs = xplore[krep]; 00232 maxdfs = xprune[krep]; 00233 00234 #ifdef CHK_DFS 00235 printf(" pop stack: krep %d,xdfs %d,maxdfs %d: ", krep,xdfs,maxdfs); 00236 for (i = xdfs; i < maxdfs; i++) printf(" %d", lsub[i]); 00237 printf("\n"); 00238 #endif 00239 } while ( kpar != EMPTY ); /* do-while - until empty stack */ 00240 00241 } /* else */ 00242 00243 } /* else */ 00244 00245 } /* for each nonz in A[*,jj] */ 00246 00247 repfnz_col += m; /* Move to next column */ 00248 dense_col += m; 00249 00250 } /* for jj ... */ 00251 00252 }