matrice2.pas

法国cromda编写的新版本MATRICE 2（矩阵和矢量运算单元）。 // ----------------------------------------------------------

      for j:=0 to pred(ColA) do
        B[j,i]:=A[i,j];
    end
end;

PROCEDURE MatUnt(N:longint;var U:RCOMat);
var i,j     : integer;
begin
  ErreurMat:=0;
  SetLength(U,N,N);
  for i:=0 to pred(N) do
    for j:=0 to pred(N) do
      if i=j then U[i,j]:=1
             else U[i,j]:=0;
end;

PROCEDURE MatZero(Lig,Col:longint;var Z:RCOMat);
var i,j     : integer;
begin
  ErreurMat:=0;
  SetLength(Z,Lig,Col);
  for i:=0 to pred(Lig) do
    for j:=0 to pred(Col) do
      Z[i,j]:=0;
end;

PROCEDURE VecZero(N:longint;var Z:RCOVec);
var i : longint;
BEGIN
  ErreurMat:=0;
  SetLength(Z,N);
  for i:=0 to pred(N) do
    Z[i]:=0;
END;

PROCEDURE MatVec(var M:RCOMat;var V,W:RCOVec);
var i, j      : longint;
    LigA,ColA : longint;
    LigV      : longint;
    Wi        : extended;
begin
  DimensionMatrice(M,LigA,ColA);
  if ErreurMat=0
    then begin
      DimensionVecteur(V,LigV);
      if ErreurMat=0
        then if ColA<>LigV
          then ErreurMat:=MatMauvDim
          else begin
            SetLength(W,LigV);
            for i:=0 to pred(LigA) do
              begin
                Wi:=0;
                for j:=0 to pred(ColA) do
                  Wi:=Wi+M[i,j]*V[j];
                W[i]:=Wi;
              end;
          end;
      end;
end;

FUNCTION  PrdScal(var U,V:RCOVec)  : extended;
var i,LigU,LigV : longint;
begin
  Result:=0;
  DimensionVecteur(U,LigU);
  if ErreurMat=0 then
    begin
      DimensionVecteur(V,LigV);
      if ErreurMat=0
        then
          if LigU<>LigV
            then ErreurMat:=MatMauvDim
            else
              for i:=0 to pred(LigU) do
                Result:=Result+U[i]*V[i];
    end;
end;

PROCEDURE ScalVec(x:extended;var U,V:RCOVec);
var i,LigU : longint;
begin
  DimensionVecteur(U,LigU);
  if ErreurMat=0
    then begin
      SetLength(V,LigU);
      for i:=0 to pred(LigU) do
        V[i]:=x*U[i];
    end;
end;

PROCEDURE ScalMat(x:extended;var A,B:RCOMat);
var i,j,LigA,ColA : integer;
BEGIN
  DimensionMatrice(A,LigA,ColA);
  if ErreurMat=0 then
    begin
      SetLength(B,LigA,ColA);
      for i:=0 to pred(LigA) do
        for j:=0 to pred(ColA) do
          B[i,j]:=x*A[i,j];
    end;
END;

PROCEDURE AddVec (Signe:extended;var U,V,W:RCOVec);
var i,LigU,LigV : integer;
begin
  DimensionVecteur(U,LigU);
  if ErreurMat=0
    then begin
      DimensionVecteur(V,LigV);
      if ErreurMat=0 then
        if LigV<>LigU
          then ErreurMat:=MatMauvDim
          else begin
            SetLength(W,LigV);
            for i:=0 to pred(LigV) do
              W[i]:=U[i]+Signe*V[i];
            end;
   end;
end;

PROCEDURE SomVec(var U,V,W:RCOVec);
begin
  AddVec(1.0,U,V,W);
end;

PROCEDURE DifVec(var U,V,W:RCOVec);
begin
  AddVec(-1.0,U,V,W);
end;

FUNCTION  Norme(var U:RCOVec)    : extended;
var i,LigU : integer;
begin
  Result:=0.0;
  DimensionVecteur(U,LigU);
  if ErreurMat=0 then
    for i:=0 to pred(LigU) do
      Result:=Result+sqr(U[i]);
  Result:=sqrt(Result);
