sparsolv.~pas

deihli写的稀疏矩阵链表存储

  end; {for row}

  for Col := 1 to Neq do if (ColCount^[Col] = 0) then begin//如果ColCount^[Col等于0
      SetErrorMsg('Empty Col', 19, Col, 0); goto Fail;
    end;

  for Row := 1 to Neq do begin
    NextActiveRow^[Row] := Row + 1;
    if RowCount^[Row] <= 1 then begin
    {get this if you have RHS but no LHS如果你只有RHS而没有LHS你将会得到这些错误信息}
      SetErrorMsg('Row without Variables', 20, Row, 0); goto Fail;
    end;
  end; {for Row:=1 to neq}
  {end setup设置完毕}
{$IFDEF TRACE}WriteLn('completed setup'); {$ENDIF}

  if Scaling then begin
{$IFDEF TRACE}WriteLn('about to scale'); {$ENDIF}
   {for each row, scale all elements, including RHS, so that the largest LHS coefficent is 1.00
   对于每一行，扫描所有的元素，包括RHS，确保最大的LHS系数是1.00}
    for Row := 1 to Neq do begin
    {find biggest element选主元高斯消去法}
      ElmPtr := FirstElm^[Row];  //为什么不直接用 FirstElm^[Row]？ 指针移动！
      Biggest := 0.0;
      while ElmPtr <> nil do begin
        if (ElmPtr^.Column > Neq) then Break;
        if (Abs(ElmPtr^.Value) > Biggest) then Biggest := Abs(ElmPtr^.Value);
        ElmPtr := ElmPtr^.Next;
      end;
      if (Biggest = 0.0) then begin //有一行全为0
        SetErrorMsg('All-Zero Row', 21, Row, 0); goto Fail;
      end;
    {divide each element by biggest用最大的数除以每个元素}
      ElmPtr := FirstElm^[Row]; //将一条记录拿出来
      while ElmPtr <> nil do begin    //如果一条记录的指针不为空
        ElmPtr^.Value := ElmPtr^.Value / Biggest;
        ElmPtr := ElmPtr^.Next;
      end;
    end; {for row行扫描完了}

    if not GetMemParaAlign(Pointer(ColScale), (1 + Neq) * SizeOf(Double)) then goto OutOfSpace;
    FillChar(ColScale^, (1 + Neq) * SizeOf(Double), 0);
   {for each LHS Col, set ColScale^[Col] to largest element in col
   对于每个LHS列，设置ColScale^[Col]到列中的最大元素}
    for Row := 1 to Neq do begin
      ElmPtr := FirstElm^[Row];
      while ElmPtr <> nil do begin
        Col := ElmPtr^.Column; if (Col > Neq) then Break;
        if (Abs(ElmPtr^.Value) > ColScale^[Col]) then
          ColScale^[Col] := Abs(ElmPtr^.Value);
        ElmPtr := ElmPtr^.Next;
      end;
    end; {for row}
    for Col := 1 to Neq do if (ColScale^[Col] = 0.0) then begin
        SetErrorMsg('All-Zero Column', 22, Col, 0); goto Fail;
      end;

   {for each LHS Col, divide all elements by largest element in col
   对于每个LHS列，用列中最大的元素除以每个元素}
    for Row := 1 to Neq do begin
      ElmPtr := FirstElm^[Row];
      while ElmPtr <> nil do begin
        Col := ElmPtr^.Column; if (Col > Neq) then Break;
        ElmPtr^.Value := ElmPtr^.Value / ColScale^[Col];
        ElmPtr := ElmPtr^.Next;
      end;
    end; {for row}

  end; {if scaling}

  {begin pivoting}  //终于开始计算了
  NextActiveRow^[0] := 1;
  NextActiveRow^[Neq] := 0;

  repeat {pivot on variables which are mentioned only once重点于只提到一次的变量，是不是只有一个变量的行就求出来}
    PivotCol := 0;
    PrevRow := 0;
    Row := NextActiveRow^[0]; //row只是当前的行，而
    while Row <> 0 do begin
      NextRow := NextActiveRow^[Row];
{$IFDEF DBUG}Assert(Row > 0, '#8033'); {$ENDIF}
{$IFDEF DBUG}Assert(Row <= Neq, '#9033'); {$ENDIF}
      SingleCount := 0;
      ElmPtr := FirstElm^[Row];// FirstElm链表中提取计算所需数据？
{$IFDEF DBUG}Assert(ElmPtr <> nil, '#34'); {$ENDIF}
{$IFDEF DBUG}Assert(ElmPtr^.Column <= Neq, '#77'); {$ENDIF}
      Col := ElmPtr^.Column;
      while (Col <= Neq) do begin
        if (ColCount^[Col] = 1) then begin
          PivotCol := Col;  //单一的列？
          Inc(SingleCount);
        end;
        ElmPtr := ElmPtr^.Next;
{$IFDEF DBUG}Assert(ElmPtr <> nil, '#35'); {$ENDIF}
        Col := ElmPtr^.Column;
      end;
      if (SingleCount > 1) then begin
        SetErrorMsg('Two Singles', 23, Row, PivotCol); goto Fail;
      end
      else if (SingleCount = 1) then begin
        Inc(PivotStep);
        PivotRow := Row;
        ElmPtr := FirstElm^[PivotRow];
{$IFDEF DBUG}Assert(ElmPtr <> nil, '#34'); {$ENDIF}
{$IFDEF DBUG}Assert(ElmPtr^.Column <= Neq, '#77'); {$ENDIF}
        Col := ElmPtr^.Column;
        while (Col <= Neq) do begin
{$IFDEF DBUG}Assert(ColCount^[Col] > 0, '#4177'); {$ENDIF}
          Dec(ColCount^[Col]);
          ElmPtr := ElmPtr^.Next;
{$IFDEF DBUG}Assert(ElmPtr <> nil, '#35'); {$ENDIF}
          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); //直到 PivotCol？？=0

  {*************main loop}
  while PivotStep < Neq do begin
    Inc(PivotStep);
{$IFDEF TRACE}WriteLn('starting step ', PivotStep); {$ENDIF}

   { Find shortest row (PivotRow) and the preceding active row (LastRow)  }
    MinRowCount := High(Integer);
    PrevRow := 0; LastRow := 0;
    Row := NextActiveRow^[0];
{$IFDEF DBUG}Assert(Row <> 0, '#33'); {$ENDIF}
    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;

