{
From: horridge@ITS-MENZ.cc.monash.edu.au (Mark Horridge)
Sparse Solver for Borland or Turbo Pascal
*****************************************

Several people asked for a copy of my Pascal Unit to solve sparse linear
systems.  Here it is.

Unit Solspar is designed to solve sparse linear systems.
These are equation systems of the form:

A X  = R

where A is a matrix of coefficients size NxN,
      R is a vector of constants size N,
and   X is a vector of variables.

The code has been in use for a number of years, being slightly adapted for
each new purpose.  This is the first attempt at a portable general purpose
version.

I expect that the code contains a few errors which are my fault - I would
be glad to hear of these.  Please check that you have not tried to solve a
singular matrix by mistake !


Below find 3 text files, separated by rows of asterisks:

Solspar.pas:
  See notes at the top for instructions.
Test1.pas:
  Simple example of using routines with medium size matrix
Test2.pas:
  Example of using routines with N = 4, showing some possible errors
{************************************************************************}
{$N+,E+}
unit SolSpar;
{ Unit Solspar is designed to solve sparse linear systems.
These are equation systems of the form:

A X  = R

where A is a matrix of coefficients size NxN,
      R is a vector of constants size N,
and   X is a vector of variables.
(remember, in matrix A, 'rows' are equations and 'cols' are variables)

Solspar computes the vector X which satisfies the above equation.
It is designed to save time and space where A is fairly large but contains
few non-zero elements.

It is designed to work in either the real, protected-mode DOS, or Windows
versions of Borland/Turbo Pascal.

Limitations of Solspar:
  It only allows one Right Hand Side (RHS) column.
  It only allows one sparse structure at a time.
  Not to be used with Mark and Release !
Points to watch:
  Normally, the A matrix grows in size (becomes less sparse)
  during the Solve phase. The biggest matrix you can store
 will be to big to solve.
}
interface

 {Unit consists of the following 5 Boolean Functions and 3 procedures.
  Each function returns True if operation was successfully completed -
  otherwise False.  If False is returned, call the procedure
  GetErrorMessage to find the reason.  Finally, call the procedure
  ReleaseStruc to return memory to the heap.}

function InitStruc(NumEq : Word) : Boolean;
 {creates and initializes sparse matrix structure - call this first.
  NumEq is number of equations/variables.}

function AddElement(ThisEqu, ThisVar : Word; ThisVal : Double) : Boolean;
 {add an element to sparse matrix for equation ThisEqu and variable ThisVar;
  if such an entry already exists, ThisVal will be added to existing value
 You can add elements in any order, but the routine will be a bit more
efficient if you add variables (columns) to any particular equation in
ascending order, and then set the RHS}

function SetRHS(ThisEqu : Word; ThisVal : Double) : Boolean;
 {Set RHS for equation ThisEqu; if RHS has already been set, ThisVal will be
  added to existing value}

function Solve1 : Boolean;
  {calculate solutions; sparse matrix is destroyed}

function GetAnswer(ThisVar : Word; var ThisVal : Double) : Boolean;
 {read solution for variable ThisVar - probably called for each variable in
  turn}

procedure ReleaseStruc;
 {releases remaining memory used by sparse matrix structure - call this
  last}

procedure GetErrorMsg(var S : String; var N1, N2, N3 : Word);
 {N1: error number; S: Error Description; N1, N2 : other possibly useful
  numbers}

procedure showmat;      { displays small sparse matrix - not for Windows}

var
  SparMemUsed    : LongInt; {no of bytes of heap used by routines}

(*
procedure GetErrorMsg is used to obtain information about any problem:

procedure GetErrorMsg(var S : String; var N1, N2, N3 : Word);
 {N1: error number; S: Error Description; N1, N2 : other possibly useful
  numbers}

List of possible error messages is as follows:
(remember, 'rows' are equations and 'cols' are variables)

N1        S                          N2        N3        Comment

 0                                 0          0      No error: default
 1  Empty Row                      Row        0
 2  Row without Variables          Row        0      Equation with RHS only
 3  Empty Col                      Col        0
 4  Two Singles                    Row      PivotCol Two variables (the 2nd
                                                     is pivotcol) are each
                                                     mentioned only in the
                                                     same equation
 5  No RHS                         Row      LastCol
 6  Numerically Singular         PivotRow PivotStep  Columns (or rows)
                                                     are not linearly
                                                     independent.
 7  Out of Space                   0          1      called from InitStruc
 7  Out of Space                   0          2      called from InitStruc
 7  Out of Space                   0          3      called from AddElement
 7  Out of Space                   PivotStep  4      called from Solve1
 8  Cols out of Order              Row        Col
25  Initialize without releasing   0          0
26  Too many equations             NumEq      MaxEq
27  Too few equations              NumEq      0
30  Row < 1                        ThisEqu  ThisVar  called from AddElement
31  Col < 1                        ThisEqu  ThisVar  called from AddElement
32  Row > Neq                      ThisEqu    Neq    called from AddElement
33  Col > Neq                      ThisVar    Neq    called from AddElement
41  VarNo < 1                      ThisVar    0      called from AddElement
42  VarNo > Neq                    ThisVar    Neq    called from AddElement
 *)

{  Hopefully, this is as far as you need to read in order to use the unit.
   There are a few comments in the code which you can look at for interest.
  To learn more about sparse matrices read:
    Direct Methods for Sparse Matrices
    Duff, Erisman and Reid         Cambridge University Press  }


implementation

const
  DBug           = False;
  {due to Borland's "constant folding" , Debug statements cost nothing with
Dbug false}
  Msg            = False;
  ParaSize       = 16;  {for paragraph alignment}
  MaxSize        = 65520; {largest variable size allowed by Turbo Pascal}
  MaxUsable      = MaxSize-ParaSize; {allows for paragraph alignment of data}
  MaxEq          = (MaxUsable div SizeOf(Double))-1; {-1 for 0-based arrays}
  {based on above, MaxEq = 8187}
  Uvalue         = 0.1; {number 0<U<1; if larger, time and memory usage are
                          increased at the gain, maybe, of accuracy}

type
  RealArr        = array[0..MaxEq] of Double;
  WordArr        = array[0..MaxEq] of Word;
  IntArr         = array[0..MaxEq] of Integer;
  PRealArr       = ^RealArr;
  PWordArr       = ^WordArr;
  PIntArr        = ^IntArr;

  { The ElmRec record type stores a single entry in the sparse matrix.  For
  efficiency reasons, it is 16 bytes long, and will be paragraph-aligned.
  The CheckRow field is padding to make up the 16 bytes.  The 'preload
  PrevPtr' trick, see below, requires that Next field be 4 bytes [= size of a
  pointer] after the start of the record.  Hence, order of fields is
  important.
    By converting the value field to a Single, and eliminating the Checkrow
  field, the record size could be reduced to 10 bytes.  The paragraph
  alignment feature would then have to be sacrificed, and the code altered in
  various places.
    The nodes in use are arranged in a series of linked lists, one for each
  equation (or row).  Within each list variables (columns) always appear in
  ascending order.  Each list is terminated by a RHS or constant entry; this
  is treated as though it corresponded to an extra variable, numbered (Neq+1).}

  PElmRec        = ^ElmRec;
  ElmRec         = record
                     Column         : Word; {variable no. of this node}
                     CheckRow       : Word; {equation no. of this node}
                     Next           : PElmRec; {pointer to next node in row}
                     Value          : Double; {coefficient value}
                   end;

  PtrArr         = array[0..MaxEq] of PElmRec;
  PPtrArr        = ^PtrArr;

  { Spare nodes, not currently in use, are linked together in one list,
  pointed to by FreePtr.  FreeCount is the number of spare nodes.  Because
  there may be very many nodes, all of the same size, they are not allocated
  directly on the heap by the System Heap Manager.  Instead a more efficient
  scheme is used.  When the list of free nodes needs to be expanded, a large
  block of memory is requested, sufficient for many nodes. Another linked list
  is used to store the addresses of these blocks, for later disposal to the
  system heap}

