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/ 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/include/epanet2.bas b/include/epanet2.bas index 0826a3b..d782c6d 100644 --- a/include/epanet2.bas +++ b/include/epanet2.bas @@ -56,6 +56,9 @@ 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_CONST_POWER = 21 +Public Const EN_SPEED = 22 Public Const EN_DURATION = 0 ' Time parameters Public Const EN_HYDSTEP = 1 @@ -154,6 +157,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 @@ -219,6 +229,7 @@ Public Const EN_CUSTOM = 2 ' user-defined custom 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 3b75915..c73d69b 100644 --- a/include/epanet2.h +++ b/include/epanet2.h @@ -115,7 +115,10 @@ typedef enum { EN_EFFICIENCY = 16, EN_HEADCURVE = 17, EN_EFFICIENCYCURVE = 18, - EN_PRICEPATTERN = 19 + EN_PRICEPATTERN = 19, + EN_STATE = 20, + EN_CONST_POWER = 21, + EN_SPEED = 22 } EN_LinkProperty; /// Time parameter codes @@ -224,8 +227,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 +234,6 @@ typedef enum { EN_RANGE = 4 } EN_StatisticType; - - typedef enum { EN_MIX1 = 0, /* Tank mixing models */ EN_MIX2 = 1, @@ -242,8 +241,6 @@ typedef enum { EN_LIFO = 3 } EN_MixingModel; - - typedef enum { EN_NOSAVE = 0, EN_SAVE = 1, @@ -251,16 +248,21 @@ typedef enum { EN_SAVE_AND_INIT = 11 } EN_SaveOption; - - typedef enum { EN_CONST_HP = 0, /* constant horsepower */ EN_POWER_FUNC = 1, /* power function */ - EN_CUSTOM = 2 /* user-defined custom curve */ + EN_CUSTOM = 2, /* user-defined custom curve */ + EN_NOCURVE = 3 /* no curve */ +} EN_PumpType; + +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) extern "C" { @@ -763,10 +765,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 @@ -1184,26 +1195,28 @@ extern "C" { int DLLEXPORT EN_getqualtype(EN_ProjectHandle ph, int *qualcode, int *tracenode); int DLLEXPORT EN_geterror(int errcode, char *errmsg, int maxLen); - int DLLEXPORT EN_getstatistic(EN_ProjectHandle ph, int code, EN_API_FLOAT_TYPE* value); - int DLLEXPORT EN_getnodeindex(EN_ProjectHandle ph, char *id, int *index); - int DLLEXPORT EN_getnodeid(EN_ProjectHandle ph, int index, char *id); - int DLLEXPORT EN_getnodetype(EN_ProjectHandle ph, int index, int *code); - int DLLEXPORT EN_getnodevalue(EN_ProjectHandle ph, int index, int code, EN_API_FLOAT_TYPE *value); - int DLLEXPORT EN_getcoord(EN_ProjectHandle ph, int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y); - int DLLEXPORT EN_setcoord(EN_ProjectHandle ph, int index, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y); - int DLLEXPORT EN_getnumdemands(EN_ProjectHandle ph, int nodeIndex, int *numDemands); - int DLLEXPORT EN_getbasedemand(EN_ProjectHandle ph, int nodeIndex, int demandIndex, EN_API_FLOAT_TYPE *baseDemand); - int DLLEXPORT EN_getdemandpattern(EN_ProjectHandle ph, int nodeIndex, int demandIndex, int *pattIndex); - int DLLEXPORT EN_getlinkindex(EN_ProjectHandle ph, char *id, int *index); - int DLLEXPORT EN_getlinkid(EN_ProjectHandle ph, int index, char *id); - int DLLEXPORT EN_getlinktype(EN_ProjectHandle ph, int index, EN_LinkType *code); - int DLLEXPORT EN_setlinktype(EN_ProjectHandle ph, char *id, EN_LinkType type); - int DLLEXPORT EN_getlinknodes(EN_ProjectHandle ph, int index, int *node1, int *node2); - int DLLEXPORT EN_getlinkvalue(EN_ProjectHandle ph, int index, EN_LinkProperty code, EN_API_FLOAT_TYPE *value); - int DLLEXPORT EN_getcurve(EN_ProjectHandle ph, int curveIndex, char* id, int *nValues, EN_API_FLOAT_TYPE **xValues, EN_API_FLOAT_TYPE **yValues); - int DLLEXPORT EN_getheadcurveindex(EN_ProjectHandle ph, int pumpIndex, int *curveIndex); - int DLLEXPORT EN_setheadcurveindex(EN_ProjectHandle ph, int pumpIndex, int curveIndex); - int DLLEXPORT EN_getpumptype(EN_ProjectHandle ph, int linkIndex, int *outType); + int DLLEXPORT EN_getstatistic(EN_Project *p, int code, EN_API_FLOAT_TYPE* value); + int DLLEXPORT EN_getnodeindex(EN_Project *p, char *id, int *index); + int DLLEXPORT EN_getnodeid(EN_Project *p, int index, char *id); + int DLLEXPORT EN_getnodetype(EN_Project *p, int index, int *code); + int DLLEXPORT EN_getnodevalue(EN_Project *p, int index, int code, EN_API_FLOAT_TYPE *value); + int DLLEXPORT EN_getcoord(EN_Project *p, int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y); + int DLLEXPORT EN_setcoord(EN_Project *p, int index, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y); + int DLLEXPORT EN_getnumdemands(EN_Project *p, int nodeIndex, int *numDemands); + int DLLEXPORT EN_getbasedemand(EN_Project *p, int nodeIndex, int demandIndex, EN_API_FLOAT_TYPE *baseDemand); + int DLLEXPORT EN_getdemandpattern(EN_Project *p, int nodeIndex, int demandIndex, int *pattIndex); + int DLLEXPORT EN_getlinkindex(EN_Project *p, char *id, int *index); + int DLLEXPORT EN_getlinkid(EN_Project *p, int index, char *id); + int DLLEXPORT EN_getlinktype(EN_Project *p, int index, EN_LinkType *code); + int DLLEXPORT EN_setlinktype(EN_Project *p, char *id, EN_LinkType type); + int DLLEXPORT EN_getlinknodes(EN_Project *p, int index, int *node1, int *node2); + int DLLEXPORT EN_getlinkvalue(EN_Project *p, int index, EN_LinkProperty code, EN_API_FLOAT_TYPE *value); + int DLLEXPORT EN_getcurve(EN_Project *p, int curveIndex, char* id, int *nValues, EN_API_FLOAT_TYPE **xValues, EN_API_FLOAT_TYPE **yValues); + 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_ProjectHandle ph, int cindex, int ctype, int lindex, EN_API_FLOAT_TYPE setting, int nindex, EN_API_FLOAT_TYPE level); int DLLEXPORT EN_setnodevalue(EN_ProjectHandle ph, int index, int code, EN_API_FLOAT_TYPE v); diff --git a/src/epanet.c b/src/epanet.c index 90b790f..7f297ee 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -436,6 +436,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); } @@ -2231,6 +2235,7 @@ int DLLEXPORT EN_getlinkvalue(EN_ProjectHandle ph, int index, EN_LinkProperty co EN_API_FLOAT_TYPE *value) { double a, h, q, v = 0.0; int returnValue = 0; + int pmp; EN_Project *p = (EN_Project*)ph; @@ -2246,8 +2251,6 @@ int DLLEXPORT EN_getlinkvalue(EN_ProjectHandle ph, int index, EN_LinkProperty co double *LinkFlows = hyd->LinkFlows; double *LinkSetting = hyd->LinkSetting; - - /* Check for valid arguments */ *value = 0.0; if (!p->Openflag) @@ -2363,6 +2366,38 @@ int DLLEXPORT EN_getlinkvalue(EN_ProjectHandle ph, int index, EN_LinkProperty co 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_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) { @@ -3187,7 +3222,7 @@ int DLLEXPORT EN_addcurve(EN_ProjectHandle ph, 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) @@ -3230,9 +3265,7 @@ int DLLEXPORT EN_setcurve(EN_ProjectHandle ph, int index, EN_API_FLOAT_TYPE *x, EN_Project *p = (EN_Project*)ph; EN_Network *net = &p->network; - Scurve *Curve = net->Curve; - - + Scurve *Curve = net->Curve; int j; /* Check for valid arguments */ @@ -3265,12 +3298,9 @@ int DLLEXPORT EN_setcurvevalue(EN_ProjectHandle ph, int index, int pnt, EN_API_F EN_Project *p = (EN_Project*)ph; EN_Network *net = &p->network; - Scurve *Curve = net->Curve; - const int Ncurves = net->Ncurves; - if (!p->Openflag) return (102); if (index <= 0 || index > Ncurves) @@ -3594,7 +3624,8 @@ int DLLEXPORT EN_setheadcurveindex(EN_ProjectHandle ph, int index, int curveinde pump->Q0 /= Ucf[FLOW]; pump->Qmax /= Ucf[FLOW]; pump->Hmax /= Ucf[HEAD]; - + + p->network.Curve[curveindex].Type = P_CURVE; return (0); } @@ -3617,6 +3648,18 @@ int DLLEXPORT EN_getpumptype(EN_ProjectHandle ph, 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 @@ -3986,7 +4029,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/funcs.h b/src/funcs.h index ef20044..b800083 100755 --- a/src/funcs.h +++ b/src/funcs.h @@ -196,24 +196,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 cd3fdf7..d7249fe 100755 --- a/src/types.h +++ b/src/types.h @@ -102,7 +102,7 @@ typedef int INT4; --------------------------------------------------------------------- */ #define MEMCHECK(x) (((x) == NULL) ? 101 : 0 ) -#define FREE(x) (free((x))) +#define FREE(x) if ((x)) free((x)) /* --------------------------------------------------------------------- @@ -153,8 +153,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 */ diff --git a/tools/before-test.cmd b/tools/before-test.cmd index 8871b44..70a7af9 100644 --- a/tools/before-test.cmd +++ b/tools/before-test.cmd @@ -25,8 +25,8 @@ set SCRIPT_HOME=%~dp0 set TEST_HOME=%~1 -set EXAMPLES_VER=1.0.2-dev -set BENCHMARK_VER=220dev-vs17 +set EXAMPLES_VER=1.0.2-dev.1 +set BENCHMARK_VER=220dev1 set TESTFILES_URL=https://github.com/OpenWaterAnalytics/epanet-example-networks/archive/v%EXAMPLES_VER%.zip diff --git a/tools/before-test.sh b/tools/before-test.sh index 8cb902b..72b1a7e 100755 --- a/tools/before-test.sh +++ b/tools/before-test.sh @@ -22,8 +22,8 @@ SCRIPT_HOME="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" TEST_HOME=$1 -EXAMPLES_VER="1.0.1" -BENCHMARK_VER="2012vs10" +EXAMPLES_VER="1.0.2-dev.1" +BENCHMARK_VER="220dev1" TEST_URL="https://github.com/OpenWaterAnalytics/epanet-example-networks/archive/v${EXAMPLES_VER}.tar.gz" BENCH_URL="https://github.com/OpenWaterAnalytics/epanet-example-networks/releases/download/v${EXAMPLES_VER}/epanet-benchmark-${BENCHMARK_VER}.tar.gz" diff --git a/tools/nrtest-epanet/main.py b/tools/nrtest-epanet/main.py index 1c0f41d..c962dab 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) @@ -120,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" @@ -130,7 +142,10 @@ 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" - result_compare(path_test, path_ref, (0.001, 0.0)) + benchmark_path = "C:\\Users\\mtryby\\Workspace\\GitRepo\\michaeltryby\\epanet-lr\\nrtestsuite\\benchmarks\\" + 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)) + 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 7bdf41b..a118695 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" @@ -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/report_diff.py b/tools/nrtest-epanet/report_diff.py new file mode 100644 index 0000000..b3f49e3 --- /dev/null +++ b/tools/nrtest-epanet/report_diff.py @@ -0,0 +1,96 @@ +# -*- 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 _binary_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. + ''' + 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): + + 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)) + + +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 diff --git a/tools/nrtest-epanet/setup.py b/tools/nrtest-epanet/setup.py index 897e621..4c7dfd7 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 ] @@ -24,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", 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" 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 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