matrice2.pas

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

// Unit MATRICE

// OBJECT : matrix and vector aritmetic (inversion, addition, multiplication ...)
// first written in Turbo-Pascal (02-05-1991)
// adapted to Delphi 4 (28-05-1999)

// NOTES

// Note 1 : these routines operate on dynamical matrix and dynamical vectors
// (RCOMat and RCOVec),whose dimension should have been be declared with
// the InitMat(Dim) procedure (error if if not initialize).

// Note 2 : the index (i,j..) are "zero based", this is a constraint of
// Delphi4 dynamical matrix structure.

// Note 3 : I checked this unit with the small program at the end of that file.
// The time required to invert a 200 x 200 extented matrix is of about 2 seconds,
// (AMD K6-2 overclocked at 400 MHz), and the accuracy compatible with extended
// accuracy (10E-17 .. 10E-18)
// ----------------------------------------------------------
// 12-01-02 : MODIFIED Matrice to Matrice2  (Delphi 6)

// All routines now operate on rectangular matrix, except (InvMat and SysLin)
// No more need to use the InitMat procedure (suppressed) :
// - the routines detect automaticaly the dimensions of matrix and vector
// - error code 'MatDimNul' is generated if zero lines or column in matrix and vector (See DimensionMatrice and DimensionVecteur)
// - error code 'MatMauvDim' is generated if the dimensions of matrix/vector don't allow valid result
// -
// The result matrix is dimensioned automaticaly
// ----------------------------------------------------------
// if you find that unit usefull, please let me know :
// remi CORMIER, [e-mail]

UNIT Matrice2;

INTERFACE
// I would'nt worry if you change "RCOMat" to "YourNameMat" ...
TYPE      RCOVec = array of extended;
          RCOMat = array of RCOVec;

VAR       ErreurMat : integer;  { Error code, see bellow : }

CONST     MatCorrect   = 0;
          MatNonInit   = 1;
          MatDebordTas = 2;
          MatDetNul    = 3;
          MatMauvDim   = 4;
          MatDimNulle  = 5;

// Returns the number of lines (Lig) and comlumn(Col) of [A]
// if length(A)=0 :
// - Lig=0 ; Col=0; and ErreurMat=MatDimNulle
// else :
// - Lig=length(A) ; Col=length(A[0]) ; ErreurMat=0
PROCEDURE DimensionMatrice(var A:RCOMat;var Lig,Col:longint);// created 12-01-02

// Returns the number of lines (Lig) (V)
// if length(V)=0 :
// - Lig=0 ; and Err=MatDimNulle
// else :
// - Lig=length(V) ; Err=0
PROCEDURE DimensionVecteur(var V:RCOVec;var Lig:longint);// created 12-01-02

// Returns [B], the invert matrix of [A] ([A]x[B]=[I])
// Check ErreurMat for valid solution (should be = 0)
PROCEDURE InvMat(var A,B:RCOMat); // modified 12-01-02

// Returns [C] = [A]x[B]
// if length(A)=0 or length(B)=0 then ErreurMat=MatDimNulle
// if Column(A) <> Lines(B)      then ErreurMat=MatMauvDim
PROCEDURE MatMat(var A,B,C:RCOMat);//modified 12-01-02

// Returns [U], a square unit matrix N x N (1 on diagonal, 0 everywhere else)
PROCEDURE MatUnt(N:longint;var U:RCOMat);//modified 12-01-02

// Returns [Z], a matrix with Lig lines and Col colmuns filled with zero
PROCEDURE MatZero(Lig,Col:longint;var Z:RCOMat); //created 12-01-02

// Returns a vector with N lines filled with zero
PROCEDURE VecZero(N:longint;var Z:RCOVec);//created 12-01-02

// Returns (W)=[M]x(V) ; [M] RCOMat, (V) and (W) RCOVec
PROCEDURE MatVec(var M:RCOMat;var V,W:RCOVec);//modified 12-01-02

// Returns the norm of vector (U)
FUNCTION  Norme(var U:RCOVec)    : extended;// modified 12-01-02

// Returns the norm of [A]
FUNCTION  NorMat(var A:RCOMat):extended;// modified 12-01-02

// Returns sqrt(sum(Aij^2)/nz) where nz is the number of non-zero elements of [A]
FUNCTION  SigMatZ(var A:RCOMat):extended;// created 12-01-02

// Returns sqrt(sum(Vi^2)/nz) where nz is the number of non-zero elements of (V)
FUNCTION  SigVecZ(var V:RCOVec):extended;// created 12-01-02

// Returns the scalar product of the vector (U) and (V)
FUNCTION  PrdScal(var U,V:RCOVec)  : extended;//modified 12-01-02

// Multiply the matrix [A] with x, the result in matrix [B]
PROCEDURE ScalMat(x:extended;var A,B:RCOMat); // modified 12-01-02

// Multiply the vector (U) with x, the result in vector (V)
PROCEDURE ScalVec(x:extended;var U,V:RCOVec); // modified 12-01-02