const
  ElmSize        = SizeOf(ElmRec);
  MaxElmsPerBlock = MaxUsable div ElmSize;
  ElmsPerBlock   = MaxElmsPerBlock; {could be set to less}
type
  TBlock         = array[1..ElmsPerBlock] of ElmRec;
  PBlock         = ^TBlock;
  PHeapRec       = ^HeapRec;
  HeapRec        = record
                     BlockPtr       : PBlock; {address of block}
                     NextRec        : PHeapRec; {address of next list item}
                   end;

var                     {this list contains variables which must have unit-wide
scope}
  OldHeapError   : Pointer; {to restore previous HeapError function}
  Reason         : String; {for error messages}
  ErrNo1, ErrNo2, ErrNo3 : Word; {for error messages}
  Answer         : PRealArr; {holds solution}
  FreePtr        : PElmRec; {points to list of free nodes}
  BlockList      : PHeapRec; {points to list of allocated node blocks}
  FreeCount      : Word; {number of free nodes}
  Neq            : Word; {number of equations}
  FirstElm       : PPtrArr; {array of pointers to first node in each equation}
  LastElm        : PPtrArr; {array of pointers to last node in each equation}
  OrigFirstElm   : PPtrArr;
  OrigLastElm    : PPtrArr;
  OrigAnswer     : PRealArr;
  OrigNextActiveRow : PIntArr;
  OrigColCount   : PIntArr;
  OrigRowCount   : PWordArr;
  OrigPivRow     : PWordArr;
