[CONTACT]

[ABOUT]

[POLICY]

[ADVERTISE]

SUBROUTINE NN,A,NWK,C

Found at: ftp.icm.edu.pl:70/packages/netlib/hompack/gmfads.f

      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

		
.


AD:

NEW PAGES:

[ODDNUGGET]

[GOPHER]