diff --git a/appveyor.yml b/appveyor.yml index d9ef393..258b330 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -17,13 +17,13 @@ matrix: environment: matrix: - - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2013 - GENERATOR: "Visual Studio 10 2010" - GROUP: "SUPPORTED" - BOOST_ROOT: "C:/Libraries/boost" - # New build on Visual Studio 15 2017 - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017 GENERATOR: "Visual Studio 15 2017" + GROUP: "SUPPORTED" + BOOST_ROOT: "C:/Libraries/boost_1_67_0" + # New build on Visual Studio 15 2017 + - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017 + GENERATOR: "Visual Studio 15 2017 Win64" GROUP: "EXPERIMENTAL" BOOST_ROOT: "C:/Libraries/boost_1_67_0" diff --git a/src/funcs.h b/src/funcs.h index b65f194..c8dfb05 100755 --- a/src/funcs.h +++ b/src/funcs.h @@ -193,24 +193,8 @@ void matrixcoeffs(EN_Project *pr); /* Finds hyd. matrix coeffs. /* ----------- SMATRIX.C ---------------*/ int createsparse(EN_Project *pr); /* Creates sparse matrix */ -int allocsparse(EN_Project *pr); /* Allocates matrix memory */ void freesparse(EN_Project *pr); /* Frees matrix memory */ -int buildlists(EN_Project *pr, int); /* Builds adjacency lists */ -int paralink(EN_Project *pr, int, int, int); /* Checks for parallel links */ -void xparalinks(EN_Project *pr); /* Removes parallel links */ -void freelists(EN_Project *pr); /* Frees adjacency lists */ -void countdegree(EN_Project *pr); /* Counts links at each node */ -int reordernodes(EN_Project *pr); /* Finds a node re-ordering */ -int mindegree(solver_t *s, int, int); /* Finds min. degree node */ -int growlist(EN_Project *pr, int); /* Augments adjacency list */ -int newlink(EN_Project *pr, Padjlist); /* Adds fill-ins for a node */ -int linked(EN_Network *net, int, int); /* Checks if 2 nodes linked */ -int addlink(EN_Network *net, int, int, int); /* Creates new fill-in */ -int storesparse(EN_Project *pr, int); /* Stores sparse matrix */ -int ordersparse(hydraulics_t *h, int); /* Orders matrix storage */ -void transpose(int,int *,int *, /* Transposes sparse matrix */ - int *,int *,int *,int *,int *); -int linsolve(solver_t *s, int); /* Solves set of inear eqns. */ +int linsolve(EN_Project *pr, int); /* Solves set of linear eqns. */ /* ----------- QUALITY.C ---------------*/ int openqual(EN_Project *pr); /* Opens WQ solver system */ diff --git a/src/genmmd.c b/src/genmmd.c new file mode 100644 index 0000000..2dccb2b --- /dev/null +++ b/src/genmmd.c @@ -0,0 +1,1001 @@ +/* MULTIPLE MINIMUM DEGREE ROW RE-ORDERING ALGORITHM + * + * Modified to work with Fortran-style arrays. + * + */ + + +#include + +int genmmd(int* neqns, int* xadj, int* adjncy, int* invp, int* perm, + int* delta, int* dhead, int* qsize, int* llist, int* marker, + int* maxint, int* nofsub); + +int mmdint_(int* neqns, int* xadj, int* adjncy, int* dhead, int* dforw, + int* dbakw, int* qsize, int* llist, int* marker); +int mmdelm_(int* mdnode, int* xadj, int* adjncy, int* dhead, int* dforw, + int* dbakw, int* qsize, int* llist, int* marker, + int* maxint, int* tag); +int mmdupd_(int* ehead, int* neqns, int* xadj, int* adjncy, int* delta, + int* mdeg, int* dhead, int* dforw, int* dbakw, int* qsize, + int* llist, int* marker, int* maxint, int* tag); +int mmdnum_(int* neqns, int* perm, int* invp, int* qsize); + +//============================================================================= + +/* genmmd.f -- translated by f2c (version of 23 April 1993 18:34:30). + You must link the resulting object file with the libraries: + -lf2c -lm (in that order) +*/ + +//#include "f2c.h" + +/* Sivan: I modified INTEGER*2 -> INTEGER*4 */ +/* *************************************************************** */ +/* *************************************************************** */ +/* **** GENMMD ..... MULTIPLE MINIMUM EXTERNAL DEGREE **** */ +/* *************************************************************** */ +/* *************************************************************** */ + +/* AUTHOR - JOSEPH W.H. LIU */ +/* DEPT OF COMPUTER SCIENCE, YORK UNIVERSITY. */ + +/* PURPOSE - THIS ROUTINE IMPLEMENTS THE MINIMUM DEGREE */ +/* ALGORITHM. IT MAKES USE OF THE IMPLICIT REPRESENTATION */ +/* OF ELIMINATION GRAPHS BY QUOTIENT GRAPHS, AND THE */ +/* NOTION OF INDISTINGUISHABLE NODES. IT ALSO IMPLEMENTS */ +/* THE MODIFICATIONS BY MULTIPLE ELIMINATION AND MINIMUM */ +/* EXTERNAL DEGREE. */ +/* --------------------------------------------- */ +/* CAUTION - THE ADJACENCY VECTOR ADJNCY WILL BE */ +/* DESTROYED. */ +/* --------------------------------------------- */ + +/* INPUT PARAMETERS - */ +/* NEQNS - NUMBER OF EQUATIONS. */ +/* (XADJ,ADJNCY) - THE ADJACENCY STRUCTURE. */ +/* DELTA - TOLERANCE VALUE FOR MULTIPLE ELIMINATION. */ +/* MAXINT - MAXIMUM MACHINE REPRESENTABLE (SHORT) INTEGER */ +/* (ANY SMALLER ESTIMATE WILL DO) FOR MARKING */ +/* NODES. */ + +/* OUTPUT PARAMETERS - */ +/* PERM - THE MINIMUM DEGREE ORDERING. */ +/* INVP - THE INVERSE OF PERM. */ +/* NOFSUB - AN UPPER BOUND ON THE NUMBER OF NONZERO */ +/* SUBSCRIPTS FOR THE COMPRESSED STORAGE SCHEME. */ + +/* WORKING PARAMETERS - */ +/* DHEAD - VECTOR FOR HEAD OF DEGREE LISTS. */ +/* INVP - USED TEMPORARILY FOR DEGREE FORWARD LINK. */ +/* PERM - USED TEMPORARILY FOR DEGREE BACKWARD LINK. */ +/* QSIZE - VECTOR FOR SIZE OF SUPERNODES. */ +/* LLIST - VECTOR FOR TEMPORARY LINKED LISTS. */ +/* MARKER - A TEMPORARY MARKER VECTOR. */ + +/* PROGRAM SUBROUTINES - */ +/* MMDELM, MMDINT, MMDNUM, MMDUPD. */ + +/* *************************************************************** */ +int genmmd(int* neqns, int* xadj, int* adjncy, int* invp, int* perm, + int* delta, int* dhead, int* qsize, int* llist, int* marker, + int* maxint, int* nofsub) +{ + /* System generated locals */ + int i__1; + + /* Local variables */ + static int mdeg, ehead, i, mdlmt, mdnode; + //extern /* Subroutine */ int mmdelm_(), mmdupd_(), mmdint_(), mmdnum_(); + static int nextmd, tag, num; + + +/* *************************************************************** */ + +/* INTEGER*2 ADJNCY(1), DHEAD(1) , INVP(1) , LLIST(1) , */ +/* 1 MARKER(1), PERM(1) , QSIZE(1) */ + +/* *************************************************************** */ + + /* Parameter adjustments */ + //--marker; --llist; --qsize; --dhead; --perm; --invp; --adjncy; --xadj; + + /* Function Body */ + if (*neqns <= 0) { + return 0; + } + +/* ------------------------------------------------ */ +/* INITIALIZATION FOR THE MINIMUM DEGREE ALGORITHM. */ +/* ------------------------------------------------ */ + *nofsub = 0; + //mmdint_(neqns, &xadj[1], &adjncy[1], &dhead[1], &invp[1], &perm[1], + // &qsize[1], &llist[1], &marker[1]); + mmdint_(neqns, xadj, adjncy, dhead, invp, perm, qsize, llist, marker); + +/* ---------------------------------------------- */ +/* NUM COUNTS THE NUMBER OF ORDERED NODES PLUS 1. */ +/* ---------------------------------------------- */ + num = 1; + +/* ----------------------------- */ +/* ELIMINATE ALL ISOLATED NODES. */ +/* ----------------------------- */ + nextmd = dhead[1]; +L100: + if (nextmd <= 0) { + goto L200; + } + mdnode = nextmd; + nextmd = invp[mdnode]; + marker[mdnode] = *maxint; + invp[mdnode] = -num; + ++num; + goto L100; + +L200: +/* ---------------------------------------- */ +/* SEARCH FOR NODE OF THE MINIMUM DEGREE. */ +/* MDEG IS THE CURRENT MINIMUM DEGREE; */ +/* TAG IS USED TO FACILITATE MARKING NODES. */ +/* ---------------------------------------- */ + if (num > *neqns) { + goto L1000; + } + tag = 1; + dhead[1] = 0; + mdeg = 2; +L300: + if (dhead[mdeg] > 0) { + goto L400; + } + ++mdeg; + goto L300; +L400: +/* ------------------------------------------------- */ +/* USE VALUE OF DELTA TO SET UP MDLMT, WHICH GOVERNS */ +/* WHEN A DEGREE UPDATE IS TO BE PERFORMED. */ +/* ------------------------------------------------- */ + mdlmt = mdeg + *delta; + ehead = 0; + +L500: + mdnode = dhead[mdeg]; + if (mdnode > 0) { + goto L600; + } + ++mdeg; + if (mdeg > mdlmt) { + goto L900; + } + goto L500; +L600: +/* ---------------------------------------- */ +/* REMOVE MDNODE FROM THE DEGREE STRUCTURE. */ +/* ---------------------------------------- */ + nextmd = invp[mdnode]; + dhead[mdeg] = nextmd; + if (nextmd > 0) { + perm[nextmd] = -mdeg; + } + invp[mdnode] = -num; + *nofsub = *nofsub + mdeg + qsize[mdnode] - 2; + if (num + qsize[mdnode] > *neqns) { + goto L1000; + } +/* ---------------------------------------------- */ +/* ELIMINATE MDNODE AND PERFORM QUOTIENT GRAPH */ +/* TRANSFORMATION. RESET TAG VALUE IF NECESSARY. */ +/* ---------------------------------------------- */ + ++tag; + if (tag < *maxint) { + goto L800; + } + tag = 1; + i__1 = *neqns; + for (i = 1; i <= i__1; ++i) { + if (marker[i] < *maxint) { + marker[i] = 0; + } +/* L700: */ + } +L800: + //mmdelm_(&mdnode, &xadj[1], &adjncy[1], &dhead[1], &invp[1], &perm[1], + // &qsize[1], &llist[1], &marker[1], maxint, &tag); + mmdelm_(&mdnode, xadj, adjncy, dhead, invp, perm, qsize, llist, + marker, maxint, &tag); + num += qsize[mdnode]; + llist[mdnode] = ehead; + ehead = mdnode; + if (*delta >= 0) { + goto L500; + } +L900: +/* ------------------------------------------- */ +/* UPDATE DEGREES OF THE NODES INVOLVED IN THE */ +/* MINIMUM DEGREE NODES ELIMINATION. */ +/* ------------------------------------------- */ + if (num > *neqns) { + goto L1000; + } + //mmdupd_(&ehead, neqns, &xadj[1], &adjncy[1], delta, &mdeg, &dhead[1], + // &invp[1], &perm[1], &qsize[1], &llist[1], &marker[1], maxint, &tag); + mmdupd_(&ehead, neqns, xadj, adjncy, delta, &mdeg, dhead, + invp, perm, qsize, llist, marker, maxint, &tag); + goto L300; + +L1000: + //mmdnum_(neqns, &perm[1], &invp[1], &qsize[1]); + mmdnum_(neqns, perm, invp, qsize); + + //++marker; ++llist; ++qsize; ++dhead; ++perm; ++invp; ++adjncy; ++xadj; + return 0; + +} /* genmmd_ */ + +/* *************************************************************** */ +/* *************************************************************** */ +/* *** MMDINT ..... MULT MINIMUM DEGREE INITIALIZATION *** */ +/* *************************************************************** */ +/* *************************************************************** */ + +/* AUTHOR - JOSEPH W.H. LIU */ +/* DEPT OF COMPUTER SCIENCE, YORK UNIVERSITY. */ + +/* PURPOSE - THIS ROUTINE PERFORMS INITIALIZATION FOR THE */ +/* MULTIPLE ELIMINATION VERSION OF THE MINIMUM DEGREE */ +/* ALGORITHM. */ + +/* INPUT PARAMETERS - */ +/* NEQNS - NUMBER OF EQUATIONS. */ +/* (XADJ,ADJNCY) - ADJACENCY STRUCTURE. */ + +/* OUTPUT PARAMETERS - */ +/* (DHEAD,DFORW,DBAKW) - DEGREE DOUBLY LINKED STRUCTURE. */ +/* QSIZE - SIZE OF SUPERNODE (INITIALIZED TO ONE). */ +/* LLIST - LINKED LIST. */ +/* MARKER - MARKER VECTOR. */ + +/* *************************************************************** */ + + +int mmdint_(int* neqns, int* xadj, int* adjncy, int* dhead, int* dforw, + int* dbakw, int* qsize, int* llist, int* marker) +{ + /* System generated locals */ + int i__1; + + /* Local variables */ + static int ndeg, node, fnode; + + +/* *************************************************************** */ + +/* INTEGER*2 ADJNCY(1), DBAKW(1) , DFORW(1) , DHEAD(1) , */ +/* 1 LLIST(1) , MARKER(1), QSIZE(1) */ + +/* *************************************************************** */ + + /* Parameter adjustments */ + //--marker; --llist; --qsize; --dbakw; --dforw; --dhead; --adjncy; --xadj; + + /* Function Body */ + i__1 = *neqns; + for (node = 1; node <= i__1; ++node) { + dhead[node] = 0; + qsize[node] = 1; + marker[node] = 0; + llist[node] = 0; +/* L100: */ + } +/* ------------------------------------------ */ +/* INITIALIZE THE DEGREE DOUBLY LINKED LISTS. */ +/* ------------------------------------------ */ + i__1 = *neqns; + for (node = 1; node <= i__1; ++node) { + ndeg = xadj[node + 1] - xadj[node] + 1; + fnode = dhead[ndeg]; + dforw[node] = fnode; + dhead[ndeg] = node; + if (fnode > 0) { + dbakw[fnode] = node; + } + dbakw[node] = -ndeg; +/* L200: */ + } + return 0; + +} /* mmdint_ */ + +/* *************************************************************** */ +/* *************************************************************** */ +/* ** MMDELM ..... MULTIPLE MINIMUM DEGREE ELIMINATION *** */ +/* *************************************************************** */ +/* *************************************************************** */ + +/* AUTHOR - JOSEPH W.H. LIU */ +/* DEPT OF COMPUTER SCIENCE, YORK UNIVERSITY. */ + +/* PURPOSE - THIS ROUTINE ELIMINATES THE NODE MDNODE OF */ +/* MINIMUM DEGREE FROM THE ADJACENCY STRUCTURE, WHICH */ +/* IS STORED IN THE QUOTIENT GRAPH FORMAT. IT ALSO */ +/* TRANSFORMS THE QUOTIENT GRAPH REPRESENTATION OF THE */ +/* ELIMINATION GRAPH. */ + +/* INPUT PARAMETERS - */ +/* MDNODE - NODE OF MINIMUM DEGREE. */ +/* MAXINT - ESTIMATE OF MAXIMUM REPRESENTABLE (SHORT) */ +/* INTEGER. */ +/* TAG - TAG VALUE. */ + +/* UPDATED PARAMETERS - */ +/* (XADJ,ADJNCY) - UPDATED ADJACENCY STRUCTURE. */ +/* (DHEAD,DFORW,DBAKW) - DEGREE DOUBLY LINKED STRUCTURE. */ +/* QSIZE - SIZE OF SUPERNODE. */ +/* MARKER - MARKER VECTOR. */ +/* LLIST - TEMPORARY LINKED LIST OF ELIMINATED NABORS. */ + +/* *************************************************************** */ + +int mmdelm_(int* mdnode, int* xadj, int* adjncy, int* dhead, int* dforw, + int* dbakw, int* qsize, int* llist, int* marker, + int* maxint, int* tag) +{ + /* System generated locals */ + int i__1, i__2; + + /* Local variables */ + static int node, link, rloc, rlmt, i, j, nabor, rnode, elmnt, xqnbr, + istop, jstop, istrt, jstrt, nxnode, pvnode, nqnbrs, npv; + + +/* *************************************************************** */ + +/* INTEGER*2 ADJNCY(1), DBAKW(1) , DFORW(1) , DHEAD(1) , */ +/* 1 LLIST(1) , MARKER(1), QSIZE(1) */ + +/* *************************************************************** */ + +/* ----------------------------------------------- */ +/* FIND REACHABLE SET AND PLACE IN DATA STRUCTURE. */ +/* ----------------------------------------------- */ + /* Parameter adjustments */ + //--marker; --llist; --qsize; --dbakw; --dforw; --dhead; + //--adjncy; --xadj; + + /* Function Body */ + marker[*mdnode] = *tag; + istrt = xadj[*mdnode]; + istop = xadj[*mdnode + 1] - 1; +/* ------------------------------------------------------- */ +/* ELMNT POINTS TO THE BEGINNING OF THE LIST OF ELIMINATED */ +/* NABORS OF MDNODE, AND RLOC GIVES THE STORAGE LOCATION */ +/* FOR THE NEXT REACHABLE NODE. */ +/* ------------------------------------------------------- */ + elmnt = 0; + rloc = istrt; + rlmt = istop; + i__1 = istop; + for (i = istrt; i <= i__1; ++i) { + nabor = adjncy[i]; + if (nabor == 0) { + goto L300; + } + if (marker[nabor] >= *tag) { + goto L200; + } + marker[nabor] = *tag; + if (dforw[nabor] < 0) { + goto L100; + } + adjncy[rloc] = nabor; + ++rloc; + goto L200; +L100: + llist[nabor] = elmnt; + elmnt = nabor; +L200: + ; + } +L300: +/* ----------------------------------------------------- */ +/* MERGE WITH REACHABLE NODES FROM GENERALIZED ELEMENTS. */ +/* ----------------------------------------------------- */ + if (elmnt <= 0) { + goto L1000; + } + adjncy[rlmt] = -elmnt; + link = elmnt; +L400: + jstrt = xadj[link]; + jstop = xadj[link + 1] - 1; + i__1 = jstop; + for (j = jstrt; j <= i__1; ++j) { + node = adjncy[j]; + link = -node; + if (node < 0) { + goto L400; + } else if (node == 0) { + goto L900; + } else { + goto L500; + } +L500: + if (marker[node] >= *tag || dforw[node] < 0) { + goto L800; + } + marker[node] = *tag; +/* --------------------------------- */ +/* USE STORAGE FROM ELIMINATED NODES */ +/* IF NECESSARY. */ +/* --------------------------------- */ +L600: + if (rloc < rlmt) { + goto L700; + } + link = -adjncy[rlmt]; + rloc = xadj[link]; + rlmt = xadj[link + 1] - 1; + goto L600; +L700: + adjncy[rloc] = node; + ++rloc; +L800: + ; + } +L900: + elmnt = llist[elmnt]; + goto L300; +L1000: + if (rloc <= rlmt) { + adjncy[rloc] = 0; + } +/* -------------------------------------------------------- */ +/* FOR EACH NODE IN THE REACHABLE SET, DO THE FOLLOWING ... */ +/* -------------------------------------------------------- */ + link = *mdnode; +L1100: + istrt = xadj[link]; + istop = xadj[link + 1] - 1; + i__1 = istop; + for (i = istrt; i <= i__1; ++i) { + rnode = adjncy[i]; + link = -rnode; + if (rnode < 0) { + goto L1100; + } else if (rnode == 0) { + goto L1800; + } else { + goto L1200; + } +L1200: +/* -------------------------------------------- */ +/* IF RNODE IS IN THE DEGREE LIST STRUCTURE ... */ +/* -------------------------------------------- */ + pvnode = dbakw[rnode]; + if (pvnode == 0 || pvnode == -(*maxint)) { + goto L1300; + } +/* ------------------------------------- */ +/* THEN REMOVE RNODE FROM THE STRUCTURE. */ +/* ------------------------------------- */ + nxnode = dforw[rnode]; + if (nxnode > 0) { + dbakw[nxnode] = pvnode; + } + if (pvnode > 0) { + dforw[pvnode] = nxnode; + } + npv = -pvnode; + if (pvnode < 0) { + dhead[npv] = nxnode; + } +L1300: +/* ---------------------------------------- */ +/* PURGE INACTIVE QUOTIENT NABORS OF RNODE. */ +/* ---------------------------------------- */ + jstrt = xadj[rnode]; + jstop = xadj[rnode + 1] - 1; + xqnbr = jstrt; + i__2 = jstop; + for (j = jstrt; j <= i__2; ++j) { + nabor = adjncy[j]; + if (nabor == 0) { + goto L1500; + } + if (marker[nabor] >= *tag) { + goto L1400; + } + adjncy[xqnbr] = nabor; + ++xqnbr; +L1400: + ; + } +L1500: +/* ---------------------------------------- */ +/* IF NO ACTIVE NABOR AFTER THE PURGING ... */ +/* ---------------------------------------- */ + nqnbrs = xqnbr - jstrt; + if (nqnbrs > 0) { + goto L1600; + } +/* ----------------------------- */ +/* THEN MERGE RNODE WITH MDNODE. */ +/* ----------------------------- */ + qsize[*mdnode] += qsize[rnode]; + qsize[rnode] = 0; + marker[rnode] = *maxint; + dforw[rnode] = -(*mdnode); + dbakw[rnode] = -(*maxint); + goto L1700; +L1600: +/* -------------------------------------- */ +/* ELSE FLAG RNODE FOR DEGREE UPDATE, AND */ +/* ADD MDNODE AS A NABOR OF RNODE. */ +/* -------------------------------------- */ + dforw[rnode] = nqnbrs + 1; + dbakw[rnode] = 0; + adjncy[xqnbr] = *mdnode; + ++xqnbr; + if (xqnbr <= jstop) { + adjncy[xqnbr] = 0; + } + +L1700: + ; + } +L1800: + return 0; + +} /* mmdelm_ */ + +/* *************************************************************** */ +/* *************************************************************** */ +/* ***** MMDUPD ..... MULTIPLE MINIMUM DEGREE UPDATE ***** */ +/* *************************************************************** */ +/* *************************************************************** */ + +/* AUTHOR - JOSEPH W.H. LIU */ +/* DEPT OF COMPUTER SCIENCE, YORK UNIVERSITY. */ + +/* PURPOSE - THIS ROUTINE UPDATES THE DEGREES OF NODES */ +/* AFTER A MULTIPLE ELIMINATION STEP. */ + +/* INPUT PARAMETERS - */ +/* EHEAD - THE BEGINNING OF THE LIST OF ELIMINATED */ +/* NODES (I.E., NEWLY FORMED ELEMENTS). */ +/* NEQNS - NUMBER OF EQUATIONS. */ +/* (XADJ,ADJNCY) - ADJACENCY STRUCTURE. */ +/* DELTA - TOLERANCE VALUE FOR MULTIPLE ELIMINATION. */ +/* MAXINT - MAXIMUM MACHINE REPRESENTABLE (SHORT) */ +/* INTEGER. */ + +/* UPDATED PARAMETERS - */ +/* MDEG - NEW MINIMUM DEGREE AFTER DEGREE UPDATE. */ +/* (DHEAD,DFORW,DBAKW) - DEGREE DOUBLY LINKED STRUCTURE. */ +/* QSIZE - SIZE OF SUPERNODE. */ +/* LLIST - WORKING LINKED LIST. */ +/* MARKER - MARKER VECTOR FOR DEGREE UPDATE. */ +/* TAG - TAG VALUE. */ + +/* *************************************************************** */ + +int mmdupd_(int* ehead, int* neqns, int* xadj, int* adjncy, int* delta, + int* mdeg, int* dhead, int* dforw, int* dbakw, int* qsize, + int* llist, int* marker, int* maxint, int* tag) +{ + /* System generated locals */ + int i__1, i__2; + + /* Local variables */ + static int node, mtag, link, mdeg0, i, j, enode, fnode, nabor, elmnt, + istop, jstop, q2head, istrt, jstrt, qxhead, iq2, deg, deg0; + + +/* *************************************************************** */ + +/* INTEGER*2 ADJNCY(1), DBAKW(1) , DFORW(1) , DHEAD(1) , */ +/* 1 LLIST(1) , MARKER(1), QSIZE(1) */ + +/* *************************************************************** */ + + /* Parameter adjustments */ + //--marker; --llist; --qsize; --dbakw; --dforw; --dhead; + //--adjncy; --xadj; + + /* Function Body */ + mdeg0 = *mdeg + *delta; + elmnt = *ehead; +L100: +/* ------------------------------------------------------- */ +/* FOR EACH OF THE NEWLY FORMED ELEMENT, DO THE FOLLOWING. */ +/* (RESET TAG VALUE IF NECESSARY.) */ +/* ------------------------------------------------------- */ + if (elmnt <= 0) { + return 0; + } + mtag = *tag + mdeg0; + if (mtag < *maxint) { + goto L300; + } + *tag = 1; + i__1 = *neqns; + for (i = 1; i <= i__1; ++i) { + if (marker[i] < *maxint) { + marker[i] = 0; + } +/* L200: */ + } + mtag = *tag + mdeg0; +L300: +/* --------------------------------------------- */ +/* CREATE TWO LINKED LISTS FROM NODES ASSOCIATED */ +/* WITH ELMNT: ONE WITH TWO NABORS (Q2HEAD) IN */ +/* ADJACENCY STRUCTURE, AND THE OTHER WITH MORE */ +/* THAN TWO NABORS (QXHEAD). ALSO COMPUTE DEG0, */ +/* NUMBER OF NODES IN THIS ELEMENT. */ +/* --------------------------------------------- */ + q2head = 0; + qxhead = 0; + deg0 = 0; + link = elmnt; +L400: + istrt = xadj[link]; + istop = xadj[link + 1] - 1; + i__1 = istop; + for (i = istrt; i <= i__1; ++i) { + enode = adjncy[i]; + link = -enode; + if (enode < 0) { + goto L400; + } else if (enode == 0) { + goto L800; + } else { + goto L500; + } + +L500: + if (qsize[enode] == 0) { + goto L700; + } + deg0 += qsize[enode]; + marker[enode] = mtag; +/* ---------------------------------- */ +/* IF ENODE REQUIRES A DEGREE UPDATE, */ +/* THEN DO THE FOLLOWING. */ +/* ---------------------------------- */ + if (dbakw[enode] != 0) { + goto L700; + } +/* ---------------------------------------*/ +/* PLACE EITHER IN QXHEAD OR Q2HEAD LISTS.*/ +/* ---------------------------------------*/ + if (dforw[enode] == 2) { + goto L600; + } + llist[enode] = qxhead; + qxhead = enode; + goto L700; +L600: + llist[enode] = q2head; + q2head = enode; +L700: + ; + } +L800: +/* -------------------------------------------- */ +/* FOR EACH ENODE IN Q2 LIST, DO THE FOLLOWING. */ +/* -------------------------------------------- */ + enode = q2head; + iq2 = 1; +L900: + if (enode <= 0) { + goto L1500; + } + if (dbakw[enode] != 0) { + goto L2200; + } + ++(*tag); + deg = deg0; +/* ------------------------------------------ */ +/* IDENTIFY THE OTHER ADJACENT ELEMENT NABOR. */ +/* ------------------------------------------ */ + istrt = xadj[enode]; + nabor = adjncy[istrt]; + if (nabor == elmnt) { + nabor = adjncy[istrt + 1]; + } +/* ------------------------------------------------ */ +/* IF NABOR IS UNELIMINATED, INCREASE DEGREE COUNT. */ +/* ------------------------------------------------ */ + link = nabor; + if (dforw[nabor] < 0) { + goto L1000; + } + deg += qsize[nabor]; + goto L2100; +L1000: +/* -------------------------------------------- */ +/* OTHERWISE, FOR EACH NODE IN THE 2ND ELEMENT, */ +/* DO THE FOLLOWING. */ +/* -------------------------------------------- */ + istrt = xadj[link]; + istop = xadj[link + 1] - 1; + i__1 = istop; + for (i = istrt; i <= i__1; ++i) { + node = adjncy[i]; + link = -node; + if (node == enode) { + goto L1400; + } + if (node < 0) { + goto L1000; + } else if (node == 0) { + goto L2100; + } else { + goto L1100; + } + +L1100: + if (qsize[node] == 0) { + goto L1400; + } + if (marker[node] >= *tag) { + goto L1200; + } +/* ------------------------------------- */ +/* CASE WHEN NODE IS NOT YET CONSIDERED. */ +/* ------------------------------------- */ + marker[node] = *tag; + deg += qsize[node]; + goto L1400; +L1200: +/* ---------------------------------------- */ +/* CASE WHEN NODE IS INDISTINGUISHABLE FROM */ +/* ENODE. MERGE THEM INTO A NEW SUPERNODE. */ +/* ---------------------------------------- */ + if (dbakw[node] != 0) { + goto L1400; + } + if (dforw[node] != 2) { + goto L1300; + } + qsize[enode] += qsize[node]; + qsize[node] = 0; + marker[node] = *maxint; + dforw[node] = -enode; + dbakw[node] = -(*maxint); + goto L1400; +L1300: +/* -------------------------------------- */ +/* CASE WHEN NODE IS OUTMATCHED BY ENODE. */ +/* -------------------------------------- */ + if (dbakw[node] == 0) { + dbakw[node] = -(*maxint); + } +L1400: + ; + } + goto L2100; +L1500: +/* ------------------------------------------------ */ +/* FOR EACH ENODE IN THE QX LIST, DO THE FOLLOWING. */ +/* ------------------------------------------------ */ + enode = qxhead; + iq2 = 0; +L1600: + if (enode <= 0) { + goto L2300; + } + if (dbakw[enode] != 0) { + goto L2200; + } + ++(*tag); + deg = deg0; +/* --------------------------------- */ +/* FOR EACH UNMARKED NABOR OF ENODE, */ +/* DO THE FOLLOWING. */ +/* --------------------------------- */ + istrt = xadj[enode]; + istop = xadj[enode + 1] - 1; + i__1 = istop; + for (i = istrt; i <= i__1; ++i) { + nabor = adjncy[i]; + if (nabor == 0) { + goto L2100; + } + if (marker[nabor] >= *tag) { + goto L2000; + } + marker[nabor] = *tag; + link = nabor; +/* ------------------------------ */ +/* IF UNELIMINATED, INCLUDE IT IN */ +/* DEG COUNT. */ +/* ------------------------------ */ + if (dforw[nabor] < 0) { + goto L1700; + } + deg += qsize[nabor]; + goto L2000; +L1700: +/* ------------------------------- */ +/* IF ELIMINATED, INCLUDE UNMARKED */ +/* NODES IN THIS ELEMENT INTO THE */ +/* DEGREE COUNT. */ +/* ------------------------------- */ + jstrt = xadj[link]; + jstop = xadj[link + 1] - 1; + i__2 = jstop; + for (j = jstrt; j <= i__2; ++j) { + node = adjncy[j]; + link = -node; + if (node < 0) { + goto L1700; + } else if (node == 0) { + goto L2000; + } else { + goto L1800; + } + +L1800: + if (marker[node] >= *tag) { + goto L1900; + } + marker[node] = *tag; + deg += qsize[node]; +L1900: + ; + } +L2000: + ; + } +L2100: +/* ------------------------------------------- */ +/* UPDATE EXTERNAL DEGREE OF ENODE IN DEGREE */ +/* STRUCTURE, AND MDEG (MIN DEG) IF NECESSARY. */ +/* ------------------------------------------- */ + deg = deg - qsize[enode] + 1; + fnode = dhead[deg]; + dforw[enode] = fnode; + dbakw[enode] = -deg; + if (fnode > 0) { + dbakw[fnode] = enode; + } + dhead[deg] = enode; + if (deg < *mdeg) { + *mdeg = deg; + } +L2200: +/* ---------------------------------- */ +/* GET NEXT ENODE IN CURRENT ELEMENT. */ +/* ---------------------------------- */ + enode = llist[enode]; + if (iq2 == 1) { + goto L900; + } + goto L1600; +L2300: +/* ----------------------------- */ +/* GET NEXT ELEMENT IN THE LIST. */ +/* ----------------------------- */ + *tag = mtag; + elmnt = llist[elmnt]; + goto L100; + +} /* mmdupd_ */ + +/* *************************************************************** */ +/* *************************************************************** */ +/* ***** MMDNUM ..... MULTI MINIMUM DEGREE NUMBERING ***** */ +/* *************************************************************** */ +/* *************************************************************** */ + +/* AUTHOR - JOSEPH W.H. LIU */ +/* DEPT OF COMPUTER SCIENCE, YORK UNIVERSITY. */ + +/* PURPOSE - THIS ROUTINE PERFORMS THE FINAL STEP IN */ +/* PRODUCING THE PERMUTATION AND INVERSE PERMUTATION */ +/* VECTORS IN THE MULTIPLE ELIMINATION VERSION OF THE */ +/* MINIMUM DEGREE ORDERING ALGORITHM. */ + +/* INPUT PARAMETERS - */ +/* NEQNS - NUMBER OF EQUATIONS. */ +/* QSIZE - SIZE OF SUPERNODES AT ELIMINATION. */ + +/* UPDATED PARAMETERS - */ +/* INVP - INVERSE PERMUTATION VECTOR. ON INPUT, */ +/* IF QSIZE(NODE)=0, THEN NODE HAS BEEN MERGED */ +/* INTO THE NODE -INVP(NODE); OTHERWISE, */ +/* -INVP(NODE) IS ITS INVERSE LABELLING. */ + +/* OUTPUT PARAMETERS - */ +/* PERM - THE PERMUTATION VECTOR. */ + +/* *************************************************************** */ + +int mmdnum_(int* neqns, int* perm, int* invp, int* qsize) +{ + /* System generated locals */ + int i__1; + + /* Local variables */ + static int node, root, nextf, father, nqsize, num; + + +/* *************************************************************** */ + +/* INTEGER*2 INVP(1) , PERM(1) , QSIZE(1) */ + +/* *************************************************************** */ + + /* Parameter adjustments */ + //--qsize; --invp; --perm; + + /* Function Body */ + i__1 = *neqns; + for (node = 1; node <= i__1; ++node) { + nqsize = qsize[node]; + if (nqsize <= 0) { + perm[node] = invp[node]; + } + if (nqsize > 0) { + perm[node] = -invp[node]; + } +/* L100: */ + } +/* ------------------------------------------------------ */ +/* FOR EACH NODE WHICH HAS BEEN MERGED, DO THE FOLLOWING. */ +/* ------------------------------------------------------ */ + i__1 = *neqns; + for (node = 1; node <= i__1; ++node) { + if (perm[node] > 0) { + goto L500; + } +/* ----------------------------------------- */ +/* TRACE THE MERGED TREE UNTIL ONE WHICH HAS */ +/* NOT BEEN MERGED, CALL IT ROOT. */ +/* ----------------------------------------- */ + father = node; +L200: + if (perm[father] > 0) { + goto L300; + } + father = -perm[father]; + goto L200; +L300: +/* ----------------------- */ +/* NUMBER NODE AFTER ROOT. */ +/* ----------------------- */ + root = father; + num = perm[root] + 1; + invp[node] = -num; + perm[root] = num; +/* ------------------------ */ +/* SHORTEN THE MERGED TREE. */ +/* ------------------------ */ + father = node; +L400: + nextf = -perm[father]; + if (nextf <= 0) { + goto L500; + } + perm[father] = -root; + father = nextf; + goto L400; +L500: + ; + } +/* ---------------------- */ +/* READY TO COMPUTE PERM. */ +/* ---------------------- */ + i__1 = *neqns; + for (node = 1; node <= i__1; ++node) { + num = -invp[node]; + invp[node] = num; + perm[num] = node; +/* L600: */ + } + return 0; + +} /* mmdnum_ */ diff --git a/src/hydsolver.c b/src/hydsolver.c index 0adf780..beeff1d 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -122,7 +122,7 @@ int hydsolve(EN_Project *pr, int *iter, double *relerr) hlosscoeff(pr, i); } matrixcoeffs(pr); - errcode = linsolve(&hyd->solver, net->Njuncs); + errcode = linsolve(pr, net->Njuncs); /* Take action depending on error code */ if (errcode < 0) { diff --git a/src/smatrix.c b/src/smatrix.c index 70ec801..8518486 100755 --- a/src/smatrix.c +++ b/src/smatrix.c @@ -15,17 +15,21 @@ module are: freesparse() -- called from closehyd() in HYDRAUL.C linsolve() -- called from netsolve() in HYDRAUL.C -Createsparse() does the following: +createsparse() does the following: 1. for each node, builds an adjacency list that identifies all links connected to the node (see buildlists()) 2. re-orders the network's nodes to minimize the number of non-zero entries in the hydraulic solution matrix - (see reorder()) - 3. converts the adjacency lists into a compact scheme + (see reorder()) + 3. symbolically factorizes the solution matrix + (see factorize()) + 4. converts the adjacency lists into a compact scheme for storing the non-zero coeffs. in the lower diagonal portion of the solution matrix (see storesparse()) -Freesparse() frees the memory used for the sparse matrix. -Linsolve() solves the linearized system of hydraulic equations. + +freesparse() frees the memory used for the sparse matrix. + +linsolve() solves the linearized system of hydraulic equations. ******************************************************************** */ @@ -38,13 +42,54 @@ Linsolve() solves the linearized system of hydraulic equations. #include #endif #include -#include "hash.h" +#include + +#include + #include "text.h" #include "types.h" -#include "epanet2.h" #include "funcs.h" -#define EXTERN extern -#include "vars.h" + +// The multiple minimum degree re-ordering routine (see genmmd.c) +extern int genmmd(int *neqns, int *xadj, int *adjncy, int *invp, int *perm, + int *delta, int *dhead, int *qsize, int *llist, int *marker, + int *maxint, int *nofsub); + + +// Local functions +static int allocsparse(EN_Project *pr); +static int buildlists(EN_Project *pr, int); +static int paralink(EN_Project *pr, int, int, int); +static void xparalinks(EN_Project *pr); +static void freelists(EN_Project *pr); +static void countdegree(EN_Project *pr); +static int reordernodes(EN_Project *pr); +static int factorize(EN_Project *pr); +static int growlist(EN_Project *pr, int); +static int newlink(EN_Project *pr, Padjlist); +static int linked(EN_Network *net, int, int); +static int addlink(EN_Network *net, int, int, int); +static int storesparse(EN_Project *pr, int); +static int sortsparse(EN_Project *pr, int); +static void transpose(int, int *, int *, int *, int *, + int *, int *, int *); + + +/************************************************************************* +* Timer macros +**************************************************************************/ + //#define cleartimer(tmr) (tmr = 0.0) + //#define starttimer(tmr) (tmr -= ((double) clock()/CLOCKS_PER_SEC)); + //#define stoptimer(tmr) (tmr += ((double) clock()/CLOCKS_PER_SEC)); + //#define gettimer(tmr) (tmr) + + +/************************************************************************* +* The following data type implements a timer +**************************************************************************/ +// typedef double timer; +// timer SmatrixTimer; + int createsparse(EN_Project *pr) /* @@ -55,53 +100,60 @@ int createsparse(EN_Project *pr) **-------------------------------------------------------------- */ { - int errcode = 0; - EN_Network *n = &pr->network; - solver_t *s = &pr->hydraulics.solver; + int errcode = 0; - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *solver = &pr->hydraulics.solver; - /* Allocate data structures */ - ERRCODE(allocsparse(pr)); +// cleartimer(SmatrixTimer); +// starttimer(SmatrixTimer); + + + /* Allocate data structures */ + ERRCODE(allocsparse(pr)); - if (errcode) { + if (errcode) { + return(errcode); + } + + /* Build node-link adjacency lists with parallel links removed. */ + solver->Degree = (int *) calloc(net->Nnodes+1, sizeof(int)); + ERRCODE(MEMCHECK(solver->Degree)); + ERRCODE(buildlists(pr, TRUE)); + if (!errcode) + { + xparalinks(pr); // Remove parallel links + countdegree(pr); // Find degree of each junction + } // (= # of adjacent links) + + // Re-order nodes to minimize number of non-zero coeffs. + // in factorized solution matrix. + hyd->Ncoeffs = net->Nlinks; + ERRCODE(reordernodes(pr)); + + // Factorize solution matrix by updating adjacency lists + // with non-zero connections due to fill-ins. + ERRCODE(factorize(pr)); + + // Allocate memory for sparse storage of positions of non-zero + // coeffs. and store these positions in vector NZSUB. + ERRCODE(storesparse(pr, net->Njuncs)); + + // Free memory used for adjacency lists and sort + // row indexes in NZSUB to optimize linsolve(). + if (!errcode) { + freelists(pr); + } + ERRCODE(sortsparse(pr, net->Njuncs)); + + // Re-build adjacency lists without removing parallel + // links for use in future connectivity checking. + ERRCODE(buildlists(pr,FALSE)); + + // Free allocated memory + free(solver->Degree); return(errcode); - } - - /* Build node-link adjacency lists with parallel links removed. */ - s->Degree = (int *) calloc(n->Nnodes+1, sizeof(int)); - ERRCODE(MEMCHECK(s->Degree)); - ERRCODE(buildlists(pr,TRUE)); - if (!errcode){ - xparalinks(pr); /* Remove parallel links */ - countdegree(pr); /* Find degree of each junction */ - } /* (= # of adjacent links) */ - - /* Re-order nodes to minimize number of non-zero coeffs. */ - /* in factorized solution matrix. At same time, adjacency */ - /* list is updated with links representing non-zero coeffs. */ - hyd->Ncoeffs = n->Nlinks; - ERRCODE(reordernodes(pr)); - - /* Allocate memory for sparse storage of positions of non-zero */ - /* coeffs. and store these positions in vector NZSUB. */ - ERRCODE(storesparse(pr,net->Njuncs)); - - /* Free memory used for adjacency lists and sort */ - /* row indexes in NZSUB to optimize linsolve(). */ - if (!errcode) { - freelists(pr); - } - ERRCODE(ordersparse(hyd,net->Njuncs)); - - /* Re-build adjacency lists without removing parallel */ - /* links for use in future connectivity checking. */ - ERRCODE(buildlists(pr,FALSE)); - - /* Free allocated memory */ - free(s->Degree); - return(errcode); } /* End of createsparse */ @@ -114,19 +166,19 @@ int allocsparse(EN_Project *pr) **-------------------------------------------------------------- */ { - EN_Network *n = &pr->network; - solver_t *s = &pr->hydraulics.solver; + EN_Network *net = &pr->network; + solver_t *solver = &pr->hydraulics.solver; - int errcode = 0; - n->Adjlist = (Padjlist *) calloc(n->Nnodes+1, sizeof(Padjlist)); - s->Order = (int *) calloc(n->Nnodes+1, sizeof(int)); - s->Row = (int *) calloc(n->Nnodes+1, sizeof(int)); - s->Ndx = (int *) calloc(n->Nlinks+1, sizeof(int)); - ERRCODE(MEMCHECK(n->Adjlist)); - ERRCODE(MEMCHECK(s->Order)); - ERRCODE(MEMCHECK(s->Row)); - ERRCODE(MEMCHECK(s->Ndx)); - return(errcode); + int errcode = 0; + net->Adjlist = (Padjlist *) calloc(net->Nnodes+1, sizeof(Padjlist)); + solver->Order = (int *) calloc(net->Nnodes+1, sizeof(int)); + solver->Row = (int *) calloc(net->Nnodes+1, sizeof(int)); + solver->Ndx = (int *) calloc(net->Nlinks+1, sizeof(int)); + ERRCODE(MEMCHECK(net->Adjlist)); + ERRCODE(MEMCHECK(solver->Order)); + ERRCODE(MEMCHECK(solver->Row)); + ERRCODE(MEMCHECK(solver->Ndx)); + return(errcode); } @@ -139,17 +191,22 @@ void freesparse(EN_Project *pr) **---------------------------------------------------------------- */ { - EN_Network *n = &pr->network; - solver_t *s = &pr->hydraulics.solver; + EN_Network *net = &pr->network; + solver_t *solver = &pr->hydraulics.solver; + +// stoptimer(SmatrixTimer); +// printf("\n"); +// printf("\n Processing Time = %7.3f s", gettimer(SmatrixTimer)); +// printf("\n"); - freelists(pr); - free(n->Adjlist); - free(s->Order); - free(s->Row); - free(s->Ndx); - free(s->XLNZ); - free(s->NZSUB); - free(s->LNZ); + freelists(pr); + FREE(net->Adjlist); + FREE(solver->Order); + FREE(solver->Row); + FREE(solver->Ndx); + FREE(solver->XLNZ); + FREE(solver->NZSUB); + FREE(solver->LNZ); } /* End of freesparse */ @@ -162,41 +219,41 @@ int buildlists(EN_Project *pr, int paraflag) **-------------------------------------------------------------- */ { - int i,j,k; - int pmark = 0; - int errcode = 0; - Padjlist alink; + int i,j,k; + int pmark = 0; + int errcode = 0; + Padjlist alink; - EN_Network *n = &pr->network; + EN_Network *net = &pr->network; - /* For each link, update adjacency lists of its end nodes */ - for (k=1; k <= n->Nlinks; k++) - { - i = n->Link[k].N1; - j = n->Link[k].N2; - if (paraflag) { - pmark = paralink(pr,i,j,k); /* Parallel link check */ - } + // For each link, update adjacency lists of its end nodes + for (k=1; k <= net->Nlinks; k++) + { + i = net->Link[k].N1; + j = net->Link[k].N2; + if (paraflag) { + pmark = paralink(pr, i, j, k); // Parallel link check + } - /* Include link in start node i's list */ - alink = (struct Sadjlist *) malloc(sizeof(struct Sadjlist)); - if (alink == NULL) return(101); - if (!pmark) alink->node = j; - else alink->node = 0; /* Parallel link marker */ - alink->link = k; - alink->next = n->Adjlist[i]; - n->Adjlist[i] = alink; + // Include link in start node i's list + alink = (struct Sadjlist *) malloc(sizeof(struct Sadjlist)); + if (alink == NULL) return(101); + if (!pmark) alink->node = j; + else alink->node = 0; // Parallel link marker + alink->link = k; + alink->next = net->Adjlist[i]; + net->Adjlist[i] = alink; - /* Include link in end node j's list */ - alink = (struct Sadjlist *) malloc(sizeof(struct Sadjlist)); - if (alink == NULL) return(101); - if (!pmark) alink->node = i; - else alink->node = 0; /* Parallel link marker */ - alink->link = k; - alink->next = n->Adjlist[j]; - n->Adjlist[j] = alink; - } - return(errcode); + // Include link in end node j's list + alink = (struct Sadjlist *) malloc(sizeof(struct Sadjlist)); + if (alink == NULL) return(101); + if (!pmark) alink->node = i; + else alink->node = 0; // Parallel link marker + alink->link = k; + alink->next = net->Adjlist[j]; + net->Adjlist[j] = alink; + } + return(errcode); } /* End of buildlists */ @@ -212,17 +269,20 @@ int paralink(EN_Project *pr, int i, int j, int k) **-------------------------------------------------------------- */ { - Padjlist alink; - for (alink = pr->network.Adjlist[i]; alink != NULL; alink = alink->next) - { - if (alink->node == j) /* Link || to k (same end nodes) */ - { - pr->hydraulics.solver.Ndx[k] = alink->link; /* Assign Ndx entry to this link */ - return(1); - } - } - pr->hydraulics.solver.Ndx[k] = k; /* Ndx entry if link not parallel */ - return(0); + Padjlist alink; + for (alink = pr->network.Adjlist[i]; alink != NULL; alink = alink->next) + { + // Link || to k (same end nodes) + if (alink->node == j) + { + // Assign Ndx entry to this link + pr->hydraulics.solver.Ndx[k] = alink->link; + return(1); + } + } + // Ndx entry if link not parallel + pr->hydraulics.solver.Ndx[k] = k; + return(0); } /* End of paralink */ @@ -235,39 +295,39 @@ void xparalinks(EN_Project *pr) **-------------------------------------------------------------- */ { - int i; - Padjlist alink, /* Current item in adjacency list */ - blink; /* Previous item in adjacency list */ - EN_Network *n = &pr->network; + int i; + Padjlist alink, // Current item in adjacency list + blink; // Previous item in adjacency list + EN_Network *net = &pr->network; - /* Scan adjacency list of each node */ - for (i=1; i <= n->Nnodes; i++) - { - alink = n->Adjlist[i]; /* First item in list */ - blink = NULL; - while (alink != NULL) - { - if (alink->node == 0) /* Parallel link marker found */ - { - if (blink == NULL) /* This holds at start of list */ + // Scan adjacency list of each node + for (i=1; i <= net->Nnodes; i++) + { + alink = net->Adjlist[i]; // First item in list + blink = NULL; + while (alink != NULL) + { + if (alink->node == 0) // Parallel link marker found { - n->Adjlist[i] = alink->next; - free(alink); /* Remove item from list */ - alink = n->Adjlist[i]; - } - else /* This holds for interior of list */ - { - blink->next = alink->next; - free(alink); /* Remove item from list */ - alink = blink->next; - } - } - else - { - blink = alink; /* Move to next item in list */ - alink = alink->next; - } - } + if (blink == NULL) // This holds at start of list + { + net->Adjlist[i] = alink->next; + free(alink); // Remove item from list + alink = net->Adjlist[i]; + } + else // This holds for interior of list + { + blink->next = alink->next; + free(alink); // Remove item from list + alink = blink->next; + } + } + else + { + blink = alink; // Move to next item in list + alink = alink->next; + } + } } } /* End of xparalinks */ @@ -281,19 +341,19 @@ void freelists(EN_Project *pr) **-------------------------------------------------------------- */ { - int i; - Padjlist alink; - EN_Network *n = &pr->network; + int i; + Padjlist alink; + EN_Network *net = &pr->network; - for (i=0; i <= n->Nnodes; i++) - { - for (alink = n->Adjlist[i]; alink != NULL; alink = n->Adjlist[i]) - { - n->Adjlist[i] = alink->next; - free(alink); - } - } + for (i=0; i <= net->Nnodes; i++) + { + for (alink = net->Adjlist[i]; alink != NULL; alink = net->Adjlist[i]) + { + net->Adjlist[i] = alink->next; + free(alink); + } + } } /* End of freelists */ @@ -306,21 +366,22 @@ void countdegree(EN_Project *pr) **---------------------------------------------------------------- */ { - int i; - Padjlist alink; - EN_Network *n = &pr->network; - memset(pr->hydraulics.solver.Degree,0,(n->Nnodes+1) * sizeof(int)); + int i; + Padjlist alink; + EN_Network *net = &pr->network; + + memset(pr->hydraulics.solver.Degree, 0, (net->Nnodes+1) * sizeof(int)); - /* NOTE: For purposes of node re-ordering, Tanks (nodes with */ - /* indexes above Njuncs) have zero degree of adjacency. */ + // NOTE: For purposes of node re-ordering, Tanks (nodes with + // indexes above Njuncs) have zero degree of adjacency. - for (i=1; i <= n->Njuncs; i++) { - for (alink = n->Adjlist[i]; alink != NULL; alink = alink->next) { - if (alink->node > 0) { - pr->hydraulics.solver.Degree[i]++; - } + for (i=1; i <= net->Njuncs; i++) { + for (alink = net->Adjlist[i]; alink != NULL; alink = alink->next) { + if (alink->node > 0) { + pr->hydraulics.solver.Degree[i]++; + } + } } - } } @@ -334,59 +395,109 @@ int reordernodes(EN_Project *pr) **-------------------------------------------------------------- */ { - int k, knode, m, n; - EN_Network *net = &pr->network; - solver_t *s = &pr->hydraulics.solver; + int k, knode, m, njuncs, nlinks; + int delta = -1; + int nofsub = 0; + int maxint = INT_MAX; //defined in limits.h + int errcode; + + EN_Network *net = &pr->network; + solver_t *solver = &pr->hydraulics.solver; + Padjlist alink; + + // Local versions of node adjacency lists + int *adjncy = NULL; + int *xadj = NULL; + + // Work arrays + int *dhead = NULL; + int *qsize = NULL; + int *llist = NULL; + int *marker = NULL; - for (k=1; k <= net->Nnodes; k++) - { - s->Row[k] = k; - s->Order[k] = k; - } - n = net->Njuncs; - for (k=1; k<=n; k++) /* Examine each junction */ - { - m = mindegree(s,k,n); /* Node with lowest degree */ - knode = s->Order[m]; /* Node's index */ - if (!growlist(pr,knode)) { - return(101); /* Augment adjacency list */ + // Default ordering + for (k=1; k <= net->Nnodes; k++) + { + solver->Row[k] = k; + solver->Order[k] = k; } - s->Order[m] = s->Order[k]; /* Switch order of nodes */ - s->Order[k] = knode; - s->Degree[knode] = 0; /* In-activate node */ - } - for (k=1; k<=n; k++) { /* Assign nodes to rows of */ - s->Row[s->Order[k]] = k; /* coeff. matrix */ - } - return(0); + njuncs = net->Njuncs; + nlinks = net->Nlinks; + + // Allocate memory + adjncy = (int *) calloc(2*nlinks+1, sizeof(int)); + xadj = (int *) calloc(njuncs+2, sizeof(int)); + dhead = (int *) calloc(njuncs+1, sizeof(int)); + qsize = (int *) calloc(njuncs + 1, sizeof(int)); + llist = (int *) calloc(njuncs + 1, sizeof(int)); + marker = (int *) calloc(njuncs + 1, sizeof(int)); + if (adjncy && xadj && dhead && qsize && llist && marker) + { + // Create local versions of node adjacency lists + xadj[1] = 1; + m = 1; + for (k = 1; k <= njuncs; k++) + { + for (alink = net->Adjlist[k]; alink != NULL; alink = alink->next) + { + knode = alink->node; + if (knode <= njuncs) + { + adjncy[m] = knode; + m++; + } + } + xadj[k+1] = m; + } + + // Generate a multiple minimum degree node re-ordering + genmmd(&njuncs, xadj, adjncy, solver->Row, solver->Order, &delta, + dhead, qsize, llist, marker, &maxint, &nofsub); + errcode = 0; + } + else errcode = 101; //insufficient memory + + // Free memory + FREE(adjncy); + FREE(xadj); + FREE(dhead); + FREE(qsize); + FREE(llist); + FREE(marker); + return errcode; } /* End of reordernodes */ -int mindegree(solver_t *s, int k, int n) +int factorize(EN_Project *pr) /* **-------------------------------------------------------------- -** Input: k = first node in list of active nodes -** n = total number of junction nodes -** Output: returns node index with fewest direct connections -** Purpose: finds active node with fewest direct connections +** Input: none +** Output: returns error code +** Purpose: symbolically factorizes the solution matrix in +** terms of its adjacency lists **-------------------------------------------------------------- */ { - int i, m; - int min = n, - imin = n; + int k, knode; + int errcode = 0; + EN_Network *net = &pr->network; + solver_t *solver = &pr->hydraulics.solver; - for (i=k; i<=n; i++) - { - m = s->Degree[s->Order[i]]; - if (m < min) - { - min = m; - imin = i; - } - } - return(imin); -} /* End of mindegree */ + // Augment each junction's adjacency list to account for + // new connections created when solution matrix is solved. + // NOTE: Only junctions (indexes <= Njuncs) appear in solution matrix. + for (k = 1; k <= net->Njuncs; k++) // Examine each junction + { + knode = solver->Order[k]; // Re-ordered index + if (!growlist(pr, knode)) // Augment adjacency list + { + errcode = 101; + break; + } + solver->Degree[knode] = 0; // In-activate node + } + return(errcode); +} /* End of factorize */ int growlist(EN_Project *pr, int knode) @@ -400,22 +511,22 @@ int growlist(EN_Project *pr, int knode) **-------------------------------------------------------------- */ { - int node; - Padjlist alink; - EN_Network *n = &pr->network; - solver_t *s = &pr->hydraulics.solver; + int node; + Padjlist alink; + EN_Network *net = &pr->network; + solver_t *solver = &pr->hydraulics.solver; - /* Iterate through all nodes connected to knode */ - for (alink = n->Adjlist[knode]; alink != NULL; alink = alink -> next) - { - node = alink->node; /* End node of connecting link */ - if (s->Degree[node] > 0) /* End node is active */ + // Iterate through all nodes connected to knode + for (alink = net->Adjlist[knode]; alink != NULL; alink = alink -> next) { - s->Degree[node]--; /* Reduce degree of adjacency */ - if (!newlink(pr,alink)) { /* Add to adjacency list */ - return(0); - } - } + node = alink->node; // End node of connecting link + if (solver->Degree[node] > 0) // End node is active + { + solver->Degree[node]--; // Reduce degree of adjacency + if (!newlink(pr, alink)) { // Add to adjacency list + return(0); + } + } } return(1); } /* End of growlist */ @@ -431,41 +542,38 @@ int newlink(EN_Project *pr, Padjlist alink) **-------------------------------------------------------------- */ { - int inode, jnode; - Padjlist blink; - EN_Network *n = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &pr->hydraulics.solver; + int inode, jnode; + Padjlist blink; + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *solver = &pr->hydraulics.solver; - /* Scan all entries in adjacency list that follow anode. */ - inode = alink->node; /* End node of connection to anode */ - for (blink = alink->next; blink != NULL; blink = blink->next) - { - jnode = blink->node; /* End node of next connection */ - - /* If jnode still active, and inode not connected to jnode, */ - /* then add a new connection between inode and jnode. */ - if (s->Degree[jnode] > 0) /* jnode still active */ + // Scan all entries in adjacency list that follow anode. + inode = alink->node; // End node of connection to anode + for (blink = alink->next; blink != NULL; blink = blink->next) { - if (!linked(n, inode,jnode)) { /* inode not linked to jnode */ - /* Since new connection represents a non-zero coeff. */ - /* in the solution matrix, update the coeff. count. */ - hyd->Ncoeffs++; + jnode = blink->node; // End node of next connection + + // If jnode still active, and inode not connected to jnode, + // then add a new connection between inode and jnode. + if (solver->Degree[jnode] > 0) // jnode still active + { + if (!linked(net, inode, jnode)) // inode not linked to jnode + { + // Since new connection represents a non-zero coeff. + // in the solution matrix, update the coeff. count. + hyd->Ncoeffs++; - /* Update adjacency lists for inode & jnode to */ - /* reflect the new connection. */ - if (!addlink(n,inode,jnode,hyd->Ncoeffs)) { - return(0); + // Update adjacency lists for inode & jnode to + // reflect the new connection. + if (!addlink(net, inode, jnode, hyd->Ncoeffs)) return(0); + if (!addlink(net, jnode, inode, hyd->Ncoeffs)) return(0); + solver->Degree[inode]++; + solver->Degree[jnode]++; + } } - if (!addlink(n,jnode,inode,hyd->Ncoeffs)) { - return(0); - } - s->Degree[inode]++; - s->Degree[jnode]++; - } } - } - return(1); + return(1); } /* End of newlink */ @@ -479,13 +587,12 @@ int linked(EN_Network *n, int i, int j) **-------------------------------------------------------------- */ { - Padjlist alink; - for (alink = n->Adjlist[i]; alink != NULL; alink = alink->next) { - if (alink->node == j) { - return(1); + Padjlist alink; + for (alink = n->Adjlist[i]; alink != NULL; alink = alink->next) + { + if (alink->node == j) return(1); } - } - return(0); + return(0); } /* End of linked */ @@ -500,14 +607,14 @@ int addlink(EN_Network *net, int i, int j, int n) **-------------------------------------------------------------- */ { - Padjlist alink; - alink = (struct Sadjlist *) malloc(sizeof(struct Sadjlist)); - if (alink == NULL) return(0); - alink->node = j; - alink->link = n; - alink->next = net->Adjlist[i]; - net->Adjlist[i] = alink; - return(1); + Padjlist alink; + alink = (struct Sadjlist *) malloc(sizeof(struct Sadjlist)); + if (alink == NULL) return(0); + alink->node = j; + alink->link = n; + alink->next = net->Adjlist[i]; + net->Adjlist[i] = alink; + return(1); } /* End of addlink */ @@ -521,49 +628,49 @@ int storesparse(EN_Project *pr, int n) **-------------------------------------------------------------- */ { - Padjlist alink; - int i, ii, j, k, l, m; - int errcode = 0; + Padjlist alink; + int i, ii, j, k, l, m; + int errcode = 0; - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &pr->hydraulics.solver; + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *solver = &pr->hydraulics.solver; - /* Allocate sparse matrix storage */ - s->XLNZ = (int *) calloc(n+2, sizeof(int)); - s->NZSUB = (int *) calloc(hyd->Ncoeffs+2, sizeof(int)); - s->LNZ = (int *) calloc(hyd->Ncoeffs+2, sizeof(int)); - ERRCODE(MEMCHECK(s->XLNZ)); - ERRCODE(MEMCHECK(s->NZSUB)); - ERRCODE(MEMCHECK(s->LNZ)); - if (errcode) { - return(errcode); - } - - /* Generate row index pointers for each column of matrix */ - k = 0; - s->XLNZ[1] = 1; - for (i=1; i<=n; i++) { /* column */ - m = 0; - ii = s->Order[i]; - for (alink = net->Adjlist[ii]; alink != NULL; alink = alink->next) + /* Allocate sparse matrix storage */ + solver->XLNZ = (int *) calloc(n+2, sizeof(int)); + solver->NZSUB = (int *) calloc(hyd->Ncoeffs+2, sizeof(int)); + solver->LNZ = (int *) calloc(hyd->Ncoeffs+2, sizeof(int)); + ERRCODE(MEMCHECK(solver->XLNZ)); + ERRCODE(MEMCHECK(solver->NZSUB)); + ERRCODE(MEMCHECK(solver->LNZ)); + if (errcode) return(errcode); + + // Generate row index pointers for each column of matrix + k = 0; + solver->XLNZ[1] = 1; + for (i=1; i<=n; i++) // column { - j = s->Row[alink->node]; /* row */ - l = alink->link; - if (j > i && j <= n) { - m++; - k++; - s->NZSUB[k] = j; - s->LNZ[k] = l; - } + m = 0; + ii = solver->Order[i]; + for (alink = net->Adjlist[ii]; alink != NULL; alink = alink->next) + { + j = solver->Row[alink->node]; // row + l = alink->link; + if (j > i && j <= n) + { + m++; + k++; + solver->NZSUB[k] = j; + solver->LNZ[k] = l; + } + } + solver->XLNZ[i+1] = solver->XLNZ[i] + m; } - s->XLNZ[i+1] = s->XLNZ[i] + m; - } - return(errcode); + return(errcode); } /* End of storesparse */ -int ordersparse(hydraulics_t *h, int n) +int sortsparse(EN_Project *pr, int n) /* **-------------------------------------------------------------- ** Input: n = number of rows in solution matrix @@ -572,43 +679,48 @@ int ordersparse(hydraulics_t *h, int n) **-------------------------------------------------------------- */ { - int i, k; - int *xlnzt, *nzsubt, *lnzt, *nzt; - int errcode = 0; - solver_t *s = &h->solver; + int i, k; + int *xlnzt, *nzsubt, *lnzt, *nzt; + int errcode = 0; + + hydraulics_t *hyd = &pr->hydraulics; + solver_t *solver = &pr->hydraulics.solver; + + int *LNZ = solver->LNZ; + int *XLNZ = solver->XLNZ; + int *NZSUB = solver->NZSUB; - xlnzt = (int *) calloc(n+2, sizeof(int)); - nzsubt = (int *) calloc(h->Ncoeffs+2, sizeof(int)); - lnzt = (int *) calloc(h->Ncoeffs+2, sizeof(int)); - nzt = (int *) calloc(n+2, sizeof(int)); - ERRCODE(MEMCHECK(xlnzt)); - ERRCODE(MEMCHECK(nzsubt)); - ERRCODE(MEMCHECK(lnzt)); - ERRCODE(MEMCHECK(nzt)); - if (!errcode) { + xlnzt = (int *) calloc(n+2, sizeof(int)); + nzsubt = (int *) calloc(hyd->Ncoeffs+2, sizeof(int)); + lnzt = (int *) calloc(hyd->Ncoeffs+2, sizeof(int)); + nzt = (int *) calloc(n+2, sizeof(int)); + ERRCODE(MEMCHECK(xlnzt)); + ERRCODE(MEMCHECK(nzsubt)); + ERRCODE(MEMCHECK(lnzt)); + ERRCODE(MEMCHECK(nzt)); + if (!errcode) + { + // Count # non-zeros in each row + for (i=1; i<=n; i++) nzt[i] = 0; + for (i=1; i<=n; i++) + { + for (k = XLNZ[i]; k < XLNZ[i+1]; k++) nzt[NZSUB[k]]++; + } + xlnzt[1] = 1; + for (i=1; i<=n; i++) xlnzt[i+1] = xlnzt[i] + nzt[i]; - /* Count # non-zeros in each row */ - for (i=1; i<=n; i++) { - nzt[i] = 0; + // Transpose matrix twice to order column indexes + transpose(n, XLNZ, NZSUB, LNZ, xlnzt, nzsubt, lnzt, nzt); + transpose(n, xlnzt, nzsubt, lnzt, XLNZ, NZSUB, LNZ, nzt); } - for (i=1; i<=n; i++) { - for (k = s->XLNZ[i]; k < s->XLNZ[i+1]; k++) nzt[s->NZSUB[k]]++; - } - xlnzt[1] = 1; - for (i=1; i<=n; i++) xlnzt[i+1] = xlnzt[i] + nzt[i]; - - /* Transpose matrix twice to order column indexes */ - transpose(n,s->XLNZ,s->NZSUB,s->LNZ,xlnzt,nzsubt,lnzt,nzt); - transpose(n,xlnzt,nzsubt,lnzt,s->XLNZ,s->NZSUB,s->LNZ,nzt); - } - /* Reclaim memory */ - free(xlnzt); - free(nzsubt); - free(lnzt); - free(nzt); - return(errcode); -} /* End of ordersparse */ + // Reclaim memory + FREE(xlnzt); + FREE(nzsubt); + FREE(lnzt); + FREE(nzt); + return(errcode); +} /* End of sortsparse */ void transpose(int n, int *il, int *jl, int *xl, int *ilt, int *jlt, @@ -623,24 +735,24 @@ void transpose(int n, int *il, int *jl, int *xl, int *ilt, int *jlt, **--------------------------------------------------------------------- */ { - int i, j, k, kk; + int i, j, k, kk; - for (i=1; i<=n; i++) nzt[i] = 0; - for (i=1; i<=n; i++) - { - for (k=il[i]; kAii; - double *Aij = s->Aij; - double *B = s->F; - int *LNZ = s->LNZ; - int *XLNZ = s->XLNZ; - int *NZSUB = s->NZSUB; + solver_t *solver = &pr->hydraulics.solver; + double *Aii = solver->Aii; + double *Aij = solver->Aij; + double *B = solver->F; + int *LNZ = solver->LNZ; + int *XLNZ = solver->XLNZ; + int *NZSUB = solver->NZSUB; int *link, *first; int i, istop, istrt, isub, j, k, kfirst, newk; diff --git a/src/types.h b/src/types.h index 3c26393..a039d8c 100755 --- a/src/types.h +++ b/src/types.h @@ -100,7 +100,7 @@ typedef int INT4; --------------------------------------------------------------------- */ #define MEMCHECK(x) (((x) == NULL) ? 101 : 0 ) -#define FREE(x) (free((x))) +#define FREE(x) if ((x)) free((x)) /* --------------------------------------------------------------------- diff --git a/tools/before-test.cmd b/tools/before-test.cmd index 8871b44..70a7af9 100644 --- a/tools/before-test.cmd +++ b/tools/before-test.cmd @@ -25,8 +25,8 @@ set SCRIPT_HOME=%~dp0 set TEST_HOME=%~1 -set EXAMPLES_VER=1.0.2-dev -set BENCHMARK_VER=220dev-vs17 +set EXAMPLES_VER=1.0.2-dev.1 +set BENCHMARK_VER=220dev1 set TESTFILES_URL=https://github.com/OpenWaterAnalytics/epanet-example-networks/archive/v%EXAMPLES_VER%.zip diff --git a/tools/before-test.sh b/tools/before-test.sh index 8cb902b..72b1a7e 100755 --- a/tools/before-test.sh +++ b/tools/before-test.sh @@ -22,8 +22,8 @@ SCRIPT_HOME="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" TEST_HOME=$1 -EXAMPLES_VER="1.0.1" -BENCHMARK_VER="2012vs10" +EXAMPLES_VER="1.0.2-dev.1" +BENCHMARK_VER="220dev1" TEST_URL="https://github.com/OpenWaterAnalytics/epanet-example-networks/archive/v${EXAMPLES_VER}.tar.gz" BENCH_URL="https://github.com/OpenWaterAnalytics/epanet-example-networks/releases/download/v${EXAMPLES_VER}/epanet-benchmark-${BENCHMARK_VER}.tar.gz" diff --git a/tools/run-nrtest.cmd b/tools/run-nrtest.cmd index f5d4c8d..e7f55a4 100644 --- a/tools/run-nrtest.cmd +++ b/tools/run-nrtest.cmd @@ -18,7 +18,7 @@ setlocal set NRTEST_SCRIPT_PATH=%~1 set TEST_SUITE_PATH=%~2 -set BENCHMARK_VER=220dev-vs17 +set BENCHMARK_VER=220dev1 set NRTEST_EXECUTE_CMD=python %NRTEST_SCRIPT_PATH%\nrtest execute diff --git a/tools/run-nrtest.sh b/tools/run-nrtest.sh index c561ad6..0c302f7 100755 --- a/tools/run-nrtest.sh +++ b/tools/run-nrtest.sh @@ -19,7 +19,7 @@ run-nrtest() return_value=0 test_suite_path=$1 -benchmark_ver="2012vs10" +benchmark_ver="220dev1" nrtest_execute_cmd="nrtest execute"