sparsolv.~pas

deihli写的稀疏矩阵链表存储

    NewBlock^[Count].Next := P;
  end;
{$ELSE}
  Ofset := PtrRec(P).OfsWord;
  for Count := 1 to ElmsPerBlock do begin
    Inc(Ofset, ElmSize);
    PtrRec(P).OfsWord := Ofset;
    NewBlock^[Count].Next := P;
  end;
{$ENDIF}

  {attach new block to free list增加新的一块到自由列表}
  NewBlock^[ElmsPerBlock].Next := FreePtr; {point end of new block to FreePtr}
  FreePtr := PElmRec(NewBlock); {point FreePtr at start of new block}
  Inc(FreeCount, ElmsPerBlock);

  TopUpFreeList := True;
end; {TopUpFreeList}

procedure SetErrorMsg(S: string; N1, N2, N3: Integer);   {这个过程干嘛的？}
begin
  Reason := S; ErrNo1 := N1; ErrNo2 := N2; ErrNo3 := N3;
end; {SetErrorMsg}

procedure GetErrorMsg(var S: string; var N1, N2, N3: Integer);   {这个过程干嘛的？}
begin
  S := Reason; N1 := ErrNo1; N2 := ErrNo2; N3 := ErrNo3;
end; {GetErrorMsg}

procedure NotInit;    {这个过程干嘛的？}
begin
  SetErrorMsg('InitStruc must be called', 1, 0, 0);
end;

function InitStruc(NumEq: Integer): Boolean; {这个函数干嘛的？初始化条件，看看满足不满足求解设定的条件}
begin
{$IFDEF TRACE}WriteLn('Entering InitStruc'); {$ENDIF}
  InitStruc := False;
  Neq := NumEq;
  Solved := False;
  SetErrorMsg('', 0, 0, 0);
  BlockList := nil;
  SparMemUsed := 0;
  FillIns := 0;
  MaxMemUsed := 0;
  if Initialized then begin
    SetErrorMsg('Initialize without releasing ', 2, 0, 0);
    Exit;
  end else Initialized := True;
  if (Neq > MaxEq) then begin
    SetErrorMsg('Too many equations ', 3, Neq, 0);
    Exit;
  end;
  if (Neq < 1) then begin
    SetErrorMsg('Too few equations ', 4, Neq, 0);
    Exit;
  end;

  if not GetMemParaAlign(Pointer(FirstElm), (1 + Neq) * SizeOf(PElmRec)) then begin
    SetErrorMsg('Out of Space', 5, 0, 0);
    Exit;
  end else FillChar(FirstElm^, (1 + Neq) * SizeOf(FreePtr), 0);

  if not GetMemParaAlign(Pointer(LastElm), (1 + Neq) * SizeOf(PElmRec)) then begin
    SetErrorMsg('Out of Space', 6, 0, 0);
    Exit;
  end else FillChar(LastElm^, (1 + Neq) * SizeOf(FreePtr), 0);

  FreePtr := nil;
  FreeCount := 0;
  InitStruc := True;
end; {InitStruc}

procedure ReleaseItem(var P: Pointer);
  {release one item from user heap从用户堆里释放一个对象}
  {note no error if P is nil or is not on user heap如果P是空的或者不再堆栈里就报告一个错误}
var NextPtr: PHeapRec;
begin
{$IFDEF TRACE}WriteLn('Entering ReleaseItem'); {$ENDIF}
  if (P = nil) then Exit;
  NextPtr := BlockList;
  while (NextPtr <> nil) do with NextPtr^ do begin
      if (NewBlockPtr = P) then begin
        FreeMem(BlockPtr, BlockSiz);
        Dec(SparMemUsed, BlockSiz);
        BlockPtr := nil;
        NewBlockPtr := nil;
        P := nil;
{$IFDEF TRACE}WriteLn('released an item'); {$ENDIF}
        Exit;
      end;
      NextPtr := NextPtr^.NextRec;
    end; {while}
{$IFDEF DBUG}Assert(False, '#3941'); {$ENDIF}
end; {ReleaseItem}

procedure ReleaseStruc;  {释放一个结构？}
var NextPtr: PHeapRec;
begin
{$IFDEF TRACE}WriteLn('Entering ReleaseStruc'); {$ENDIF}
  {get rid of user heap}
  while (BlockList <> nil) do with BlockList^ do begin
      if (BlockPtr <> nil) then begin
        FreeMem(BlockPtr, BlockSiz);
        Dec(SparMemUsed, BlockSiz);
{$IFDEF TRACE}WriteLn('releasing a block'); {$ENDIF}
      end;
      NextPtr := BlockList^.NextRec;
      Dispose(BlockList); Dec(SparMemUsed, SizeOf(HeapRec));
      BlockList := NextPtr;
    end; {while}
  Initialized := False;
