[Back to MATH SWAG index] [Back to Main SWAG index] [Original]

    {
    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.

[Back to MATH SWAG index] [Back to Main SWAG index] [Original]