const
  Initialized    : Boolean = False;

  procedure Assert(P : Boolean; S : String);
    {used for debugging}
  begin
    if P then Exit;
    WriteLn('assertion failed: ', S);
    Halt(1);
  end;                  {Assert}

  function HeapFunc(Size : Word) : Integer; far; {allows trapping GetMem
errors}
  begin HeapFunc := 1; end;


  function GetMemParaAlign(var P, OrigP : Pointer; Size : Word) : Boolean;
  { Returns a pointer P which is paragraph aligned (a multiple of ParaSize)
  and which points to a new block of memory of which at least Size bytes can
  be used.  The 486 cache is well suited to paragraph aligned data.
    In more detail, GetMemParaAlign obtains a block of memory from the system
  heap which is ParaSize bytes larger than Size.  OrigP points to this
  original block.  P is the first address after OrigP which is a multiple of
  ParaSize.  OrigP must be saved for passing (with Size) to a later
  FreeMemParaAlign call.}
  type                  {used only for typecasting}
    PtrRec         = record OfsWord, SegWord : Word; end;
  var
    Ofset          : Word;
  begin
    if Msg then WriteLn('Entering GetMemParaAlign');
    P := nil; OrigP := nil;
    GetMemParaAlign := False;
    OldHeapError := HeapError; HeapError := @HeapFunc;
    GetMem(OrigP, Size+ParaSize);
    HeapError := OldHeapError;
    if (OrigP = nil) then Exit;
    Inc(SparMemUsed, Size+ParaSize);
    P := OrigP;         {to load segment}
    Ofset := PtrRec(P).OfsWord;
    {adjust offset to paragraph boundary}
    PtrRec(P).OfsWord := ParaSize+ParaSize*(Ofset div ParaSize);
    GetMemParaAlign := True;
  end;                  {GetMemParaAlign}

  procedure FreeMemParaAlign(var OrigP : Pointer; Size : Word);
 {If OrigP<>Nil, returns a block at OrigP of (Size+ParaSize) bytes to the
  system heap.  Failure will cause a runtime error.  If this occurs, it is
  most likely due to an error in this unit!  }
  begin
    if Msg then WriteLn('Entering FreeMemParaAlign');
    if (OrigP = nil) then Exit; {already freed or never allocated, we presume}
    FreeMem(OrigP, Size+ParaSize); {failure will cause runtime error}
    Dec(SparMemUsed, Size+ParaSize);
    OrigP := nil;       {to indicate now freed}
  end;                  {FreeMemParaAlign}

  function FrePrime : Boolean;
  type                  {used only for typecasting}
    PtrRec         = record OfsWord, SegWord : Word; end;
  var
    NewBlock, OrigBlock : PBlock;
    NewRec         : PHeapRec;
    Ofset, Count   : Word;
    P              : PElmRec;
  begin
    if Msg then WriteLn('Entering FrePrime');
    FrePrime := False;

    {add new node to list of blocks}
    OldHeapError := HeapError; HeapError := @HeapFunc;
    NewRec := nil;
    New(NewRec); Inc(SparMemUsed, SizeOf(HeapRec));
    HeapError := OldHeapError;
    if (NewRec = nil) then Exit;
    NewRec^.NextRec := BlockList;
    NewRec^.BlockPtr := nil;
    BlockList := NewRec;


    {get new block}
    if not GetMemParaAlign(Pointer(NewBlock), Pointer(OrigBlock),
SizeOf(TBlock)) then Exit;
    NewRec^.BlockPtr := OrigBlock;
    {fill new block with linked nodes}
    P := PElmRec(NewBlock);
    Ofset := PtrRec(P).OfsWord;
    for Count := 1 to ElmsPerBlock do begin
      Inc(Ofset, ElmSize);
      PtrRec(P).OfsWord := Ofset;
      NewBlock^[Count].Next := P;
    end;
    Inc(FreeCount, ElmsPerBlock);

    NewBlock^[ElmsPerBlock].Next := FreePtr; {point end of new block to
FreePtr}
    FreePtr := PElmRec(NewBlock); {point FreePtr at start of new block}

    FrePrime := True;
  end;                  {FrePrime}

  procedure SetErrorMsg(S : String; N1, N2, N3 : Word);
  begin
    Reason := S; ErrNo1 := N1; ErrNo2 := N2; ErrNo3 := N3;
  end;                  {SetErrorMsg}

  procedure GetErrorMsg(var S : String; var N1, N2, N3 : Word);
  begin
    S := Reason; N1 := ErrNo1; N2 := ErrNo2; N3 := ErrNo3;
  end;                  {GetErrorMsg}

  procedure FillStruct(var Dest; Count : Word; var Filler; FillerSize : Word);
    {-Fill memory starting at Dest with Count instances of Filler}
  inline(
    $58/$5B/$5A/$59/$5F/$07/$E3/$11/$FC/$1E/$8E/$DA/
    $89/$CA/$89/$DE/$89/$C1/$F2/$A4/$89/$D1/$E2/$F4/$1F);

  function InitStruc(NumEq : Word) : Boolean;
  var
    Col            : Word;
  label Fail;
  begin
    if Msg then WriteLn('Entering InitStruc');
    InitStruc := False;
    SetErrorMsg('', 0, 0, 0);
    SparMemUsed := 0;

    OrigAnswer := nil;
    OrigColCount := nil;
    OrigFirstElm := nil;
    OrigLastElm := nil;
    OrigNextActiveRow := nil;
    OrigPivRow := nil;
    OrigRowCount := nil;

    if Initialized then begin
      SetErrorMsg('Initialize without releasing ', 25, 0, 0);
      Exit;
    end;
    Initialized := True;
    if (NumEq > MaxEq) then begin
      SetErrorMsg('Too many equations ', 26, NumEq, MaxEq);
      Exit;
    end;
    if (NumEq < 1) then begin
      SetErrorMsg('Too few equations ', 27, NumEq, 0);
      Exit;
    end;

    Neq := NumEq;
    if not GetMemParaAlign(Pointer(FirstElm), Pointer(OrigFirstElm),
(1+NumEq)*SizeOf(PElmRec)) then begin
      SetErrorMsg('Out of Space', 7, 0, 1);
      Exit;
    end;
    if not GetMemParaAlign(Pointer(LastElm), Pointer(OrigLastElm),
(1+NumEq)*SizeOf(PElmRec)) then begin
      SetErrorMsg('Out of Space', 7, 0, 2);
      Exit;
    end;

    BlockList := nil;
    FreePtr := nil;
    Col := 0;
    FreeCount := 0;
    FillStruct(FirstElm^, 1+Neq, FreePtr, SizeOf(FreePtr));
    FillStruct(LastElm^, 1+Neq, FreePtr, SizeOf(FreePtr));
    InitStruc := True;
  end;                  {InitStruc}

  procedure ReleaseStruc;
  var
    OrigBlock      : PBlock;
    NextPtr        : PHeapRec;
  begin
    if Msg then WriteLn('Entering ReleaseStruc');
    Initialized := False;

    {some of these may have been released already}
    FreeMemParaAlign(Pointer(OrigAnswer), (1+Neq)*SizeOf(Double));
    FreeMemParaAlign(Pointer(OrigColCount), (1+Neq)*SizeOf(Integer));
    FreeMemParaAlign(Pointer(OrigFirstElm), (1+Neq)*SizeOf(PElmRec));
    FreeMemParaAlign(Pointer(OrigLastElm), (1+Neq)*SizeOf(PElmRec));
    FreeMemParaAlign(Pointer(OrigNextActiveRow), (1+Neq)*SizeOf(Integer));
    FreeMemParaAlign(Pointer(OrigPivRow), (1+Neq)*SizeOf(Word));
    FreeMemParaAlign(Pointer(OrigRowCount), (1+Neq)*SizeOf(Word));


    {get rid of user heap}
    while (BlockList <> nil) do begin
      OrigBlock := BlockList^.BlockPtr;
      if (OrigBlock <> nil) then begin
        FreeMemParaAlign(Pointer(OrigBlock), SizeOf(TBlock));
        if Msg then WriteLn('releasing a block');
      end;
      NextPtr := BlockList^.NextRec;
      Dispose(BlockList); Dec(SparMemUsed, SizeOf(HeapRec));
      BlockList := NextPtr;
    end;                {while}

    Initialized := False;
  end;                  {ReleaseStruc}

  function AddElement(ThisEqu, ThisVar : Word; ThisVal : Double) : Boolean;
  var
    PrevPtr, ElmPtr, NewPtr : PElmRec;
  begin                 {PivRow[Row] points to last element }
    AddElement := False;
    if (ThisEqu < 1) then begin SetErrorMsg('Row < 1', 30, ThisEqu, ThisVar);
