#pragma once
 
#include "stdafx.h"
 
extern FILE     *InFile,               /* Input file pointer           */
                *OutFile,              /* Output file pointer          */
                *RptFile,              /* Report file pointer          */
                *HydFile,              /* Hydraulics file pointer      */
                *TmpOutFile;           /* Temporary file handle        */
extern long     HydOffset,             /* Hydraulics file byte offset  */
                OutOffset1,            /* 1st output file byte offset  */
                OutOffset2;            /* 2nd output file byte offset  */
extern char     Msg[MAXMSG+1],         /* Text of output message       */
                InpFname[MAXFNAME+1],  /* Input file name              */
                Rpt1Fname[MAXFNAME+1], /* Primary report file name     */
                Rpt2Fname[MAXFNAME+1], /* Secondary report file name   */
                HydFname[MAXFNAME+1],  /* Hydraulics file name         */
                OutFname[MAXFNAME+1],  /* Binary output file name      */
                MapFname[MAXFNAME+1],  /* Map file name                */
                Title[MAXTITLE][MAXMSG+1], /* Problem title            */
                ChemName[MAXID+1],     /* Name of chemical             */
                ChemUnits[MAXID+1],    /* Units of chemical            */
                DefPatID[MAXID+1],     /* Default demand pattern ID    */
                Atime[13],             /* Clock time (hrs:min:sec)     */
                Hydflag,               /* Hydraulics flag              */
                Qualflag,              /* Water quality flag           */
                Unitsflag,             /* Unit system flag             */
                Flowflag,              /* Flow units flag              */
                Pressflag,             /* Pressure units flag          */
                Formflag,              /* Hydraulic formula flag       */
                Rptflag,               /* Report flag                  */
                Summaryflag,           /* Report summary flag          */
                Messageflag,           /* Error/warning message flag   */
                Statflag,              /* Status report flag           */
                Energyflag,            /* Energy report flag           */
                Nodeflag,              /* Node report flag             */
                Linkflag,              /* Link report flag             */
                Tstatflag,             /* Time statistics flag         */
                Warnflag,              /* Warning flag                 */
                Openflag,              /* Input processed flag         */
                OpenHflag,             /* Hydraul. system opened flag  */
                SaveHflag,             /* Hydraul. results saved flag  */
                OpenQflag,             /* Quality system opened flag   */
                SaveQflag,             /* Quality results saved flag   */
                Saveflag;              /* General purpose save flag    */
extern int      MaxNodes,              /* Node count from input file   */
                MaxLinks,              /* Link count from input file   */
                MaxJuncs,              /* Junction count               */
                MaxPipes,              /* Pipe count                   */
                MaxTanks,              /* Tank count                   */
                MaxPumps,              /* Pump count                   */
                MaxValves,             /* Valve count                  */
                MaxControls,           /* Control count                */
                MaxRules,              /* Rule count                   */
                MaxPats,               /* Pattern count                */
                MaxCurves,             /* Curve count                  */
                Nnodes,                /* Number of network nodes      */
                Ntanks,                /* Number of tanks              */
                Njuncs,                /* Number of junction nodes     */
                Nlinks,                /* Number of network links      */
                Npipes,                /* Number of pipes              */
                Npumps,                /* Number of pumps              */
                Nvalves,               /* Number of valves             */
                Ncontrols,             /* Number of simple controls    */
                Nrules,                /* Number of control rules      */
                Npats,                 /* Number of time patterns      */
                Ncurves,               /* Number of data curves        */
                Nperiods,              /* Number of reporting periods  */
                Ncoeffs,               /* Number of non-0 matrix coeffs*/
                DefPat,                /* Default demand pattern       */
                Epat,                  /* Energy cost time pattern     */
                MaxIter,               /* Max. hydraulic trials        */
                ExtraIter,             /* Extra hydraulic trials       */
                TraceNode,             /* Source node for flow tracing */
                PageSize,              /* Lines/page in output report  */
                CheckFreq,             /* Hydraulics solver parameter  */
                MaxCheck;              /* Hydraulics solver parameter  */
