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;