UNIT TrnMatrix;

INTERFACE


TYPE    
       RMatrixLine = ARRAY[1..4096] OF DOUBLE;
       RMatrix     = RECORD
                        XM,YM : WORD;
                        RM    : ^RMatrixLine;
                     END;


{
Иницииирует матрицу рахмером X на Y
}
                     
PROCEDURE InitMatrix(VAR M : RMatrix; X,Y: WORD);

{
Удаляет матрицу
}

PROCEDURE CloseMatrix(VAR M : RMatrix);


{
Заносит значение R в матрицу M в ячейку X,Y
}
PROCEDURE ToMat(VAR M: RMatrix; X,Y : WORD; R : DOUBLE);


{
Возвращает значение в ячейке X,Y матрицы M
}
FUNCTION FromMat(VAR M: RMatrix; X,Y : WORD) : DOUBLE;



{
==================================

Реализация вышеуказанных процедур и функций 
не приводится за ненадобностью - 
каждый может написать свои. 
И, возможно, даже лучше, тем те, что были у меня.

==================================
}


{
Обращает матрицу M и помещает результат в матрицу REZ. 
Матрица REZ не должна быть инициализирована.
}
PROCEDURE TurnMatrix(VAR M : RMatrix;VAR REZ : RMatrix);


IMPLEMENTATION

{

....

....

....

}


PROCEDURE TurnMatrix;
TYPE STRSTR = ARRAY[1..60] OF DOUBLE;

VAR BPL     : RMATRIX;
    SA,SAP,
    C,D,B,
    B1,AST  : STRSTR;
    NO,BK   : DOUBLE;
    I,J,K,L : WORD;

PROCEDURE PSV(SP : STRSTR;VAR ST:STRSTR);
VAR S : DOUBLE;
    ii: WORD;
BEGIN
 S := 0.0;
 FOR II := 1 TO 60 DO ST[II] := 0.0;
 FOR II := 1 TO M.XM DO
    BEGIN
      S := S +SQR(SP[II]);
    END;
 S := 1.0 /S;
 FOR II := 1 TO M.XM DO ST[II] := S*SP[II];
END;

PROCEDURE NORM(SN : STRSTR; VAR NRM : DOUBLE);
VAR II : WORD;

BEGIN
 NRM :=0.0;
 FOR II:=1 TO M.XM DO NRM :=NRM + SN[II]*SN[II];
END;

BEGIN
  IF M.XM <> M.YM THEN EXIT;
  INITMATRIX(REZ,M.XM,M.XM);
  INITMATRIX(BPL,M.XM,M.XM);
  
  FOR I := 1 TO M.XM DO
     BEGIN
       FOR J := 1 TO M.XM DO ToMat(REZ,I,J,0.0);
     END;
  FOR I := 1 TO M.XM DO SA[I] := FromMat(M,I,1);
  NO := 0.0;
  FOR I := 1 TO M.XM DO NO := NO + FromMat(M,I,1);
  IF NO = 0.0 THEN
    BEGIN
       FOR I := 1 TO M.XM DO ToMat(REZ,1,I,0.0);
    END
              ELSE
    BEGIN
       PSV(SA,SAP);
       FOR J := 1 TO M.XM DO ToMat(REZ,1,J,SAP[J]);
    END;
  FOR K := 2 TO M.XM DO
    BEGIN
      FOR I := 1 TO M.XM DO 
        BEGIN
          D[I] := 0.0;
          C[I] := 0.0;
          B1[I] := 0.0;
          B[I] := 0.0;
        END;


      FOR I := 1 TO M.XM DO AST[I] := FromMat(M,I,K);

      FOR I := 1 TO K-1 DO
        BEGIN
           FOR J := 1 TO M.XM DO
              BEGIN
                 D[I] := D[I] + FromMat(REZ,I,J)*AST[J];
              END;
        END;

      FOR I := 1 TO M.XM DO
        BEGIN
          J := 1;

          WHILE J <= K-1 DO
            BEGIN
              C[I] := C[I] + FromMat(M,I,J)*D[J];
              INC(J);
            END;

        END;

      FOR J := 1 TO M.XM DO C[J] := AST[J]-C[J];
      NORM(C,NO);
      IF NO < 1E-10 THEN
         BEGIN
           BK := 0.0;
           FOR I := 1 TO K-1 DO BK := BK +D[I]*D[I];
           BK := 1.0 /(1.0 + BK);
           FOR I := 1 TO K-1 DO B1[I] := BK*D[I];
           FOR J := 1 TO M.XM DO 
              BEGIN
                I := 1;
                WHILE I <= K-1 DO
                  BEGIN
                    B[J] := B[J] +B1[I]* FromMat(REZ,I,J);
                    INC(I);
                  END;
              END;
         END
                    ELSE
         BEGIN
            PSV(C,B);
         END;
      FOR I := 1 TO M.XM DO
        BEGIN
           FOR J := 1 TO M.XM DO ToMat(BPL,I,J,0.0);
        END;

      FOR I := 1 TO K-1 DO
         BEGIN
            FOR J := 1 TO M.XM DO ToMat(BPL,I,J,FromMat(BPL,I,J)+D[I]*B[J]);
         END;

      FOR I := 1 TO K-1 DO
         BEGIN
           FOR J := 1 TO M.XM DO ToMat(BPL,I,J,FromMat(REZ,I,J)-FromMat(BPL,I,J));
         END;

      FOR I := 1 TO K-1 DO
        BEGIN
          FOR J := 1 TO M.XM DO ToMat(REZ,I,J,FromMat(BPL,I,J));
        END;

      FOR J := 1 TO M.XM DO ToMat(REZ,K,J,B[J]);
    END;
  CloseMatrix(BPL);
END;
    
end.