// Addition of two square matrix [A] and [B], result in [C]
PROCEDURE SomMat(var A,B,C:RCOMat); // modified 12-01-02

// Addition of two vectors, (U) and (V), result in (W)
PROCEDURE SomVec(var U,V,W:RCOVec); // modified 12-01-02

// Returns (W) = (U)-(V) ; U,V,W : Vectors
PROCEDURE DifVec(var U,V,W:RCOVec); // modified 12-01-02

// Solution of the linear system [A]x(X)=(B)
// vector (X) is the solution
// check ErreurMat for valid solution (should be = 0)
PROCEDURE SysLin(var A:RCOMat;var X:RCOVec;var B:RCOVec); // modified 12-01-02

// Returns the transposed matrix of [A] into [B] (Bij = Aji)
PROCEDURE TransMat(var A,B:RCOMat); // modified 12-01-02

// copy [A] into [B]
PROCEDURE CopieMat(var A,B:RCOMat); // created 12-01-02

// copy [U] into [V]
PROCEDURE CopieVec(var U,V:RCOVec); // created 12-01-02

// V[j] = A[Ligne,j]
PROCEDURE CopieLigneMat(var A:RCOMat;Ligne:longint;var U:RCOVec); // created 12-01-02

// V[i] = A[i,Colonne]
PROCEDURE CopieColonneMat(var A:RCOMat;Colonne:longint;var U:RCOVec); // created 12-01-02

//returns the max of A[i,j]
FUNCTION MatMax(var A:RCOMat):extended; // created 12-01-02

//returns the min of A[i,j]
FUNCTION MatMin(var A:RCOMat):extended; // created 12-01-02

//Returns [C] so as the linear system [A]x(X)=(B) if transform into :
//[C]x(Y)=(D)
//i=p : Dp = Xq ; i<>p : Di=Bi
//i=q : Yq = Bp ; i<>q : Yi=Xi
//(as "known value", Bp is replaced by Xq
// as "unknown value", Xq is replaced by Bp)
PROCEDURE PermVar(A:RCOMat;p,q:longint;var D:RCOMat;var Resu : boolean);// created 10-11-01 ; modified 12-01-02

// Returns the text (in french) associated with the error # no
FUNCTION  MsgErreurMat(no:integer):string;

IMPLEMENTATION

TYPE  RCoVecE   = array of integer;

FUNCTION Signe(x:extended):longint;
BEGIN
  if x>0
    then Result:=1
    else
      if x=0
        then Result:=0
        else Result:=-1;
END;

PROCEDURE DimensionMatrice(var A:RCOMat;var Lig,Col:longint);
BEGIN
  ErreurMat:=0;Col:=0;
  Lig:=Length(A);
  if Lig=0
    then ErreurMat:=MatDimNulle
    else Col:=length(A[0]);
END;

PROCEDURE DimensionVecteur(var V:RCOVec;var Lig:longint);
BEGIN
  Lig:=length(V);
  if Lig=0
    then ErreurMat:=MatDimNulle
    else ErreurMat:=0;
END;

PROCEDURE CopieMat(var A,B:RCOMat);
var i,j,LigA,ColA : longint;
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]:=A[i,j];
    end;
END;

PROCEDURE CopieVec(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]:=U[i];
    end;
END;

PROCEDURE CopieLigneMat(var A:RCOMat;Ligne:longint;var U:RCOVec);
var j,LigA,ColA : longint;
BEGIN
  DimensionMatrice(A,LigA,ColA);
  if ErreurMat=0 then begin
    if (ligne>=0) and (ligne<=pred(LigA))
      then begin
        SetLength(U,ColA);
        for j:=0 to pred(ColA) do
          U[j]:=A[Ligne,j];
        end
      else ErreurMat:=MatMauvDim;
    end;
END;

PROCEDURE CopieColonneMat(var A:RCOMat;Colonne:longint;var U:RCOVec);
var i,LigA,ColA : longint;
BEGIN
  DimensionMatrice(A,LigA,ColA);
  if ErreurMat=0 then begin
    if (Colonne>=0) and (Colonne<=pred(ColA))
      then begin
        SetLength(U,LigA);
        for i:=0 to pred(LigA) do
          U[i]:=A[i,Colonne];
        end
      else ErreurMat:=MatMauvDim;
    end;
END;

FUNCTION  SigMatZ(var A:RCOMat):extended;
var i,j,nz,LigA,ColA : longint;
    Aij              : extended;
BEGIN
  Result:=0;nz:=0;
  DimensionMatrice(A,LigA,ColA);
  if ErreurMat=0 then
    for i:=0 to pred(LigA) do
      for j:=0 to pred(ColA) do
        begin
          Aij:=A[i,j];
          if Aij<>0.0 then begin
            inc(nz);
            Result:=Result+sqr(Aij);
            end;
    end;
  if nz<>0 then Result:=sqrt(Result/nz);
END;

FUNCTION  SigVecZ(var V:RCOVec):extended;
var i,LigV,nz : longint;
    Vi        : extended;