Exit; end;
    if (ThisVar < 1) then begin SetErrorMsg('Col < 1', 31, ThisEqu, ThisVar);
Exit; end;
    if (ThisEqu > Neq) then begin SetErrorMsg('Row > Neq', 32, ThisEqu, Neq);
Exit; end;
    if (ThisVar > (1+Neq)) then begin SetErrorMsg('Col > Neq', 33, ThisVar,
Neq); Exit; end;
    if ThisVal = 0.0 then begin AddElement := True; Exit; end;
    if (FreeCount < Neq) then begin
      if not FrePrime then begin
        SetErrorMsg('Out of Space', 7, 0, 3); Exit;
      end;
    end;
    NewPtr := FreePtr; FreePtr := FreePtr^.Next; Dec(FreeCount);

    NewPtr^.Value := ThisVal;
    NewPtr^.Column := ThisVar;
    NewPtr^.Next := nil;

    if FirstElm^[ThisEqu] = nil then begin
      FirstElm^[ThisEqu] := NewPtr;
      LastElm^[ThisEqu] := NewPtr;
    end
    else if (LastElm^[ThisEqu]^.Column < ThisVar) then begin
      LastElm^[ThisEqu]^.Next := NewPtr;
      LastElm^[ThisEqu] := NewPtr;
    end
    else begin          {insertion sort}
      ElmPtr := FirstElm^[ThisEqu];
      PrevPtr := nil;
      while (ElmPtr^.Column < ThisVar) do begin
        PrevPtr := ElmPtr;
        ElmPtr := ElmPtr^.Next;
      end;
      if (ElmPtr^.Column = ThisVar) then begin
        ElmPtr^.Value := ElmPtr^.Value+ThisVal;
        {new node not needed: return to free list}
        NewPtr^.Next := FreePtr; FreePtr := NewPtr; Inc(FreeCount);
      end
      else {ElmPtr^.Column > ThisVar} begin
        if PrevPtr = nil then FirstElm^[ThisEqu] := NewPtr
        else PrevPtr^.Next := NewPtr;
        NewPtr^.Next := ElmPtr;
      end;
    end;
    AddElement := True;
  end;                  {AddElement}

  function SetRHS(ThisEqu : Word; ThisVal : Double) : Boolean;
  begin
    SetRHS := AddElement(ThisEqu, 1+Neq, ThisVal);
  end;                  {SetRHS}

  function GetAnswer(ThisVar : Word; var ThisVal : Double) : Boolean;
    {should fail if solve not called}
  begin
    if ((ThisVar > 0) and (ThisVar <= Neq)) then begin
      GetAnswer := True;
      ThisVal := Answer^[ThisVar];
    end
    else begin
      GetAnswer := False;
      if (ThisVar < 1) then SetErrorMsg('VarNo < 1', 41, ThisVar, 0);
      if (ThisVar > Neq) then SetErrorMsg('VarNo > Neq', 42, ThisVar, Neq);
    end;                {else}
  end;                  {GetAnswer}

  procedure showmat;
  var
    Row, Col, c, LastCol : Word;
    ElmPtr         : PElmRec;
  begin
    for Row := 1 to Neq do begin
      ElmPtr := FirstElm^[Row];
      LastCol := 0;
      while ElmPtr <> nil do begin
        Col := ElmPtr^.Column;
        for c := (LastCol+1) to (Col-1) do Write('nil':6);
        Write(ElmPtr^.Value:6:2);
        LastCol := Col;
        ElmPtr := ElmPtr^.Next;
      end;
      for c := (LastCol+1) to (Neq+1) do Write('nil':6);
      WriteLn;
    end;                {for row}
  end;                  {showmat}


