JIANG/EPANET
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
7e6ceb436634a0ab5c81e200170ba7bfd472af86
EPANET/src/smatrix.c
Lew Rossman 86ffbcf614 Improved node re-ordering method (Issue #162)
2018-07-02 13:55:41 -04:00

920 lines
30 KiB
C
Executable File
Raw Blame History

/*
*******************************************************************
SMATRIX.C -- Sparse matrix routines for EPANET program.
VERSION: 2.00
DATE: 5/8/00
AUTHOR: L. Rossman
US EPA - NRMRL
This module contains the sparse matrix routines used to solve
a network's hydraulic equations. The entry points into this
module are:
createsparse() -- called from openhyd() in HYDRAUL.C
freesparse() -- called from closehyd() in HYDRAUL.C
linsolve() -- called from netsolve() in HYDRAUL.C
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. 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.
********************************************************************
*/
#include <stdio.h>
#include <string.h>
#ifndef __APPLE__
#include <malloc.h>
#else
#include <stdlib.h>
#endif
#include <math.h>
#include <limits.h>
#include <time.h>
#include "text.h"
#include "types.h"
#include "funcs.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)
/*
**--------------------------------------------------------------
** Input: none
** Output: returns error code
** Purpose: creates sparse representation of coeff. matrix
**--------------------------------------------------------------
*/
{
int errcode = 0;
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));
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);
} /* End of createsparse */
int allocsparse(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: returns error code
** Purpose: allocates memory for indexing the solution matrix
**--------------------------------------------------------------
*/
{
EN_Network *net = &pr->network;
solver_t *solver = &pr->hydraulics.solver;
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);
}
void freesparse(EN_Project *pr)
/*
**----------------------------------------------------------------
** Input: None
** Output: None
** Purpose: Frees memory used for sparse matrix storage
**----------------------------------------------------------------
*/
{
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(net->Adjlist);
FREE(solver->Order);
FREE(solver->Row);
FREE(solver->Ndx);
FREE(solver->XLNZ);
FREE(solver->NZSUB);
FREE(solver->LNZ);
} /* End of freesparse */
int buildlists(EN_Project *pr, int paraflag)
/*
**--------------------------------------------------------------
** Input: paraflag = TRUE if list marks parallel links
** Output: returns error code
** Purpose: builds linked list of links adjacent to each node
**--------------------------------------------------------------
*/
{
int i,j,k;
int pmark = 0;
int errcode = 0;
Padjlist alink;
EN_Network *net = &pr->network;
// 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 = 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 = net->Adjlist[j];
net->Adjlist[j] = alink;
}
return(errcode);
} /* End of buildlists */
int paralink(EN_Project *pr, int i, int j, int k)
/*
**--------------------------------------------------------------
** Input: i = index of start node of link
** j = index of end node of link
** k = link index
** Output: returns 1 if link k parallels another link, else 0
** Purpose: checks for parallel links between nodes i and j
**
**--------------------------------------------------------------
*/
{
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 */
void xparalinks(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: none
** Purpose: removes parallel links from nodal adjacency lists
**--------------------------------------------------------------
*/
{
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 <= net->Nnodes; i++)
{
alink = net->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
{
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 */
void freelists(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: none
** Purpose: frees memory used for nodal adjacency lists
**--------------------------------------------------------------
*/
{
int i;
Padjlist alink;
EN_Network *net = &pr->network;
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 */
void countdegree(EN_Project *pr)
/*
**----------------------------------------------------------------
** Input: none
** Output: none
** Purpose: counts number of nodes directly connected to each node
**----------------------------------------------------------------
*/
{
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.
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]++;
}
}
}
}
int reordernodes(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: returns 1 if successful, 0 if not
** Purpose: re-orders nodes to minimize # of non-zeros that
** will appear in factorized solution matrix
**--------------------------------------------------------------
*/
{
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++)
{
solver->Row[k] = k;
solver->Order[k] = k;
}
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 factorize(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: returns error code
** Purpose: symbolically factorizes the solution matrix in
** terms of its adjacency lists
**--------------------------------------------------------------
*/
{
int k, knode;
int errcode = 0;
EN_Network *net = &pr->network;
solver_t *solver = &pr->hydraulics.solver;
// 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)
/*
**--------------------------------------------------------------
** Input: knode = node index
** Output: returns 1 if successful, 0 if not
** Purpose: creates new entries in knode's adjacency list for
** all unlinked pairs of active nodes that are
** adjacent to knode
**--------------------------------------------------------------
*/
{
int node;
Padjlist alink;
EN_Network *net = &pr->network;
solver_t *solver = &pr->hydraulics.solver;
// 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 (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 */
int newlink(EN_Project *pr, Padjlist alink)
/*
**--------------------------------------------------------------
** Input: alink = element of node's adjacency list
** Output: returns 1 if successful, 0 if not
** Purpose: links end of current adjacent link to end nodes of
** all links that follow it on adjacency list
**--------------------------------------------------------------
*/
{
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 (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(net, inode, jnode, hyd->Ncoeffs)) return(0);
if (!addlink(net, jnode, inode, hyd->Ncoeffs)) return(0);
solver->Degree[inode]++;
solver->Degree[jnode]++;
}
}
}
return(1);
} /* End of newlink */
int linked(EN_Network *n, int i, int j)
/*
**--------------------------------------------------------------
** Input: i = node index
** j = node index
** Output: returns 1 if nodes i and j are linked, 0 if not
** Purpose: checks if nodes i and j are already linked.
**--------------------------------------------------------------
*/
{
Padjlist alink;
for (alink = n->Adjlist[i]; alink != NULL; alink = alink->next)
{
if (alink->node == j) return(1);
}
return(0);
} /* End of linked */
int addlink(EN_Network *net, int i, int j, int n)
/*
**--------------------------------------------------------------
** Input: i = node index
** j = node index
** n = link index
** Output: returns 1 if successful, 0 if not
** Purpose: augments node i's adjacency list with node j
**--------------------------------------------------------------
*/
{
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 */
int storesparse(EN_Project *pr, int n)
/*
**--------------------------------------------------------------
** Input: n = number of rows in solution matrix
** Output: returns error code
** Purpose: stores row indexes of non-zeros of each column of
** lower triangular portion of factorized matrix
**--------------------------------------------------------------
*/
{
Padjlist alink;
int i, ii, j, k, l, m;
int errcode = 0;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
solver_t *solver = &pr->hydraulics.solver;
/* 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
{
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;
}
return(errcode);
} /* End of storesparse */
int sortsparse(EN_Project *pr, int n)
/*
**--------------------------------------------------------------
** Input: n = number of rows in solution matrix
** Output: returns eror code
** Purpose: puts row indexes in ascending order in NZSUB
**--------------------------------------------------------------
*/
{
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(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];
// 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);
return(errcode);
} /* End of sortsparse */
void transpose(int n, int *il, int *jl, int *xl, int *ilt, int *jlt,
int *xlt, int *nzt)
/*
**---------------------------------------------------------------------
** Input: n = matrix order
** il,jl,xl = sparse storage scheme for original matrix
** nzt = work array
** Output: ilt,jlt,xlt = sparse storage scheme for transposed matrix
** Purpose: Determines sparse storage scheme for transpose of a matrix
**---------------------------------------------------------------------
*/
{
int i, j, k, kk;
for (i=1; i<=n; i++) nzt[i] = 0;
for (i=1; i<=n; i++)
{
for (k=il[i]; k<il[i+1]; k++)
{
j = jl[k];
kk = ilt[j] + nzt[j];
jlt[kk] = i;
xlt[kk] = xl[k];
nzt[j]++;
}
}
} /* End of transpose */
int linsolve(EN_Project *pr, int n)
/*
**--------------------------------------------------------------
** Input: s = solver struct
n = number of equations
** Output: s->F = solution values
** returns 0 if solution found, or index of
** equation causing system to be ill-conditioned
** Purpose: solves sparse symmetric system of linear
** equations using Cholesky factorization
**
** NOTE: This procedure assumes that the solution matrix has
** been symbolically factorized with the positions of
** the lower triangular, off-diagonal, non-zero coeffs.
** stored in the following integer arrays:
** XLNZ (start position of each column in NZSUB)
** NZSUB (row index of each non-zero in each column)
** LNZ (position of each NZSUB entry in Aij array)
**
** This procedure has been adapted from subroutines GSFCT and
** GSSLV in the book "Computer Solution of Large Sparse
** Positive Definite Systems" by A. George and J. W-H Liu
** (Prentice-Hall, 1981).
**--------------------------------------------------------------
*/
{
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;
int errcode = 0;
double bj, diagj, ljk;
double *temp;
temp = (double *) calloc(n+1, sizeof(double));
link = (int *) calloc(n+1,sizeof(int));
first = (int *) calloc(n+1,sizeof(int));
ERRCODE(MEMCHECK(temp));
ERRCODE(MEMCHECK(link));
ERRCODE(MEMCHECK(first));
if (errcode)
{
errcode = -errcode;
goto ENDLINSOLVE;
}
memset(temp,0,(n+1)*sizeof(double));
memset(link,0,(n+1)*sizeof(int));
/* Begin numerical factorization of matrix A into L */
/* Compute column L(*,j) for j = 1,...n */
for (j=1; j<=n; j++)
{
/* For each column L(*,k) that affects L(*,j): */
diagj = 0.0;
newk = link[j];
k = newk;
while (k != 0)
{
/* Outer product modification of L(*,j) by */
/* L(*,k) starting at first[k] of L(*,k). */
newk = link[k];
kfirst = first[k];
ljk = Aij[LNZ[kfirst]];
diagj += ljk*ljk;
istrt = kfirst + 1;
istop = XLNZ[k+1] - 1;
if (istop >= istrt)
{
/* Before modification, update vectors 'first' */
/* and 'link' for future modification steps. */
first[k] = istrt;
isub = NZSUB[istrt];
link[k] = link[isub];
link[isub] = k;
/* The actual mod is saved in vector 'temp'. */
for (i=istrt; i<=istop; i++)
{
isub = NZSUB[i];
temp[isub] += Aij[LNZ[i]]*ljk;
}
}
k = newk;
}
/* Apply the modifications accumulated */
/* in 'temp' to column L(*,j). */
diagj = Aii[j] - diagj;
if (diagj <= 0.0) /* Check for ill-conditioning */
{
errcode = j;
goto ENDLINSOLVE;
}
diagj = sqrt(diagj);
Aii[j] = diagj;
istrt = XLNZ[j];
istop = XLNZ[j+1] - 1;
if (istop >= istrt)
{
first[j] = istrt;
isub = NZSUB[istrt];
link[j] = link[isub];
link[isub] = j;
for (i=istrt; i<=istop; i++)
{
isub = NZSUB[i];
bj = (Aij[LNZ[i]] - temp[isub])/diagj;
Aij[LNZ[i]] = bj;
temp[isub] = 0.0;
}
}
} /* next j */
/* Foward substitution */
for (j=1; j<=n; j++)
{
bj = B[j]/Aii[j];
B[j] = bj;
istrt = XLNZ[j];
istop = XLNZ[j+1] - 1;
if (istop >= istrt)
{
for (i=istrt; i<=istop; i++)
{
isub = NZSUB[i];
B[isub] -= Aij[LNZ[i]]*bj;
}
}
}
/* Backward substitution */
for (j=n; j>=1; j--)
{
bj = B[j];
istrt = XLNZ[j];
istop = XLNZ[j+1] - 1;
if (istop >= istrt)
{
for (i=istrt; i<=istop; i++)
{
isub = NZSUB[i];
bj -= Aij[LNZ[i]]*B[isub];
}
}
B[j] = bj/Aii[j];
}
ENDLINSOLVE:
free(temp);
free(link);
free(first);
return(errcode);
} /* End of linsolve */
/************************ END OF SMATRIX.C ************************/
Reference in New Issue View Git Blame Copy Permalink