BEGIN
  Result:=0;nz:=0;
  DimensionVecteur(V,LigV);
  if ErreurMat=0 then begin
    for i:=0 to pred(LigV) do begin
      Vi:=V[i];
      if Vi<>0 then begin
        inc(nz);
        Result:=Result+sqr(Vi);
        end;
      end;
    end;
  if nz<>0 then result:=sqrt(Result/nz);  
END;

PROCEDURE PermVar(A:RCOMat;p,q:longint;var D:RCOMat;var Resu : boolean);
VAR Apq,Aiq : extended;
    i,j,LigA,ColA : longint;
BEGIN
  DimensionMatrice(A,LigA,ColA);
  if ErreurMat=0 then begin
    if LigA<>ColA
      then ErreurMat:=MatMauvDim
      else begin
        Apq:=A[p,q];
        if Apq=0
          then Resu:=false
          else
            begin
              for j:=0 to pred(LigA) do
                if j=q
                  then D[p,j]:=1/Apq
                  else D[p,j]:=-A[p,j]/Apq;
              for i:=0 to pred(LigA) do
                if i<>p then
                  begin
                    Aiq:=A[i,q];
                    for j:=0 to pred(LigA) do
                      if j=q
                        then D[i,j]:=Aiq*D[p,j]
                        else D[i,j]:=A[i,j]+Aiq*D[p,j];
                  end;
              Resu:=true;
            end;
          end;
    end;
END;

FUNCTION  MsgErreurMat(no:integer):string;
BEGIN
  MsgErreurMat:='Matrice2:';
  case no of
    MatCorrect   : Result:=Result+'OK';
    MatNonInit   : Result:=Result+'unité non initialisée';
    MatDebordTas : Result:=Result+'pas assez de mémoire';
    MatDetNul    : Result:=Result+'déterminant nul';
    MatMauvDim   : Result:=Result+'dimensions inadaptées';
    MatDimNulle  : Result:=Result+'dimensions nulles';
    end;
END;

FUNCTION  NorMat(var A:RCOMat):extended;
VAR i,j,ColA,LigA : integer;
BEGIN
  Result:=0;
  DimensionMatrice(A,LigA,ColA);
  if ErreurMat=0 then begin
    for i:=0 to pred(LigA) do
      for j:=0 to pred(ColA) do
        Result:=Result+sqr(A[i,j]);
    Result:=sqrt(Result/LigA/ColA);
    end;
END;

FUNCTION MatMax(var A:RCOMat):extended;
var i,j,LigA,ColA : longint;
BEGIN
  Result:=0;
  DimensionMatrice(A,LigA,ColA);
  if ErreurMat=0 then begin
    Result:=A[0,0];
    for i:=0 to pred(LigA) do
      for j:=0 to pred(Cola) do
        if A[i,j]>Result Then Result:=A[i,j];
    end;
END;

FUNCTION MatMin(var A:RCOMat):extended;
var i,j,LigA,ColA : longint;
BEGIN
  Result:=0;
  DimensionMatrice(A,LigA,ColA);
  if ErreurMat=0 then begin
    Result:=A[0,0];
    for i:=0 to pred(LigA) do
      for j:=0 to pred(Cola) do
        if A[i,j]<Result Then Result:=A[i,j];
    end;
END;

PROCEDURE Init;
begin
end;

//PROCEDURE MatMat(A,B:RCOMat;var C:RCOMat);
PROCEDURE MatMat(var A,B,C:RCOMat);
var i,j,k,
    LigA,ColA,LigB,ColB : longint;
    Cij      : extended;
begin
  DimensionMatrice(A,LigA,ColA);
  if ErreurMat=0 then begin
    DimensionMatrice(B,LigB,ColB);
    if ErreurMat=0 then
      if ColA<>LigB
        then ErreurMat:=MatMauvDim
        else begin
          SetLength(C,LigA,ColB);
          for i:=0 to pred(LigA)  do
            for j:=0 to pred(ColB) do
              begin
                Cij:=0;
                for k:=0 to pred(ColA) do
                  Cij:=Cij+A[i,k]*B[k,j];
                C[i,j]:=Cij;
              end;
          end;
    end;
end;

PROCEDURE SomMat(var A,B,C:RCOMat);
var i,j,LigA,ColA,LigB,ColB : longint;
BEGIN
  DimensionMatrice(A,LigA,ColA);
  if ErreurMat=0 then
    begin
      DimensionMatrice(B,LigB,ColB);
      if ErreurMat=0 then
        if (LigA<>LigB) or (ColA<>ColB)
          then ErreurMat:=MatMauvDim
          else begin
            SetLength(C,LigA,ColA);
            for i:=0 to pred(LigA) do
              for j:=0 to pred(ColA) do
                C[i,j]:=A[i,j]+B[i,j];
            end;
   end;
END;

PROCEDURE TransMat(var A,B:RCOMat);
var i, j      : longint;
    LigA,ColA : longint;
begin
  DimensionMatrice(A,LigA,ColA);
  if ErreurMat=0
  then begin
    SetLength(B,ColA,LigA);
    for i:=0 to pred(LigA)  do