var
  SumTerm        : Extended;
  Factor         : Extended;
  RHS            : Extended;
  Biggest        : Extended;
  Coeff          : Extended;
  PivotValue     : Extended;
  BestPtr        : PElmRec;
  PrevPtr        : PElmRec;
  Next_Pivot     : PElmRec;
  Next_Tar       : PElmRec;
  NewPtr         : PElmRec;
  ElmPtr         : PElmRec;
  MinRowCount    : Word;
  Best_Addelm    : Word;
  NextTarCol     : Word;
  NumToFind      : Word;
  AddElm         : Word;
  Col, LastCol   : Word;
  SingleCount    : Word;
  PivotStep      : Word;
  PivotCol       : Word;
  NextPivotCol   : Word;
  LastRow        : Word;
  PrevRow, NextRow : Word;
  PivotRow       : Word;
  Row            : Word;
  NextActiveRow  : PIntArr;
  ColCount       : PIntArr;
  RowCount       : PWordArr;
  PivRow         : PWordArr;
  function Solve1 : Boolean;

  label Go_on, AsBigAs, NextOne;
  label TwoSingles, NumericallySingular, NoVars, EmptyRow, EmptyCol,
OutOfSpace, ColsOutofOrder, NoRHS;
  label InsertElm, AdjustElm, OutCase;
  begin                 {solve1}
    if Msg then WriteLn('Entering Solve1');
    PivotStep := 0;
    FreeMemParaAlign(Pointer(OrigLastElm), (1+Neq)*SizeOf(PElmRec));


    if not GetMemParaAlign(Pointer(NextActiveRow), Pointer(OrigNextActiveRow),
(1+Neq)*SizeOf(Integer)) then goto OutOfSpace;
    if not GetMemParaAlign(Pointer(ColCount), Pointer(OrigColCount),
(1+Neq)*SizeOf(Integer)) then goto OutOfSpace;
    if not GetMemParaAlign(Pointer(RowCount), Pointer(OrigRowCount),
(1+Neq)*SizeOf(Word)) then goto OutOfSpace;
    if not GetMemParaAlign(Pointer(PivRow), Pointer(OrigPivRow),
(1+Neq)*SizeOf(Word)) then goto OutOfSpace;


    Col := 0;           {set vectors to zero}
    FillStruct(RowCount^, 1+Neq, Col, SizeOf(Word));
    FillStruct(PivRow^, 1+Neq, Col, SizeOf(Word));
    FillStruct(ColCount^, 1+Neq, Col, SizeOf(Word));

    Solve1 := False;
    if Msg then WriteLn('about to set up row and column counts');
    for Row := 1 to Neq do begin
      ElmPtr := FirstElm^[Row];
      if (ElmPtr = nil) then goto EmptyRow;
      LastCol := 0;
      while ElmPtr <> nil do begin
        Col := ElmPtr^.Column;
        if DBug then ElmPtr^.CheckRow := Row;
        if (Col <= LastCol) then goto ColsOutofOrder
        else LastCol := Col;
        Inc(RowCount^[Row]);
        if DBug then Assert(Col > 0, '#2133');
        if DBug then Assert(Col <= (1+Neq), '#2134');
        if (Col <= Neq) then Inc(ColCount^[Col]);
        ElmPtr := ElmPtr^.Next;
      end;
      if (LastCol <> (1+Neq)) then goto NoRHS;
    end;                {for row}

    if Msg then WriteLn('about to complete setup');
    Row := 0;
    while (Row < Neq) do begin
      NextActiveRow^[Row] := Row+1;
      Inc(Row);
      if ColCount^[Row] = 0 then goto EmptyCol;
      if RowCount^[Row] = 1 then goto NoVars;
    end;                {for Row:=1 to neq}
    NextActiveRow^[Neq] := 0;
    {end setup}

    repeat              {pivot on variables which are mentioned only once}
      PivotCol := 0;
      PrevRow := 0;
      Row := NextActiveRow^[0];
      while Row <> 0 do begin
        NextRow := NextActiveRow^[Row];
        if DBug then Assert(Row > 0, '#8033');
        if DBug then Assert(Row <= Neq, '#9033');
        SingleCount := 0;
        ElmPtr := FirstElm^[Row]; if DBug then Assert(ElmPtr <> nil, '#34');
        if DBug then Assert(ElmPtr^.Column <= Neq, '#77');
        Col := ElmPtr^.Column;
        while (Col <= Neq) do begin
          if (ColCount^[Col] = 1) then begin
            PivotCol := Col;
            Inc(SingleCount);
          end;
          ElmPtr := ElmPtr^.Next; if DBug then Assert(ElmPtr <> nil, '#35');
          Col := ElmPtr^.Column;
        end;
        if (SingleCount > 1) then goto TwoSingles
        else if (SingleCount = 1) then begin
          Inc(PivotStep);
          PivotRow := Row;
          ElmPtr := FirstElm^[PivotRow]; if DBug then Assert(ElmPtr <> nil,
'#34');
          if DBug then Assert(ElmPtr^.Column <= Neq, '#77');
          Col := ElmPtr^.Column;
          while (Col <= Neq) do begin
            if DBug then Assert(ColCount^[Col] > 0, '#4177');
            Dec(ColCount^[Col]);
            ElmPtr := ElmPtr^.Next; if DBug then Assert(ElmPtr <> nil, '#35');
            Col := ElmPtr^.Column;
          end;
          PivRow^[PivotStep] := PivotRow;
          NextActiveRow^[PrevRow] := NextActiveRow^[PivotRow];
          NextActiveRow^[PivotRow] := -1; {useful ?}
          RowCount^[PivotRow] := PivotCol; {change of meaning}
          ColCount^[PivotCol] := -1; {mark as done}
        end             {if (SingleCount=1) }
        else {no Singles} PrevRow := Row;
        Row := NextRow;
      end;

    until (PivotCol = 0);


    {*************main loop}
    while PivotStep < Neq do begin
      Inc(PivotStep);
      if Msg then WriteLn('starting step ', PivotStep);

      { Find shortest row (PivotRow) and the preceding active row (LastRow)  }
      MinRowCount := MaxInt;
      PrevRow := 0; LastRow := 0;
      Row := NextActiveRow^[0];
      if DBug then Assert(Row <> 0, '#33');
      while Row <> 0 do begin
        if RowCount^[Row] < MinRowCount then begin
          MinRowCount := RowCount^[Row];
          PivotRow := Row;

          LastRow := PrevRow;
        end;
        PrevRow := Row;
        Row := NextActiveRow^[Row];
      end;

      if Msg then WriteLn('Pivotrow: ', PivotRow, ' Rowcount ', MinRowCount);

      Biggest := -1;
      ElmPtr := FirstElm^[PivotRow]; if DBug then Assert(ElmPtr <> nil, '#34');
      if DBug then Assert(ElmPtr^.Column <= Neq, '#77');
      while (ElmPtr^.Column <= Neq) do begin
        if (Abs(ElmPtr^.Value) > Biggest) then Biggest := Abs(ElmPtr^.Value);
        ElmPtr := ElmPtr^.Next; if DBug then Assert(ElmPtr <> nil, '#35');
      end;
      if DBug then Assert(Biggest >= 0, '#45');
      if Biggest = 0 then goto NumericallySingular;

      if Msg then WriteLn('Biggest was :', Biggest);

      Biggest := Biggest*Uvalue;
      BestPtr := nil;
      Best_Addelm := MaxInt;
      ElmPtr := FirstElm^[PivotRow]; if DBug then Assert(ElmPtr <> nil, '#36');
      while (ElmPtr^.Column <= Neq) do begin
        Col := ElmPtr^.Column;
        Dec(ColCount^[Col]);
        if (Abs(ElmPtr^.Value) >= Biggest) then begin
          AddElm := ColCount^[Col];
          {addelm is the number of additional nonzeros which would be added}
          {if this Pivot were chosen}
          if AddElm < Best_Addelm then begin
            BestPtr := ElmPtr;
            Best_Addelm := AddElm;
          end;
        end;            {if (....>=UValue)}
        ElmPtr := ElmPtr^.Next;
        if DBug then Assert(ElmPtr <> nil, '#37');
      end;
      if DBug then Assert(BestPtr <> nil, '#38');

      PivotCol := BestPtr^.Column;
      PivotValue := BestPtr^.Value;
      {Mark Pivot Row as inactive}
      NextActiveRow^[LastRow] := NextActiveRow^[PivotRow];

      {Answer Values for use by backsub}
      PivRow^[PivotStep] := PivotRow;
      {note change of meaning}
      NextActiveRow^[PivotRow] := -1; {useful ?}
      RowCount^[PivotRow] := PivotCol; {change of meaning}
      NumToFind := ColCount^[PivotCol];
      ColCount^[PivotCol] := -1; {mark as done}


      if Msg then Write('Start Pivot, ');
      if Msg then Write(' Pivotrow: ', PivotRow, ' Rowcount ', MinRowCount);
      if Msg then WriteLn(' PivotCol: ', PivotCol, ' ColCount ', NumToFind);

      Row := NextActiveRow^[0];
      while ((Row <> 0) and (NumToFind > 0)) do begin
        if (FreeCount < Neq) then if not FrePrime then goto OutOfSpace;

        {preload PrevPtr so that: PrevPtr^.Next :=  FirstElm^[Row]}
        PrevPtr := Addr(FirstElm^[Row-1]);
        ElmPtr := FirstElm^[Row];

        {the goto's in this section are intended to ease transition to assembly
language}

