Page MenuHomec4science
Files F121678152

symfct.c
No OneTemporary
Actions

File Metadata

Created
Sun, Jul 13, 02:15

symfct.c
View Options

/* symfct.f -- translated by f2c (version 19951025).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "mex.h"
typedef mwSignedIndex integer; /* removed "long" */
#if !defined(max)
#define max(A, B) ((A) > (B) ? (A) : (B))
#endif
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: January 12, 1995 */
/* Authors: Esmond G. Ng and Barry W. Peyton */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* ************** SFINIT ..... SET UP FOR SYMB. FACT. ************ */
/* *********************************************************************** */
/* *********************************************************************** */
/* PURPOSE: */
/* THIS SUBROUTINE COMPUTES THE STORAGE REQUIREMENTS AND SETS UP */
/* PRELIMINARY DATA STRUCTURES FOR THE SYMBOLIC FACTORIZATION. */
/* NOTE: */
/* THIS VERSION PRODUCES THE MAXIMAL SUPERNODE PARTITION (I.E., */
/* THE ONE WITH THE FEWEST POSSIBLE SUPERNODES). */
/* INPUT PARAMETERS: */
/* NEQNS - NUMBER OF EQUATIONS. */
/* NNZA - LENGTH OF ADJACENCY STRUCTURE. */
/* XADJ(*) - ARRAY OF LENGTH NEQNS+1, CONTAINING POINTERS */
/* TO THE ADJACENCY STRUCTURE. */
/* ADJNCY(*) - ARRAY OF LENGTH XADJ(NEQNS+1)-1, CONTAINING */
/* THE ADJACENCY STRUCTURE. */
/* IWSIZ - SIZE OF INTEGER WORKING STORAGE. */
/* UPDATED PARAMETERS: */
/* (* JFS Sept 1, 1998: moved "PERM" and "INVP" to updated section *) */
/* (PERM,INVP) - ON INPUT, THE GIVEN PERM AND INVERSE PERM */
/* VECTORS. ON OUTPUT, THE NEW PERM AND */
/* INVERSE PERM VECTORS OF THE EQUIVALENT */
/* ORDERING. */
/* OUTPUT PARAMETERS: */
/* COLCNT(*) - ARRAY OF LENGTH NEQNS, CONTAINING THE NUMBER */
/* OF NONZEROS IN EACH COLUMN OF THE FACTOR, */
/* INCLUDING THE DIAGONAL ENTRY. */
/* NNZL - NUMBER OF NONZEROS IN THE FACTOR, INCLUDING */
/* THE DIAGONAL ENTRIES. */
/* NSUB - NUMBER OF SUBSCRIPTS. */
/* NSUPER - NUMBER OF SUPERNODES (<= NEQNS). */
/* SNODE(*) - ARRAY OF LENGTH NEQNS FOR RECORDING */
/* SUPERNODE MEMBERSHIP. */
/* XSUPER(*) - ARRAY OF LENGTH NEQNS+1, CONTAINING THE */
/* SUPERNODE PARTITIONING. */
/* IFLAG(*) - ERROR FLAG. */
/* 0: SUCCESSFUL SF INITIALIZATION. */
/* -1: INSUFFICENT WORKING STORAGE */
/* [IWORK(*)]. */
/* WORK PARAMETERS: */
/* IWORK(*) - INTEGER WORK ARRAY OF LENGTH 7*NEQNS+3. */
/* FIRST CREATED ON NOVEMBER 14, 1994. */
/* LAST UPDATED ON January 12, 1995. */
/* *********************************************************************** */
/* Subroutine */ int sfinit_(neqns, nnza, xadj, adjncy, perm, invp, colcnt,
nnzl, nsub, nsuper, snode, xsuper, iwsiz, iwork, iflag)
integer *neqns, *nnza, *xadj, *adjncy, *perm, *invp, *colcnt, *nnzl, *nsub, *
nsuper, *snode, *xsuper, *iwsiz, *iwork, *iflag;
{
/* System generated locals */
integer i__1;
/* Local variables */
extern /* Subroutine */ int fsup1_(), fsup2_();
static integer i__;
extern /* Subroutine */ int fcnthn_(), chordr_(), etordr_();
/* ----------- */
/* PARAMETERS. */
/* ----------- */
/* ***********************************************************************
*/
/* -------------------------------------------------------- */
/* RETURN IF THERE IS INSUFFICIENT INTEGER WORKING STORAGE. */
/* -------------------------------------------------------- */
/* Parameter adjustments */
--iwork;
--xsuper;
--snode;
--colcnt;
--invp;
--perm;
--xadj;
--adjncy;
/* Function Body */
*iflag = 0;
if (*iwsiz < *neqns * 7 + 3) {
*iflag = -1;
return 0;
}
/* (* JFS: */
/* ------------------------------------------------------------ */
/* Handle the case of diagonal matrices separately, to avoid */
/* an unsolved BUG in FCNTHN. */
/* This patch is due to Jos F. Sturm, Sept 1, 1998. */
/* ------------------------------------------------------------ */
if (xadj[*neqns + 1] - 1 == 0) {
i__1 = *neqns;
for (i__ = 1; i__ <= i__1; ++i__) {
colcnt[i__] = 1;
snode[i__] = i__;
xsuper[i__] = i__;
/* L10: */
}
xsuper[*neqns + 1] = *neqns + 1;
*nnzl = *neqns;
*nsub = *neqns;
*nsuper = *neqns;
*iflag = 0;
return 0;
}
/* ------------------------------------------------------------ */
/* end of patch JFS *) */
/* ------------------------------------------ */
/* COMPUTE ELIMINATION TREE AND POSTORDERING. */
/* ------------------------------------------ */
etordr_(neqns, &xadj[1], &adjncy[1], &perm[1], &invp[1], &iwork[1], &
iwork[*neqns + 1], &iwork[(*neqns << 1) + 1], &iwork[*neqns * 3 +
1]);
/* --------------------------------------------- */
/* COMPUTE ROW AND COLUMN FACTOR NONZERO COUNTS. */
/* --------------------------------------------- */
fcnthn_(neqns, nnza, &xadj[1], &adjncy[1], &perm[1], &invp[1], &iwork[1],
&snode[1], &colcnt[1], nnzl, &iwork[*neqns + 1], &iwork[(*neqns <<
1) + 1], &xsuper[1], &iwork[*neqns * 3 + 1], &iwork[(*neqns << 2)
+ 2], &iwork[*neqns * 5 + 3], &iwork[*neqns * 6 + 4]);
/* --------------------------------------------------------- */
/* REARRANGE CHILDREN SO THAT THE LAST CHILD HAS THE MAXIMUM */
/* NUMBER OF NONZEROS IN ITS COLUMN OF L. */
/* --------------------------------------------------------- */
chordr_(neqns, &xadj[1], &adjncy[1], &perm[1], &invp[1], &colcnt[1], &
iwork[1], &iwork[*neqns + 1], &iwork[(*neqns << 1) + 1], &iwork[*
neqns * 3 + 1]);
/* ---------------- */
/* FIND SUPERNODES. */
/* ---------------- */
fsup1_(neqns, &iwork[1], &colcnt[1], nsub, nsuper, &snode[1]);
fsup2_(neqns, nsuper, &iwork[1], &snode[1], &xsuper[1]);
return 0;
} /* sfinit_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: December 27, 1994 */
/* Authors: Joseph W.H. Liu */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* ********** ETORDR ..... ELIMINATION TREE REORDERING *********** */
/* *********************************************************************** */
/* *********************************************************************** */
/* WRITTEN BY JOSEPH LIU (JUL 17, 1985) */
/* PURPOSE: */
/* TO DETERMINE AN EQUIVALENT REORDERING BASED ON THE STRUCTURE OF */
/* THE ELIMINATION TREE. A POSTORDERING OF THE GIVEN ELIMINATION */
/* TREE IS RETURNED. */
/* INPUT PARAMETERS: */
/* NEQNS - NUMBER OF EQUATIONS. */
/* (XADJ,ADJNCY) - THE ADJACENCY STRUCTURE. */
/* UPDATED PARAMETERS: */
/* (PERM,INVP) - ON INPUT, THE GIVEN PERM AND INVERSE PERM */
/* VECTORS. ON OUTPUT, THE NEW PERM AND */
/* INVERSE PERM VECTORS OF THE EQUIVALENT */
/* ORDERING. */
/* OUTPUT PARAMETERS: */
/* PARENT - THE PARENT VECTOR OF THE ELIMINATION TREE */
/* ASSOCIATED WITH THE NEW ORDERING. */
/* WORKING PARAMETERS: */
/* FSON - THE FIRST SON VECTOR. */
/* BROTHR - THE BROTHER VECTOR. */
/* INVPOS - THE INVERSE PERM VECTOR FOR THE */
/* POSTORDERING. */
/* PROGRAM SUBROUTINES: */
/* BETREE, ETPOST, ETREE , INVINV. */
/* *********************************************************************** */
/* Subroutine */ int etordr_(neqns, xadj, adjncy, perm, invp, parent, fson,
brothr, invpos)
integer *neqns, *xadj, *adjncy, *perm, *invp, *parent, *fson, *brothr, *
invpos;
{
extern /* Subroutine */ int etree_(), betree_(), invinv_(), etpost_();
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* ----------------------------- */
/* COMPUTE THE ELIMINATION TREE. */
/* ----------------------------- */
/* Parameter adjustments */
--invpos;
--brothr;
--fson;
--parent;
--invp;
--perm;
--adjncy;
--xadj;
/* Function Body */
etree_(neqns, &xadj[1], &adjncy[1], &perm[1], &invp[1], &parent[1], &
invpos[1]);
/* -------------------------------------------------------- */
/* COMPUTE A BINARY REPRESENTATION OF THE ELIMINATION TREE. */
/* -------------------------------------------------------- */
betree_(neqns, &parent[1], &fson[1], &brothr[1]);
/* ------------------------------- */
/* POSTORDER THE ELIMINATION TREE. */
/* ------------------------------- */
etpost_(neqns, &fson[1], &brothr[1], &invpos[1], &parent[1], &perm[1]);
/* -------------------------------------------------------- */
/* COMPOSE THE ORIGINAL ORDERING WITH THE NEW POSTORDERING. */
/* -------------------------------------------------------- */
invinv_(neqns, &invp[1], &invpos[1], &perm[1]);
return 0;
} /* etordr_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: December 27, 1994 */
/* Authors: Joseph W.H. Liu */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* **************** ETREE ..... ELIMINATION TREE ***************** */
/* *********************************************************************** */
/* *********************************************************************** */
/* WRITTEN BY JOSEPH LIU (JUL 17, 1985) */
/* PURPOSE: */
/* TO DETERMINE THE ELIMINATION TREE FROM A GIVEN ORDERING AND */
/* THE ADJACENCY STRUCTURE. THE PARENT VECTOR IS RETURNED. */
/* INPUT PARAMETERS: */
/* NEQNS - NUMBER OF EQUATIONS. */
/* (XADJ,ADJNCY) - THE ADJACENCY STRUCTURE. */
/* (PERM,INVP) - PERMUTATION AND INVERSE PERMUTATION VECTORS */
/* OUTPUT PARAMETERS: */
/* PARENT - THE PARENT VECTOR OF THE ELIMINATION TREE. */
/* WORKING PARAMETERS: */
/* ANCSTR - THE ANCESTOR VECTOR. */
/* *********************************************************************** */
/* Subroutine */ int etree_(neqns, xadj, adjncy, perm, invp, parent, ancstr)
integer *neqns, *xadj, *adjncy, *perm, *invp, *parent, *ancstr;
{
/* System generated locals */
integer i__1, i__2;
/* Local variables */
static integer node, next, i__, j, jstop, jstrt, nbr;
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* Parameter adjustments */
--ancstr;
--parent;
--invp;
--perm;
--adjncy;
--xadj;
/* Function Body */
if (*neqns <= 0) {
return 0;
}
i__1 = *neqns;
for (i__ = 1; i__ <= i__1; ++i__) {
parent[i__] = 0;
ancstr[i__] = 0;
node = perm[i__];
jstrt = xadj[node];
jstop = xadj[node + 1] - 1;
if (jstrt <= jstop) {
i__2 = jstop;
for (j = jstrt; j <= i__2; ++j) {
nbr = adjncy[j];
nbr = invp[nbr];
if (nbr < i__) {
/* --------------------------------
----------- */
/* FOR EACH NBR, FIND THE ROOT OF IT
S CURRENT */
/* ELIMINATION TREE. PERFORM PATH C
OMPRESSION */
/* AS THE SUBTREE IS TRAVERSED. */
/* --------------------------------
----------- */
L100:
if (ancstr[nbr] == i__) {
goto L300;
}
if (ancstr[nbr] > 0) {
next = ancstr[nbr];
ancstr[nbr] = i__;
nbr = next;
goto L100;
}
/* --------------------------------
------------ */
/* NOW, NBR IS THE ROOT OF THE SUBTR
EE. MAKE I */
/* THE PARENT NODE OF THIS ROOT. */
/* --------------------------------
------------ */
parent[nbr] = i__;
ancstr[nbr] = i__;
}
L300:
;
}
}
/* L400: */
}
return 0;
} /* etree_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: December 27, 1994 */
/* Authors: Joseph W.H. Liu */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* ****** BETREE ..... BINARY TREE REPRESENTATION OF ETREE ******* */
/* *********************************************************************** */
/* *********************************************************************** */
/* WRITTEN BY JOSEPH LIU (JUL 17, 1985) */
/* PURPOSE: */
/* TO DETERMINE THE BINARY TREE REPRESENTATION OF THE ELIMINATION */
/* TREE GIVEN BY THE PARENT VECTOR. THE RETURNED REPRESENTATION */
/* WILL BE GIVEN BY THE FIRST-SON AND BROTHER VECTORS. THE ROOT */
/* OF THE BINARY TREE IS ALWAYS NEQNS. */
/* INPUT PARAMETERS: */
/* NEQNS - NUMBER OF EQUATIONS. */
/* PARENT - THE PARENT VECTOR OF THE ELIMINATION TREE. */
/* IT IS ASSUMED THAT PARENT(I) > I EXCEPT OF */
/* THE ROOTS. */
/* OUTPUT PARAMETERS: */
/* FSON - THE FIRST SON VECTOR. */
/* BROTHR - THE BROTHER VECTOR. */
/* *********************************************************************** */
/* Subroutine */ int betree_(neqns, parent, fson, brothr)
integer *neqns, *parent, *fson, *brothr;
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer node, ndpar, lroot;
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* Parameter adjustments */
--brothr;
--fson;
--parent;
/* Function Body */
if (*neqns <= 0) {
return 0;
}
i__1 = *neqns;
for (node = 1; node <= i__1; ++node) {
fson[node] = 0;
brothr[node] = 0;
/* L100: */
}
lroot = *neqns;
/* ------------------------------------------------------------ */
/* FOR EACH NODE := NEQNS-1 STEP -1 DOWNTO 1, DO THE FOLLOWING. */
/* ------------------------------------------------------------ */
if (*neqns <= 1) {
return 0;
}
for (node = *neqns - 1; node >= 1; --node) {
ndpar = parent[node];
if (ndpar <= 0 || ndpar == node) {
/* -------------------------------------------------
*/
/* NODE HAS NO PARENT. GIVEN STRUCTURE IS A FOREST.
*/
/* SET NODE TO BE ONE OF THE ROOTS OF THE TREES. */
/* -------------------------------------------------
*/
brothr[lroot] = node;
lroot = node;
} else {
/* ------------------------------------------- */
/* OTHERWISE, BECOMES FIRST SON OF ITS PARENT. */
/* ------------------------------------------- */
brothr[node] = fson[ndpar];
fson[ndpar] = node;
}
/* L300: */
}
brothr[lroot] = 0;
return 0;
} /* betree_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: December 27, 1994 */
/* Authors: Joseph W.H. Liu */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* *************** ETPOST ..... ETREE POSTORDERING *************** */
/* *********************************************************************** */
/* *********************************************************************** */
/* WRITTEN BY JOSEPH LIU (SEPT 17, 1986) */
/* PURPOSE: */
/* BASED ON THE BINARY REPRESENTATION (FIRST-SON,BROTHER) OF */
/* THE ELIMINATION TREE, A POSTORDERING IS DETERMINED. THE */
/* CORRESPONDING PARENT VECTOR IS ALSO MODIFIED TO REFLECT */
/* THE REORDERING. */
/* INPUT PARAMETERS: */
/* ROOT - ROOT OF THE ELIMINATION TREE (USUALLY IT */
/* IS NEQNS). */
/* FSON - THE FIRST SON VECTOR. */
/* BROTHR - THE BROTHR VECTOR. */
/* UPDATED PARAMETERS: */
/* PARENT - THE PARENT VECTOR. */
/* OUTPUT PARAMETERS: */
/* INVPOS - INVERSE PERMUTATION FOR THE POSTORDERING. */
/* WORKING PARAMETERS: */
/* STACK - THE STACK FOR POSTORDER TRAVERSAL OF THE */
/* TREE. */
/* *********************************************************************** */
/* Subroutine */ int etpost_(root, fson, brothr, invpos, parent, stack)
integer *root, *fson, *brothr, *invpos, *parent, *stack;
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer node, itop, ndpar, nunode, num;
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* Parameter adjustments */
--stack;
--parent;
--invpos;
--brothr;
--fson;
/* Function Body */
num = 0;
itop = 0;
node = *root;
/* ------------------------------------------------------------- */
/* TRAVERSE ALONG THE FIRST SONS POINTER AND PUSH THE TREE NODES */
/* ALONG THE TRAVERSAL INTO THE STACK. */
/* ------------------------------------------------------------- */
L100:
++itop;
stack[itop] = node;
node = fson[node];
if (node > 0) {
goto L100;
}
/* ---------------------------------------------------------- */
/* IF POSSIBLE, POP A TREE NODE FROM THE STACK AND NUMBER IT. */
/* ---------------------------------------------------------- */
L200:
if (itop <= 0) {
goto L300;
}
node = stack[itop];
--itop;
++num;
invpos[node] = num;
/* ---------------------------------------------------- */
/* THEN, TRAVERSE TO ITS YOUNGER BROTHER IF IT HAS ONE. */
/* ---------------------------------------------------- */
node = brothr[node];
if (node <= 0) {
goto L200;
}
goto L100;
L300:
/* ------------------------------------------------------------ */
/* DETERMINE THE NEW PARENT VECTOR OF THE POSTORDERING. BROTHR */
/* IS USED TEMPORARILY FOR THE NEW PARENT VECTOR. */
/* ------------------------------------------------------------ */
i__1 = num;
for (node = 1; node <= i__1; ++node) {
nunode = invpos[node];
ndpar = parent[node];
if (ndpar > 0) {
ndpar = invpos[ndpar];
}
brothr[nunode] = ndpar;
/* L400: */
}
i__1 = num;
for (nunode = 1; nunode <= i__1; ++nunode) {
parent[nunode] = brothr[nunode];
/* L500: */
}
return 0;
} /* etpost_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: December 27, 1994 */
/* Authors: Joseph W.H. Liu */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* *********** INVINV ..... CONCATENATION OF TWO INVP ************ */
/* *********************************************************************** */
/* *********************************************************************** */
/* WRITTEN BY JOSEPH LIU (JUL 17, 1985) */
/* PURPOSE: */
/* TO PERFORM THE MAPPING OF */
/* ORIGINAL-INVP --> INTERMEDIATE-INVP --> NEW INVP */
/* AND THE RESULTING ORDERING REPLACES INVP. THE NEW PERMUTATION */
/* VECTOR PERM IS ALSO COMPUTED. */
/* INPUT PARAMETERS: */
/* NEQNS - NUMBER OF EQUATIONS. */
/* INVP2 - THE SECOND INVERSE PERMUTATION VECTOR. */
/* UPDATED PARAMETERS: */
/* INVP - THE FIRST INVERSE PERMUTATION VECTOR. ON */
/* OUTPUT, IT CONTAINS THE NEW INVERSE */
/* PERMUTATION. */
/* OUTPUT PARAMETER: */
/* PERM - NEW PERMUTATION VECTOR (CAN BE THE SAME AS */
/* INVP2). */
/* *********************************************************************** */
/* Subroutine */ int invinv_(neqns, invp, invp2, perm)
integer *neqns, *invp, *invp2, *perm;
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer node, i__, interm;
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* Parameter adjustments */
--perm;
--invp2;
--invp;
/* Function Body */
i__1 = *neqns;
for (i__ = 1; i__ <= i__1; ++i__) {
interm = invp[i__];
invp[i__] = invp2[interm];
/* L100: */
}
i__1 = *neqns;
for (i__ = 1; i__ <= i__1; ++i__) {
node = invp[i__];
perm[node] = i__;
/* L200: */
}
return 0;
} /* invinv_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: December 27, 1994 */
/* Authors: Esmond G. Ng and Barry W. Peyton */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* ********** CHORDR ..... CHILD REORDERING *********** */
/* *********************************************************************** */
/* *********************************************************************** */
/* PURPOSE: */
/* REARRANGE THE CHILDREN OF EACH VERTEX SO THAT THE LAST ONE */
/* MAXIMIZES (AMONG THE CHILDREN) THE NUMBER OF NONZEROS IN THE */
/* CORRESPONDING COLUMN OF L. ALSO DETERMINE AN NEW POSTORDERING */
/* BASED ON THE STRUCTURE OF THE MODIFIED ELIMINATION TREE. */
/* INPUT PARAMETERS: */
/* NEQNS - NUMBER OF EQUATIONS. */
/* (XADJ,ADJNCY) - THE ADJACENCY STRUCTURE. */
/* UPDATED PARAMETERS: */
/* (PERM,INVP) - ON INPUT, THE GIVEN PERM AND INVERSE PERM */
/* VECTORS. ON OUTPUT, THE NEW PERM AND */
/* INVERSE PERM VECTORS OF THE NEW */
/* POSTORDERING. */
/* COLCNT - COLUMN COUNTS IN L UNDER INITIAL ORDERING; */
/* MODIFIED TO REFLECT THE NEW ORDERING. */
/* OUTPUT PARAMETERS: */
/* PARENT - THE PARENT VECTOR OF THE ELIMINATION TREE */
/* ASSOCIATED WITH THE NEW ORDERING. */
/* WORKING PARAMETERS: */
/* FSON - THE FIRST SON VECTOR. */
/* BROTHR - THE BROTHER VECTOR. */
/* INVPOS - THE INVERSE PERM VECTOR FOR THE */
/* POSTORDERING. */
/* PROGRAM SUBROUTINES: */
/* BTREE2, EPOST2, INVINV. */
/* *********************************************************************** */
/* Subroutine */ int chordr_(neqns, xadj, adjncy, perm, invp, colcnt, parent,
fson, brothr, invpos)
integer *neqns, *xadj, *adjncy, *perm, *invp, *colcnt, *parent, *fson, *
brothr, *invpos;
{
extern /* Subroutine */ int btree2_(), epost2_(), invinv_();
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* ---------------------------------------------------------- */
/* COMPUTE A BINARY REPRESENTATION OF THE ELIMINATION TREE, */
/* SO THAT EACH "LAST CHILD" MAXIMIZES AMONG ITS SIBLINGS THE */
/* NUMBER OF NONZEROS IN THE CORRESPONDING COLUMNS OF L. */
/* ---------------------------------------------------------- */
/* Parameter adjustments */
--invpos;
--brothr;
--fson;
--parent;
--colcnt;
--invp;
--perm;
--adjncy;
--xadj;
/* Function Body */
btree2_(neqns, &parent[1], &colcnt[1], &fson[1], &brothr[1], &invpos[1]);
/* ---------------------------------------------------- */
/* POSTORDER THE ELIMINATION TREE (USING THE NEW BINARY */
/* REPRESENTATION. */
/* ---------------------------------------------------- */
epost2_(neqns, &fson[1], &brothr[1], &invpos[1], &parent[1], &colcnt[1], &
perm[1]);
/* -------------------------------------------------------- */
/* COMPOSE THE ORIGINAL ORDERING WITH THE NEW POSTORDERING. */
/* -------------------------------------------------------- */
invinv_(neqns, &invp[1], &invpos[1], &perm[1]);
return 0;
} /* chordr_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: January 12, 1995 */
/* Authors: Esmond G. Ng and Barry W. Peyton */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* ****** BTREE2 ..... BINARY TREE REPRESENTATION OF ETREE ******* */
/* *********************************************************************** */
/* *********************************************************************** */
/* PURPOSE: */
/* TO DETERMINE A BINARY TREE REPRESENTATION OF THE ELIMINATION */
/* TREE, FOR WHICH EVERY "LAST CHILD" HAS THE MAXIMUM POSSIBLE */
/* COLUMN NONZERO COUNT IN THE FACTOR. THE RETURNED REPRESENTATION */
/* WILL BE GIVEN BY THE FIRST-SON AND BROTHER VECTORS. THE ROOT OF */
/* THE BINARY TREE IS ALWAYS NEQNS. */
/* INPUT PARAMETERS: */
/* NEQNS - NUMBER OF EQUATIONS. */
/* PARENT - THE PARENT VECTOR OF THE ELIMINATION TREE. */
/* IT IS ASSUMED THAT PARENT(I) > I EXCEPT OF */
/* THE ROOTS. */
/* COLCNT - COLUMN NONZERO COUNTS OF THE FACTOR. */
/* OUTPUT PARAMETERS: */
/* FSON - THE FIRST SON VECTOR. */
/* BROTHR - THE BROTHER VECTOR. */
/* WORKING PARAMETERS: */
/* LSON - LAST SON VECTOR. */
/* *********************************************************************** */
/* Subroutine */ int btree2_(neqns, parent, colcnt, fson, brothr, lson)
integer *neqns, *parent, *colcnt, *fson, *brothr, *lson;
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer node, ndpar, lroot, ndlson;
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* Parameter adjustments */
--lson;
--brothr;
--fson;
--colcnt;
--parent;
/* Function Body */
if (*neqns <= 0) {
return 0;
}
i__1 = *neqns;
for (node = 1; node <= i__1; ++node) {
fson[node] = 0;
brothr[node] = 0;
lson[node] = 0;
/* L100: */
}
lroot = *neqns;
/* ------------------------------------------------------------ */
/* FOR EACH NODE := NEQNS-1 STEP -1 DOWNTO 1, DO THE FOLLOWING. */
/* ------------------------------------------------------------ */
if (*neqns <= 1) {
return 0;
}
for (node = *neqns - 1; node >= 1; --node) {
ndpar = parent[node];
if (ndpar <= 0 || ndpar == node) {
/* -------------------------------------------------
*/
/* NODE HAS NO PARENT. GIVEN STRUCTURE IS A FOREST.
*/
/* SET NODE TO BE ONE OF THE ROOTS OF THE TREES. */
/* -------------------------------------------------
*/
brothr[lroot] = node;
lroot = node;
} else {
/* ------------------------------------------- */
/* OTHERWISE, BECOMES FIRST SON OF ITS PARENT. */
/* ------------------------------------------- */
ndlson = lson[ndpar];
if (ndlson != 0) {
if (colcnt[node] >= colcnt[ndlson]) {
brothr[node] = fson[ndpar];
fson[ndpar] = node;
} else {
brothr[ndlson] = node;
lson[ndpar] = node;
}
} else {
fson[ndpar] = node;
lson[ndpar] = node;
}
}
/* L300: */
}
brothr[lroot] = 0;
return 0;
} /* btree2_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: December 27, 1994 */
/* Authors: Esmond G. Ng and Barry W. Peyton */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* *************** EPOST2 ..... ETREE POSTORDERING #2 *************** */
/* *********************************************************************** */
/* *********************************************************************** */
/* PURPOSE: */
/* BASED ON THE BINARY REPRESENTATION (FIRST-SON,BROTHER) OF THE */
/* ELIMINATION TREE, A POSTORDERING IS DETERMINED. THE */
/* CORRESPONDING PARENT AND COLCNT VECTORS ARE ALSO MODIFIED TO */
/* REFLECT THE REORDERING. */
/* INPUT PARAMETERS: */
/* ROOT - ROOT OF THE ELIMINATION TREE (USUALLY IT */
/* IS NEQNS). */
/* FSON - THE FIRST SON VECTOR. */
/* BROTHR - THE BROTHR VECTOR. */
/* UPDATED PARAMETERS: */
/* PARENT - THE PARENT VECTOR. */
/* COLCNT - COLUMN NONZERO COUNTS OF THE FACTOR. */
/* OUTPUT PARAMETERS: */
/* INVPOS - INVERSE PERMUTATION FOR THE POSTORDERING. */
/* WORKING PARAMETERS: */
/* STACK - THE STACK FOR POSTORDER TRAVERSAL OF THE */
/* TREE. */
/* *********************************************************************** */
/* Subroutine */ int epost2_(root, fson, brothr, invpos, parent, colcnt,
stack)
integer *root, *fson, *brothr, *invpos, *parent, *colcnt, *stack;
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer node, itop, ndpar, nunode, num;
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* Parameter adjustments */
--stack;
--colcnt;
--parent;
--invpos;
--brothr;
--fson;
/* Function Body */
num = 0;
itop = 0;
node = *root;
/* ------------------------------------------------------------- */
/* TRAVERSE ALONG THE FIRST SONS POINTER AND PUSH THE TREE NODES */
/* ALONG THE TRAVERSAL INTO THE STACK. */
/* ------------------------------------------------------------- */
L100:
++itop;
stack[itop] = node;
node = fson[node];
if (node > 0) {
goto L100;
}
/* ---------------------------------------------------------- */
/* IF POSSIBLE, POP A TREE NODE FROM THE STACK AND NUMBER IT. */
/* ---------------------------------------------------------- */
L200:
if (itop <= 0) {
goto L300;
}
node = stack[itop];
--itop;
++num;
invpos[node] = num;
/* ---------------------------------------------------- */
/* THEN, TRAVERSE TO ITS YOUNGER BROTHER IF IT HAS ONE. */
/* ---------------------------------------------------- */
node = brothr[node];
if (node <= 0) {
goto L200;
}
goto L100;
L300:
/* ------------------------------------------------------------ */
/* DETERMINE THE NEW PARENT VECTOR OF THE POSTORDERING. BROTHR */
/* IS USED TEMPORARILY FOR THE NEW PARENT VECTOR. */
/* ------------------------------------------------------------ */
i__1 = num;
for (node = 1; node <= i__1; ++node) {
nunode = invpos[node];
ndpar = parent[node];
if (ndpar > 0) {
ndpar = invpos[ndpar];
}
brothr[nunode] = ndpar;
/* L400: */
}
i__1 = num;
for (nunode = 1; nunode <= i__1; ++nunode) {
parent[nunode] = brothr[nunode];
/* L500: */
}
/* ---------------------------------------------- */
/* PERMUTE COLCNT(*) TO REFLECT THE NEW ORDERING. */
/* ---------------------------------------------- */
i__1 = num;
for (node = 1; node <= i__1; ++node) {
nunode = invpos[node];
stack[nunode] = colcnt[node];
/* L600: */
}
i__1 = num;
for (node = 1; node <= i__1; ++node) {
colcnt[node] = stack[node];
/* L700: */
}
return 0;
} /* epost2_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: January 12, 1995 */
/* Authors: Esmond G. Ng and Barry W. Peyton */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* ************** FCNTHN ..... FIND NONZERO COUNTS *************** */
/* *********************************************************************** */
/* *********************************************************************** */
/* PURPOSE: */
/* THIS SUBROUTINE DETERMINES THE ROW COUNTS AND COLUMN COUNTS IN */
/* THE CHOLESKY FACTOR. IT USES A DISJOINT SET UNION ALGORITHM. */
/* TECHNIQUES: */
/* 1) SUPERNODE DETECTION. */
/* 2) PATH HALVING. */
/* 3) NO UNION BY RANK. */
/* NOTES: */
/* 1) ASSUMES A POSTORDERING OF THE ELIMINATION TREE. */
/* INPUT PARAMETERS: */
/* (I) NEQNS - NUMBER OF EQUATIONS. */
/* (I) ADJLEN - LENGTH OF ADJACENCY STRUCTURE. */
/* (I) XADJ(*) - ARRAY OF LENGTH NEQNS+1, CONTAINING POINTERS */
/* TO THE ADJACENCY STRUCTURE. */
/* (I) ADJNCY(*) - ARRAY OF LENGTH XADJ(NEQNS+1)-1, CONTAINING */
/* THE ADJACENCY STRUCTURE. */
/* (I) PERM(*) - ARRAY OF LENGTH NEQNS, CONTAINING THE */
/* POSTORDERING. */
/* (I) INVP(*) - ARRAY OF LENGTH NEQNS, CONTAINING THE */
/* INVERSE OF THE POSTORDERING. */
/* (I) ETPAR(*) - ARRAY OF LENGTH NEQNS, CONTAINING THE */
/* ELIMINATION TREE OF THE POSTORDERED MATRIX. */
/* OUTPUT PARAMETERS: */
/* (I) ROWCNT(*) - ARRAY OF LENGTH NEQNS, CONTAINING THE NUMBER */
/* OF NONZEROS IN EACH ROW OF THE FACTOR, */
/* INCLUDING THE DIAGONAL ENTRY. */
/* (I) COLCNT(*) - ARRAY OF LENGTH NEQNS, CONTAINING THE NUMBER */
/* OF NONZEROS IN EACH COLUMN OF THE FACTOR, */
/* INCLUDING THE DIAGONAL ENTRY. */
/* (I) NLNZ - NUMBER OF NONZEROS IN THE FACTOR, INCLUDING */
/* THE DIAGONAL ENTRIES. */
/* WORK PARAMETERS: */
/* (I) SET(*) - ARRAY OF LENGTH NEQNS USED TO MAINTAIN THE */
/* DISJOINT SETS (I.E., SUBTREES). */
/* (I) PRVLF(*) - ARRAY OF LENGTH NEQNS USED TO RECORD THE */
/* PREVIOUS LEAF OF EACH ROW SUBTREE. */
/* (I) LEVEL(*) - ARRAY OF LENGTH NEQNS+1 CONTAINING THE LEVEL */
/* (DISTANCE FROM THE ROOT). */
/* (I) WEIGHT(*) - ARRAY OF LENGTH NEQNS+1 CONTAINING WEIGHTS */
/* USED TO COMPUTE COLUMN COUNTS. */
/* (I) FDESC(*) - ARRAY OF LENGTH NEQNS+1 CONTAINING THE */
/* FIRST (I.E., LOWEST-NUMBERED) DESCENDANT. */
/* (I) NCHILD(*) - ARRAY OF LENGTH NEQNS+1 CONTAINING THE */
/* NUMBER OF CHILDREN. */
/* (I) PRVNBR(*) - ARRAY OF LENGTH NEQNS USED TO RECORD THE */
/* PREVIOUS ``LOWER NEIGHBOR'' OF EACH NODE. */
/* FIRST CREATED ON APRIL 12, 1990. */
/* LAST UPDATED ON JANUARY 12, 1995. */
/* (*JFS Sept 1, 1998: there is a BUG in fcnthn: if adjlen = 0, i.e. */
/* the matrix is purely diagonal, then "segment violation" *) */
/* *********************************************************************** */
/* Subroutine */ int fcnthn_(neqns, adjlen, xadj, adjncy, perm, invp, etpar,
rowcnt, colcnt, nlnz, set, prvlf, level, weight, fdesc, nchild,
prvnbr)
integer *neqns, *adjlen, *xadj, *adjncy, *perm, *invp, *etpar, *rowcnt, *
colcnt, *nlnz, *set, *prvlf, *level, *weight, *fdesc, *nchild, *
prvnbr;
{
/* System generated locals */
integer i__1, i__2;
/* Local variables */
static integer temp, xsup, last1, last2, j, k, lflag, pleaf, hinbr, jstop,
jstrt, ifdesc, oldnbr, parent, lownbr, lca;
/* ----------- */
/* PARAMETERS. */
/* ----------- */
/* ---------------- */
/* LOCAL VARIABLES. */
/* ---------------- */
/* ***********************************************************************
*/
/* -------------------------------------------------- */
/* COMPUTE LEVEL(*), FDESC(*), NCHILD(*). */
/* INITIALIZE ROWCNT(*), COLCNT(*), */
/* SET(*), PRVLF(*), WEIGHT(*), PRVNBR(*). */
/* -------------------------------------------------- */
/* Parameter adjustments */
--prvnbr;
--prvlf;
--set;
--colcnt;
--rowcnt;
--etpar;
--invp;
--perm;
--adjncy;
--xadj;
/* Function Body */
level[0] = 0;
for (k = *neqns; k >= 1; --k) {
rowcnt[k] = 1;
colcnt[k] = 0;
set[k] = k;
prvlf[k] = 0;
level[k] = level[etpar[k]] + 1;
weight[k] = 1;
fdesc[k] = k;
nchild[k] = 0;
prvnbr[k] = 0;
/* L100: */
}
nchild[0] = 0;
fdesc[0] = 0;
i__1 = *neqns;
for (k = 1; k <= i__1; ++k) {
parent = etpar[k];
weight[parent] = 0;
++nchild[parent];
ifdesc = fdesc[k];
if (ifdesc < fdesc[parent]) {
fdesc[parent] = ifdesc;
}
/* L200: */
}
/* ------------------------------------ */
/* FOR EACH ``LOW NEIGHBOR'' LOWNBR ... */
/* ------------------------------------ */
i__1 = *neqns;
for (lownbr = 1; lownbr <= i__1; ++lownbr) {
lflag = 0;
ifdesc = fdesc[lownbr];
oldnbr = perm[lownbr];
jstrt = xadj[oldnbr];
jstop = xadj[oldnbr + 1] - 1;
/* ----------------------------------------------- */
/* FOR EACH ``HIGH NEIGHBOR'', HINBR OF LOWNBR ... */
/* ----------------------------------------------- */
i__2 = jstop;
for (j = jstrt; j <= i__2; ++j) {
hinbr = invp[adjncy[j]];
if (hinbr > lownbr) {
if (ifdesc > prvnbr[hinbr]) {
/* ------------------------- */
/* INCREMENT WEIGHT(LOWNBR). */
/* ------------------------- */
++weight[lownbr];
pleaf = prvlf[hinbr];
/* --------------------------------
--------- */
/* IF HINBR HAS NO PREVIOUS ``LOW NE
IGHBOR'' */
/* THEN ... */
/* --------------------------------
--------- */
if (pleaf == 0) {
/* ------------------------
----------------- */
/* ... ACCUMULATE LOWNBR-->H
INBR PATH LENGTH */
/* IN ROWCNT(HINBR). */
/* ------------------------
----------------- */
rowcnt[hinbr] = rowcnt[hinbr] + level[lownbr] - level[
hinbr];
} else {
/* ------------------------
----------------- */
/* ... OTHERWISE, LCA <-- FI
ND(PLEAF), WHICH */
/* IS THE LEAST COMMON A
NCESTOR OF PLEAF */
/* AND LOWNBR. */
/* (PATH HALVING.) */
/* ------------------------
----------------- */
last1 = pleaf;
last2 = set[last1];
lca = set[last2];
L300:
if (lca != last2) {
set[last1] = lca;
last1 = lca;
last2 = set[last1];
lca = set[last2];
goto L300;
}
/* ------------------------
------------- */
/* ACCUMULATE PLEAF-->LCA PA
TH LENGTH IN */
/* ROWCNT(HINBR). */
/* DECREMENT WEIGHT(LCA). */
/* ------------------------
------------- */
rowcnt[hinbr] = rowcnt[hinbr] + level[lownbr] - level[
lca];
--weight[lca];
}
/* --------------------------------
-------------- */
/* LOWNBR NOW BECOMES ``PREVIOUS LEA
F'' OF HINBR. */
/* --------------------------------
-------------- */
prvlf[hinbr] = lownbr;
lflag = 1;
}
/* ----------------------------------------
---------- */
/* LOWNBR NOW BECOMES ``PREVIOUS NEIGHBOR''
OF HINBR. */
/* ----------------------------------------
---------- */
prvnbr[hinbr] = lownbr;
}
/* L500: */
}
/* ---------------------------------------------------- */
/* DECREMENT WEIGHT ( PARENT(LOWNBR) ). */
/* SET ( P(LOWNBR) ) <-- SET ( P(LOWNBR) ) + SET(XSUP). */
/* ---------------------------------------------------- */
parent = etpar[lownbr];
--weight[parent];
if (lflag == 1 || nchild[lownbr] >= 2) {
xsup = lownbr;
}
set[xsup] = parent;
/* L600: */
}
/* --------------------------------------------------------- */
/* USE WEIGHTS TO COMPUTE COLUMN (AND TOTAL) NONZERO COUNTS. */
/* --------------------------------------------------------- */
*nlnz = 0;
i__1 = *neqns;
for (k = 1; k <= i__1; ++k) {
temp = colcnt[k] + weight[k];
colcnt[k] = temp;
*nlnz += temp;
parent = etpar[k];
if (parent != 0) {
colcnt[parent] += temp;
}
/* L700: */
}
return 0;
} /* fcnthn_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: December 27, 1994 */
/* Authors: Esmond G. Ng and Barry W. Peyton */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* **************** FSUP1 ..... FIND SUPERNODES #1 ***************** */
/* *********************************************************************** */
/* *********************************************************************** */
/* PURPOSE: */
/* THIS SUBROUTINE IS THE FIRST OF TWO ROUTINES FOR FINDING A */
/* MAXIMAL SUPERNODE PARTITION. IT RETURNS ONLY THE NUMBER OF */
/* SUPERNODES NSUPER AND THE SUPERNODE MEMBERSHIP VECTOR SNODE(*), */
/* WHICH IS OF LENGTH NEQNS. THE VECTORS OF LENGTH NSUPER ARE */
/* COMPUTED SUBSEQUENTLY BY THE COMPANION ROUTINE FSUP2. */
/* METHOD AND ASSUMPTIONS: */
/* THIS ROUTINE USES THE ELIMINATION TREE AND THE FACTOR COLUMN */
/* COUNTS TO COMPUTE THE SUPERNODE PARTITION; IT ALSO ASSUMES A */
/* POSTORDERING OF THE ELIMINATION TREE. */
/* INPUT PARAMETERS: */
/* (I) NEQNS - NUMBER OF EQUATIONS. */
/* (I) ETPAR(*) - ARRAY OF LENGTH NEQNS, CONTAINING THE */
/* ELIMINATION TREE OF THE POSTORDERED MATRIX. */
/* (I) COLCNT(*) - ARRAY OF LENGTH NEQNS, CONTAINING THE */
/* FACTOR COLUMN COUNTS: I.E., THE NUMBER OF */
/* NONZERO ENTRIES IN EACH COLUMN OF L */
/* (INCLUDING THE DIAGONAL ENTRY). */
/* OUTPUT PARAMETERS: */
/* (I) NOFSUB - NUMBER OF SUBSCRIPTS. */
/* (I) NSUPER - NUMBER OF SUPERNODES (<= NEQNS). */
/* (I) SNODE(*) - ARRAY OF LENGTH NEQNS FOR RECORDING */
/* SUPERNODE MEMBERSHIP. */
/* FIRST CREATED ON JANUARY 18, 1992. */
/* LAST UPDATED ON NOVEMBER 11, 1994. */
/* *********************************************************************** */
/* Subroutine */ int fsup1_(neqns, etpar, colcnt, nofsub, nsuper, snode)
integer *neqns, *etpar, *colcnt, *nofsub, *nsuper, *snode;
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer kcol;
/* ***********************************************************************
*/
/* ----------- */
/* PARAMETERS. */
/* ----------- */
/* ---------------- */
/* LOCAL VARIABLES. */
/* ---------------- */
/* ***********************************************************************
*/
/* -------------------------------------------- */
/* COMPUTE THE FUNDAMENTAL SUPERNODE PARTITION. */
/* -------------------------------------------- */
/* Parameter adjustments */
--snode;
--colcnt;
--etpar;
/* Function Body */
*nsuper = 1;
snode[1] = 1;
*nofsub = colcnt[1];
i__1 = *neqns;
for (kcol = 2; kcol <= i__1; ++kcol) {
if (etpar[kcol - 1] == kcol) {
if (colcnt[kcol - 1] == colcnt[kcol] + 1) {
snode[kcol] = *nsuper;
goto L300;
}
}
++(*nsuper);
snode[kcol] = *nsuper;
*nofsub += colcnt[kcol];
L300:
;
}
return 0;
} /* fsup1_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: December 27, 1994 */
/* Authors: Esmond G. Ng and Barry W. Peyton */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* **************** FSUP2 ..... FIND SUPERNODES #2 ***************** */
/* *********************************************************************** */
/* *********************************************************************** */
/* PURPOSE: */
/* THIS SUBROUTINE IS THE SECOND OF TWO ROUTINES FOR FINDING A */
/* MAXIMAL SUPERNODE PARTITION. IT'S SOLE PURPOSE IS TO */
/* CONSTRUCT THE NEEDED VECTOR OF LENGTH NSUPER: XSUPER(*). THE */
/* FIRST ROUTINE FSUP1 COMPUTES THE NUMBER OF SUPERNODES AND THE */
/* SUPERNODE MEMBERSHIP VECTOR SNODE(*), WHICH IS OF LENGTH NEQNS. */
/* ASSUMPTIONS: */
/* THIS ROUTINE ASSUMES A POSTORDERING OF THE ELIMINATION TREE. IT */
/* ALSO ASSUMES THAT THE OUTPUT FROM FSUP1 IS AVAILABLE. */
/* INPUT PARAMETERS: */
/* (I) NEQNS - NUMBER OF EQUATIONS. */
/* (I) NSUPER - NUMBER OF SUPERNODES (<= NEQNS). */
/* (I) ETPAR(*) - ARRAY OF LENGTH NEQNS, CONTAINING THE */
/* ELIMINATION TREE OF THE POSTORDERED MATRIX. */
/* (I) SNODE(*) - ARRAY OF LENGTH NEQNS FOR RECORDING */
/* SUPERNODE MEMBERSHIP. */
/* OUTPUT PARAMETERS: */
/* (I) XSUPER(*) - ARRAY OF LENGTH NSUPER+1, CONTAINING THE */
/* SUPERNODE PARTITIONING. */
/* FIRST CREATED ON JANUARY 18, 1992. */
/* LAST UPDATED ON NOVEMEBER 22, 1994. */
/* *********************************************************************** */
/* Subroutine */ int fsup2_(neqns, nsuper, etpar, snode, xsuper)
integer *neqns, *nsuper, *etpar, *snode, *xsuper;
{
static integer kcol, ksup, lstsup;
/* ***********************************************************************
*/
/* ----------- */
/* PARAMETERS. */
/* ----------- */
/* ---------------- */
/* LOCAL VARIABLES. */
/* ---------------- */
/* ***********************************************************************
*/
/* ------------------------------------------------- */
/* COMPUTE THE SUPERNODE PARTITION VECTOR XSUPER(*). */
/* ------------------------------------------------- */
/* Parameter adjustments */
--xsuper;
--snode;
--etpar;
/* Function Body */
lstsup = *nsuper + 1;
for (kcol = *neqns; kcol >= 1; --kcol) {
ksup = snode[kcol];
if (ksup != lstsup) {
xsuper[lstsup] = kcol + 1;
}
lstsup = ksup;
/* L100: */
}
xsuper[1] = 1;
return 0;
} /* fsup2_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: February 13, 1995 */
/* Authors: Esmond G. Ng and Barry W. Peyton */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* ************* SYMFCT ..... SYMBOLIC FACTORIZATION ************** */
/* *********************************************************************** */
/* *********************************************************************** */
/* PURPOSE: */
/* THIS ROUTINE CALLS SYMFC2 WHICH PERFORMS SUPERNODAL SYMBOLIC */
/* FACTORIZATION ON A REORDERED LINEAR SYSTEM. */
/* INPUT PARAMETERS: */
/* (I) NEQNS - NUMBER OF EQUATIONS */
/* (I) ADJLEN - LENGTH OF THE ADJACENCY LIST. */
/* (I) XADJ(*) - ARRAY OF LENGTH NEQNS+1 CONTAINING POINTERS */
/* TO THE ADJACENCY STRUCTURE. */
/* (I) ADJNCY(*) - ARRAY OF LENGTH XADJ(NEQNS+1)-1 CONTAINING */
/* THE ADJACENCY STRUCTURE. */
/* (I) PERM(*) - ARRAY OF LENGTH NEQNS CONTAINING THE */
/* POSTORDERING. */
/* (I) INVP(*) - ARRAY OF LENGTH NEQNS CONTAINING THE */
/* INVERSE OF THE POSTORDERING. */
/* (I) COLCNT(*) - ARRAY OF LENGTH NEQNS, CONTAINING THE NUMBER */
/* OF NONZEROS IN EACH COLUMN OF THE FACTOR, */
/* INCLUDING THE DIAGONAL ENTRY. */
/* (I) NSUPER - NUMBER OF SUPERNODES. */
/* (I) XSUPER(*) - ARRAY OF LENGTH NSUPER+1, CONTAINING THE */
/* FIRST COLUMN OF EACH SUPERNODE. */
/* (I) SNODE(*) - ARRAY OF LENGTH NEQNS FOR RECORDING */
/* SUPERNODE MEMBERSHIP. */
/* (I) NOFSUB - NUMBER OF SUBSCRIPTS TO BE STORED IN */
/* LINDX(*). */
/* (I) IWSIZ - SIZE OF INTEGER WORKING STORAGE. */
/* OUTPUT PARAMETERS: */
/* (I) XLINDX - ARRAY OF LENGTH NEQNS+1, CONTAINING POINTERS */
/* INTO THE SUBSCRIPT VECTOR. */
/* (I) LINDX - ARRAY OF LENGTH MAXSUB, CONTAINING THE */
/* COMPRESSED SUBSCRIPTS. */
/* (I) XLNZ - COLUMN POINTERS FOR L. */
/* (I) FLAG - ERROR FLAG: */
/* 0 - NO ERROR. */
/* -1 - INSUFFICIENT INTEGER WORKING SPACE. */
/* -2 - INCONSISTANCY IN THE INPUT. */
/* WORKING PARAMETERS: */
/* (I) IWORK - WORKING ARRAY OF LENGTH NSUPER+2*NEQNS. */
/* *********************************************************************** */
/* Subroutine */ int symfct_(neqns, adjlen, xadj, adjncy, perm, invp, colcnt,
nsuper, xsuper, snode, nofsub, xlindx, lindx, xlnz, iwsiz, iwork,
flag__)
integer *neqns, *adjlen, *xadj, *adjncy, *perm, *invp, *colcnt, *nsuper, *
xsuper, *snode, *nofsub, *xlindx, *lindx, *xlnz, *iwsiz, *iwork, *
flag__;
{
extern /* Subroutine */ int symfc2_();
/* ***********************************************************************
*/
/* ----------- */
/* PARAMETERS. */
/* ----------- */
/* ***********************************************************************
*/
/* Parameter adjustments */
--xlnz;
--snode;
--colcnt;
--invp;
--perm;
--xadj;
--adjncy;
--iwork;
--xlindx;
--xsuper;
--lindx;
/* Function Body */
*flag__ = 0;
if (*iwsiz < *nsuper + (*neqns << 1) + 1) {
*flag__ = -1;
return 0;
}
symfc2_(neqns, adjlen, &xadj[1], &adjncy[1], &perm[1], &invp[1], &colcnt[
1], nsuper, &xsuper[1], &snode[1], nofsub, &xlindx[1], &lindx[1],
&xlnz[1], &iwork[1], &iwork[*nsuper + 1], &iwork[*nsuper + *neqns
+ 2], flag__);
return 0;
} /* symfct_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: February 13, 1995 */
/* Authors: Esmond G. Ng and Barry W. Peyton */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* ************* SYMFC2 ..... SYMBOLIC FACTORIZATION ************** */
/* *********************************************************************** */
/* *********************************************************************** */
/* PURPOSE: */
/* THIS ROUTINE PERFORMS SUPERNODAL SYMBOLIC FACTORIZATION ON A */
/* REORDERED LINEAR SYSTEM. IT ASSUMES ACCESS TO THE COLUMNS */
/* COUNTS, SUPERNODE PARTITION, AND SUPERNODAL ELIMINATION TREE */
/* ASSOCIATED WITH THE FACTOR MATRIX L. */
/* INPUT PARAMETERS: */
/* (I) NEQNS - NUMBER OF EQUATIONS */
/* (I) ADJLEN - LENGTH OF THE ADJACENCY LIST. */
/* (I) XADJ(*) - ARRAY OF LENGTH NEQNS+1 CONTAINING POINTERS */
/* TO THE ADJACENCY STRUCTURE. */
/* (I) ADJNCY(*) - ARRAY OF LENGTH XADJ(NEQNS+1)-1 CONTAINING */
/* THE ADJACENCY STRUCTURE. */
/* (I) PERM(*) - ARRAY OF LENGTH NEQNS CONTAINING THE */
/* POSTORDERING. */
/* (I) INVP(*) - ARRAY OF LENGTH NEQNS CONTAINING THE */
/* INVERSE OF THE POSTORDERING. */
/* (I) COLCNT(*) - ARRAY OF LENGTH NEQNS, CONTAINING THE NUMBER */
/* OF NONZEROS IN EACH COLUMN OF THE FACTOR, */
/* INCLUDING THE DIAGONAL ENTRY. */
/* (I) NSUPER - NUMBER OF SUPERNODES. */
/* (I) XSUPER(*) - ARRAY OF LENGTH NSUPER+1, CONTAINING THE */
/* FIRST COLUMN OF EACH SUPERNODE. */
/* (I) SNODE(*) - ARRAY OF LENGTH NEQNS FOR RECORDING */
/* SUPERNODE MEMBERSHIP. */
/* (I) NOFSUB - NUMBER OF SUBSCRIPTS TO BE STORED IN */
/* LINDX(*). */
/* OUTPUT PARAMETERS: */
/* (I) XLINDX - ARRAY OF LENGTH NEQNS+1, CONTAINING POINTERS */
/* INTO THE SUBSCRIPT VECTOR. */
/* (I) LINDX - ARRAY OF LENGTH MAXSUB, CONTAINING THE */
/* COMPRESSED SUBSCRIPTS. */
/* (I) XLNZ - COLUMN POINTERS FOR L. */
/* (I) FLAG - ERROR FLAG: */
/* 0 - NO ERROR. */
/* 1 - INCONSISTANCY IN THE INPUT. */
/* WORKING PARAMETERS: */
/* (I) MRGLNK - ARRAY OF LENGTH NSUPER, CONTAINING THE */
/* CHILDREN OF EACH SUPERNODE AS A LINKED LIST. */
/* (I) RCHLNK - ARRAY OF LENGTH NEQNS+1, CONTAINING THE */
/* CURRENT LINKED LIST OF MERGED INDICES (THE */
/* "REACH" SET). */
/* (I) MARKER - ARRAY OF LENGTH NEQNS USED TO MARK INDICES */
/* AS THEY ARE INTRODUCED INTO EACH SUPERNODE'S */
/* INDEX SET. */
/* *********************************************************************** */
/* Subroutine */ int symfc2_(neqns, adjlen, xadj, adjncy, perm, invp, colcnt,
nsuper, xsuper, snode, nofsub, xlindx, lindx, xlnz, mrglnk, rchlnk,
marker, flag__)
integer *neqns, *adjlen, *xadj, *adjncy, *perm, *invp, *colcnt, *nsuper, *
xsuper, *snode, *nofsub, *xlindx, *lindx, *xlnz, *mrglnk, *rchlnk, *
marker, *flag__;
{
/* System generated locals */
integer i__1, i__2;
/* Local variables */
static integer head, node, tail, pcol, newi, jptr, kptr, jsup, ksup, psup,
i__, nzbeg, nzend, width, nexti, point, jnzbeg, knzbeg, length,
jnzend, jwidth, fstcol, knzend, lstcol, knz;
/* ***********************************************************************
*/
/* ----------- */
/* PARAMETERS. */
/* ----------- */
/* ---------------- */
/* LOCAL VARIABLES. */
/* ---------------- */
/* ***********************************************************************
*/
/* Parameter adjustments */
--marker;
--xlnz;
--snode;
--colcnt;
--invp;
--perm;
--xadj;
--adjncy;
--mrglnk;
--xlindx;
--xsuper;
--lindx;
/* Function Body */
*flag__ = 0;
if (*neqns <= 0) {
return 0;
}
/* --------------------------------------------------- */
/* INITIALIZATIONS ... */
/* NZEND : POINTS TO THE LAST USED SLOT IN LINDX. */
/* TAIL : END OF LIST INDICATOR */
/* (IN RCHLNK(*), NOT MRGLNK(*)). */
/* MRGLNK : CREATE EMPTY LISTS. */
/* MARKER : "UNMARK" THE INDICES. */
/* --------------------------------------------------- */
nzend = 0;
head = 0;
tail = *neqns + 1;
point = 1;
i__1 = *neqns;
for (i__ = 1; i__ <= i__1; ++i__) {
marker[i__] = 0;
xlnz[i__] = point;
point += colcnt[i__];
/* L50: */
}
xlnz[*neqns + 1] = point;
point = 1;
i__1 = *nsuper;
for (ksup = 1; ksup <= i__1; ++ksup) {
mrglnk[ksup] = 0;
fstcol = xsuper[ksup];
xlindx[ksup] = point;
point += colcnt[fstcol];
/* L100: */
}
xlindx[*nsuper + 1] = point;
/* --------------------------- */
/* FOR EACH SUPERNODE KSUP ... */
/* --------------------------- */
i__1 = *nsuper;
for (ksup = 1; ksup <= i__1; ++ksup) {
/* ---------------------------------------------------------
*/
/* INITIALIZATIONS ... */
/* FSTCOL : FIRST COLUMN OF SUPERNODE KSUP. */
/* LSTCOL : LAST COLUMN OF SUPERNODE KSUP. */
/* KNZ : WILL COUNT THE NONZEROS OF L IN COLUMN KCOL.
*/
/* RCHLNK : INITIALIZE EMPTY INDEX LIST FOR KCOL. */
/* ---------------------------------------------------------
*/
fstcol = xsuper[ksup];
lstcol = xsuper[ksup + 1] - 1;
width = lstcol - fstcol + 1;
length = colcnt[fstcol];
knz = 0;
rchlnk[head] = tail;
jsup = mrglnk[ksup];
/* ------------------------------------------------- */
/* IF KSUP HAS CHILDREN IN THE SUPERNODAL E-TREE ... */
/* ------------------------------------------------- */
if (jsup > 0) {
/* --------------------------------------------- */
/* COPY THE INDICES OF THE FIRST CHILD JSUP INTO */
/* THE LINKED LIST, AND MARK EACH WITH THE VALUE */
/* KSUP. */
/* --------------------------------------------- */
jwidth = xsuper[jsup + 1] - xsuper[jsup];
jnzbeg = xlindx[jsup] + jwidth;
jnzend = xlindx[jsup + 1] - 1;
i__2 = jnzbeg;
for (jptr = jnzend; jptr >= i__2; --jptr) {
newi = lindx[jptr];
++knz;
marker[newi] = ksup;
rchlnk[newi] = rchlnk[head];
rchlnk[head] = newi;
/* L200: */
}
/* ------------------------------------------ */
/* FOR EACH SUBSEQUENT CHILD JSUP OF KSUP ... */
/* ------------------------------------------ */
jsup = mrglnk[jsup];
L300:
if (jsup != 0 && knz < length) {
/* ----------------------------------------
*/
/* MERGE THE INDICES OF JSUP INTO THE LIST,
*/
/* AND MARK NEW INDICES WITH VALUE KSUP. */
/* ----------------------------------------
*/
jwidth = xsuper[jsup + 1] - xsuper[jsup];
jnzbeg = xlindx[jsup] + jwidth;
jnzend = xlindx[jsup + 1] - 1;
nexti = head;
i__2 = jnzend;
for (jptr = jnzbeg; jptr <= i__2; ++jptr) {
newi = lindx[jptr];
L400:
i__ = nexti;
nexti = rchlnk[i__];
if (newi > nexti) {
goto L400;
}
if (newi < nexti) {
++knz;
rchlnk[i__] = newi;
rchlnk[newi] = nexti;
marker[newi] = ksup;
nexti = newi;
}
/* L500: */
}
jsup = mrglnk[jsup];
goto L300;
}
}
/* --------------------------------------------------- */
/* STRUCTURE OF A(*,FSTCOL) HAS NOT BEEN EXAMINED YET. */
/* "SORT" ITS STRUCTURE INTO THE LINKED LIST, */
/* INSERTING ONLY THOSE INDICES NOT ALREADY IN THE */
/* LIST. */
/* --------------------------------------------------- */
if (knz < length) {
node = perm[fstcol];
knzbeg = xadj[node];
knzend = xadj[node + 1] - 1;
i__2 = knzend;
for (kptr = knzbeg; kptr <= i__2; ++kptr) {
newi = adjncy[kptr];
newi = invp[newi];
if (newi > fstcol && marker[newi] != ksup) {
/* --------------------------------
*/
/* POSITION AND INSERT NEWI IN LIST
*/
/* AND MARK IT WITH KCOL. */
/* --------------------------------
*/
nexti = head;
L600:
i__ = nexti;
nexti = rchlnk[i__];
if (newi > nexti) {
goto L600;
}
++knz;
rchlnk[i__] = newi;
rchlnk[newi] = nexti;
marker[newi] = ksup;
}
/* L700: */
}
}
/* --------------------------------------------------------
---- */
/* IF KSUP HAS NO CHILDREN, INSERT FSTCOL INTO THE LINKED LI
ST. */
/* --------------------------------------------------------
---- */
if (rchlnk[head] != fstcol) {
rchlnk[fstcol] = rchlnk[head];
rchlnk[head] = fstcol;
++knz;
}
/* -------------------------------------------- */
/* COPY INDICES FROM LINKED LIST INTO LINDX(*). */
/* -------------------------------------------- */
nzbeg = nzend + 1;
nzend += knz;
if (nzend + 1 != xlindx[ksup + 1]) {
goto L8000;
}
i__ = head;
i__2 = nzend;
for (kptr = nzbeg; kptr <= i__2; ++kptr) {
i__ = rchlnk[i__];
lindx[kptr] = i__;
/* L800: */
}
/* --------------------------------------------------- */
/* IF KSUP HAS A PARENT, INSERT KSUP INTO ITS PARENT'S */
/* "MERGE" LIST. */
/* --------------------------------------------------- */
if (length > width) {
pcol = lindx[xlindx[ksup] + width];
psup = snode[pcol];
mrglnk[ksup] = mrglnk[psup];
mrglnk[psup] = ksup;
}
/* L1000: */
}
return 0;
/* ----------------------------------------------- */
/* INCONSISTENCY IN DATA STRUCTURE WAS DISCOVERED. */
/* ----------------------------------------------- */
L8000:
*flag__ = -2;
return 0;
} /* symfc2_ */
#ifdef DO_BFINIT
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: December 27, 1994 */
/* Authors: Esmond G. Ng and Barry W. Peyton */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* ****** BFINIT ..... INITIALIZATION FOR BLOCK FACTORIZATION ****** */
/* *********************************************************************** */
/* *********************************************************************** */
/* PURPOSE: */
/* THIS SUBROUTINE COMPUTES ITEMS NEEDED BY THE LEFT-LOOKING */
/* BLOCK-TO-BLOCK CHOLESKY FACTORITZATION ROUTINE BLKFCT. */
/* INPUT PARAMETERS: */
/* NEQNS - NUMBER OF EQUATIONS. */
/* NSUPER - NUMBER OF SUPERNODES. */
/* XSUPER - INTEGER ARRAY OF SIZE (NSUPER+1) CONTAINING */
/* THE SUPERNODE PARTITIONING. */
/* SNODE - SUPERNODE MEMBERSHIP. */
/* (XLINDX,LINDX) - ARRAYS DESCRIBING THE SUPERNODAL STRUCTURE. */
/* CACHSZ - CACHE SIZE (IN KBYTES). */
/* OUTPUT PARAMETERS: */
/* TMPSIZ - SIZE OF WORKING STORAGE REQUIRED BY BLKFCT. */
/* SPLIT - SPLITTING OF SUPERNODES SO THAT THEY FIT */
/* INTO CACHE. */
/* *********************************************************************** */
/* Subroutine */ int bfinit_(neqns, nsuper, xsuper, snode, xlindx, lindx,
cachsz, tmpsiz, split)
integer *neqns, *nsuper, *xsuper, *snode, *xlindx, *lindx, *cachsz, *tmpsiz, *
split;
{
extern /* Subroutine */ int fnsplt_(), fntsiz_();
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* --------------------------------------------------- */
/* DETERMINE FLOATING POINT WORKING SPACE REQUIREMENT. */
/* --------------------------------------------------- */
/* Parameter adjustments */
--split;
--lindx;
--xlindx;
--snode;
--xsuper;
/* Function Body */
fntsiz_(nsuper, &xsuper[1], &snode[1], &xlindx[1], &lindx[1], tmpsiz);
/* ------------------------------- */
/* PARTITION SUPERNODES FOR CACHE. */
/* ------------------------------- */
fnsplt_(neqns, nsuper, &xsuper[1], &xlindx[1], cachsz, &split[1]);
return 0;
} /* bfinit_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: December 27, 1994 */
/* Authors: Esmond G. Ng and Barry W. Peyton */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* ****** FNTSIZ ..... COMPUTE WORK STORAGE SIZE FOR BLKFCT ****** */
/* *********************************************************************** */
/* *********************************************************************** */
/* PURPOSE: */
/* THIS SUBROUTINE DETERMINES THE SIZE OF THE WORKING STORAGE */
/* REQUIRED BY BLKFCT. */
/* INPUT PARAMETERS: */
/* NSUPER - NUMBER OF SUPERNODES. */
/* XSUPER - INTEGER ARRAY OF SIZE (NSUPER+1) CONTAINING */
/* THE SUPERNODE PARTITIONING. */
/* SNODE - SUPERNODE MEMBERSHIP. */
/* (XLINDX,LINDX) - ARRAYS DESCRIBING THE SUPERNODAL STRUCTURE. */
/* OUTPUT PARAMETERS: */
/* TMPSIZ - SIZE OF WORKING STORAGE REQUIRED BY BLKFCT. */
/* *********************************************************************** */
/* Subroutine */ int fntsiz_(nsuper, xsuper, snode, xlindx, lindx, tmpsiz)
integer *nsuper, *xsuper, *snode, *xlindx, *lindx, *tmpsiz;
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer iend, clen, ksup, i__, bound, ncols, width, tsize, ibegin,
length, cursup, nxtsup;
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* RETURNS SIZE OF TEMP ARRAY USED BY BLKFCT FACTORIZATION ROUTINE.
*/
/* NOTE THAT THE VALUE RETURNED IS AN ESTIMATE, THOUGH IT IS USUALLY
*/
/* TIGHT. */
/* ---------------------------------------- */
/* COMPUTE SIZE OF TEMPORARY STORAGE VECTOR */
/* NEEDED BY BLKFCT. */
/* ---------------------------------------- */
/* Parameter adjustments */
--lindx;
--xlindx;
--snode;
--xsuper;
/* Function Body */
*tmpsiz = 0;
for (ksup = *nsuper; ksup >= 1; --ksup) {
ncols = xsuper[ksup + 1] - xsuper[ksup];
ibegin = xlindx[ksup] + ncols;
iend = xlindx[ksup + 1] - 1;
length = iend - ibegin + 1;
bound = length * (length + 1) / 2;
if (bound > *tmpsiz) {
cursup = snode[lindx[ibegin]];
clen = xlindx[cursup + 1] - xlindx[cursup];
width = 0;
i__1 = iend;
for (i__ = ibegin; i__ <= i__1; ++i__) {
nxtsup = snode[lindx[i__]];
if (nxtsup == cursup) {
++width;
if (i__ == iend) {
if (clen > length) {
tsize = length * width - (width - 1) * width / 2;
*tmpsiz = max(tsize,*tmpsiz);
}
}
} else {
if (clen > length) {
tsize = length * width - (width - 1) * width / 2;
*tmpsiz = max(tsize,*tmpsiz);
}
length -= width;
bound = length * (length + 1) / 2;
if (bound <= *tmpsiz) {
goto L500;
}
width = 1;
cursup = nxtsup;
clen = xlindx[cursup + 1] - xlindx[cursup];
}
/* L400: */
}
}
L500:
;
}
return 0;
} /* fntsiz_ */
/* *********************************************************************** */
/* *********************************************************************** */
/* Version: 0.3 */
/* Last modified: December 27, 1994 */
/* Authors: Esmond G. Ng and Barry W. Peyton */
/* Mathematical Sciences Section, Oak Ridge National Laboratory */
/* *********************************************************************** */
/* *********************************************************************** */
/* **** FNSPLT ..... COMPUTE FINE PARTITIONING OF SUPERNODES ***** */
/* *********************************************************************** */
/* *********************************************************************** */
/* PURPOSE: */
/* THIS SUBROUTINE DETERMINES A FINE PARTITIONING OF SUPERNODES */
/* WHEN THERE IS A CACHE AVAILABLE ON THE MACHINE. THE FINE */
/* PARTITIONING IS CHOSEN SO THAT DATA RE-USE IS MAXIMIZED. */
/* INPUT PARAMETERS: */
/* NEQNS - NUMBER OF EQUATIONS. */
/* NSUPER - NUMBER OF SUPERNODES. */
/* XSUPER - INTEGER ARRAY OF SIZE (NSUPER+1) CONTAINING */
/* THE SUPERNODE PARTITIONING. */
/* XLINDX - INTEGER ARRAY OF SIZE (NSUPER+1) CONTAINING */
/* POINTERS IN THE SUPERNODE INDICES. */
/* CACHSZ - CACHE SIZE IN KILO BYTES. */
/* IF THERE IS NO CACHE, SET CACHSZ = 0. */
/* OUTPUT PARAMETERS: */
/* SPLIT - INTEGER ARRAY OF SIZE NEQNS CONTAINING THE */
/* FINE PARTITIONING. */
/* *********************************************************************** */
/* Subroutine */ int fnsplt_(neqns, nsuper, xsuper, xlindx, cachsz, split)
integer *neqns, *nsuper, *xsuper, *xlindx, *cachsz, *split;
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer kcol, used, ksup, cache, ncols, width, height, curcol,
fstcol, lstcol, nxtblk;
/* ***********************************************************************
*/
/* ----------- */
/* PARAMETERS. */
/* ----------- */
/* ---------------- */
/* LOCAL VARIABLES. */
/* ---------------- */
/* ******************************************************************* */
/* -------------------------------------------- */
/* COMPUTE THE NUMBER OF 8-BYTE WORDS IN CACHE. */
/* -------------------------------------------- */
/* Parameter adjustments */
--split;
--xlindx;
--xsuper;
/* Function Body */
if (*cachsz <= 0) {
cache = 2000000000;
} else {
cache = (float) (*cachsz) * (float)1024. / (float)8. * (float).9;
}
/* --------------- */
/* INITIALIZATION. */
/* --------------- */
i__1 = *neqns;
for (kcol = 1; kcol <= i__1; ++kcol) {
split[kcol] = 0;
/* L100: */
}
/* --------------------------- */
/* FOR EACH SUPERNODE KSUP ... */
/* --------------------------- */
i__1 = *nsuper;
for (ksup = 1; ksup <= i__1; ++ksup) {
/* ----------------------- */
/* ... GET SUPERNODE INFO. */
/* ----------------------- */
height = xlindx[ksup + 1] - xlindx[ksup];
fstcol = xsuper[ksup];
lstcol = xsuper[ksup + 1] - 1;
width = lstcol - fstcol + 1;
nxtblk = fstcol;
/* -------------------------------------- */
/* ... UNTIL ALL COLUMNS OF THE SUPERNODE */
/* HAVE BEEN PROCESSED ... */
/* -------------------------------------- */
curcol = fstcol - 1;
L200:
/* ------------------------------------------- */
/* ... PLACE THE FIRST COLUMN(S) IN THE CACHE. */
/* ------------------------------------------- */
++curcol;
if (curcol < lstcol) {
++curcol;
ncols = 2;
used = height * 3 - 1;
height += -2;
} else {
ncols = 1;
used = height << 1;
--height;
}
/* -------------------------------------- */
/* ... WHILE THE CACHE IS NOT FILLED AND */
/* THERE ARE COLUMNS OF THE SUPERNODE */
/* REMAINING TO BE PROCESSED ... */
/* -------------------------------------- */
L300:
if (used + height < cache && curcol < lstcol) {
/* -------------------------------- */
/* ... ADD ANOTHER COLUMN TO CACHE. */
/* -------------------------------- */
++curcol;
++ncols;
used += height;
--height;
goto L300;
}
/* ------------------------------------- */
/* ... RECORD THE NUMBER OF COLUMNS THAT */
/* FILLED THE CACHE. */
/* ------------------------------------- */
split[nxtblk] = ncols;
++nxtblk;
/* -------------------------- */
/* ... GO PROCESS NEXT BLOCK. */
/* -------------------------- */
if (curcol < lstcol) {
goto L200;
}
/* L1000: */
}
return 0;
} /* fnsplt_ */
#endif

Event Timeline

c4science · Help