JIANG/EPANET
1
0
Fork 0
Code Issues Pull Requests Actions 4 Packages Projects Releases Wiki Activity

Improved node re-ordering method (Issue #162)

Browse Source
This commit is contained in:
Lew Rossman
2018-07-02 13:55:41 -04:00
parent 55a7179ca6
commit 86ffbcf614
5 changed files with 1476 additions and 379 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
src/funcs.h
View File

@@ -193,24 +193,8 @@ void matrixcoeffs(EN_Project *pr); /* Finds hyd. matrix coeffs.
/* ----------- SMATRIX.C ---------------*/
int createsparse(EN_Project *pr); /* Creates sparse matrix */
int allocsparse(EN_Project *pr); /* Allocates matrix memory */
void freesparse(EN_Project *pr); /* Frees matrix memory */
int buildlists(EN_Project *pr, int); /* Builds adjacency lists */
int paralink(EN_Project *pr, int, int, int); /* Checks for parallel links */
void xparalinks(EN_Project *pr); /* Removes parallel links */
void freelists(EN_Project *pr); /* Frees adjacency lists */
void countdegree(EN_Project *pr); /* Counts links at each node */
int reordernodes(EN_Project *pr); /* Finds a node re-ordering */
int mindegree(solver_t *s, int, int); /* Finds min. degree node */
int growlist(EN_Project *pr, int); /* Augments adjacency list */
int newlink(EN_Project *pr, Padjlist); /* Adds fill-ins for a node */
int linked(EN_Network *net, int, int); /* Checks if 2 nodes linked */
int addlink(EN_Network *net, int, int, int); /* Creates new fill-in */
int storesparse(EN_Project *pr, int); /* Stores sparse matrix */
int ordersparse(hydraulics_t *h, int); /* Orders matrix storage */
void transpose(int,int *,int *, /* Transposes sparse matrix */
int *,int *,int *,int *,int *);
int linsolve(solver_t *s, int); /* Solves set of inear eqns. */
int linsolve(EN_Project *pr, int); /* Solves set of linear eqns. */
/* ----------- QUALITY.C ---------------*/
int openqual(EN_Project *pr); /* Opens WQ solver system */

1001
src/genmmd.c Normal file
View File

File diff suppressed because it is too large Load Diff

2
src/hydsolver.c
View File

@@ -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) {

520
src/smatrix.c
View File

@@ -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
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 <stdlib.h>
#endif
#include <math.h>
#include "hash.h"
#include <limits.h>
#include <time.h>
#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)
/*
@@ -56,11 +101,14 @@ int createsparse(EN_Project *pr)
*/
{
int errcode = 0;
EN_Network *n = &pr->network;
solver_t *s = &pr->hydraulics.solver;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
solver_t *solver = &pr->hydraulics.solver;
// cleartimer(SmatrixTimer);
// starttimer(SmatrixTimer);
/* Allocate data structures */
ERRCODE(allocsparse(pr));
@@ -70,37 +118,41 @@ int createsparse(EN_Project *pr)
}
/* 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) */
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. At same time, adjacency */
/* list is updated with links representing non-zero coeffs. */
hyd->Ncoeffs = n->Nlinks;
// Re-order nodes to minimize number of non-zero coeffs.
// in factorized solution matrix.
hyd->Ncoeffs = net->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));
// Factorize solution matrix by updating adjacency lists
// with non-zero connections due to fill-ins.
ERRCODE(factorize(pr));
/* Free memory used for adjacency lists and sort */
/* row indexes in NZSUB to optimize linsolve(). */
// 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));
ERRCODE(sortsparse(pr, net->Njuncs));
/* Re-build adjacency lists without removing parallel */
/* links for use in future connectivity checking. */
// Re-build adjacency lists without removing parallel
// links for use in future connectivity checking.
ERRCODE(buildlists(pr,FALSE));
/* Free allocated memory */
free(s->Degree);
// Free allocated memory
free(solver->Degree);
return(errcode);
} /* End of createsparse */
@@ -114,18 +166,18 @@ 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));
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);
FREE(net->Adjlist);
FREE(solver->Order);
FREE(solver->Row);
FREE(solver->Ndx);
FREE(solver->XLNZ);
FREE(solver->NZSUB);
FREE(solver->LNZ);
} /* End of freesparse */
@@ -167,34 +224,34 @@ int buildlists(EN_Project *pr, int paraflag)
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++)
// For each link, update adjacency lists of its end nodes
for (k=1; k <= net->Nlinks; k++)
{
i = n->Link[k].N1;
j = n->Link[k].N2;
i = net->Link[k].N1;
j = net->Link[k].N2;
if (paraflag) {
pmark = paralink(pr,i,j,k); /* Parallel link check */
pmark = paralink(pr, i, j, k); // Parallel link check
}
/* Include link in start node i's list */
// 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 */
else alink->node = 0; // Parallel link marker
alink->link = k;
alink->next = n->Adjlist[i];
n->Adjlist[i] = alink;
alink->next = net->Adjlist[i];
net->Adjlist[i] = alink;
/* Include link in end node j's list */
// 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 */
else alink->node = 0; // Parallel link marker
alink->link = k;
alink->next = n->Adjlist[j];
n->Adjlist[j] = alink;
alink->next = net->Adjlist[j];
net->Adjlist[j] = alink;
}
return(errcode);
} /* End of buildlists */
@@ -215,13 +272,16 @@ 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) */
// Link || to k (same end nodes)
if (alink->node == j)
{
pr->hydraulics.solver.Ndx[k] = alink->link; /* Assign Ndx entry to this link */
// Assign Ndx entry to this link
pr->hydraulics.solver.Ndx[k] = alink->link;
return(1);
}
}
pr->hydraulics.solver.Ndx[k] = k; /* Ndx entry if link not parallel */
// Ndx entry if link not parallel
pr->hydraulics.solver.Ndx[k] = k;
return(0);
} /* End of paralink */
@@ -236,35 +296,35 @@ 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;
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++)
// Scan adjacency list of each node
for (i=1; i <= net->Nnodes; i++)
{
alink = n->Adjlist[i]; /* First item in list */
alink = net->Adjlist[i]; // First item in list
blink = NULL;
while (alink != NULL)
{
if (alink->node == 0) /* Parallel link marker found */
if (alink->node == 0) // Parallel link marker found
{
if (blink == NULL) /* This holds at start of list */
if (blink == NULL) // This holds at start of list
{
n->Adjlist[i] = alink->next;
free(alink); /* Remove item from list */
alink = n->Adjlist[i];
net->Adjlist[i] = alink->next;
free(alink); // Remove item from list
alink = net->Adjlist[i];
}
else /* This holds for interior of list */
else // This holds for interior of list
{
blink->next = alink->next;
free(alink); /* Remove item from list */
free(alink); // Remove item from list
alink = blink->next;
}
}
else
{
blink = alink; /* Move to next item in list */
blink = alink; // Move to next item in list
alink = alink->next;
}
}
@@ -283,14 +343,14 @@ void freelists(EN_Project *pr)
{
int i;
Padjlist alink;
EN_Network *n = &pr->network;
EN_Network *net = &pr->network;
for (i=0; i <= n->Nnodes; i++)
for (i=0; i <= net->Nnodes; i++)
{
for (alink = n->Adjlist[i]; alink != NULL; alink = n->Adjlist[i])
for (alink = net->Adjlist[i]; alink != NULL; alink = net->Adjlist[i])
{
n->Adjlist[i] = alink->next;
net->Adjlist[i] = alink->next;
free(alink);
}
}
@@ -308,14 +368,15 @@ 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));
EN_Network *net = &pr->network;
/* NOTE: For purposes of node re-ordering, Tanks (nodes with */
/* indexes above Njuncs) have zero degree of adjacency. */
memset(pr->hydraulics.solver.Degree, 0, (net->Nnodes+1) * sizeof(int));
for (i=1; i <= n->Njuncs; i++) {
for (alink = n->Adjlist[i]; alink != NULL; alink = alink->next) {
// NOTE: For purposes of node re-ordering, Tanks (nodes with
// indexes above Njuncs) have zero degree of adjacency.
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;
// Default ordering
for (k=1; k <= net->Nnodes; k++)
{
s->Row[k] = k;
s->Order[k] = k;
solver->Row[k] = k;
solver->Order[k] = k;
}
n = net->Njuncs;
for (k=1; k<=n; k++) /* Examine each junction */
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)
{
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 */
// 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++;
}
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 */
xadj[k+1] = m;
}
return(0);
// 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++)
// 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
{
m = s->Degree[s->Order[i]];
if (m < min)
knode = solver->Order[k]; // Re-ordered index
if (!growlist(pr, knode)) // Augment adjacency list
{
min = m;
imin = i;
errcode = 101;
break;
}
solver->Degree[knode] = 0; // In-activate node
}
return(imin);
} /* End of mindegree */
return(errcode);
} /* End of factorize */
int growlist(EN_Project *pr, int knode)
@@ -402,17 +513,17 @@ int growlist(EN_Project *pr, int knode)
{
int node;
Padjlist alink;
EN_Network *n = &pr->network;
solver_t *s = &pr->hydraulics.solver;
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)
// Iterate through all nodes connected to knode
for (alink = net->Adjlist[knode]; alink != NULL; alink = alink -> next)
{
node = alink->node; /* End node of connecting link */
if (s->Degree[node] > 0) /* End node is active */
node = alink->node; // End node of connecting link
if (solver->Degree[node] > 0) // End node is active
{
s->Degree[node]--; /* Reduce degree of adjacency */
if (!newlink(pr,alink)) { /* Add to adjacency list */
solver->Degree[node]--; // Reduce degree of adjacency
if (!newlink(pr, alink)) { // Add to adjacency list
return(0);
}
}
@@ -433,35 +544,32 @@ int newlink(EN_Project *pr, Padjlist alink)
{
int inode, jnode;
Padjlist blink;
EN_Network *n = &pr->network;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
solver_t *s = &pr->hydraulics.solver;
solver_t *solver = &pr->hydraulics.solver;
/* Scan all entries in adjacency list that follow anode. */
inode = alink->node; /* End node of connection to anode */
// 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 */
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 */
// 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(n, inode,jnode)) { /* inode not linked to jnode */
/* Since new connection represents a non-zero coeff. */
/* in the solution matrix, update the coeff. count. */
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);
}
if (!addlink(n,jnode,inode,hyd->Ncoeffs)) {
return(0);
}
s->Degree[inode]++;
s->Degree[jnode]++;
// 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]++;
}
}
}
@@ -480,10 +588,9 @@ 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);
}
for (alink = n->Adjlist[i]; alink != NULL; alink = alink->next)
{
if (alink->node == j) return(1);
}
return(0);
} /* End of linked */
@@ -527,43 +634,43 @@ int storesparse(EN_Project *pr, int n)
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
solver_t *s = &pr->hydraulics.solver;
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);
}
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 */
// Generate row index pointers for each column of matrix
k = 0;
s->XLNZ[1] = 1;
for (i=1; i<=n; i++) { /* column */
solver->XLNZ[1] = 1;
for (i=1; i<=n; i++) // column
{
m = 0;
ii = s->Order[i];
ii = solver->Order[i];
for (alink = net->Adjlist[ii]; alink != NULL; alink = alink->next)
{
j = s->Row[alink->node]; /* row */
j = solver->Row[alink->node]; // row
l = alink->link;
if (j > i && j <= n) {
if (j > i && j <= n)
{
m++;
k++;
s->NZSUB[k] = j;
s->LNZ[k] = l;
solver->NZSUB[k] = j;
solver->LNZ[k] = l;
}
}
s->XLNZ[i+1] = s->XLNZ[i] + m;
solver->XLNZ[i+1] = solver->XLNZ[i] + m;
}
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
@@ -575,40 +682,45 @@ int ordersparse(hydraulics_t *h, int n)
int i, k;
int *xlnzt, *nzsubt, *lnzt, *nzt;
int errcode = 0;
solver_t *s = &h->solver;
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));
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 = s->XLNZ[i]; k < s->XLNZ[i+1]; k++) nzt[s->NZSUB[k]]++;
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];
/* 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);
// 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);
}
/* Reclaim memory */
free(xlnzt);
free(nzsubt);
free(lnzt);
free(nzt);
// Reclaim memory
FREE(xlnzt);
FREE(nzsubt);
FREE(lnzt);
FREE(nzt);
return(errcode);
} /* End of ordersparse */
} /* End of sortsparse */
void transpose(int n, int *il, int *jl, int *xl, int *ilt, int *jlt,
@@ -640,7 +752,7 @@ void transpose(int n, int *il, int *jl, int *xl, int *ilt, int *jlt,
} /* End of transpose */
int linsolve(solver_t *s, int n)
int linsolve(EN_Project *pr, int n)
/*
**--------------------------------------------------------------
** Input: s = solver struct
@@ -666,13 +778,13 @@ int linsolve(solver_t *s, int n)
**--------------------------------------------------------------
*/
{
double *Aii = s->Aii;
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;

2
src/types.h
View File

@@ -100,7 +100,7 @@ typedef int INT4;
---------------------------------------------------------------------
*/
#define MEMCHECK(x) (((x) == NULL) ? 101 : 0 )
#define FREE(x) (free((x)))
#define FREE(x) if ((x)) free((x))
/*
---------------------------------------------------------------------