extern double   Ucf[MAXVAR],           /* Unit conversion factors      */
                Ctol,                  /* Water quality tolerance      */
                Htol,                  /* Hydraulic head tolerance     */
                Qtol,                  /* Flow rate tolerance          */
                RQtol,                 /* Flow resistance tolerance    */
                Hexp,                  /* Exponent in headloss formula */
                Qexp,                  /* Exponent in orifice formula  */
                Dmult,                 /* Demand multiplier            */
                Hacc,                  /* Hydraulics solution accuracy */
                BulkOrder,             /* Bulk flow reaction order     */
                WallOrder,             /* Pipe wall reaction order     */
                TankOrder,             /* Tank reaction order          */
                Kbulk,                 /* Global bulk reaction coeff.  */
                Kwall,                 /* Global wall reaction coeff.  */
                Climit,                /* Limiting potential quality   */
                Rfactor,               /* Roughness-reaction factor    */
                Diffus,                /* Diffusivity (sq ft/sec)      */
                Viscos,                /* Kin. viscosity (sq ft/sec)   */
                SpGrav,                /* Specific gravity             */
                Ecost,                 /* Base energy cost per kwh     */
                Dcost,                 /* Energy demand charge/kw/day  */
                Epump,                 /* Global pump efficiency       */
                Emax,                  /* Peak energy usage            */
                Dsystem,               /* Total system demand          */
                Wbulk,                 /* Avg. bulk reaction rate      */
                Wwall,                 /* Avg. wall reaction rate      */
                Wtank,                 /* Avg. tank reaction rate      */
                Wsource;               /* Avg. mass inflow             */
extern long     Tstart,                /* Starting time of day (sec)   */
                Hstep,                 /* Nominal hyd. time step (sec) */
                Qstep,                 /* Quality time step (sec)      */
                Pstep,                 /* Time pattern time step (sec) */
                Pstart,                /* Starting pattern time (sec)  */
                Rstep,                 /* Reporting time step (sec)    */
                Rstart,                /* Time when reporting starts   */
                Rtime,                 /* Next reporting time          */
                Htime,                 /* Current hyd. time (sec)      */
                Qtime,                 /* Current quality time (sec)   */
                Hydstep,               /* Actual hydraulic time step   */
                Rulestep,              /* Rule evaluation time step    */
                Dur;                   /* Duration of simulation (sec) */
extern SField   Field[MAXVAR];         /* Output reporting fields      */
extern char     *S;                    /* Link status                  */
extern char     *OldStat;              /* Previous link/tank status    */
extern double   *D,                    /* Node actual demand           */
                *C,                    /* Node actual quality          */
                *E,                    /* Emitter flows                */
                *K,                    /* Link settings                */
                *Q,                    /* Link flows                   */
                *R,                    /* Pipe reaction rate           */
                *X;                    /* General purpose array        */
extern double   *H;                    /* Node heads                   */
extern STmplist *Patlist;              /* Temporary time pattern list  */ 
extern STmplist *Curvelist;            /* Temporary list of curves     */
extern Spattern *Pattern;              /* Time patterns                */
extern Scurve   *Curve;                /* Curve data                   */
extern Snode    *Node;                 /* Node data                    */
extern Slink    *Link;                 /* Link data                    */
extern Stank    *Tank;                 /* Tank data                    */
extern Spump    *Pump;                 /* Pump data                    */
extern Svalve   *Valve;                /* Valve data                   */
extern Scontrol *Control;              /* Control data                 */
 
 
typedef struct HTentry *HTtable;
 
HTtable *HTcreate(void);
int     HTinsert(HTtable *, char *, int);
int 	HTfind(HTtable *, char *);
char    *HTfindKey(HTtable *, char *);
void	HTfree(HTtable *);
 
extern HTtable  *Nht, *Lht;            /* Hash tables for ID labels    */
extern Padjlist *Adjlist;              /* Node adjacency lists         */
extern double   *Aii,        /* Diagonal coeffs. of A               */
                *Aij,        /* Non-zero, off-diagonal coeffs. of A */
                *F;          /* Right hand side coeffs.             */
extern double   *P,          /* Inverse headloss derivatives        */
                *Y;          /* Flow correction factors             */
extern int      *Order,      /* Node-to-row of A                    */
                *Row,        /* Row-to-node of A                    */
                *Ndx;        /* Index of link's coeff. in Aij       */
/*
** The following arrays store the positions of the non-zero coeffs.    
** of the lower triangular portion of A whose values are stored in Aij:
*/
extern int      *XLNZ,       /* Start position of each column in NZSUB  */
                *NZSUB,      /* Row index of each coeff. in each column */
                *LNZ;        /* Position of each coeff. in Aij array    */
 
 
 