end; {ReleaseStruc}

function AddElement(ThisEqu, ThisVar: Integer; ThisVal: Double): Boolean;//
  {AddElement will run quicker if: for each row, you first set LHS elements in ascending
  order, THEN set the RHS 增加单元会运行得更快，如果每一行你首先按升序设置LHS元素，然后设置RHS元素}
  //用几个链表来存储呢   两个左边矩阵，右边值向量，两个都用一个元素增加函数，怎么区别？
  //值是怎么存的  通过单链表，在单链表中增加元素的做法
  // 存储后怎么计算呢？
var
  PrevPtr, ElmPtr, NewPtr: PElmRec;   // PrevPtr,以前的链表， ElmPtr元素链表，NewPtr新的链表
begin {PivRow[Row] points to last element  PivRow[Row]指向最后一个节点 }
  AddElement := False;
  if (ThisEqu < 1) then begin SetErrorMsg('Row < 1', 7, ThisEqu, ThisVar); Exit; end;
  if (ThisEqu > Neq) then begin SetErrorMsg('Row > Neq', 8, ThisEqu, Neq); Exit; end;
  if (FreeCount < Neq) then begin
    if not TopUpFreeList then begin
      SetErrorMsg('Out of Space', 9, 0, 0); Exit;
    end;
  end;

  {get node from free list and set its values从自由列表得到节点并设置值，数据先放到新节点然后将节点插入到链表中}
  NewPtr := FreePtr; FreePtr := FreePtr^.Next; Dec(FreeCount);
  NewPtr^.Value := ThisVal;
  NewPtr^.Column := ThisVar;
  NewPtr^.Next := nil;    //这好像是用链表节点记录数据

  {insert node in proper place in linked list在链表中合适的地方增加节点}
  if FirstElm^[ThisEqu] = nil then begin {this is first entry这是第一个实体}
    FirstElm^[ThisEqu] := NewPtr;
    LastElm^[ThisEqu] := NewPtr;
       end
       else if (LastElm^[ThisEqu]^.Column < ThisVar) then begin //这是在新的列大于前面一个元素的列
      {attach at end of row加在行链表的末端}
       LastElm^[ThisEqu]^.Next := NewPtr;
       LastElm^[ThisEqu] := NewPtr;
           end
          else begin {traverse the row till we find the right place遍历行知道我们找到正确的位置}
        ElmPtr := FirstElm^[ThisEqu];  //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;
          {return node to free list}
           NewPtr^.Next := FreePtr; FreePtr := NewPtr; Inc(FreeCount);
                end
         else begin {ElmPtr^.Column > ThisVar}
        {insert new node before elmptr}
            if PrevPtr = nil then FirstElm^[ThisEqu] := NewPtr
            else PrevPtr^.Next := NewPtr;
         NewPtr^.Next := ElmPtr;
           end;
     end;
  AddElement := True;
end; {AddElement}

function AddLHS(ThisEqu, ThisVar: Integer; ThisVal: Double): Boolean;   {增加LHS左边矩阵值就存在这里了！}
//判断是否是0，然后决定存储与否，返回判断
//ThisEqu是行 ThisVar是列
begin
  if (ThisVal = 0.0) then begin AddLHS := True; Exit; end;
  AddLHS := False;
  if not Initialized then begin NotInit; Exit; end;
  if (ThisVar < 1) then begin SetErrorMsg('Col < 1', 10, ThisEqu, ThisVar); Exit; end;
  if (ThisVar > Neq) then begin SetErrorMsg('Col > Neq', 11, ThisEqu, ThisVar); Exit; end;
  AddLHS := AddElement(ThisEqu, ThisVar, ThisVal)//都给链表增加一个值但是RHS？？;
end; {AddLHS}

function AddRHS(ThisEqu: Integer; ThisVal: Double): Boolean;
begin
  AddRHS := False;
  if not Initialized then begin NotInit; Exit; end;
  AddRHS := AddElement(ThisEqu, 1 + Neq, ThisVal);
end; {AddRHS}

function GetAnswer(ThisVar: Integer; var ThisVal: Double): Boolean;
  {should fail if solve not called 如果solve命令没有执行，这个函数就会失败}
begin
  GetAnswer := False;
  if not Initialized then begin NotInit; Exit; end;
  if not Solved then begin SetErrorMsg('System not solved', 12, 0, 0); Exit; end;
  if (ThisVar < 1) then begin SetErrorMsg('VarNo < 1', 13, ThisVar, 0); Exit; end;
  if (ThisVar > Neq) then begin SetErrorMsg('VarNo >Neq', 14, ThisVar, Neq); Exit; end;
  ThisVal := Answer^[ThisVar];
  GetAnswer := True;
end; {GetAnswer}

