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.