sparsolv.~pas

deihli写的稀疏矩阵链表存储

{ DEFINE DBUG}
{$IFDEF DBUG}
{$R+,S+,Q+,X+}
{$ELSE}
{$R-,S+,Q-,X+}
{$ENDIF}
{$X+}

{you may have to adjust these defines .........}
{$DEFINE DELPHI}   {if using Delphi}
{$DEFINE DELPHI32} {if using Delphi 2 or later or any 32-bit pascal}
{$DEFINE BIT32}    {if using any 32-bit pascal}
{x $DEFINE BIT16}   {if using 16-bit pascal; just ONE of BIT32 or BIT16 must be defined}


{ $DEFINE ALIGNEDBLOCKS} {probably only useful for 16bit pascal}


unit SparSolv;
interface
   { 包含下面5个逻辑函数和3个过程，如果运算成功地执行每个函数返回到真否则是假。如果返回的是假，调
用出错信息对话框来寻找原因。最后调用过程释放结构来给还l内存空间给堆栈。    }
 {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: Integer): Boolean;
{ 初始化 ，        产生一个初始化的稀疏矩阵程序-第一步。
   那个NumEq是方程或者变量的个数}
 {creates and initializes sparse matrix structure - call this first.
  NumEq is number of equations/variables.}

function AddLHS(ThisEqu, ThisVar: Integer; ThisVal: Double): Boolean;
//叠加元素到 LHS = Left Hand Side 左手边
{给方程 ThisEqu和变量 ThisVar 增加一个元素到稀疏矩阵
如果那样的实体真的存在， ThisVal 将会被加到存在的值}
 {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}

function AddRHS(ThisEqu: Integer; ThisVal: Double): Boolean;
//叠加元素到  RHS = Right Hand Side 右手边
{给方程ThisEqu 设置右手边的数据，如果已经存在，就加到上面 }
 {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: Integer; 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: Integer);
//异常
 {N1: error number; S: Error Description; N2, N3 : other possibly useful
  numbers}

procedure ShowMat;
//显示小规模稀疏矩阵
{ displays small sparse matrix }

procedure WriteMatrixAsMTX(const filename: string);
//以MTX格式写矩阵
 {writes matrix in MTX format}

var
  SparMemUsed: LongInt; {no of bytes of heap currently used by routines}  {}
  MaxMemUsed: LongInt; {highest value reached by SparMemUsed}    {被 SparMemUsed用到的最大的数值}
  FillIns: LongInt; {no of elements added during solve}     {在求解的过程中被加的元素的数量}
const
{$IFDEF BIT16}
  MaxMemToUse: LongInt = 4 * 16 * High(Word); {upper limit to heap use: default 4MB} {堆的上界：默认的是4MB}
{$ELSE}
  MaxMemToUse: LongInt = High(LongInt); {upper limit to heap use: default no limit}   {堆的上界：默认的是无穷}
{$ENDIF}

implementation // see also line 166

const
  Scaling = True;
  //是否已经预处理
  {if true, matrix is preconditioned: see below}
  Msg = False;
{$IFDEF ALIGNEDBLOCKS}
  ParaAlign = 16; {suggest for paragraph alignment - try 1 to see slowdown} {堆栈段的建议值-试一下1让其停下}
{$ENDIF}
{$IFDEF BIT32}
  MaxEq = 500000; {too big doesnt matter}  {太大也没问题}
{$ELSE} {16-bit, so 64k structure limit applies}
  MaxSize = 65520; {largest variable size allowed by Turbo Pascal}   {Turbo Pascal最大的允许变量规模}
  MaxUsable = MaxSize - ParaAlign; {allows for paragraph alignment of data}  {堆栈段的允许值}
  MaxEq = (MaxUsable div SizeOf(Double)) - 1; {allows 1 extra for 0-based arrays}   {}
 {based on above, MaxEq = 8187} {在前面的基础上最大的值是8178}
{$ENDIF}
  UValue = 0.1; {number 0<U<1; if larger, time and memory usage are
                          increased at the gain, maybe, of accuracy}      {0<u<1,如果太大，同样计算准确度的时间和内存将会增加}