procedure ShowMat;//  展示一下较小的矩阵
var
  Row, Col, c, LastCol: Integer;
  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}

procedure WriteMatrixAsMTX(const filename: string);
var
  Row, Col, Count: Integer;
  ElmPtr: PElmRec;
  Outfile: text;
begin
  Count := 0;
  for Row := 1 to Neq do begin {count elements}
    ElmPtr := FirstElm^[Row];
    while ElmPtr <> nil do begin
      Col := ElmPtr^.Column;
      if (Col <= Neq) then Inc(Count);
      ElmPtr := ElmPtr^.Next;
    end;
  end; {for row}

  Assign(outfile, filename); Rewrite(Outfile);
  Writeln(Outfile, '%%MatrixMarket matrix coordinate real general');
  Writeln(Outfile, Neq, ' ', Neq, ' ', Count); //
  for Row := 1 to Neq do begin
    ElmPtr := FirstElm^[Row];
    while ElmPtr <> nil do begin
      Col := ElmPtr^.Column;
      if (Col <= Neq) then Writeln(Outfile, Row, ' ', Col, ' ', ElmPtr^.Value: 0);
      ElmPtr := ElmPtr^.Next;
    end;
  end; {for row}
  Close(OutFile);
end; {showmat}

{All these variables are used only by Solve1, but have been declared
here so that their address is known at compile time.  This can add to
speed 所有这些变量都在Solve1程序中被用到，但是在这里定义了，
这样在比编译的时候做他们的地址就可以知道，这样可以加快速度}
var
  PrevPtr: PElmRec; //
  ElmPtr: PElmRec;  //当前记录的指针
  SumTerm: Extended;
  Factor: Extended;
  RHS: Extended;
  Biggest: Extended;
  Coeff: Extended;
  PivotValue: Extended;
  BestPtr: PElmRec;
  Next_Pivot: PElmRec;
  Next_Tar: PElmRec;
  NewPtr: PElmRec;
  NextActiveRow: PIntArr;
  ColCount: PIntArr; //不能出现为空的列啊
  RowCount: PIntArr; {for active rows=no of LHS elements; 活动的行就等于是  LHS元素的序号
                      for solved rows=the variable solved for解出的行就等于 求解的变量}
  ColScale: PRealArr;//这是用来干嘛的？
  PivRow: PIntArr; {used to remember in which order pivot rows were chosen用来记录 选择的支点行 在哪里}
  MinRowCount: Integer;
  Best_Addelm: Integer;
  NextTarCol: Integer;
  NumToFind: Integer;
  AddElm: Integer;
  LastCol: Integer;
  SingleCount: Integer;
  PivotStep: Integer;
  PivotCol: Integer;
  NextPivotCol: Integer;
  LastRow: Integer;
  PrevRow, NextRow: Integer;
  PivotRow: Integer;

function Solve1: Boolean;
var Row, Col: Integer; {to avoid delphi 2 warnings}
label Fail, OutOfSpace;
begin {solve1}
{$IFDEF TRACE}WriteLn('Entering Solve1'); {$ENDIF}
  Solve1 := False;
  PivotStep := 0;
  if not Initialized then begin NotInit; goto Fail; end;
  if Solved then begin SetErrorMsg('System already solved', 15, 0, 0); goto Fail; end;

  ReleaseItem(Pointer(LastElm));

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

  {set vectors to zero将向量设置为0}
  FillChar(RowCount^, (1 + Neq) * SizeOf(Integer), 0);
  FillChar(PivRow^, (1 + Neq) * SizeOf(Integer), 0);
  FillChar(ColCount^, (1 + Neq) * SizeOf(Integer), 0);

{$IFDEF TRACE}WriteLn('about to set up row and column counts现在成立行和列计数'); {$ENDIF}
  for Row := 1 to Neq do begin
    ElmPtr := FirstElm^[Row];
    if (ElmPtr = nil) then begin
      SetErrorMsg('Empty Row', 16, Row, 0); goto Fail;
    end;
    LastCol := 0;
    while ElmPtr <> nil do begin
      Col := ElmPtr^.Column;
      if (Col <= LastCol) then begin
        SetErrorMsg('Cols out of Order', 17, Row, Col); goto Fail;
      end
      else LastCol := Col;
      Inc(RowCount^[Row]); //取整加一
{$IFDEF DBUG}Assert(Col > 0, '#2133'); {$ENDIF}
{$IFDEF DBUG}Assert(Col <= (1 + Neq), '#2134'); {$ENDIF}
      if (Col <= Neq) then Inc(ColCount^[Col]);
      ElmPtr := ElmPtr^.Next;//将当前记录的指针移到下面一个记录
    end;
    if (LastCol <> (1 + Neq)) then begin
      SetErrorMsg('No RHS', 18, LastCol, 0); goto Fail;
    end;