From aa81bff09efd7a3c6ca2e9929ad86dd21d0b790e Mon Sep 17 00:00:00 2001 From: Michael Tryby Date: Tue, 19 Jun 2018 17:09:18 -0400 Subject: [PATCH 01/17] Adding build and testing instructions. --- .github/BuildAndTest.md | 183 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 183 insertions(+) create mode 100644 .github/BuildAndTest.md diff --git a/.github/BuildAndTest.md b/.github/BuildAndTest.md new file mode 100644 index 0000000..a37121e --- /dev/null +++ b/.github/BuildAndTest.md @@ -0,0 +1,183 @@ +## Building OWA EPANET From Source on Windows +by Michael E. Tryby + +Created on: March 13, 2018 + + +### Introduction + +Building OWA's fork of EPANET from source is a basic skill that all developers +interested in contributing to the project should know how to perform. This +document describes the build process step-by-step. You will learn 1) how to +configure your machine to build the project locally; 2) how to obtain the +project files using git; 3) how to use cmake to generate build files and build +the project; and 4) how to use ctest and nrtest to perform unit and regression +testing on the build artifacts produced. Be advised, you will need local admin +privileges on your machine to follow this tutorial. Let’s begin! + +### Dependencies + +Before the project can be built the required tools must be installed. The OWA +EPANET project adheres to a platform compiler policy - for each platform there +is a designated compiler. The platform compiler for Windows is Visual +Studio cl, for Linux gcc, and for Mac OS clang. These instructions describe how +to build EPANET on Windows. CMake is a cross platform build, testing, and packaging +tool that is used to automate the EPANET build workflow. Boost is a free portable +peer-reviewed C++ library. Unit tests are linked with Boost unit test libraries. +Lastly, git is a free and open source distributed version control system. Git must +be installed to obtain the project source code from the OWA EPANET repository +found on GitHub. + +### Summary of Build Dependencies + - Platform Compiler + - Windows: Visual Studio 10.0 32-bit cl (version 16.00.40219.01 for 80x86) + - CMake (version 3.0.0 or greater) + - Boost Libraries (version 1.58 or greater) + - git (version 2.6.0 or greater) + +### Build Procedure +1. Install Dependencies + * A. Install Visual Studio 2010 Express and SP1 + Our current benchmark platform and compiler is Windows 32-bit Visual Studio 10 + 2010. Older versions of Visual Studio are available for download here: + + https://www.visualstudio.com/vs/older-downloads/ + + A service pack for Visual Studio 10 2010 is available here: + + https://www.microsoft.com/en-us/download/details.aspx?id=34677 + + * B. Install Boost + Boost binaries for Windows offer a convenient installation solution. Be sure to + select for download the boost installer exe that corresponds to the version of Visual Studio you have installed. + + https://sourceforge.net/projects/boost/files/boost-binaries/1.58.0/ + + Although newer version of Boost are available, a link to Boost 1.58 is provided. This is the library version that the unit tests have been written against. Older versions of Boost may not work. The default install location for the Boost + Libraries is C:\local\boost_1_58_0 + + * C. Install Chocolatey, CMake, and git + Chocolatey is a Windows package manager that makes installing some of these + dependencies a little easier. When working with Chocolatey it is useful to have + local admin privileges. Chocolatey is available here: + + https://chocolatey.org/install + + Once Chocolately is installed, from a command prompt running with admin privileges + issue the following commands + ``` + \>choco install -y cmake --installargs 'ADD_CMAKE_TO_PATH=User' + \>choco install -y git --installargs /GitOnlyOnPath + \>refreshenv + ``` + + * D. Common Problems + Using chocolatey requires a command prompt with admin privileges. + Check to make sure installed applications are on the command path. + Make note of the Boost Library install location. + + +2. Build The Project + As administrator open a Visual Studio 2010 Command Prompt. Change directories to + the location where you wish to build the EPANET project. Now we will issue a series + of commands to create a parent directory for the project root and clone the project + from OWA's GitHub repository. + + * A. Clone the EPANET Repository + ``` + \>mkdir OWA + \>cd OWA + \>git clone --branch=dev https://github.com/OpenWaterAnalytics/EPANET.git + \>cd EPANET + ``` + The present working directory is now the project root EPANET. The directory contains + the same files that are visibly present in the GitHub Repo by browsing to the URL + https://github.com/OpenWaterAnalytics/EPANET/tree/dev. + + Now we will create a build products directory and generate the platform build + file using cmake. + + * B. Generate the build files + ``` + \>mkdir buildprod + \>cd buildprod + \>set BOOST_ROOT=C:\local\boost_1_58_0 + \>cmake -G "Visual Studio 10 2010" -DBOOST_ROOT="%BOOST_ROOT%" -DBoost_USE_STATIC_LIBS="ON" .. + ``` + + Now that the dependencies have been installed and the build system has been + generated, building EPANET is a simple CMake command. + + * C. Build EPANET + \>cmake --build . --config Debug + + * D. Common Problems + CMake may not be able to find the project CMakeLists.txt file or the Boost + library install location. + + +3. Testing + Unit Testing uses Boost Unit Test library and CMake ctest as the test runner. + Cmake has been configured to register tests with ctest as part of the build process. + + * A. Unit Testing + ``` + \>cd tests + \>ctest -C Debug + ``` + The unit tests run quietly. Ctest redirects stdout to a log file which can be + found in the "tests\Testing\Temporary" folder. This is useful when a test fails. + + Regression testing is somewhat more complicated because it relies on Python + to execute EPANET for each test and compare the binary files and report files. + To run regression tests first python and any required packages must be installed. + If Python is already installed on your local machine the installation of + miniconda can be skipped. + + * B. Installing Regression Testing Dependencies + ``` + cd ..\.. + \>choco install -y miniconda --installargs '/AddToPath=1' + \>choco install -y curl + \>choco install -y 7zip + \>refreshenv + \>pip install -r tools/requirements-appveyor.txt + ``` + + With Python and the necessary dependencies installed, regression testing can be run + using the before-test and run-nrtest helper scripts found in the tools folder. The script + before-test stages the test and benchmark files for nrtest. The script run-nrtest calls + nrtest execute and nrtest compare to perform the regression test. + + To run the executable under test, nrtest needs the absolute path to it and a + unique identifier for it such as the version number. The project cmake build places build + artifacts in the buildprod\bin\ folder. On Windows the build configuration "Debug" or + "Release" must also be indicated. On Windows it is also necessary to specify the path to + the Python Scripts folder so the nrtest execute and compare commands can be found. You + need to substitute bracketed fields below like "" with the values for + your setup. + + * C. Regression Testing + ``` + \>tools\before-test.cmd + \>tools\run-nrtest.cmd + ``` + + * D. Common Problems + The batch file before-test.cmd needs to run with admin privileges. The nrtest script complains when it can't find manifest files. + +That concludes this tutorial on building OWA EPANET from source on Windows. +You have learned how to configure your machine satisfying project dependencies +and how to acquire, build, and test EPANET on your local machine. To be sure, +there is a lot more to learn, but this is a good start! Learn more about project +build and testing dependencies by following the links provided below. + +### Further Reading + * Visual Studio - https://msdn.microsoft.com/en-us/library/dd831853(v=vs.100).aspx + * CMake - https://cmake.org/documentation/ + * Boost - http://www.boost.org/doc/ + * git - https://git-scm.com/doc + * Miniconda - https://conda.io/docs/user-guide/index.html + * curl - https://curl.haxx.se/ + * 7zip - https://www.7-zip.org/ + * nrtest - https://nrtest.readthedocs.io/en/latest/ From 86ffbcf614440d82ff648fe0ab8c655a3dc2c5dd Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Mon, 2 Jul 2018 13:55:41 -0400 Subject: [PATCH 02/17] Improved node re-ordering method (Issue #162) --- src/funcs.h | 18 +- src/genmmd.c | 1001 +++++++++++++++++++++++++++++++++++++++++++++++ src/hydsolver.c | 2 +- src/smatrix.c | 832 ++++++++++++++++++++++----------------- src/types.h | 2 +- 5 files changed, 1476 insertions(+), 379 deletions(-) create mode 100644 src/genmmd.c 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)) /* --------------------------------------------------------------------- From f38a637aaf1dc6af1724969b5ded392940f3924b Mon Sep 17 00:00:00 2001 From: Michael Tryby Date: Tue, 10 Jul 2018 18:50:33 -0400 Subject: [PATCH 03/17] Adding mincdd comparison and better debug information --- tools/nrtest-epanet/main.py | 39 +++++++----- tools/nrtest-epanet/nrtest_epanet/__init__.py | 63 ++++++++++++++++--- .../nrtest_epanet/output_reader.py | 6 +- tools/nrtest-epanet/setup.py | 1 + 4 files changed, 84 insertions(+), 25 deletions(-) diff --git a/tools/nrtest-epanet/main.py b/tools/nrtest-epanet/main.py index 1c0f41d..d51c44e 100644 --- a/tools/nrtest-epanet/main.py +++ b/tools/nrtest-epanet/main.py @@ -4,7 +4,9 @@ import time import cStringIO import itertools as it -import epanet_reader as er +# project import +import nrtest_epanet.output_reader as er + def result_compare(path_test, path_ref, comp_args): @@ -15,29 +17,36 @@ def result_compare(path_test, path_ref, comp_args): total = 0 output = cStringIO.StringIO() eps = np.finfo(float).eps + min_cdd = 100.0 start = time.time() - test_reader = er.reader(path_test) - ref_reader = er.reader(path_ref) + test_reader = er.output_generator(path_test) + ref_reader = er.output_generator(path_ref) for test, ref in it.izip(test_reader, ref_reader): total += 1 if total%100000 == 0: print(total) - if test.size != ref.size: + if len(test[0]) != len(ref[0]): raise ValueError('Inconsistent lengths') # Skip results if they are zero or equal - if np.array_equal(test, ref): - equal += 1 - continue + #if np.array_equal(test, ref): + # equal += 1 + # continue else: try: - np.testing.assert_allclose(test, ref, 1.0e-06, 2*eps) - close += 1 + diff = np.fabs(np.subtract(test[0], ref[0])) + idx = np.unravel_index(np.argmax(diff), diff.shape) + if diff[idx] != 0.0: + tmp = - np.log10(diff[idx]) + + if tmp < min_cdd: + min_cdd = tmp; + except AssertionError as ae: notclose += 1 output.write(str(ae)) @@ -49,8 +58,9 @@ def result_compare(path_test, path_ref, comp_args): print(output.getvalue()) output.close() - print('equal: %d close: %d notclose: %d total: %d in %f (sec)\n' % - (equal, close, notclose, total, (stop - start))) + print('mincdd: %d in %f (sec)' % (np.floor(min_cdd), (stop - start))) + #print('equal: %d close: %d notclose: %d total: %d in %f (sec)\n' % + # (equal, close, notclose, total, (stop - start))) if notclose > 0: print('%d differences found\n' % notclose) @@ -130,7 +140,8 @@ if __name__ == "__main__": # ref_path = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\epanet-testsuite\\benchmarks\\v2012" # print(nrtest_compare(test_path, ref_path, (0.001, 0.0))) - - path_test = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\epanet-testsuite\\benchmarks\\v2011a\\Example_3\\example3.out" - path_ref = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\epanet-testsuite\\benchmarks\\v2012\\Example_3\\example3.out" + benchmark_path = "C:\\Users\\mtryby\\Workspace\\GitRepo\\michaeltryby\\epanet-lr\\nrtestsuite\\benchmarks\\" + path_test = benchmark_path + "epanet-220dev\\example1\\example1.out" + path_ref = benchmark_path + "epanet-2012\\example1\\example1.out" + result_compare(path_test, path_ref, (0.001, 0.0)) diff --git a/tools/nrtest-epanet/nrtest_epanet/__init__.py b/tools/nrtest-epanet/nrtest_epanet/__init__.py index 7bdf41b..6135873 100644 --- a/tools/nrtest-epanet/nrtest_epanet/__init__.py +++ b/tools/nrtest-epanet/nrtest_epanet/__init__.py @@ -38,7 +38,7 @@ __status = "Development" def epanet_allclose_compare(path_test, path_ref, rtol, atol): ''' - Compares results in two EPANET binary files. Using the comparison criteria + Compares results in two EPANET binary files using the comparison criteria described in the numpy assert_allclose documentation. (test_value - ref_value) <= atol + rtol * abs(ref_value) @@ -67,22 +67,67 @@ def epanet_allclose_compare(path_test, path_ref, rtol, atol): for (test, ref) in it.izip(ordr.output_generator(path_test), ordr.output_generator(path_ref)): - if len(test) != len(ref): + if len(test[0]) != len(ref[0]): raise ValueError('Inconsistent lengths') # Skip over arrays that are equal - if np.array_equal(test, ref): + if np.array_equal(test[0], ref[0]): continue else: - np.testing.assert_allclose(test, ref, rtol, atol) + np.testing.assert_allclose(test[0], ref[0], rtol, atol) return True -# def epanet_better_compare(path_test, path_ref, rtol, atol): -# ''' -# If you don't like assert_allclose you can add another function here. -# ''' -# pass + +def epanet_mincdd_compare(path_test, path_ref, rtol, atol): + ''' + Compares the results of two EPANET binary files using a correct decimal + digits (cdd) comparison criteria: + + min cdd(test, ref) >= atol + + Returns true if min cdd in the file is greater than or equal to atol, + otherwise an AssertionError is thrown. + + Arguments: + path_test - path to result file being testedgit + path_ref - path to reference result file + rtol - ignored + atol - minimum allowable cdd value (i.e. 3) + + Returns: + True + + Raises: + ValueError() + AssertionError() + ''' + min_cdd = 100.0 + + for (test, ref) in it.izip(ordr.output_generator(path_test), + ordr.output_generator(path_ref)): + + if len(test[0]) != len(ref[0]): + raise ValueError('Inconsistent lengths') + + # Skip over arrays that are equal + if np.array_equal(test[0], ref[0]): + continue + else: + diff = np.fabs(np.subtract(test[0], ref[0])) + idx = np.unravel_index(np.argmax(diff), diff.shape) + + if diff[idx] != 0.0: + tmp = - np.log10(diff[idx]) + + if tmp < min_cdd: + min_cdd = tmp; + + if np.floor(min_cdd) >= atol: + return True + else: + raise AssertionError('min_cdd=%d less than atol=%g' % (min_cdd, atol)) + def epanet_report_compare(path_test, path_ref, rtol, atol): ''' diff --git a/tools/nrtest-epanet/nrtest_epanet/output_reader.py b/tools/nrtest-epanet/nrtest_epanet/output_reader.py index 0bf640f..6877d91 100644 --- a/tools/nrtest-epanet/nrtest_epanet/output_reader.py +++ b/tools/nrtest-epanet/nrtest_epanet/output_reader.py @@ -23,7 +23,8 @@ def output_generator(path_ref): yield element attributes. It is useful for comparing contents of binary files for numerical regression testing. - The generator yields a Python list containing element attributes. + The generator yields a Python tuple containing an array of element + attributes and a tuple containing the element type, period, and attribute. Arguments: path_ref - path to result file @@ -38,7 +39,8 @@ def output_generator(path_ref): for element_type in oapi.ElementType: for attribute in br.elementAttributes[element_type]: - yield br.element_attribute(element_type, period_index, attribute) + yield (br.element_attribute(element_type, period_index, attribute), + (element_type, period_index, attribute)) class OutputReader(): diff --git a/tools/nrtest-epanet/setup.py b/tools/nrtest-epanet/setup.py index 897e621..1a18763 100644 --- a/tools/nrtest-epanet/setup.py +++ b/tools/nrtest-epanet/setup.py @@ -17,6 +17,7 @@ except ImportError: entry_points = { 'nrtest.compare': [ 'epanet allclose = nrtest_epanet:epanet_allclose_compare', + 'epanet mincdd = nrtest_epanet:epanet_mincdd_compare', 'epanet report = nrtest_epanet:epanet_report_compare', # Add entry point for new comparison functions here ] From e3cd4ae41aa40c39f79a1f3d3b4d5e08762f4eee Mon Sep 17 00:00:00 2001 From: Michael Tryby Date: Wed, 11 Jul 2018 10:04:10 -0400 Subject: [PATCH 04/17] Bumping version --- tools/nrtest-epanet/nrtest_epanet/__init__.py | 2 +- tools/nrtest-epanet/setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/nrtest-epanet/nrtest_epanet/__init__.py b/tools/nrtest-epanet/nrtest_epanet/__init__.py index 6135873..d6ff9bd 100644 --- a/tools/nrtest-epanet/nrtest_epanet/__init__.py +++ b/tools/nrtest-epanet/nrtest_epanet/__init__.py @@ -28,7 +28,7 @@ __copyright__ = "None" __credits__ = "Colleen Barr, Maurizio Cingi, Mark Gray, David Hall, Bryant McDonnell" __license__ = "CC0 1.0 Universal" -__version__ = "0.4.0" +__version__ = "0.5.0" __date__ = "September 6, 2017" __maintainer__ = "Michael Tryby" diff --git a/tools/nrtest-epanet/setup.py b/tools/nrtest-epanet/setup.py index 1a18763..4c7dfd7 100644 --- a/tools/nrtest-epanet/setup.py +++ b/tools/nrtest-epanet/setup.py @@ -25,7 +25,7 @@ entry_points = { setup( name='nrtest-epanet', - version='0.4.0', + version='0.5.0', description="EPANET extension for nrtest", author="Michael E. Tryby", From c482f38639bdbbd1fab110d06b87167396ddccb8 Mon Sep 17 00:00:00 2001 From: Michael Tryby Date: Wed, 11 Jul 2018 17:19:18 -0400 Subject: [PATCH 05/17] Initial commit for report_diff --- tools/nrtest-epanet/main.py | 10 ++- tools/nrtest-epanet/nrtest_epanet/__init__.py | 2 +- tools/nrtest-epanet/report_diff.py | 69 +++++++++++++++++++ 3 files changed, 77 insertions(+), 4 deletions(-) create mode 100644 tools/nrtest-epanet/report_diff.py diff --git a/tools/nrtest-epanet/main.py b/tools/nrtest-epanet/main.py index d51c44e..c962dab 100644 --- a/tools/nrtest-epanet/main.py +++ b/tools/nrtest-epanet/main.py @@ -130,6 +130,8 @@ def nrtest_execute(app_path, test_path, output_path): exit(not success) +import report_diff as rd + if __name__ == "__main__": # app_path = "apps\\swmm-5.1.11.json" # test_path = "tests\\examples\\example1.json" @@ -141,7 +143,9 @@ if __name__ == "__main__": # print(nrtest_compare(test_path, ref_path, (0.001, 0.0))) benchmark_path = "C:\\Users\\mtryby\\Workspace\\GitRepo\\michaeltryby\\epanet-lr\\nrtestsuite\\benchmarks\\" - path_test = benchmark_path + "epanet-220dev\\example1\\example1.out" - path_ref = benchmark_path + "epanet-2012\\example1\\example1.out" + path_test = benchmark_path + "epanet-220dev\\example2\\example2.out" + path_ref = benchmark_path + "epanet-2012\\example2\\example2.out" - result_compare(path_test, path_ref, (0.001, 0.0)) + #result_compare(path_test, path_ref, (0.001, 0.0)) + rd.report_diff(path_test, path_ref, 2) + \ No newline at end of file diff --git a/tools/nrtest-epanet/nrtest_epanet/__init__.py b/tools/nrtest-epanet/nrtest_epanet/__init__.py index d6ff9bd..a118695 100644 --- a/tools/nrtest-epanet/nrtest_epanet/__init__.py +++ b/tools/nrtest-epanet/nrtest_epanet/__init__.py @@ -118,7 +118,7 @@ def epanet_mincdd_compare(path_test, path_ref, rtol, atol): idx = np.unravel_index(np.argmax(diff), diff.shape) if diff[idx] != 0.0: - tmp = - np.log10(diff[idx]) + tmp = - np.log10(diff[idx]) if tmp < min_cdd: min_cdd = tmp; diff --git a/tools/nrtest-epanet/report_diff.py b/tools/nrtest-epanet/report_diff.py new file mode 100644 index 0000000..689d7aa --- /dev/null +++ b/tools/nrtest-epanet/report_diff.py @@ -0,0 +1,69 @@ +# -*- coding: utf-8 -*- + +# +# report_diff.py +# +# Date Created: July 11, 2018 +# +# Author: Michael E. Tryby +# US EPA - ORD/NRMRL +# + +# system imports +import itertools as it + +# third party imports +import numpy as np + +# project imports +import nrtest_epanet.output_reader as ordr + + +def report_diff(path_test, path_ref, min_cdd): + for (test, ref) in it.izip(ordr.output_generator(path_test), + ordr.output_generator(path_ref)): + + if len(test[0]) != len(ref[0]): + raise ValueError('Inconsistent lengths') + + # Skip over arrays that are equal + if np.array_equal(test[0], ref[0]): + continue + else: + lre = log_relative_error(test[0], ref[0]) + idx = np.unravel_index(np.argmin(lre), lre.shape) + + if lre[idx] < min_cdd: + print_diff(idx, lre, test, ref) + + return + + +def log_relative_error(q, c): + ''' + Computes log relative error, a measure of numerical accuracy. + + Single precision machine epsilon is between 2^-24 and 2^-23. + + Reference: + McCullough, B. D. "Assessing the Reliability of Statistical Software: Part I." + The American Statistician, vol. 52, no. 4, 1998, pp. 358�366. + JSTOR, JSTOR, www.jstor.org/stable/2685442. + ''' + diff = np.subtract(q, c) + tmp_c = np.copy(c) + # If ref value is small compute absolute error + tmp_c[np.fabs(tmp_c) <= 1.0e-6] = 1.0 + + re = np.fabs(diff)/np.fabs(tmp_c) + # If re is tiny set lre to number of digits + re[re < 1.0e-7] = 1.0e-7 + # If re is very large set lre to zero + re[re > 2.0] = 1.0 + + return np.negative(np.log10(re)) + + +def print_diff(idx, lre, test, ref): + print("LRE: %f\nIdx: %s\nSut: %f\nRef: %f\n" + % ((lre[idx]),(idx[0], ref[1]),(test[0][idx[0]]),(ref[0][idx[0]]))) \ No newline at end of file From 86c55a237328708de3a3b06484922ec22490f610 Mon Sep 17 00:00:00 2001 From: Michael Tryby Date: Wed, 11 Jul 2018 17:22:26 -0400 Subject: [PATCH 06/17] Fixing typo --- tools/nrtest-epanet/report_diff.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tools/nrtest-epanet/report_diff.py b/tools/nrtest-epanet/report_diff.py index 689d7aa..05f5e26 100644 --- a/tools/nrtest-epanet/report_diff.py +++ b/tools/nrtest-epanet/report_diff.py @@ -47,8 +47,7 @@ def log_relative_error(q, c): Reference: McCullough, B. D. "Assessing the Reliability of Statistical Software: Part I." - The American Statistician, vol. 52, no. 4, 1998, pp. 358�366. - JSTOR, JSTOR, www.jstor.org/stable/2685442. + The American Statistician, vol. 52, no. 4, 1998, pp. 358-366. ''' diff = np.subtract(q, c) tmp_c = np.copy(c) From db40fa7bbe324fa49879cf94fdec851fc7f54676 Mon Sep 17 00:00:00 2001 From: Michael Tryby Date: Thu, 12 Jul 2018 10:28:25 -0400 Subject: [PATCH 07/17] Improving diff report format --- tools/nrtest-epanet/report_diff.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/tools/nrtest-epanet/report_diff.py b/tools/nrtest-epanet/report_diff.py index 05f5e26..00570fe 100644 --- a/tools/nrtest-epanet/report_diff.py +++ b/tools/nrtest-epanet/report_diff.py @@ -17,6 +17,8 @@ import numpy as np # project imports import nrtest_epanet.output_reader as ordr +from numpy.f2py.auxfuncs import hasinitvalue +from numpy.tests.test_numpy_version import test_valid_numpy_version def report_diff(path_test, path_ref, min_cdd): @@ -52,7 +54,7 @@ def log_relative_error(q, c): diff = np.subtract(q, c) tmp_c = np.copy(c) # If ref value is small compute absolute error - tmp_c[np.fabs(tmp_c) <= 1.0e-6] = 1.0 + tmp_c[np.fabs(tmp_c) < 1.0e-6] = 1.0 re = np.fabs(diff)/np.fabs(tmp_c) # If re is tiny set lre to number of digits @@ -64,5 +66,13 @@ def log_relative_error(q, c): def print_diff(idx, lre, test, ref): - print("LRE: %f\nIdx: %s\nSut: %f\nRef: %f\n" - % ((lre[idx]),(idx[0], ref[1]),(test[0][idx[0]]),(ref[0][idx[0]]))) \ No newline at end of file + + idx_val = (idx[0], ref[1]) + test_val = (test[0][idx[0]]) + ref_val = (ref[0][idx[0]]) + diff_val = (test_val - ref_val) + lre_val = (lre[idx[0]]) + + print("Idx: %s\nSut: %f Ref: %f Diff: %f LRE: %f\n" + % (idx_val, test_val, ref_val, diff_val, lre_val)) + \ No newline at end of file From 10c5608210f2977c6e18832c0c67d3b427dc0582 Mon Sep 17 00:00:00 2001 From: Michael Tryby Date: Mon, 16 Jul 2018 14:25:53 -0400 Subject: [PATCH 08/17] Adding command line interface to report_diff --- tools/nrtest-epanet/report_diff.py | 32 +++++++++++++++++++++++------- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/tools/nrtest-epanet/report_diff.py b/tools/nrtest-epanet/report_diff.py index 00570fe..b3f49e3 100644 --- a/tools/nrtest-epanet/report_diff.py +++ b/tools/nrtest-epanet/report_diff.py @@ -17,11 +17,9 @@ import numpy as np # project imports import nrtest_epanet.output_reader as ordr -from numpy.f2py.auxfuncs import hasinitvalue -from numpy.tests.test_numpy_version import test_valid_numpy_version -def report_diff(path_test, path_ref, min_cdd): +def _binary_diff(path_test, path_ref, min_cdd): for (test, ref) in it.izip(ordr.output_generator(path_test), ordr.output_generator(path_ref)): @@ -32,16 +30,16 @@ def report_diff(path_test, path_ref, min_cdd): if np.array_equal(test[0], ref[0]): continue else: - lre = log_relative_error(test[0], ref[0]) + lre = _log_relative_error(test[0], ref[0]) idx = np.unravel_index(np.argmin(lre), lre.shape) if lre[idx] < min_cdd: - print_diff(idx, lre, test, ref) + _print_diff(idx, lre, test, ref) return -def log_relative_error(q, c): +def _log_relative_error(q, c): ''' Computes log relative error, a measure of numerical accuracy. @@ -65,7 +63,7 @@ def log_relative_error(q, c): return np.negative(np.log10(re)) -def print_diff(idx, lre, test, ref): +def _print_diff(idx, lre, test, ref): idx_val = (idx[0], ref[1]) test_val = (test[0][idx[0]]) @@ -75,4 +73,24 @@ def print_diff(idx, lre, test, ref): print("Idx: %s\nSut: %f Ref: %f Diff: %f LRE: %f\n" % (idx_val, test_val, ref_val, diff_val, lre_val)) + + +def report(args): + _binary_diff(args.test, args.ref, args.mincdd) + + +if __name__ == '__main__': + from argparse import ArgumentParser + + parser = ArgumentParser(description='EPANET benchmark difference reporting') + parser.set_defaults(func=report) + parser.add_argument('-t', '--test', default=None, + help='Path to test benchmark') + parser.add_argument('-r', '--ref', default=None, + help='Path to reference benchmark') + parser.add_argument('-mc', '--mincdd', type=int, default=3, + help='Minimum correct decimal digits') + + args = parser.parse_args() + args.func(args) \ No newline at end of file From 711bde01d147ab95ea3bb0efa02d9bcdd510f423 Mon Sep 17 00:00:00 2001 From: Michael Tryby Date: Mon, 16 Jul 2018 17:02:35 -0400 Subject: [PATCH 09/17] Switching build to VS 2017 and bumping benchmark version --- appveyor.yml | 10 +++++----- tools/before-test.cmd | 4 ++-- tools/before-test.sh | 4 ++-- 3 files changed, 9 insertions(+), 9 deletions(-) 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/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..afd543a 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="2020dev1" 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" From e38d7ece505efdea2937cf5d70fa4294e26cffa2 Mon Sep 17 00:00:00 2001 From: Michael Tryby Date: Mon, 16 Jul 2018 17:09:02 -0400 Subject: [PATCH 10/17] Bumping benchmark version --- tools/run-nrtest.cmd | 2 +- tools/run-nrtest.sh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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" From e6e794258505099e305cb226a0259526f9dafb2e Mon Sep 17 00:00:00 2001 From: Michael Tryby Date: Mon, 16 Jul 2018 17:18:09 -0400 Subject: [PATCH 11/17] Fixing typo --- tools/before-test.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/before-test.sh b/tools/before-test.sh index afd543a..72b1a7e 100755 --- a/tools/before-test.sh +++ b/tools/before-test.sh @@ -23,7 +23,7 @@ SCRIPT_HOME="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" TEST_HOME=$1 EXAMPLES_VER="1.0.2-dev.1" -BENCHMARK_VER="2020dev1" +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" From 528e99a791c04a19ec838633f3ac25e4bfa6386b Mon Sep 17 00:00:00 2001 From: Elad Salomons Date: Tue, 17 Jul 2018 07:17:33 +0300 Subject: [PATCH 12/17] Update win make file to include genmmd.c --- win_build/WinSDK/Makefile.bat | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/win_build/WinSDK/Makefile.bat b/win_build/WinSDK/Makefile.bat index d0a00c3..22fb38f 100644 --- a/win_build/WinSDK/Makefile.bat +++ b/win_build/WinSDK/Makefile.bat @@ -21,9 +21,9 @@ Find /i "x86" < checkOS.tmp > StringCheck.tmp If %ERRORLEVEL% == 1 ( CALL "%SDK_PATH%bin\"SetEnv.cmd /x64 /release rem : create EPANET2.DLL - cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL + cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL rem : create EPANET2.EXE - cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link + cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /I ..\src /link md "%Build_PATH%"\64bit move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\64bit move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\64bit @@ -35,9 +35,9 @@ rem : 32 bit with DEF CALL "%SDK_PATH%bin\"SetEnv.cmd /x86 /release echo "32 bit with epanet2.def mapping" rem : create EPANET2.DLL -cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP +cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP rem : create EPANET2.EXE -cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link +cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /I ..\src /link md "%Build_PATH%"\32bit move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\32bit move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\32bit From 4e48c3114794dbdd7e8cffd562bcb396fcc0e154 Mon Sep 17 00:00:00 2001 From: Elad Salomons Date: Wed, 15 Aug 2018 07:37:19 +0300 Subject: [PATCH 13/17] Fix PumpType and CurveType enums Also adds a general \default curve type. Fixes #208. --- include/epanet2.bas | 7 +++++++ include/epanet2.h | 15 +++++++-------- src/epanet.c | 11 +++-------- src/types.h | 5 +++-- 4 files changed, 20 insertions(+), 18 deletions(-) diff --git a/include/epanet2.bas b/include/epanet2.bas index 0826a3b..615f045 100644 --- a/include/epanet2.bas +++ b/include/epanet2.bas @@ -154,6 +154,13 @@ Public Const EN_INITFLOW = 10 ' Re-initialize flow flag Public Const EN_CONST_HP = 0 ' constant horsepower Public Const EN_POWER_FUNC = 1 ' power function Public Const EN_CUSTOM = 2 ' user-defined custom curve +Public Const EN_NOCURVE = 3 ' no curve + +Public Const EN_V_CURVE = 0 ' volume curve +Public Const EN_P_CURVE = 1 ' pump curve +Public Const EN_E_CURVE = 2 ' efficiency curve +Public Const EN_H_CURVE = 3 ' head loss curve +Public Const EN_G_CURVE = 4 ' General\default curve 'These are the external functions that comprise the DLL diff --git a/include/epanet2.h b/include/epanet2.h index 81b469e..b41ae0f 100644 --- a/include/epanet2.h +++ b/include/epanet2.h @@ -224,8 +224,6 @@ typedef enum { EN_TIMEOFDAY = 3 } EN_ControlType; - - typedef enum { EN_AVERAGE = 1, /* Time statistic types. */ EN_MINIMUM = 2, /* See TstatType in TYPES.H */ @@ -233,8 +231,6 @@ typedef enum { EN_RANGE = 4 } EN_StatisticType; - - typedef enum { EN_MIX1 = 0, /* Tank mixing models */ EN_MIX2 = 1, @@ -242,8 +238,6 @@ typedef enum { EN_LIFO = 3 } EN_MixingModel; - - typedef enum { EN_NOSAVE = 0, EN_SAVE = 1, @@ -251,8 +245,6 @@ typedef enum { EN_SAVE_AND_INIT = 11 } EN_SaveOption; - - typedef enum { EN_CONST_HP = 0, /* constant horsepower */ EN_POWER_FUNC = 1, /* power function */ @@ -260,6 +252,13 @@ typedef enum { } EN_CurveType; +typedef enum { + EN_V_CURVE = 0, /* volume curve */ + EN_P_CURVE = 1, /* pump curve */ + EN_E_CURVE = 2, /* efficiency curve */ + EN_H_CURVE = 3, /* head loss curve */ + EN_G_CURVE = 4 /* General\default curve */ +} EN_CurveType; // --- Declare the EPANET toolkit functions #if defined(__cplusplus) diff --git a/src/epanet.c b/src/epanet.c index f0f0ff8..72f59d8 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -2933,7 +2933,7 @@ int DLLEXPORT EN_addcurve(EN_Project *p, char *id) { strcpy(tmpCur[n].ID, id); tmpCur[n].Npts = 1; - tmpCur[n].Type = -1; + tmpCur[n].Type = G_CURVE; tmpCur[n].X = (double *)calloc(tmpCur[n].Npts, sizeof(double)); tmpCur[n].Y = (double *)calloc(tmpCur[n].Npts, sizeof(double)); if (tmpCur[n].X == NULL) @@ -2974,9 +2974,7 @@ int DLLEXPORT EN_addcurve(EN_Project *p, char *id) { int DLLEXPORT EN_setcurve(EN_Project *p, int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y, int n) { EN_Network *net = &p->network; - Scurve *Curve = net->Curve; - - + Scurve *Curve = net->Curve; int j; /* Check for valid arguments */ @@ -3007,12 +3005,9 @@ int DLLEXPORT EN_setcurve(EN_Project *p, int index, EN_API_FLOAT_TYPE *x, EN_API int DLLEXPORT EN_setcurvevalue(EN_Project *p, int index, int pnt, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y) { EN_Network *net = &p->network; - Scurve *Curve = net->Curve; - const int Ncurves = net->Ncurves; - if (!p->Openflag) return (102); if (index <= 0 || index > Ncurves) @@ -3713,7 +3708,7 @@ int allocdata(EN_Project *p) } for (n = 0; n <= par->MaxCurves; n++) { net->Curve[n].Npts = 0; - net->Curve[n].Type = -1; + net->Curve[n].Type = G_CURVE; net->Curve[n].X = NULL; net->Curve[n].Y = NULL; } diff --git a/src/types.h b/src/types.h index a039d8c..be6e6fa 100755 --- a/src/types.h +++ b/src/types.h @@ -151,8 +151,9 @@ typedef enum { V_CURVE, /* volume curve */ P_CURVE, /* pump curve */ E_CURVE, /* efficiency curve */ - H_CURVE -} CurveType; /* head loss curve */ + H_CURVE, /* head loss curve */ + G_CURVE /* General\default curve */ +} CurveType; typedef enum { CONST_HP, /* constant horsepower */ From e06368eca0fe51a9de867bb79e2a35cd8914dc4e Mon Sep 17 00:00:00 2001 From: Elad Salomons Date: Wed, 15 Aug 2018 16:45:51 +0300 Subject: [PATCH 14/17] Adtional chnage for #208 --- include/epanet2.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/epanet2.h b/include/epanet2.h index b41ae0f..0470d9e 100644 --- a/include/epanet2.h +++ b/include/epanet2.h @@ -248,9 +248,9 @@ typedef enum { typedef enum { EN_CONST_HP = 0, /* constant horsepower */ EN_POWER_FUNC = 1, /* power function */ - EN_CUSTOM = 2 /* user-defined custom curve */ -} EN_CurveType; - + EN_CUSTOM = 2, /* user-defined custom curve */ + EN_NOCURVE = 3 /* no curve */ +} EN_PumpType; typedef enum { EN_V_CURVE = 0, /* volume curve */ From 99e71511726271fffde360a8f526e143e054a23d Mon Sep 17 00:00:00 2001 From: Elad Salomons Date: Fri, 17 Aug 2018 15:43:53 +0300 Subject: [PATCH 15/17] Add option to get link state Add EN_STATE to ENgetlinkvalue. Related to #218 --- include/epanet2.bas | 1 + include/epanet2.h | 3 ++- src/epanet.c | 17 +++++++++++++++-- 3 files changed, 18 insertions(+), 3 deletions(-) diff --git a/include/epanet2.bas b/include/epanet2.bas index 615f045..585d556 100644 --- a/include/epanet2.bas +++ b/include/epanet2.bas @@ -56,6 +56,7 @@ Public Const EN_EFFICIENCY = 16 Public Const EN_HEADCURVE = 17 Public Const EN_EFFICIENCYCURVE = 18 Public Const EN_PRICEPATTERN = 19 +Public Const EN_STATE = 20 Public Const EN_DURATION = 0 ' Time parameters Public Const EN_HYDSTEP = 1 diff --git a/include/epanet2.h b/include/epanet2.h index 0470d9e..3b26116 100644 --- a/include/epanet2.h +++ b/include/epanet2.h @@ -115,7 +115,8 @@ typedef enum { EN_EFFICIENCY = 16, EN_HEADCURVE = 17, EN_EFFICIENCYCURVE = 18, - EN_PRICEPATTERN = 19 + EN_PRICEPATTERN = 19, + EN_STATE = 20 } EN_LinkProperty; /// Time parameter codes diff --git a/src/epanet.c b/src/epanet.c index 72f59d8..0fe9b72 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -1994,6 +1994,7 @@ int DLLEXPORT EN_getlinknodes(EN_Project *p, int index, int *node1, int DLLEXPORT EN_getlinkvalue(EN_Project *p, int index, EN_LinkProperty code, EN_API_FLOAT_TYPE *value) { double a, h, q, v = 0.0; int returnValue = 0; + int pmp; EN_Network *net = &p->network; hydraulics_t *hyd = &p->hydraulics; @@ -2007,8 +2008,6 @@ int DLLEXPORT EN_getlinkvalue(EN_Project *p, int index, EN_LinkProperty code, EN double *LinkFlows = hyd->LinkFlows; double *LinkSetting = hyd->LinkSetting; - - /* Check for valid arguments */ *value = 0.0; if (!p->Openflag) @@ -2124,6 +2123,20 @@ int DLLEXPORT EN_getlinkvalue(EN_Project *p, int index, EN_LinkProperty code, EN else v = 1.0; break; + + case EN_STATE: + v = hyd->LinkStatus[index]; + + if (Link[index].Type == EN_PUMP) { + pmp = findpump(net, index); + if (hyd->LinkStatus[index] >= OPEN) { + if (hyd->LinkFlows[index] > hyd->LinkSetting[index] * Pump[pmp].Qmax) + v = XFLOW; + if (hyd->LinkFlows[index] < 0.0) + v = XHEAD; + } + } + break; case EN_SETTING: if (Link[index].Type == EN_PIPE || Link[index].Type == EN_CVPIPE) { From 634651d4e6c79e3ebd556f29447a6235a304a8a4 Mon Sep 17 00:00:00 2001 From: Elad Salomons Date: Sat, 18 Aug 2018 21:54:52 +0300 Subject: [PATCH 16/17] Add option to get pump's constant power and speed Details in #215 --- include/epanet2.bas | 2 ++ include/epanet2.h | 4 +++- src/epanet.c | 18 ++++++++++++++++++ 3 files changed, 23 insertions(+), 1 deletion(-) diff --git a/include/epanet2.bas b/include/epanet2.bas index 585d556..497ca1c 100644 --- a/include/epanet2.bas +++ b/include/epanet2.bas @@ -57,6 +57,8 @@ Public Const EN_HEADCURVE = 17 Public Const EN_EFFICIENCYCURVE = 18 Public Const EN_PRICEPATTERN = 19 Public Const EN_STATE = 20 +Public Const EN_CONST_POWER = 21 +Public Const EN_SPEED = 22 Public Const EN_DURATION = 0 ' Time parameters Public Const EN_HYDSTEP = 1 diff --git a/include/epanet2.h b/include/epanet2.h index 3b26116..feb8a37 100644 --- a/include/epanet2.h +++ b/include/epanet2.h @@ -116,7 +116,9 @@ typedef enum { EN_HEADCURVE = 17, EN_EFFICIENCYCURVE = 18, EN_PRICEPATTERN = 19, - EN_STATE = 20 + EN_STATE = 20, + EN_CONST_POWER = 21, + EN_SPEED = 22 } EN_LinkProperty; /// Time parameter codes diff --git a/src/epanet.c b/src/epanet.c index 0fe9b72..aec91b9 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -2137,6 +2137,24 @@ int DLLEXPORT EN_getlinkvalue(EN_Project *p, int index, EN_LinkProperty code, EN } } break; + + case EN_CONST_POWER: + v = 0; + if (Link[index].Type == EN_PUMP) { + pmp = findpump(net, index); + if (Pump[pmp].Ptype == CONST_HP) { + v = Link[index].Km; // Power in HP + } + } + break; + + case EN_SPEED: + v = 0; + if (Link[index].Type == EN_PUMP) { + pmp = findpump(net, index); + v = Link[index].Kc; + } + break; case EN_SETTING: if (Link[index].Type == EN_PIPE || Link[index].Type == EN_CVPIPE) { From 761199b08e8b885d14a8800256d80890d901af96 Mon Sep 17 00:00:00 2001 From: Elad Salomons Date: Mon, 20 Aug 2018 21:10:49 +0300 Subject: [PATCH 17/17] Added ENgetcurvetype API Also update curve type on EN_setheadcurveindex. Co-Authored-By: milad ghiami --- include/epanet2.bas | 1 + include/epanet2.h | 14 ++++++++++++-- src/epanet.c | 19 ++++++++++++++++++- win_build/WinSDK/epanet2.def | 3 ++- 4 files changed, 33 insertions(+), 4 deletions(-) diff --git a/include/epanet2.bas b/include/epanet2.bas index 497ca1c..d782c6d 100644 --- a/include/epanet2.bas +++ b/include/epanet2.bas @@ -229,6 +229,7 @@ Public Const EN_G_CURVE = 4 ' General\default curve Declare Function ENgetcurve Lib "epanet2.dll" (ByVal curveIndex As Long, ByVal CurveID As String, nValues As Long, xValues As Any, yValues As Any) As Long Declare Function ENgetheadcurveindex Lib "epanet2.dll" (ByVal pumpIndex As Long, curveIndex As Long) As Long Declare Function ENgetpumptype Lib "epanet2.dll" (ByVal index As Long, PumpType As Long) As Long + Declare Function ENgetcurvetype Lib "epanet2.dll" (ByVal curveindex As Long, CurveType As Long) As Long Declare Function ENgetversion Lib "epanet2.dll" (value As Long) As Long diff --git a/include/epanet2.h b/include/epanet2.h index feb8a37..e070f37 100644 --- a/include/epanet2.h +++ b/include/epanet2.h @@ -763,10 +763,19 @@ extern "C" { @param linkIndex The index of the pump element @param[out] outType The integer-typed pump curve type signifier (output parameter) @return Error code - @see EN_CurveType + @see EN_PumpType */ int DLLEXPORT ENgetpumptype(int linkIndex, int *outType); - + + /** + @brief Get the type of a curve + @param curveIndex The index of the curve element + @param[out] outType The integer-typed curve curve type signifier (output parameter) + @return Error code + @see EN_CurveType + */ + int DLLEXPORT ENgetcurvetype(int curveIndex, int *outType); + /** @brief Get the version number. This number is to be interpreted with implied decimals, i.e., "20100" == "2(.)01(.)00" @param[out] version The version of EPANET @@ -1189,6 +1198,7 @@ extern "C" { int DLLEXPORT EN_getheadcurveindex(EN_Project *p, int pumpIndex, int *curveIndex); int DLLEXPORT EN_setheadcurveindex(EN_Project *p, int pumpIndex, int curveIndex); int DLLEXPORT EN_getpumptype(EN_Project *p, int linkIndex, int *outType); + int DLLEXPORT EN_getcurvetype(EN_Project *p, int curveIndex, int *outType); int DLLEXPORT EN_getversion(int *version); int DLLEXPORT EN_setcontrol(EN_Project *p, int cindex, int ctype, int lindex, EN_API_FLOAT_TYPE setting, int nindex, EN_API_FLOAT_TYPE level); int DLLEXPORT EN_setnodevalue(EN_Project *p, int index, int code, EN_API_FLOAT_TYPE v); diff --git a/src/epanet.c b/src/epanet.c index aec91b9..03d1300 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -435,6 +435,10 @@ int DLLEXPORT ENgetpumptype(int index, int *type) { return EN_getpumptype(_defaultModel, index, type); } +int DLLEXPORT ENgetcurvetype(int curveindex, int *type) { + return EN_getcurvetype(_defaultModel, curveindex, type); +} + int DLLEXPORT ENgetnumdemands(int nodeIndex, int *numDemands) { return EN_getnumdemands(_defaultModel, nodeIndex, numDemands); } @@ -3349,7 +3353,8 @@ int DLLEXPORT EN_setheadcurveindex(EN_Project *p, int index, int curveindex) { pump->Q0 /= Ucf[FLOW]; pump->Qmax /= Ucf[FLOW]; pump->Hmax /= Ucf[HEAD]; - + + p->network.Curve[curveindex].Type = P_CURVE; return (0); } @@ -3370,6 +3375,18 @@ int DLLEXPORT EN_getpumptype(EN_Project *p, int index, int *type) { return (0); } +int DLLEXPORT EN_getcurvetype(EN_Project *p, int curveindex, int *type) { + + EN_Network *net = &p->network; + + if (!p->Openflag) + return (102); + if (curveindex < 1 || curveindex > net->Ncurves) + return (206); + *type = net->Curve[curveindex].Type; + return (0); +} + /* ---------------------------------------------------------------- Functions for opening files diff --git a/win_build/WinSDK/epanet2.def b/win_build/WinSDK/epanet2.def index 09acdb7..e6a674a 100644 --- a/win_build/WinSDK/epanet2.def +++ b/win_build/WinSDK/epanet2.def @@ -92,4 +92,5 @@ EXPORTS ENaddlink = _ENaddlink@16 ENdeletelink = _ENdeletelink@4 ENdeletenode = _ENdeletenode@4 - ENsetlinktype = _ENsetlinktype@8 \ No newline at end of file + ENsetlinktype = _ENsetlinktype@8 + ENgetcurvetype = _ENgetcurvetype@8 \ No newline at end of file