SUBROUTINE GMFADS(NN,A,NWK,MAXA) C C This subroutine computes the LDU decomposition of a symmetric positive C definite matrix B where only the upper triangular skyline structure C is stored. The decomposition is done by the Gill-Murray C strategy from P.E. Gill and W. Murray, Newton type Methods C for Unconstrained and Linearly Constrained Optimization, C Mathematical Programming, 7, 311-350 (1974) and gives an C approximate decomposition in the case of a nonpositive C definite or ill-conditioned matrix. C C Input variables: C C NN -- dimension of B. C C A -- one dimensional real array containing the upper C triangular skyline portion of a symmetric matrix B in C packed skyline storage format. C C NWK -- number of elements in A. C C MAXA -- an integer array of dimension NN+1 containing the C locations of the diagonal elements of B in A. C By convention, MAXA(NN+1)=NWK+1. C C Output variables: C C A -- the upper triangular skyline portion of the LDU C decomposition of the symmetric matrix B (or B + E if B C was not sufficiently positive definite). C C C No working storage is required by this routine. C C Subroutines called: D1MACH C INTEGER I,I0,I1,I2,I3,I4,J,J1,K,K1,K2,KH,KL,KN,KU,KZ,L,L1, & L2,L3,M,M1,MAXA(NN+1),N1,NN,NNN,NWK DOUBLE PRECISION A(NWK),BET,DEL,DJ,D1MACH,G,GAM,GAM1,PHI, & THE,THE1,XT1,XT2,ZET,ZET1 C LOGICAL GMALT C GMALT=.FALSE. G=0.0 GAM=0.0 DO 1 I=1,NN K=MAXA(I) G=G+A(K)*A(K) GAM1=ABS(A(K)) IF(GAM1.GT.GAM)GAM=GAM1 1 CONTINUE ZET=0.0 DO 3 I=1,NN K=MAXA(I) K1=MAXA(I+1)-1 K2=K1-K IF(K2.EQ.0)GO TO 3 L=K+1 DO 2 J=L,K1 G=G+2.0*A(J)*A(J) ZET1=ABS(A(J)) IF(ZET1.GT.ZET)ZET=ZET1 2 CONTINUE 3 CONTINUE ZET=ZET/NN DEL=D1MACH(4) BET=DEL IF(ZET.GT.BET)BET=ZET IF(GAM.GT.BET)BET=GAM G=SQRT(G) IF(G.GT.1.0)DEL=DEL*G DO 4 I=1,NN N1=I-1 KN=MAXA(I) KL=KN+1 KU=MAXA(I+1)-1 KH=KU-KL PHI=A(KN) IF(KH.LT.0)GO TO 10 K1=KN+1 K2=I DO 5 J=K1,KU K2=K2-1 KZ=MAXA(K2) PHI=PHI-A(J)*A(J)*A(KZ) 5 CONTINUE C10 IF(PHI.LE.0.0)GMALT=.TRUE. 10 PHI=ABS(PHI) L=I+1 THE=0.0 NNN=NN+1 IF(L.EQ.NNN)GO TO 11 DO 6 J=L,NN L1=MAXA(J) L2=MAXA(J+1) L3=L2-L1-1 M=J-I IF(L3.LT.M)GO TO 6 M1=L1+M IF(N1.EQ.0)GO TO 7 DO 8 J1=1,N1 I0=MAXA(J1) I1=MAXA(L) I2=I-J1 I3=I1-KN-1 I4=J-J1 IF(I3.LT.I2)GO TO 8 IF(L3.LT.I4)GO TO 8 XT1=A(KN+I2) XT2=A(L1+I4) A(M1)=A(M1)-XT1*XT2*A(I0) 8 CONTINUE 7 THE1=ABS(A(M1)) IF(THE.LT.THE1)THE=THE1 6 CONTINUE 11 THE=THE*THE/BET DJ=DEL IF(PHI.GT.DJ)DJ=PHI IF(THE.GT.DJ)DJ=THE C IF(ABS(DJ).NE.PHI)GMALT=.TRUE. A(KN)=DJ IF(L.EQ.NNN)GO TO 4 DO 9 J=L,NN L1=MAXA(J) L2=MAXA(J+1) L3=L2-L1-1 M=J-I IF(L3.LT.M)GO TO 9 M1=L1+M A(M1)=A(M1)/A(KN) 9 CONTINUE 4 CONTINUE RETURN END