Go_on:  if (ElmPtr^.Column >= PivotCol) then goto AsBigAs;
        PrevPtr := ElmPtr;
        ElmPtr := ElmPtr^.Next; {not got to it yet}
        if DBug then Assert(ElmPtr <> nil, '#55');
        if (ElmPtr^.Column >= PivotCol) then goto AsBigAs;
        PrevPtr := ElmPtr;
        ElmPtr := ElmPtr^.Next; {not got to it yet}
        if DBug then Assert(ElmPtr <> nil, '#55');

        goto Go_on;
AsBigAs: if (ElmPtr^.Column <> PivotCol) then goto NextOne;

        {current row contains pivot col}
        if Msg then WriteLn('Altering Row ', Row);
        Factor := ElmPtr^.Value/PivotValue;
        Dec(NumToFind);


        { DELETE pivot col in current row}
        Dec(RowCount^[Row]);
        PrevPtr^.Next := ElmPtr^.Next;
        ElmPtr^.Next := FreePtr; FreePtr := ElmPtr; Inc(FreeCount);

        Next_Pivot := FirstElm^[PivotRow];
        if DBug then Assert(Next_Pivot <> nil, '#333');
        PrevPtr := Addr(FirstElm^[Row-1]); {PrevPtr^.Next :=  FirstElm^[Row]}
        Next_Tar := FirstElm^[Row];
        if DBug then Assert(Next_Tar <> nil, '#334');
        NextTarCol := Next_Tar^.Column;
        while Next_Pivot <> nil do begin
          NextPivotCol := Next_Pivot^.Column;
          if (NextPivotCol = PivotCol) then goto OutCase;

          while NextTarCol < NextPivotCol do begin
            PrevPtr := Next_Tar;
            Next_Tar := Next_Tar^.Next; if DBug then Assert(Next_Tar <> nil,
'#99');
            NextTarCol := Next_Tar^.Column
          end;
          if NextTarCol = NextPivotCol then goto AdjustElm;