end;

PROCEDURE TrianguleO({A:RCOMat;}Dim:longint;var AT:RCOMat;var Jx,Iy:RCOVecE;
                     var SignDet:longint;var LnDet:extended);

var i, k     : integer;
    im , jm  : integer;
    ATkk     : extended;
    C        : extended;
    Max      : extended;
//    LnDet    : extended;
//    SignDet  : longint;

  PROCEDURE InitJxIy;
  var i:integer;
  begin
    for i:=0 to pred(Dim) do
    begin
      Jx[i]:=i;
      Iy[i]:=i;
    end;
  end;
{
  PROCEDURE CopiAAT;
  var i,j:integer;
  begin
    for i:=0 to pred(Dim) do
    for j:=0 to pred(Dim) do
      AT[i,j]:=A[i,j];
  end;
}
  PROCEDURE Reduit;
  var i,j:integer;
  begin
    ATkk:=AT[k,k];
      for i:=k+1 to pred(Dim) do
        begin
          C   :=AT[i,k]/ATkk;
          for j:=k+1 to pred(Dim) do
            AT[i,j]:=AT[i,j]-C*AT[k,j];
          for j:=0 to k do
            if j<>k
              then AT[i,j]:=AT[i,j]-C*AT[k,j]
              else AT[i,j]:=-C;
        end;
  end;

  PROCEDURE ChercheMax(var Max:extended;var im,jm:integer);
  var i,j : integer;
      M   : extended;
  begin
    Max:=abs(AT[k,k]);
    im:=k;jm:=k;
    for i:=k to pred(Dim) do
      for j:=k to pred(Dim) do
        begin
          M:=abs(AT[i,j]);
          if M>Max then
            begin
              Max:=M;
              im:=i;jm:=j;
            end;
        end;
  end;

  PROCEDURE PermuteCol;
  var Ind : integer;
      i   : integer;
      x   : extended;
  begin
    for i:=0 to pred(Dim) do
      begin
        x:=AT[i,k];
        AT[i,k]:=AT[i,jm];
        AT[i,jm]:=x;
      end;
    Ind   :=Jx[k];
    Jx[k] :=Jx[jm];
    Jx[jm]:=Ind;
    SignDet:=-SignDet;
  end;

  PROCEDURE PermuteLig;
  var Ind : integer;
      j   : integer;
      x   : extended;
  begin
    for j:=0 to pred(Dim) do
      begin
        x:=AT[k,j];
        AT[k,j]:=AT[im,j];
        AT[im,j]:=x;
      end;
    Ind:=Iy[k];
    Iy[k]:=Iy[im];
    Iy[im]:=ind;
    SignDet:=-SignDet;
  end;

begin
  InitJxIy;
{
  CopiAAT;
}
  k:=0;
  SignDet:=1;
  while (k<=pred(Dim)-1) and (ErreurMat=0) do
    begin
      ChercheMax(Max,im,jm);
      if Max<>0
        then
          begin
            if k<>jm then PermuteCol;
            if k<>im then PermuteLig;
            Reduit;
          end
        else ErreurMat:=MatDetNul;
      k:=k+1;
    end;
{
  Det:=SignDet;
  for i:=0 to pred(Dim) do
    Det:=Det*AT[i,i];
  if Det=0 then ErreurMat:=MatDetNul;
}
//
  LnDet:= 0.0;
  i    :=-1;
  while (i<=pred(Dim)-1) and (ErreurMat=0) do
    begin
      inc(i);
      if AT[i,i]=0.0
        then ErreurMat:=MatDetNul
        else begin
          LnDet  :=LnDet+ln(abs(AT[i,i]));
          SignDet:=SignDet*Signe(AT[i,i]);
          end;
    end;
  
end;

PROCEDURE ResouTrO(var A:RCOMat;Dim:longint;var X,Y:RCOVec;Jx,Iy:RCOVecE);
var i, j     : integer;
    PVX      : RCOVec;
    Xi       : extended;