File : smatrix.c
 
#include <stdio.h>
#include <string.h>
#include <malloc.h>
#include <math.h>
#include "stdafx.h"
#include "hash.h"
#include "text.h"
#include "types.h"
#include "funcs.h"
#include "vars.h"
#define  EXTERN  extern
 
int      *Degree;     /* Number of links adjacent to each node  */
 
 
int  createsparse()
{
   int errcode = 0;
 
   /* Allocate data structures */
   ERRCODE(allocsparse());
   if (errcode) return(errcode);
   Degree = (int *) calloc(Nnodes+1, sizeof(int));
   ERRCODE(MEMCHECK(Degree));
   ERRCODE(buildlists(TRUE));
   if (!errcode)
   {
      xparalinks();    /* Remove parallel links */
      countdegree();   /* Find degree of each junction */
   }                   /* (= # of adjacent links)  */
 
    Ncoeffs = Nlinks;
   ERRCODE(reordernodes());
 
   ERRCODE(storesparse(Njuncs));
   if (!errcode) freelists();
   ERRCODE(ordersparse(Njuncs));
 
    ERRCODE(buildlists(FALSE));
 
    free(Degree);
   return(errcode);
}                      
 
int  allocsparse()
{
   int errcode = 0;
   Adjlist = (Padjlist *) calloc(Nnodes+1,  sizeof(Padjlist));
   Order  = (int *)   calloc(Nnodes+1,  sizeof(int));
   Row    = (int *)   calloc(Nnodes+1,  sizeof(int));
   Ndx    = (int *)   calloc(Nlinks+1,  sizeof(int));
   ERRCODE(MEMCHECK(Adjlist));
   ERRCODE(MEMCHECK(Order));
   ERRCODE(MEMCHECK(Row));
   ERRCODE(MEMCHECK(Ndx));
   return(errcode);
}
 
 
void  freesparse()
{
   freelists();
   free(Adjlist);
   free(Order);
   free(Row);
   free(Ndx);
   free(XLNZ);
   free(NZSUB);
   free(LNZ); 
}                        /* End of freesparse */
 
 
int  buildlists(int paraflag)
{
   int    i,j,k;
   int    pmark = 0;
   int    errcode = 0;
   Padjlist  alink;
 
   /* For each link, update adjacency lists of its end nodes */
   for (k=1; k<=Nlinks; k++)
   {
      i = Link[k].N1;
      j = Link[k].N2;
      if (paraflag) pmark = paralink(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 = Adjlist[i];
      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 = Adjlist[j];
      Adjlist[j] = alink;
   }
   return(errcode);
}                        /* End of buildlists */
 
 
int  paralink(int i, int j, int k)
{
   Padjlist alink;
   for (alink = Adjlist[i]; alink != NULL; alink = alink->next)
   {
      if (alink->node == j)     /* Link || to k (same end nodes) */
      {
         Ndx[k] = alink->link;  /* Assign Ndx entry to this link */
         return(1);
      }
   }
   Ndx[k] = k;                  /* Ndx entry if link not parallel */
   return(0);
}                        /* End of paralink */
 
 
void  xparalinks()
{
   int    i;
   Padjlist  alink,       /* Current item in adjacency list */
             blink;       /* Previous item in adjacency list */
 
   /* Scan adjacency list of each node */
   for (i=1; i<=Nnodes; i++)
   {
      alink = 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 */
            {
               Adjlist[i] = alink->next;
               free(alink);             /* Remove item from list */
               alink = 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()
{
   int   i;
   Padjlist alink;
 
   for (i=0; i<=Nnodes; i++)
   {
      for (alink = Adjlist[i]; alink != NULL; alink = Adjlist[i])
      {
         Adjlist[i] = alink->next;
         free(alink);
      }
   }
}                        /* End of freelists */
 
 
void  countdegree()
{
    int   i;
    Padjlist alink;
    memset(Degree,0,(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<=Njuncs; i++)
        for (alink = Adjlist[i]; alink != NULL; alink = alink->next)
            if (alink->node > 0) Degree[i]++;
}
 
 
int   reordernodes()
{
   int k, knode, m, n;
   for (k=1; k<=Nnodes; k++)
   {
      Row[k] = k;
      Order[k] = k;
   }
   n = Njuncs;
   for (k=1; k<=n; k++)                   /* Examine each junction    */
   {
      m = mindegree(k,n);                 /* Node with lowest degree  */
      knode = Order[m];                   /* Node's index             */
      if (!growlist(knode)) return(101);  /* Augment adjacency list   */
      Order[m] = Order[k];                /* Switch order of nodes    */
      Order[k] = knode;
      Degree[knode] = 0;                  /* In-activate node         */
   }
   for (k=1; k<=n; k++)                   /* Assign nodes to rows of  */
     Row[Order[k]] = k;                   /*   coeff. matrix          */
   return(0);
}                        /* End of reordernodes */
 
 
int  mindegree(int k, int n)
{
   int i, m;
   int min = n,
       imin = n;
 
   for (i=k; i<=n; i++)
   {
      m = Degree[Order[i]];
      if (m < min)
      {
         min = m;
         imin = i;
      }
   }
   return(imin);
}                        /* End of mindegree */
 
 
int  growlist(int knode)
{
   int   node;
   Padjlist alink;
 
   /* Iterate through all nodes connected to knode */
   for (alink = Adjlist[knode]; alink != NULL; alink = alink -> next)
   {
      node = alink->node;       /* End node of connecting link  */
      if (Degree[node] > 0)     /* End node is active           */
      {
         Degree[node]--;        /* Reduce degree of adjacency   */
         if (!newlink(alink))   /* Add to adjacency list        */
            return(0);
      }
   }
   return(1);
}                        /* End of growlist */
 
 
int  newlink(Padjlist alink)
{
   int   inode, jnode;
   Padjlist blink;
 
   /* 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 (Degree[jnode] > 0)        /* jnode still active */
      {
         if (!linked(inode,jnode))  /* inode not linked to jnode */
         {
 
            /* Since new connection represents a non-zero coeff. */
	    /* in the solution matrix, update the coeff. count.  */
            Ncoeffs++;
 
	    /* Update adjacency lists for inode & jnode to */
	    /* reflect the new connection.                 */
            if (!addlink(inode,jnode,Ncoeffs)) return(0);
            if (!addlink(jnode,inode,Ncoeffs)) return(0);
            Degree[inode]++;
            Degree[jnode]++;
         }
      }
   }
   return(1);
}                        /* End of newlink */
 
 
int  linked(int i, int j)
 
{
   Padjlist alink;
   for (alink = Adjlist[i]; alink != NULL; alink = alink->next)
      if (alink->node == j) return(1);
   return(0);
}                        /* End of linked */
 
 
int  addlink(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 = Adjlist[i];
   Adjlist[i] = alink;
   return(1);
}                        /* End of addlink */
 
 
int  storesparse(int n)
{
   Padjlist alink;
   int   i, ii, j, k, l, m;
   int   errcode = 0;
 
   /* Allocate sparse matrix storage */
   XLNZ  = (int *) calloc(n+2, sizeof(int));
   NZSUB = (int *) calloc(Ncoeffs+2, sizeof(int));
   LNZ   = (int *) calloc(Ncoeffs+2, sizeof(int));
   ERRCODE(MEMCHECK(XLNZ));
   ERRCODE(MEMCHECK(NZSUB));
   ERRCODE(MEMCHECK(LNZ));
   if (errcode) return(errcode);
 
   /* Generate row index pointers for each column of matrix */
   k = 0;
   XLNZ[1] = 1;
   for (i=1; i<=n; i++)             /* column */
   {
       m = 0;
       ii = Order[i];
       for (alink = Adjlist[ii]; alink != NULL; alink = alink->next)
       {
          j = Row[alink->node];    /* row */
          l = alink->link;
          if (j > i && j <= n)
          {
             m++;
             k++;
             NZSUB[k] = j;
             LNZ[k] = l;
          }
       }
       XLNZ[i+1] = XLNZ[i] + m;
   }
   return(errcode);
}                        /* End of storesparse */
 
 
int  ordersparse(int n)
 
{
   int  i, k;
   int  *xlnzt, *nzsubt, *lnzt, *nzt;
   int  errcode = 0;
 
   xlnzt  = (int *) calloc(n+2, sizeof(int));
   nzsubt = (int *) calloc(Ncoeffs+2, sizeof(int));
   lnzt   = (int *) calloc(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 ordersparse */
 
 
void  transpose(int n, int *il, int *jl, int *xl, int *ilt, int *jlt,
                int *xlt, int *nzt)
{
   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(int n, double *Aii, double *Aij, double *B)
{
   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);
}