1      SUBROUTINE CPODI(A,LDA,N,DET,JOB)
2C***BEGIN PROLOGUE  CPODI
3C***DATE WRITTEN   780814   (YYMMDD)
4C***REVISION DATE  820801   (YYMMDD)
5C***REVISION HISTORY  (YYMMDD)
6C   000330  Modified array declarations.  (JEC)
7C***CATEGORY NO.  D2D1B,D3D1B
8C***KEYWORDS  COMPLEX,DETERMINANT,FACTOR,INVERSE,LINEAR ALGEBRA,LINPACK,
9C             MATRIX,POSITIVE DEFINITE
10C***AUTHOR  MOLER, C. B., (U. OF NEW MEXICO)
11C***PURPOSE  Computes the determinant and inverse of a certain COMPLEX
12C            HERMITIAN POSITIVE DEFINITE matrix using factors of CPOCO,
13C            CPOFA, or CQRDC.
14C***DESCRIPTION
15C
16C     CPODI computes the determinant and inverse of a certain
17C     complex Hermitian positive definite matrix (see below)
18C     using the factors computed by CPOCO, CPOFA or CQRDC.
19C
20C     On Entry
21C
22C        A       COMPLEX(LDA, N)
23C                the output  A  from CPOCO or CPOFA
24C                or the output  X  from CQRDC.
25C
26C        LDA     INTEGER
27C                the leading dimension of the array  A .
28C
29C        N       INTEGER
30C                the order of the matrix  A .
31C
32C        JOB     INTEGER
33C                = 11   both determinant and inverse.
34C                = 01   inverse only.
35C                = 10   determinant only.
36C
37C     On Return
38C
39C        A       If CPOCO or CPOFA was used to factor  A  then
40C                CPODI produces the upper half of INVERSE(A) .
41C                If CQRDC was used to decompose  X  then
42C                CPODI produces the upper half of INVERSE(CTRANS(X)*X)
43C                where CTRANS(X) is the conjugate transpose.
44C                Elements of  A  below the diagonal are unchanged.
45C                If the units digit of JOB is zero,  A  is unchanged.
46C
47C        DET     REAL(2)
48C                determinant of  A  or of  CTRANS(X)*X  if requested.
49C                Otherwise not referenced.
50C                Determinant = DET(1) * 10.0**DET(2)
51C                with  1.0 .LE. DET(1) .LT. 10.0
52C                or  DET(1) .EQ. 0.0 .
53C
54C     Error Condition
55C
56C        a division by zero will occur if the input factor contains
57C        a zero on the diagonal and the inverse is requested.
58C        It will not occur if the subroutines are called correctly
59C        and if CPOCO or CPOFA has set INFO .EQ. 0 .
60C
61C     LINPACK.  This version dated 08/14/78 .
62C     Cleve Moler, University of New Mexico, Argonne National Lab.
63C
64C     Subroutines and Functions
65C
66C     BLAS CAXPY,CSCAL
67C     Fortran CONJG,MOD,REAL
68C***REFERENCES  DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W.,
69C                 *LINPACK USERS  GUIDE*, SIAM, 1979.
70C***ROUTINES CALLED  CAXPY,CSCAL
71C***END PROLOGUE  CPODI
72      INTEGER LDA,N,JOB
73      COMPLEX A(LDA,*)
74      REAL DET(2)
75C
76      COMPLEX T
77      REAL S
78      INTEGER I,J,JM1,K,KP1
79C
80C     COMPUTE DETERMINANT
81C
82C***FIRST EXECUTABLE STATEMENT  CPODI
83      IF (JOB/10 .EQ. 0) GO TO 70
84         DET(1) = 1.0E0
85         DET(2) = 0.0E0
86         S = 10.0E0
87         DO 50 I = 1, N
88            DET(1) = REAL(A(I,I))**2*DET(1)
89C        ...EXIT
90            IF (DET(1) .EQ. 0.0E0) GO TO 60
91   10       IF (DET(1) .GE. 1.0E0) GO TO 20
92               DET(1) = S*DET(1)
93               DET(2) = DET(2) - 1.0E0
94            GO TO 10
95   20       CONTINUE
96   30       IF (DET(1) .LT. S) GO TO 40
97               DET(1) = DET(1)/S
98               DET(2) = DET(2) + 1.0E0
99            GO TO 30
100   40       CONTINUE
101   50    CONTINUE
102   60    CONTINUE
103   70 CONTINUE
104C
105C     COMPUTE INVERSE(R)
106C
107      IF (MOD(JOB,10) .EQ. 0) GO TO 140
108         DO 100 K = 1, N
109            A(K,K) = (1.0E0,0.0E0)/A(K,K)
110            T = -A(K,K)
111            CALL CSCAL(K-1,T,A(1,K),1)
112            KP1 = K + 1
113            IF (N .LT. KP1) GO TO 90
114            DO 80 J = KP1, N
115               T = A(K,J)
116               A(K,J) = (0.0E0,0.0E0)
117               CALL CAXPY(K,T,A(1,K),1,A(1,J),1)
118   80       CONTINUE
119   90       CONTINUE
120  100    CONTINUE
121C
122C        FORM  INVERSE(R) * CTRANS(INVERSE(R))
123C
124         DO 130 J = 1, N
125            JM1 = J - 1
126            IF (JM1 .LT. 1) GO TO 120
127            DO 110 K = 1, JM1
128               T = CONJG(A(K,J))
129               CALL CAXPY(K,T,A(1,J),1,A(1,K),1)
130  110       CONTINUE
131  120       CONTINUE
132            T = CONJG(A(J,J))
133            CALL CSCAL(J,T,A(1,J),1)
134  130    CONTINUE
135  140 CONTINUE
136      RETURN
137      END
138