type
  RealArr = array[0..MaxEq] of Double;
  IntArr = array[0..MaxEq] of Integer;
  PRealArr = ^RealArr; //双精度数组 RealArr的指针？
  PIntArr = ^IntArr; //
    {ElmRec 记录类型存储了一个单链条的在稀疏矩阵中的实体。基于效率的考虑它是16字节长，同时可以被段对齐。
    一个填充字段构成的16个字节，“载入”“”技巧要求下面的字段在记录的开始。因此，排列好字段尤为重要。
    【注：在DELPHI2或者更前的版本中整数要求4个bit，其他的2个bit就够了】，
    通过转换 数值字段为单精度，定位字段为字符同时消除填充字段，记录大小可以减小到10bit.
    但是精确度就会极大的降低。用到的节点按一系列的连接的列表排列，每一个代表一个方程，或者一行。
    在每一个清单变量（栏）总是出现在升序排列。每一个列表都被RHS或者常实体终止；这个可视为相当于额外变量编号（数字+1） }
  { The ElmRec record type stores a single entry in the sparse matrix.
  For efficiency reasons, it is 16 bytes long, and should be paragraph-
  aligned.  A Padding field makes up the 16 bytes.  The 'preload
  PrevPtr' trick, see below, requires that Next field be at the start
  of the record.  Hence, order of fields is important.  [Note:
  Integers are 4 bytes in Delphi 2 and above, otherwise 2 bytes.]
    By converting the Value field to a Single, the Column Field to a
  Word, and eliminating the Padding field, the record size could be
  reduced to 10 bytes.  However, accuracy would be greatly reduced.
    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记录的指针

  ElmRec = record   {定义记录}
    Next: PElmRec; {pointer to next node in row}// 行上的下一个指针
    Column: Integer; {variable no. of this node} //这个节点的定位数字 为什么只要列的定位呢？
{$IFDEF BIT16}
    Padding: Integer; {only needed if Integer is 2 bytes long}  {只有在整数是2bit时才需要}
{$ENDIF}
    Value: Double; {coefficient value}  {系数值}
  end;

  PtrArr = array[0..MaxEq] of PElmRec; // 指针数组，他的每个元素都是一个指针
  PPtrArr = ^PtrArr;// 指针数组记录集PtrArr的指针
    {被指针FreePtr 指向的不是立即要用的备用节点，被连在一个列表。FreeCount 是这些节点的数目。
    因为有许多同样大小的节点没有被系统堆栈管理器直接分配地址。而是一个更加有效率计划被用了。
    当自由节点需要扩展足够的空间给许多节点时，需要一大块的内存。
    为了后面系统堆栈的处理，另外一个列表被用来存储这些块的地址。 }
  { 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);//
{$IFDEF ALIGNEDBLOCKS}    {see GetMemParaAlign below}
  ElmsPerBlock = (65520 - ParaAlign) div ElmSize;//
{$ELSE}
  ElmsPerBlock = 65520 div ElmSize; //
{$ENDIF}
{这个集块的大小都小于64K-最大的16-bit Pascal允许大小。因此每个块拥有16000个节点。这个对32-bit Pascal同样适用}
  {This sets block size at just under 64k - the largest size allowed
   by 16-bit Pascal. Hence each block holds around 16000 nodes. This
   should be fine for 32-bit Pascal too}

type
  TNodeBlock = array[1..ElmsPerBlock] of ElmRec; //这个类定义干什么的？ 分块存储？
  PNodeBlock = ^TNodeBlock;

  PHeapRec = ^HeapRec;
  HeapRec = record
    BlockPtr: PNodeBlock; {original address of block块的原始地址}
    NewBlockPtr: PNodeBlock; {adjusted (aligned) address of bloc块的调整的地址}
    BlockSiz: LongInt; {size of block块的大小}
    NextRec: PHeapRec; {address of next list item下一个列表对象的地址}
  end;

var {this list contains variables which must have unit-wide scope这个列表包含了单元范围作用的变量}
  Reason: string; {for error messages为了错误信息}
  ErrNo1, ErrNo2, ErrNo3: Integer; {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: Integer; {number of free nodes自由节点的数量}
  Neq: Integer; {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每个方程的最后依个节点的指针数组}

{$IFDEF DELPHI32}{$J+}{$ENDIF} {allow writeable constants}

const
  Initialized: Boolean = False; {marks whether unit is in use 记录单元是否在用}
  Solved: Boolean = False; {marks whether solution was successful记录方程否求解顺利}

//implementation  // with "implementation" here, all is exposed !

procedure Assert(P: Boolean; S: string);
  {used for debugging用来调试}
begin
  if P then Exit;
  WriteLn('assertion failed: ', S);
  WriteLn('hit enter to finish');
  ReadLn;
  Halt(1);
end; {Assert}

procedure RecordAllocation(const OrigP, NewP: Pointer; const Siz: LongInt);
//预定分配堆
  {prepend a new heap record to the list}
var NewRec: PHeapRec;
begin
  New(NewRec);//分配内存，初始化
  NewRec^.NextRec := BlockList;
  NewRec^.BlockSiz := Siz;
  NewRec^.BlockPtr := OrigP;
  NewRec^.NewBlockPtr := NewP;
  BlockList := NewRec;
  Inc(SparMemUsed, Siz + SizeOf(HeapRec)); {tally of memory used计算所用的内存}
  if (SparMemUsed > MaxMemUsed) then MaxMemUsed := SparMemUsed;
end;

{$IFDEF ALIGNEDBLOCKS}
      {如果所有的节点记录开始于段边界，程序会快一些。遗憾的是delphi16-bit堆管理不是总是发行按段排列的地址。
      GetMemParaAlign程序纠正这个问题。

      32-bit Delphi 和免费的Pascal堆管理器会发行段排列的地址。因此 GetMemParaAlign被更简单的版本代替。。如下}

      
{routine will be faster if all node records start on a paragraph boundary.
 Unfortunately the heap manager in Borland 16-bit Pascals did not
 always issue paragraph-aligned addresses. The GetMemParaAlign function
 corrects this problem.

 It appears that the heap manager in 32-bit Delphi and Free Pascal does
 issue paragraph-aligned addresses. So GetMemParaAlign is replaced
 by a simpler version...see below. }

function GetMemParaAlign(var P: Pointer; Size: LongInt): Boolean;
{返回一个段排列的指针点（多重的段排列），这个指针同时也是指向一个至少几个字节的新的记忆块。
486高速缓冲存储器对段排列的数据很适用。更详细的说GetMemParaAlign从系统栈中获得一块内存，系统栈
是比尺寸大的段排列字节。OrigP指向这些原始块。P是在OrigP的第一个地址，OrigP是很多的段排列。
为了后面的释放OrigP和分派的大小尺寸被保存在连在一起的列表}
  { Returns a pointer P which is paragraph aligned (a multiple of ParaAlign)
  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 ParaAlign bytes larger than Size.  OrigP points to this
  original block.  P is the first address after OrigP which is a multiple of
  ParaAlign.  OrigP and the allocated size are saved in a linked list for
  later release.}
var
  OrigP: Pointer;
  AllocatedSize: LongInt;

{$IFDEF BIT16}
type {used only for typecasting给字符常量用}
  PtrRec = record OfsWord, SegWord: Word; end;
var
  Ofset: Word;
{$ENDIF}
begin
{$IFDEF TRACE}WriteLn('Entering GetMemParaAlign'); {$ENDIF}
  P := nil;
  GetMemParaAlign := False;
  AllocatedSize := Size + ParaAlign;
{$IFNDEF DELPHI32}
  if (AllocatedSize > MaxAvail) then Exit;
{$ENDIF}
  if ((MaxMemToUse <> -1) and
    ((SparMemUsed + AllocatedSize) > MaxMemToUse)) then Exit;
  GetMem(OrigP, AllocatedSize);
  if (OrigP = nil) then Exit;
{$IFDEF BIT32}
  P := Pointer(ParaAlign + ParaAlign * (LongInt(OrigP) div ParaAlign));
{$ELSE}
  P := OrigP; {to load segment}
  Ofset := PtrRec(P).OfsWord;
  {adjust offset to paragraph boundary}
  PtrRec(P).OfsWord := ParaAlign + ParaAlign * (Ofset div ParaAlign);
{$ENDIF}
  RecordAllocation(OrigP, P, AllocatedSize);
  GetMemParaAlign := True;
end; {GetMemParaAlign}

{$ELSE}

function GetMemParaAlign(var P: Pointer; Size: LongInt): Boolean;
  { replacement GetMemParaAlign function if original gives problems.
    Returns a pointer P which points to a new block of memory.
  be used.  P and the allocated size are saved in a linked list for
  later release.
  如果原始数据出现问题就替代GetMemParaAlign，返回一个指向新内存块的P指针，P和
  为了释放的方便指针P和相应分配的尺寸保存在链表里}
var
  OrigP: Pointer;
begin
  P := nil;
  GetMemParaAlign := False;
  if ((MaxMemToUse <> -1) and
    ((SparMemUsed + Size) > MaxMemToUse)) then Exit;
  GetMem(OrigP, Size);
  if (OrigP = nil) then Exit;
  P := OrigP;
  RecordAllocation(OrigP, P, Size);
  GetMemParaAlign := True;
end; {GetMemParaAlign}

{$ENDIF}

function TopUpFreeList: Boolean;
type {used only for typecasting}
  PtrRec = record OfsWord, SegWord: Word; end;
var
  NewBlock: PNodeBlock;
{$IFDEF BIT16}
  Ofset: Word;
{$ENDIF}
  Count: Integer;
  P: PElmRec;
begin
{$IFDEF TRACE}WriteLn('Entering TopUpFreeList'); {$ENDIF}
  TopUpFreeList := False;
  {get new block}
  if not GetMemParaAlign(Pointer(NewBlock), SizeOf(TNodeBlock)) then Exit;
  {fill new block with linked nodes}
  P := PElmRec(NewBlock);
{$IFDEF BIT32}
  for Count := 1 to ElmsPerBlock do begin
    Inc(P); {increments P by size of ElmRec}