begin
  SetLength(X,Dim);
  SetLength(PVX,Dim);
  for i:=pred(Dim) downto 0 do
    begin
      Xi:=Y[Iy[i]];
      for j:=0 to i-1 do
        Xi:=Xi+A[i,j]*Y[Iy[j]];
      for j:=i+1 to pred(Dim) do
        Xi:=Xi-A[i,j]*PVX[j];
      PVX[i]:=Xi/A[i,i];
    end;
  for i:=0 to pred(Dim) do
    X[Jx[i]]:=PVX[i];
  Setlength(PVX,0);
  ErreurMat:=0;
end;

PROCEDURE SysLin(var A:RCOMat;var X:RCOVec;var B:RCOVec);
var {AT            : RCOMat;}
    Jx,Iy         : RCOVecE;
    LigA,ColA,Dim : longint;
    SignDet       : longint;
    LnDet         : extended;

  PROCEDURE Alloc;
  begin {SetLength(AT,Dim,Dim);}SetLength(Jx,Dim);SetLength(Iy,Dim);end;

  PROCEDURE DesAlloc;
  begin SetLength(Iy,0);SetLength(Jx,0);{SetLength(At,0,0)};end;

begin
  DimensionMatrice(A,LigA,ColA);
  if ErreurMat=0 then
    if (LigA<>ColA)
      then ErreurMat:=MatMauvDim
      else begin
        Dim:=LigA;
        Alloc;
        TrianguleO({A,}Dim,A{T},Jx,Iy,SignDet,LnDet);
        if ErreurMat=0 then ResouTrO(A{T},Dim,X,B,Jx,Iy);
        DesAlloc;
        end;
end;

PROCEDURE InvMat(var A,B:RCOMat);
var j,LigA,ColA,Dim  : longint;
{
    AT               : RCOMAt;
}
    PVX              : RCOVec;
    Jx,Iy            : RCOVecE;
    SignDet          : longint;
    LnDet            : extended;

  PROCEDURE Alloc;
  begin {Setlength(AT,Dim,Dim);}SetLength(PVX,Dim);SetLength(Jx,Dim);SetLength(Iy,Dim);end;

  PROCEDURE DesAlloc;
  begin SetLength(Iy,0);SetLength(Jx,0);SetLength(PVX,0);{SetLength(AT,0,0);}end;

  PROCEDURE initVX;
  var i : integer;
  begin
    for i:=0 to pred(Dim) do
      PVX[i]:=0;
    PVX[j]:=1;
  end;

  PROCEDURE TransfertCol;
  var i : integer;
  begin
    for i:=0 to pred(Dim) do
      B[i,j]:=PVX[i];
  end;

begin
  DimensionMatrice(A,LigA,ColA);
  if ErreurMat=0 then
    if LigA<>ColA
      then ErreurMat:=MatMauvDim
      else begin
        Dim:=LigA;
        Alloc;
        TrianguleO({A,}Dim,A{T},Jx,Iy,SignDet,LnDet);
        if ErreurMat=0
          then
            begin
              SetLength(B,Dim,Dim);
              for j:=0 to pred(Dim) do
              begin
                InitVX;
                ResouTrO(A{T},Dim,PVX,PVX,Jx,Iy);
                TransfertCol;
              end;
            end;
        DesAlloc;
        end;
end;

END.


// Test program : (I adapted it to Matrice2, but I didn't checked it)  11-01-02

procedure TForm1.FormCreate(Sender: TObject);
Const Dim = 200;
var   A,B,C : RCOMat;
      V,W,U : RCOVec;
      i,j   : integer;
      Det   : extended;
begin
  SetLength(A,Dim,Dim);
  SetLength(B,Dim,Dim);
  SetLength(C,Dim,Dim);
  Randomize;
  for i:=0 to pred(Dim) do
    begin
      V[i]:=random*10;
      for j:=0 to pred(Dim) do
      A[i,j]:=random*10;
    end;
  InvMat(A,B,Det);
  C:=MatMat(A,B);
end;

// [C] should be the unit matrix