/*
* 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.1 1998/11/27 17:59:40 karypis Exp $
*
*/
#include <metis.h>
/*************************************************************************
* This function sets up data structures for fill-in computations
**************************************************************************/
void ComputeFillIn(GraphType *graph, idxtype *iperm)
{
long i, j, k, nvtxs, maxlnz, maxsub;
idxtype *xadj, *adjncy;
idxtype *perm, *xlnz, *xnzsub, *nzsub;
double opc;
/*
printf("\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)) {
free(nzsub);
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!");
}
GKfree(&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)
{
long 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)) {
free(nzsub);
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]);
GKfree(&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.
*
*******************************************************************/
long smbfct(long neqns, idxtype *xadj, idxtype *adjncy, idxtype *perm, idxtype *invp, idxtype *xlnz,
long *maxlnz, idxtype *xnzsub, idxtype *nzsub, long *maxsub)
{
/* Local variables */
long node, rchm, mrgk, lmax, i, j, k, m, nabor, nzbeg, nzend;
long kxsub, jstop, jstrt, mrkflg, inz, knz, flag,lgind;
idxtype *mrglnk, *marker, *rchlnk;
/* jyb */
long snode, lnode, xnzi, xnzip1,nbi, nbip1;
long aux, minsn, maxsn;
double opc;
long taille1, k1, taille2, k2, taille3, k3,nd;
idxtype *parent, *nliens, *supnd, *tbnode, *parsnode;
FILE *fpout;
/* extern void ecri11_; */
/* jyb fin */
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];
}
parent = idxsmalloc(neqns, 0, "smbfct: parent");
parent--;
/* Calculate the elimination tree */
for (i=1; i<=neqns; i++) {
if (xlnz[i] < xlnz[i+1])
parent[i] = nzsub[xnzsub[i]];
else
parent[i] = -1; /* Indicates no parent */
}
supnd = idxsmalloc(neqns+1, 0, "smbfct: supnd");
supnd--;
/* super-noeuds */
nbip1 = xlnz[2] -xlnz[1] ;
xnzip1 = xnzsub[1];
i=1 ;
snode = 1 ; supnd[snode] = 1 ;
lgind = nbip1;
marker[i]= snode;
nbi = nbip1;
xnzi = xnzip1;
marker[i]= snode;
i = i+ 1;
/* C.Rose ligne suivante : correction de bug 13/03/02 neqns+1 remplace neqns
dans l'instruction suivante
Il y avait un pb pour ssll102a avec le dernier SN qui ne comprend qu'un noeud*/
while ( i < neqns+1)
{
nbip1 = xlnz[i+1] -xlnz[i] ;
xnzip1 = xnzsub[i];
if ( nbip1 != (nbi-1))
{
snode ++;supnd[snode] = i;
lgind += nbip1;
}
else if ( xnzip1 != (xnzi+1))
{
snode ++;supnd[snode] = i;
lgind += nbip1;
}
nbi = nbip1;
xnzi = xnzip1;
marker[i]= snode;
i = i+ 1;
}
supnd[snode+1] = neqns +1;
lnode =neqns;
marker[neqns] = snode ;
/* C.Rose ajout de la ligne precedante : correction de bug 11/03/02 */
/* printf(" nbsn : %d \n",snode);
for(i=1;i<=neqns;i++) printf(" smarker %8ld %8ld\n ",i,marker[i]);
for(i=1;i<=snode+1;i++) printf(" snode %8ld %8ld\n ",i,supnd[i]);
*/
/* fin de cr*/
minsn = neqns;
maxsn = 0;
k1 = 0;
k2 = 0;
k3 = 0;
taille1 = 50;
taille2 = 20;
taille3 = 05;
for (i = 1; i<=snode; i++)
{
aux = supnd[i+1]-supnd[i];
if(aux < minsn) minsn = aux;
if(aux > maxsn) maxsn = aux;
if(aux < taille1) k1++;
if(aux < taille2) k2++;
if(aux < taille3) k3++;
}
printf("\n ");
printf(" taille du plus petit super-noeud %8ld\n ",minsn);
printf(" taille du plus grand super-noeud %8ld\n ",maxsn);
printf(" taille moyenne d un super-noeud %8ld\n ",lnode/snode);
printf(" nb de super-noeuds de taille < %4ld : %8ld\n ",taille1,k1);
printf(" nb de super-noeuds de taille < %4ld : %8ld\n ",taille2,k2);
printf(" nb de super-noeuds de taille < %4ld : %8ld\n ",taille3,k3);
/* for (i=1; i<=snode+1; i++) {printf(" i pt. %4ld%4d\n ",i,supnd[i]);} */
/* parent des super-noeuds */
/* principe : on a stocke dans marker pour chaque noeud son super-noeud */
/* maintenant pour chaque super-noeud on prend son dernier noeud, */
/* on cherche le parent de ce noeud, puis le supernoeud auquel */
/* appartient ce parent et on l'affecte comme parent au supernoeud */
/* quand il n'y a pas de parent on met 0 (convention C. Rose) */
parsnode = idxsmalloc(snode+1, 0, "smbfct: parsnode");
parsnode--;
tbnode = idxsmalloc(neqns, 0, "smbfct: tbnode");
tbnode--;
/*sauvegarde de parent nodal */
for(i = 1; i<=neqns; i++) {
tbnode[i] = parent[i];
/* printf(" i %4ld parent %8ld \n", i, tbnode[i]);*/
}
for(i = 1; i<=snode; i++)
{
nd = supnd[i+1] - 1;
if(tbnode[nd] != -1)
{
parsnode[i] = marker[tbnode[nd]];
}
else
{
parsnode[i] = 0;
}
/* printf(" noeud %4ld parent %8ld \n ",i, parsnode[i]); */
}
opc = 0;
for (i=1; i<=neqns; i++)
xlnz[i]--;
for (i=1; i<=neqns; i++) {
opc += (xlnz[i+1]-xlnz[i])*(xlnz[i+1]-xlnz[i]) - (xlnz[i+1]-xlnz[i]);
}
printf(" Nonzeros: %ld, \tOperation Count: %6.4le \n", *maxlnz, opc);
/* ecriture sur le fichier 85 */
ecri11_(invp,perm,supnd,parsnode,&neqns,&snode,&opc,maxlnz,&lgind);
marker++;
mrglnk++;
rchlnk++;
nzsub++;
xnzsub++;
xlnz++;
invp++;
perm++;
adjncy++;
xadj++;
GKfree(&rchlnk, &mrglnk, &marker, LTERM);
return flag;
}
syntax highlighted by Code2HTML, v. 0.9.1