{$IFDEF TRACE}WriteLn('Pivotrow: ', PivotRow, ' Rowcount ', MinRowCount); {$ENDIF}
   {find Biggest Element in Pivot Row}
    Biggest := -1;
    ElmPtr := FirstElm^[PivotRow];
    while (ElmPtr^.Column <= Neq) do begin
      if (Abs(ElmPtr^.Value) > Biggest) then Biggest := Abs(ElmPtr^.Value);
      ElmPtr := ElmPtr^.Next;
    end;
    if (Biggest < 0.0) then begin {row had no elements}
      SetErrorMsg('Structurally Singular', 26, PivotRow, 0); goto Fail;
    end;
    if (Biggest = 0.0) then begin
      SetErrorMsg('Numerically Singular', 24, PivotRow, 0); goto Fail;
    end;

{$IFDEF TRACE}WriteLn('Biggest was :', Biggest); {$ENDIF}

      {find element in pivotrow with sparsest column,
      as long as it has absolute value at least UValue*biggest element}
    Biggest := Biggest * UValue;
    BestPtr := nil;
    Best_Addelm := High(Integer);
    ElmPtr := FirstElm^[PivotRow];
{$IFDEF DBUG}Assert(ElmPtr <> nil, '#36'); {$ENDIF}
    while (ElmPtr^.Column <= Neq) do begin
      Col := ElmPtr^.Column;
      Dec(ColCount^[Col]);
      if (Abs(ElmPtr^.Value) >= Biggest) then begin
        AddElm := ColCount^[Col];
     {addelm*rowcount 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;
    end;
{$IFDEF DBUG}Assert(BestPtr <> nil, '#38'); {$ENDIF}

    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}