InsertElm: {NextTarCol > NextPivotCol}
          {element in pivot row but not in current row: add in new element}
          if DBug then Assert(NextTarCol > NextPivotCol, '#69');
          if DBug then Assert(FreePtr <> nil, '#89');
          NewPtr := FreePtr;
          NewPtr^.Value := -Factor*Next_Pivot^.Value; {up for copro}
          FreePtr := FreePtr^.Next; Dec(FreeCount);
          PrevPtr^.Next := NewPtr;
          PrevPtr := NewPtr;
          NewPtr^.Column := NextPivotCol;
          if DBug then NewPtr^.CheckRow := Row;
          NewPtr^.Next := Next_Tar;
          Inc(ColCount^[NextPivotCol]);
          Inc(RowCount^[Row]);
          goto OutCase;

AdjustElm: if DBug then Assert(NextTarCol = NextPivotCol, '#67');
          Next_Tar^.Value := Next_Tar^.Value-Factor*Next_Pivot^.Value;
          PrevPtr := Next_Tar;
          Next_Tar := Next_Tar^.Next;
          if (Next_Tar <> nil) then
            NextTarCol := Next_Tar^.Column
          else NextTarCol := 2+Neq; {sentinel value}

OutCase:
          Next_Pivot := Next_Pivot^.Next; {move along pivot row}
        end;            {while Next_Pivot <>  Nil}

NextOne:
        Row := NextActiveRow^[Row];
      end;              {while row}
      if DBug then Assert(NumToFind = 0, '#66');
    end;                {main loop}

    {release un-needed vectors}
    FreeMemParaAlign(Pointer(OrigColCount), (1+Neq)*SizeOf(Integer));
    FreeMemParaAlign(Pointer(OrigNextActiveRow), (1+Neq)*SizeOf(Integer));


    {create Answer vector}
    if not GetMemParaAlign(Pointer(Answer), Pointer(OrigAnswer),
(1+Neq)*SizeOf(Double))
    then goto OutOfSpace;
    if DBug then for Row := 1 to Neq do Answer^[Row] := -99;


    PivotStep := 1+Neq;
    while (PivotStep > 1) do begin
      Dec(PivotStep);
      Row := PivRow^[PivotStep];
      PivotCol := RowCount^[Row]; {note change of meaning}
      SumTerm := 0.0;
      Coeff := 0.0;
      ElmPtr := FirstElm^[Row];
      if DBug then Assert(ElmPtr <> nil, '#188');
      while ElmPtr <> nil do begin
        Col := ElmPtr^.Column;
        if (Col = PivotCol) then
          Coeff := ElmPtr^.Value
        else if (Col <= Neq) then begin
          if DBug then Assert(Answer^[Col] <> -99, '#177');
          SumTerm := SumTerm+Answer^[Col]*ElmPtr^.Value;
        end
        else
          RHS := ElmPtr^.Value;

        ElmPtr := ElmPtr^.Next;
      end;              {until (elmptr=Nil)}
      if DBug then Assert(Answer^[PivotCol] = -99, '#77');
      Answer^[PivotCol] := (RHS-SumTerm)/Coeff;
    end;                {for PivotRow:=neq downto 1}


    FreeMemParaAlign(Pointer(OrigRowCount), (1+Neq)*SizeOf(Word));
    FreeMemParaAlign(Pointer(OrigPivRow), (1+Neq)*SizeOf(Word));


    Solve1 := True; Exit; {normal}


EmptyRow: SetErrorMsg('Empty Row', 1, Row, 0); Exit;
NoVars: SetErrorMsg('Row without Variables', 2, Row, 0); Exit;
EmptyCol: SetErrorMsg('Empty Col', 3, Col, 0); Exit;
TwoSingles: SetErrorMsg('Two Singles', 4, Row, PivotCol); Exit;
NoRHS: SetErrorMsg('No RHS', 5, Row, LastCol); Exit;
NumericallySingular: SetErrorMsg('Numerically Singular', 6, PivotRow,
PivotStep); Exit;
OutOfSpace: SetErrorMsg('Out of Space', 7, PivotStep, 4); Exit;
ColsOutofOrder: SetErrorMsg('Cols out of Order', 8, Row, Col); Exit;

  end;                  {Solve1}

end.
{************************************************************************}
{$N+,E+}
uses SolSpar;
var
  Reason         : String;
  ErrNo1, ErrNo2, ErrNo3 : Word;
  Col, Row, N    : Integer;
  Total, Value   : Double;
label Fail;
begin
  WriteLn('Initial Heap Size ', MemAvail);
  N := 1500;            {size of matrix}
  if not InitStruc(N) then goto Fail;

