NUMAL Section 3.1.1.1.1.3.1

BEGIN SECTION : 3.1.1.1.1.3.1 (, december 1978 )

AUTHORS : J.R.BUNCH,L,KAUFMAN,B.N.PARLETT.

CONTRIBUTOR : C.H.CONVENT.

INSTITUTE : UNIVERSITY OF AMSTERDAM.

RECEIVED : 770712.

BRIEF DESCRIPTION :

    DECSYM2 CALCULATES THE LDL' DECOMPOSITION OF A SYMMETRIC MATRIX.
    THE MATRIX MAY BE INDEFINITE AND/OR SINGULAR;

KEYWORDS :

    GENERAL SYMMETRIC MATRIX,
    LDL' DECOMPOSITION,
    BLOCK DIAGONAL PIVOTING;

CALLING SEQUENCE :

    THE HEADING OF THE PROCEDURE READS :
    "PROCEDURE" DECSYM2(A,N,TOL,AUX,P,DETAUX);
    "VALUE" N;"INTEGER" N;"REAL" TOL;
    "ARRAY" A,DETAUX;"INTEGER" "ARRAY" AUX,P;
    "CODE" 34291;

    THE MEANING OF THE FORMAL PARAMETERS IS :
    A :     <ARRAY IDENTIFIER>;
           "ARRAY" A[1:N,1:N];
           ENTRY : THE SYMMETRIC COEFFICIENT MATRIX;
           EXIT  : THE ELEMENTS OF THE LDL' DECOMPOSITION OF A ARE
                   STORED IN THE UPPER TRIANGULAR PART OF A. HERE D
                   IS A BLOCK DIAGONAL MATRIX WITH BLOCKS OF ORDER 1
                   OR 2. FOR A BLOCK OF ORDER 2 WE ALWAYS HAVE
                   D[I,I+1]^=0 AND L[I+1,I]=0,SO THAT D AND L' FIT
                   IN THE UPPER TRIANGULAR PART OF A.
                   THE STRICTLY LOWER TRIANGULAR PART OF A IS LEFT
                   UNDISTURBED.
    N :     <ARITHMETIC EXPRESSION>;
           THE ORDER OF THE MATRIX;

    TOL :  <ARITHMETIC EXPRESSION>;
          A RELATIVE TOLERANCE, USED TO CONTROL THE
          CALCULATION OF THE BLOCK DIAGONAL ELEMENTS;
    AUX :   <ARRAY IDENTIFIER>;
           "INTEGER" "ARRAY" AUX[2:5];
           EXIT :
           AUX[2] : IF THE MATRIX IS SYMMETRIC THEN 1, OTHER-
                    WISE 0;IN THE LAST CASE NO DECOMPOSITION
                    IS PERFORMED;
           AUX[3] : IF THE MATRIX IS SYMMETRIC THEN THE
                    NUMBER OF ITS POSITIVE EIGENVALUES,
                    OTHERWISE 0. IF AUX[3]=N THEN THE
                    MATRIX IS POSITIVE DEFINITE;
           AUX[4] : IF THE MATRIX IS SYMMETRIC THEN THE
                    NUMBER OF ITS NEGATIVE EIGENVALUES,
                    OTHERWISE 0. IF AUX[4]=N THEN THE
                    MATRIX IS NEGATIVE DEFINITE;
           AUX[5] : IF THE MATRIX IS SYMMETRIC THEN THE
                    NUMBER OF ITS ZERO EIGENVALUES, OTHERWISE
                    N; SO, IF AUX[5]=0 THEN THE MATRIX IS
                    SYMMETRIC AND NON-SINGULAR;
    P :     <ARRAY IDENTIFIER>;
           "INTEGER" "ARRAY" P[1:N];
           EXIT :
           A VECTOR RECORDING
           1) THE INTERCHANGES PERFORMED ON A DURING THE
              COMPUTATION OF THE DECOMPOSITION AND
           2) THE BLOCK STRUCTURE OF D.
           IF P[I]>0 AND P[I+1]=0 A 2*2 BLOCK HAS BEEN
           FOUND I.E. D[I,I+1]^=0 AND L[I+1,I]=0;
    DETAUX : <ARRAY IDENTIFIER>;
            "ARRAY" DETAUX[1:N];
            EXIT :
            IF P[I]>0 AND P[I+1]>0 THEN DETAUX[I] EQUALS
            THE EXIT-VALUE OF A[I,I].
            IF P[I]>0 AND P[I+1]=0 THEN DETAUX[I]=1 AND
            DETAUX[I+1] EQUALS THE VALUE OF THE DETERMINANT
            OF THE CORRESPONDING 2*2 DIAGONAL BLOCK AS
            DETERMINED BY DECSYM2;

PROCEDURES USED :

    ELMROW=CP34030.
    ICHROW=CP34032.
    ICHROWCOL=CP34033.

RUNNING TIME : ROUGHLY PROPORTIONAL TO N**3.

METHOD AND PERFORMANCE :

    THE PROCEDURE DECSYM2 COMPUTES THE LDL' DECOMPOSITION OF A
    SYMMETRIC MATRIX,ACCORDING TO A METHOD DUE TO BUNCH,KAUFMAN AND
    PARLETT (SEE [1],[2]). IT USES BLOCK DIAGONAL PIVOTING.
    THE BLOCK DIAGONAL MATRIX D IS DELIVERED IN THE BLOCK DIAGONAL
    OF A. IF P[I]>0 AND P[I+1]=0 A 2*2 BLOCK HAS BEEN FOUND AND
    FURTHERMORE : L[I+1,I]=0 WHEN D[I,I+1]^=0.
    THE STRICTLY UPPER TRIANGULAR PART OF L' IS DELIVERED IN THE
    STRICTLY UPPER TRIANGULAR PART OF A. FOR THE INERTIA PROBLEM
    IT IS IMPORTANT THAT DECSYM2 CAN ACCEPT SINGULAR MATRICES.
    NOTE, HOWEVER, THAT IN ORDER TO FIND THE NUMBER OF ZERO
    EIGENVALUES OF SINGULAR MATRICES, THE SINGULAR VALUE
    DECOMPOSITION MIGHT BE PREFERRED.
    BEFORE THE DECOMPOSITION IS PERFORMED A CHECK IS MADE TO SEE
    WHETHER THE MATRIX IS SYMMETRIC. IF THE MATRIX IS ASYMMETRIC
    THEN NO DECOMPOSITION IS PERFORMED;

REFERENCES :

    1) J.R.BUNCH,L.KAUFMAN.
      SOME STABLE METHODS FOR CALCULATING INERTIA AND SOLVING
      SYMMETRIC LINEAR SYSTEMS.
      MATHEMATICS OF COMPUTATION 31,P 163-180,1977.
    2) J.R.BUNCH,L.KAUFMAN,B.N.PARLETT.
      DECOMPOSITION OF A SYMMETRIC MATRIX.
      NUMERISCHE MATHEMATIK 27,P 95-109,1976.

SOURCE TEXT :
"CODE" 34291;