{$IFDEF TRACE}
    Writeln('Start Pivot, Pivotrow: ', PivotRow, ' Rowcount ', MinRowCount);
    WriteLn(' PivotCol: ', PivotCol, ' ColCount ', NumToFind);
{$ENDIF}

    Row := NextActiveRow^[0];
    while ((Row <> 0) and (NumToFind > 0)) do begin
    {check that FreeList has enough items for the maximum possible number of insertions}
      if (FreeCount < (Neq - PivotStep)) then if not TopUpFreeList then goto OutOfSpace;

        {preload PrevPtr so that: PrevPtr^.Next :=  FirstElm^[Row]
        this works because the Next field of an ElmRec is at the start of the record}
      PrevPtr := Addr(FirstElm^[Row]);
      ElmPtr := FirstElm^[Row];

        {search along row looking for pivot column:
         if we find a bigger column we have gone far enough}
      while (ElmPtr^.Column < PivotCol) do begin
        PrevPtr := ElmPtr;
        ElmPtr := ElmPtr^.Next;
      end;

      if (ElmPtr^.Column = PivotCol) then begin

          {current row  (called "Tar", after target)contains pivot col,
           so we must add to it a multiple of pivot row}
{$IFDEF TRACE}WriteLn('Altering Row ', Row); {$ENDIF}

        Factor := ElmPtr^.Value / PivotValue;
        Dec(NumToFind);

     {unlink pivot col from current row}
        PrevPtr^.Next := ElmPtr^.Next; Dec(RowCount^[Row]);
     {prepend discarded node to Free List}
        ElmPtr^.Next := FreePtr; FreePtr := ElmPtr; Inc(FreeCount);

        Next_Pivot := FirstElm^[PivotRow];
{$IFDEF DBUG}Assert(Next_Pivot <> nil, '#333'); {$ENDIF}
        PrevPtr := Addr(FirstElm^[Row]); {so that PrevPtr^.Next :=  FirstElm^[Row]}
        Next_Tar := FirstElm^[Row];
{$IFDEF DBUG}Assert(Next_Tar <> nil, '#334'); {$ENDIF}
        NextTarCol := Next_Tar^.Column;
        while Next_Pivot <> nil do begin
          NextPivotCol := Next_Pivot^.Column;
          if (NextPivotCol <> PivotCol) then begin

       {Skip along Tar Row until we find a column as big as NextPivotCol}
            while NextTarCol < NextPivotCol do begin
              PrevPtr := Next_Tar;
              Next_Tar := Next_Tar^.Next;
{$IFDEF DBUG}Assert(Next_Tar <> nil, '#99'); {$ENDIF}
              NextTarCol := Next_Tar^.Column
            end;

            if (NextTarCol > NextPivotCol) then begin
        {element in pivot row but not in current row: add in new element}
{$IFDEF DBUG}Assert(NextTarCol > NextPivotCol, '#69'); {$ENDIF}
        {get element from free list}
              NewPtr := FreePtr; FreePtr := FreePtr^.Next; Dec(FreeCount); Inc(Fillins);
        {set values for new item}
              NewPtr^.Value := -Factor * Next_Pivot^.Value;
              NewPtr^.Column := NextPivotCol;
              NewPtr^.Next := Next_Tar;
        {connect previous item in list to new item}
              PrevPtr^.Next := NewPtr;
        {update PrevPtr}
              PrevPtr := NewPtr;
        {update column and row counts}
              Inc(ColCount^[NextPivotCol]);
              Inc(RowCount^[Row]);
            end {if (NextTarCol > NextPivotCol)}
            else begin
        {element in pivot row and also in current row: adjust value}
{$IFDEF DBUG}Assert(NextTarCol = NextPivotCol, '#67'); {$ENDIF}
              Next_Tar^.Value := Next_Tar^.Value - Factor * Next_Pivot^.Value;
            end; {else begin}
          end; {if (NextPivotCol <> PivotCol)}
          Next_Pivot := Next_Pivot^.Next; {move along pivot row}
        end; {while Next_Pivot <>  Nil}
      end; {if (ElmPtr^.Column = PivotCol)}
      Row := NextActiveRow^[Row];
    end; {while row}
{$IFDEF DBUG}Assert(NumToFind = 0, '#66'); {$ENDIF}
  end; {main loop}

  {release un-needed vectors}
  ReleaseItem(Pointer(ColCount));
  ReleaseItem(Pointer(NextActiveRow));

  {create Answer vector}
  if not GetMemParaAlign(Pointer(Answer), (1 + Neq) * SizeOf(Double))
    then goto OutOfSpace;
{$IFDEF DBUG}for Row := 1 to Neq do Answer^[Row] := -99; {$ENDIF}

  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];
{$IFDEF DBUG}Assert(ElmPtr <> nil, '#188'); {$ENDIF}
    while ElmPtr <> nil do begin
      Col := ElmPtr^.Column;
      if (Col = PivotCol) then
        Coeff := ElmPtr^.Value
      else if (Col <= Neq) then begin
{$IFDEF DBUG}Assert(Answer^[Col] <> -99, '#177'); {$ENDIF}
        SumTerm := SumTerm + Answer^[Col] * ElmPtr^.Value;
      end
      else RHS := ElmPtr^.Value;
      ElmPtr := ElmPtr^.Next;
    end; {until (elmptr=Nil)}
{$IFDEF DBUG}Assert(Answer^[PivotCol] = -99, '#77'); {$ENDIF}
    Answer^[PivotCol] := (RHS - SumTerm) / Coeff;
  end; {for PivotRow:=neq downto 1}

  if Scaling then for Col := 1 to Neq do Answer^[Col] := Answer^[Col] / ColScale^[Col];

  Solved := True;
  Solve1 := True;
  Exit; {normal exit}

  {exit for error conditions}
  OutOfSpace: SetErrorMsg('Out of Space', 25, PivotStep, 0); Exit;

  Fail:

end; {Solve1}

end.