{constructs a 'band' matrix with column border and solves it}

{for testing purposes, the RHS is set up so that variable(n)
has value n. The 'Total' variable is part of this. In normal use
the RHS would be set to some other value}

  for Row := 1 to N do begin
    Total := 0.0;

    Value := 0.1+Random;
    Col := Row;
    if not AddElement(Row, Col, Value) then goto Fail;
    Total := Total+Col*Value;

    Col := Row-1;
    if (Col > 0) then begin
      Value := 0.1+Random;
      if not AddElement(Row, Col, Value) then goto Fail;
      Total := Total+Col*Value;
    end;

    Col := Row+1;
    if (Col <= N) then begin
      Value := 0.1+Random;
      if not AddElement(Row, Col, Value) then goto Fail;
      Total := Total+Col*Value;
    end;


    Col := N;
    Value := 0.1+Random;
    if not AddElement(Row, Col, Value) then goto Fail;
    Total := Total+Col*Value;

    if not SetRHS(Row, Total) then goto Fail;
  end;                  {for Row:= 1 TO N}
  WriteLn(' Memory Used To Store Matrix ', SparMemUsed);

  if not Solve1 then goto Fail;
  for Row := 1 to N do if ((Row mod 100) = 0) then begin
    if not GetAnswer(Row, Value) then goto Fail;
    WriteLn(Row:3, Value:15:5);
  end;
  WriteLn(' Max Memory Used To Solve Matrix ', SparMemUsed);
  ReleaseStruc;
  WriteLn(' Memory Used after ReleaseStruc ', SparMemUsed);
  WriteLn('  Final Heap Size ', MemAvail);

  Halt;

Fail:
  GetErrorMsg(Reason, ErrNo1, ErrNo2, ErrNo3);
  WriteLn('Failed:  Error ', ErrNo1:0, ' ', Reason, ' ', ' ', ErrNo2:3, ' ',
ErrNo3:3);
  ReleaseStruc;
  WriteLn(' Memory Used after ReleaseStruc ', SparMemUsed);
  WriteLn('  Final Heap Size ', MemAvail);
end.
{************************************************************************}
{$N+,E+}
uses SolSpar;
var
  Reason         : String;
  ErrNo1, ErrNo2, ErrNo3 : Word;
  Example, Col, Row, N : Integer;
  Total, Value   : Double;
label Fail, EndProg;
begin
  N := 4;               {size of matrix}

{Shows how error messages work. Compile and run
   test2 >tem
 off DOS command line.  Then compare file 'tem' with code below}

{for testing purposes, the RHS is set up so that variable(n)
has value n. The 'Total' variable is part of this. In normal use
the RHS would be set to some other value}

  for Example := 1 to 8 do begin
    WriteLn('Example ', Example);
    if not InitStruc(N) then goto Fail;

    if not AddElement(1, 1, 0.50) then goto Fail;
    if not AddElement(1, 3, 2.00) then goto Fail;
    if not SetRHS(1, 6.50) then goto Fail;

    if (Example <> 7) then
      if not AddElement(2, 2, 0.20) then goto Fail;
    if not AddElement(2, 4, 5.00) then goto Fail;
    if not SetRHS(2, 20.40) then goto Fail;

    if not AddElement(3, 3, 1.00) then goto Fail;
    if not AddElement(3, 1, 0.75) then goto Fail;
    if not SetRHS(3, 3.75) then goto Fail;

    if (Example = 1) or (Example = 7) then begin
      if not AddElement(4, 4, 6.00) then goto Fail;
      if not AddElement(4, 1, 2.00) then goto Fail;
      if not SetRHS(4, 26.00) then goto Fail;
    end
    else if (Example = 2) then begin {4th row is multiple of second}
      if not AddElement(4, 2, 0.40) then goto Fail;
      if not AddElement(4, 4, 10.0) then goto Fail;
      if not SetRHS(4, 26.00) then goto Fail;
    end
    else if (Example = 3) then begin {variables 2 and 4 each appear only in row
2}
      if not AddElement(4, 3, 1.00) then goto Fail;
      if not AddElement(4, 1, 0.75) then goto Fail;
      if not SetRHS(4, 26.00) then goto Fail;
    end
    else if (Example = 4) then begin {no vars in row 4}
      if not SetRHS(4, 26.00) then goto Fail;
    end
    else if (Example = 5) then begin    {no rhs}
      if not AddElement(4, 3, 1.00) then goto Fail;
      if not AddElement(4, 1, 0.75) then goto Fail;
    end
    else if (Example = 6) then begin   {no row 4}
    end;
    showmat;
    if not Solve1 then begin
      GetErrorMsg(Reason, ErrNo1, ErrNo2, ErrNo3);
      WriteLn('Failed:  Error ', ErrNo1:0, ' ', Reason, ' ', ' ', ErrNo2:3, '
', ErrNo3:3);
    end
    else begin
      Write('Solved: ');
      for Row := 1 to N do if GetAnswer(Row, Value)
        then Write('  X(', Row, ')=', Value:0:3)
        else goto Fail;
    end;
    WriteLn;

    if (Example < 7) then ReleaseStruc;
  end;                  {for example}

  Halt;

Fail:
  GetErrorMsg(Reason, ErrNo1, ErrNo2, ErrNo3);
  WriteLn('Failed:  Error ', ErrNo1:0, ' ', Reason, ' ', ' ', ErrNo2:3, ' ',
ErrNo3:3);
end.