|
|
/* |
|
* Copyright 1997, Regents of the University of Minnesota |
|
* |
|
* smbfactor.c |
|
* |
|
* This file performs the symbolic factorization of a matrix |
|
* |
|
* Started 8/1/97 |
|
* George |
|
* |
|
* $Id: smbfactor.c,v 1.2 2002/08/10 06:02:55 karypis Exp $ |
|
* |
|
*/ |
|
|
|
#include <metisbin.h> |
|
|
|
|
|
/************************************************************************* |
|
* This function sets up data structures for fill-in computations |
|
**************************************************************************/ |
|
void ComputeFillIn(GraphType *graph, idxtype *iperm) |
|
{ |
|
idxtype i, j, k, nvtxs, maxlnz, maxsub; |
|
idxtype *xadj, *adjncy; |
|
idxtype *perm, *xlnz, *xnzsub, *nzsub; |
|
double opc; |
|
|
|
/* |
|
mprintf("\nSymbolic factorization... --------------------------------------------\n"); |
|
*/ |
|
|
|
nvtxs = graph->nvtxs; |
|
xadj = graph->xadj; |
|
adjncy = graph->adjncy; |
|
|
|
maxsub = 4*xadj[nvtxs]; |
|
|
|
/* Relabel the vertices so that it starts from 1 */ |
|
k = xadj[nvtxs]; |
|
for (i=0; i<k; i++) |
|
adjncy[i]++; |
|
for (i=0; i<nvtxs+1; i++) |
|
xadj[i]++; |
|
|
|
/* Allocate the required memory */ |
|
perm = idxmalloc(nvtxs+1, "ComputeFillIn: perm"); |
|
xlnz = idxmalloc(nvtxs+1, "ComputeFillIn: xlnz"); |
|
xnzsub = idxmalloc(nvtxs+1, "ComputeFillIn: xnzsub"); |
|
nzsub = idxmalloc(maxsub, "ComputeFillIn: nzsub"); |
|
|
|
/* Construct perm from iperm and change the numbering of iperm */ |
|
for (i=0; i<nvtxs; i++) |
|
perm[iperm[i]] = i; |
|
for (i=0; i<nvtxs; i++) { |
|
iperm[i]++; |
|
perm[i]++; |
|
} |
|
|
|
/* |
|
* Call sparspak routine. |
|
*/ |
|
if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub)) { |
|
gk_free((void **)&nzsub, LTERM); |
|
|
|
maxsub = 4*maxsub; |
|
nzsub = idxmalloc(maxsub, "ComputeFillIn: nzsub"); |
|
if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub)) |
|
errexit("MAXSUB is too small!"); |
|
} |
|
|
|
opc = 0; |
|
for (i=0; i<nvtxs; i++) |
|
xlnz[i]--; |
|
for (i=0; i<nvtxs; i++) |
|
opc += (xlnz[i+1]-xlnz[i])*(xlnz[i+1]-xlnz[i]) - (xlnz[i+1]-xlnz[i]); |
|
|
|
mprintf(" Nonzeros: %D, \tOperation Count: %6.4le\n", maxlnz, opc); |
|
|
|
|
|
gk_free((void **)&perm, &xlnz, &xnzsub, &nzsub, LTERM); |
|
|
|
|
|
/* Relabel the vertices so that it starts from 0 */ |
|
for (i=0; i<nvtxs; i++) |
|
iperm[i]--; |
|
for (i=0; i<nvtxs+1; i++) |
|
xadj[i]--; |
|
k = xadj[nvtxs]; |
|
for (i=0; i<k; i++) |
|
adjncy[i]--; |
|
|
|
} |
|
|
|
|
|
|
|
/************************************************************************* |
|
* This function sets up data structures for fill-in computations |
|
**************************************************************************/ |
|
idxtype ComputeFillIn2(GraphType *graph, idxtype *iperm) |
|
{ |
|
idxtype i, j, k, nvtxs, maxlnz, maxsub; |
|
idxtype *xadj, *adjncy; |
|
idxtype *perm, *xlnz, *xnzsub, *nzsub; |
|
double opc; |
|
|
|
nvtxs = graph->nvtxs; |
|
xadj = graph->xadj; |
|
adjncy = graph->adjncy; |
|
|
|
maxsub = 4*xadj[nvtxs]; |
|
|
|
/* Relabel the vertices so that it starts from 1 */ |
|
k = xadj[nvtxs]; |
|
for (i=0; i<k; i++) |
|
adjncy[i]++; |
|
for (i=0; i<nvtxs+1; i++) |
|
xadj[i]++; |
|
|
|
/* Allocate the required memory */ |
|
perm = idxmalloc(nvtxs+1, "ComputeFillIn: perm"); |
|
xlnz = idxmalloc(nvtxs+1, "ComputeFillIn: xlnz"); |
|
xnzsub = idxmalloc(nvtxs+1, "ComputeFillIn: xnzsub"); |
|
nzsub = idxmalloc(maxsub, "ComputeFillIn: nzsub"); |
|
|
|
/* Construct perm from iperm and change the numbering of iperm */ |
|
for (i=0; i<nvtxs; i++) |
|
perm[iperm[i]] = i; |
|
for (i=0; i<nvtxs; i++) { |
|
iperm[i]++; |
|
perm[i]++; |
|
} |
|
|
|
/* |
|
* Call sparspak routine. |
|
*/ |
|
if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub)) { |
|
gk_free((void **)&nzsub, LTERM); |
|
|
|
maxsub = 4*maxsub; |
|
nzsub = idxmalloc(maxsub, "ComputeFillIn: nzsub"); |
|
if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub)) |
|
errexit("MAXSUB is too small!"); |
|
} |
|
|
|
opc = 0; |
|
for (i=0; i<nvtxs; i++) |
|
xlnz[i]--; |
|
for (i=0; i<nvtxs; i++) |
|
opc += (xlnz[i+1]-xlnz[i])*(xlnz[i+1]-xlnz[i]) - (xlnz[i+1]-xlnz[i]); |
|
|
|
|
|
gk_free((void **)&perm, &xlnz, &xnzsub, &nzsub, LTERM); |
|
|
|
|
|
/* Relabel the vertices so that it starts from 0 */ |
|
for (i=0; i<nvtxs; i++) |
|
iperm[i]--; |
|
for (i=0; i<nvtxs+1; i++) |
|
xadj[i]--; |
|
k = xadj[nvtxs]; |
|
for (i=0; i<k; i++) |
|
adjncy[i]--; |
|
|
|
return maxlnz; |
|
|
|
} |
|
|
|
|
|
/***************************************************************** |
|
********** SMBFCT ..... SYMBOLIC FACTORIZATION ********* |
|
****************************************************************** |
|
* PURPOSE - THIS ROUTINE PERFORMS SYMBOLIC FACTORIZATION |
|
* ON A PERMUTED LINEAR SYSTEM AND IT ALSO SETS UP THE |
|
* COMPRESSED DATA STRUCTURE FOR THE SYSTEM. |
|
* |
|
* INPUT PARAMETERS - |
|
* NEQNS - NUMBER OF EQUATIONS. |
|
* (XADJ, ADJNCY) - THE ADJACENCY STRUCTURE. |
|
* (PERM, INVP) - THE PERMUTATION VECTOR AND ITS INVERSE. |
|
* |
|
* UPDATED PARAMETERS - |
|
* MAXSUB - SIZE OF THE SUBSCRIPT ARRAY NZSUB. ON RETURN, |
|
* IT CONTAINS THE NUMBER OF SUBSCRIPTS USED |
|
* |
|
* OUTPUT PARAMETERS - |
|
* XLNZ - INDEX INTO THE NONZERO STORAGE VECTOR LNZ. |
|
* (XNZSUB, NZSUB) - THE COMPRESSED SUBSCRIPT VECTORS. |
|
* MAXLNZ - THE NUMBER OF NONZEROS FOUND. |
|
* |
|
*******************************************************************/ |
|
idxtype smbfct(idxtype neqns, idxtype *xadj, idxtype *adjncy, idxtype *perm, idxtype *invp, |
|
idxtype *xlnz, idxtype *maxlnz, idxtype *xnzsub, idxtype *nzsub, idxtype *maxsub) |
|
{ |
|
/* Local variables */ |
|
idxtype node, rchm, mrgk, lmax, i, j, k, m, nabor, nzbeg, nzend; |
|
idxtype kxsub, jstop, jstrt, mrkflg, inz, knz, flag; |
|
idxtype *mrglnk, *marker, *rchlnk; |
|
|
|
rchlnk = idxmalloc(neqns+1, "smbfct: rchlnk"); |
|
marker = idxsmalloc(neqns+1, 0, "smbfct: marker"); |
|
mrglnk = idxsmalloc(neqns+1, 0, "smbfct: mgrlnk"); |
|
|
|
/* Parameter adjustments */ |
|
--marker; |
|
--mrglnk; |
|
--rchlnk; |
|
--nzsub; |
|
--xnzsub; |
|
--xlnz; |
|
--invp; |
|
--perm; |
|
--adjncy; |
|
--xadj; |
|
|
|
/* Function Body */ |
|
flag = 0; |
|
nzbeg = 1; |
|
nzend = 0; |
|
xlnz[1] = 1; |
|
|
|
/* FOR EACH COLUMN KNZ COUNTS THE NUMBER OF NONZEROS IN COLUMN K ACCUMULATED IN RCHLNK. */ |
|
for (k = 1; k <= neqns; ++k) { |
|
knz = 0; |
|
mrgk = mrglnk[k]; |
|
mrkflg = 0; |
|
marker[k] = k; |
|
if (mrgk != 0) |
|
marker[k] = marker[mrgk]; |
|
xnzsub[k] = nzend; |
|
node = perm[k]; |
|
|
|
if (xadj[node] >= xadj[node+1]) { |
|
xlnz[k+1] = xlnz[k]; |
|
continue; |
|
} |
|
|
|
/* USE RCHLNK TO LINK THROUGH THE STRUCTURE OF A(*,K) BELOW DIAGONAL */ |
|
rchlnk[k] = neqns+1; |
|
for (j=xadj[node]; j<xadj[node+1]; j++) { |
|
nabor = invp[adjncy[j]]; |
|
if (nabor <= k) |
|
continue; |
|
rchm = k; |
|
|
|
do { |
|
m = rchm; |
|
rchm = rchlnk[m]; |
|
} while (rchm <= nabor); |
|
|
|
knz++; |
|
rchlnk[m] = nabor; |
|
rchlnk[nabor] = rchm; |
|
if (marker[nabor] != marker[k]) |
|
mrkflg = 1; |
|
} |
|
|
|
/* TEST FOR MASS SYMBOLIC ELIMINATION */ |
|
lmax = 0; |
|
if (mrkflg != 0 || mrgk == 0 || mrglnk[mrgk] != 0) |
|
goto L350; |
|
xnzsub[k] = xnzsub[mrgk] + 1; |
|
knz = xlnz[mrgk + 1] - (xlnz[mrgk] + 1); |
|
goto L1400; |
|
|
|
|
|
/* LINK THROUGH EACH COLUMN I THAT AFFECTS L(*,K) */ |
|
L350: |
|
i = k; |
|
while ((i = mrglnk[i]) != 0) { |
|
inz = xlnz[i+1] - (xlnz[i]+1); |
|
jstrt = xnzsub[i] + 1; |
|
jstop = xnzsub[i] + inz; |
|
|
|
if (inz > lmax) { |
|
lmax = inz; |
|
xnzsub[k] = jstrt; |
|
} |
|
|
|
/* MERGE STRUCTURE OF L(*,I) IN NZSUB INTO RCHLNK. */ |
|
rchm = k; |
|
for (j = jstrt; j <= jstop; ++j) { |
|
nabor = nzsub[j]; |
|
do { |
|
m = rchm; |
|
rchm = rchlnk[m]; |
|
} while (rchm < nabor); |
|
|
|
if (rchm != nabor) { |
|
knz++; |
|
rchlnk[m] = nabor; |
|
rchlnk[nabor] = rchm; |
|
rchm = nabor; |
|
} |
|
} |
|
} |
|
|
|
/* CHECK IF SUBSCRIPTS DUPLICATE THOSE OF ANOTHER COLUMN */ |
|
if (knz == lmax) |
|
goto L1400; |
|
|
|
/* OR IF TAIL OF K-1ST COLUMN MATCHES HEAD OF KTH */ |
|
if (nzbeg > nzend) |
|
goto L1200; |
|
|
|
i = rchlnk[k]; |
|
for (jstrt = nzbeg; jstrt <= nzend; ++jstrt) { |
|
if (nzsub[jstrt] < i) |
|
continue; |
|
|
|
if (nzsub[jstrt] == i) |
|
goto L1000; |
|
else |
|
goto L1200; |
|
} |
|
goto L1200; |
|
|
|
L1000: |
|
xnzsub[k] = jstrt; |
|
for (j = jstrt; j <= nzend; ++j) { |
|
if (nzsub[j] != i) |
|
goto L1200; |
|
|
|
i = rchlnk[i]; |
|
if (i > neqns) |
|
goto L1400; |
|
} |
|
nzend = jstrt - 1; |
|
|
|
/* COPY THE STRUCTURE OF L(*,K) FROM RCHLNK TO THE DATA STRUCTURE (XNZSUB, NZSUB) */ |
|
L1200: |
|
nzbeg = nzend + 1; |
|
nzend += knz; |
|
|
|
if (nzend > *maxsub) { |
|
flag = 1; /* Out of memory */ |
|
break; |
|
} |
|
|
|
i = k; |
|
for (j=nzbeg; j<=nzend; ++j) { |
|
i = rchlnk[i]; |
|
nzsub[j] = i; |
|
marker[i] = k; |
|
} |
|
xnzsub[k] = nzbeg; |
|
marker[k] = k; |
|
|
|
/* |
|
* UPDATE THE VECTOR MRGLNK. NOTE COLUMN L(*,K) JUST FOUND |
|
* IS REQUIRED TO DETERMINE COLUMN L(*,J), WHERE |
|
* L(J,K) IS THE FIRST NONZERO IN L(*,K) BELOW DIAGONAL. |
|
*/ |
|
L1400: |
|
if (knz > 1) { |
|
kxsub = xnzsub[k]; |
|
i = nzsub[kxsub]; |
|
mrglnk[k] = mrglnk[i]; |
|
mrglnk[i] = k; |
|
} |
|
|
|
xlnz[k + 1] = xlnz[k] + knz; |
|
} |
|
|
|
if (flag == 0) { |
|
*maxlnz = xlnz[neqns] - 1; |
|
*maxsub = xnzsub[neqns]; |
|
xnzsub[neqns + 1] = xnzsub[neqns]; |
|
} |
|
|
|
marker++; |
|
mrglnk++; |
|
rchlnk++; |
|
nzsub++; |
|
xnzsub++; |
|
xlnz++; |
|
invp++; |
|
perm++; |
|
adjncy++; |
|
xadj++; |
|
gk_free((void **)&rchlnk, &mrglnk, &marker, LTERM); |
|
|
|
return flag; |
|
|
|
} |
|
|
|