Contributor: MARK HORRIDGE
{
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 0Nil, 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.