1*DECK CPODI
2      SUBROUTINE CPODI (A, LDA, N, DET, JOB)
3C***BEGIN PROLOGUE  CPODI
4C***PURPOSE  Compute the determinant and inverse of a certain complex
5C            Hermitian positive definite matrix using the factors
6C            computed by CPOCO, CPOFA, or CQRDC.
7C***LIBRARY   SLATEC (LINPACK)
8C***CATEGORY  D2D1B, D3D1B
9C***TYPE      COMPLEX (SPODI-S, DPODI-D, CPODI-C)
10C***KEYWORDS  DETERMINANT, INVERSE, LINEAR ALGEBRA, LINPACK, MATRIX,
11C             POSITIVE DEFINITE
12C***AUTHOR  Moler, C. B., (U. of New Mexico)
13C***DESCRIPTION
14C
15C     CPODI computes the determinant and inverse of a certain
16C     complex Hermitian positive definite matrix (see below)
17C     using the factors computed by CPOCO, CPOFA or CQRDC.
18C
19C     On Entry
20C
21C        A       COMPLEX(LDA, N)
22C                the output  A  from CPOCO or CPOFA
23C                or the output  X  from CQRDC.
24C
25C        LDA     INTEGER
26C                the leading dimension of the array  A .
27C
28C        N       INTEGER
29C                the order of the matrix  A .
30C
31C        JOB     INTEGER
32C                = 11   both determinant and inverse.
33C                = 01   inverse only.
34C                = 10   determinant only.
35C
36C     On Return
37C
38C        A       If CPOCO or CPOFA was used to factor  A  then
39C                CPODI produces the upper half of INVERSE(A) .
40C                If CQRDC was used to decompose  X  then
41C                CPODI produces the upper half of INVERSE(CTRANS(X)*X)
42C                where CTRANS(X) is the conjugate transpose.
43C                Elements of  A  below the diagonal are unchanged.
44C                If the units digit of JOB is zero,  A  is unchanged.
45C
46C        DET     REAL(2)
47C                determinant of  A  or of  CTRANS(X)*X  if requested.
48C                Otherwise not referenced.
49C                Determinant = DET(1) * 10.0**DET(2)
50C                with  1.0 .LE. DET(1) .LT. 10.0
51C                or  DET(1) .EQ. 0.0 .
52C
53C     Error Condition
54C
55C        a division by zero will occur if the input factor contains
56C        a zero on the diagonal and the inverse is requested.
57C        It will not occur if the subroutines are called correctly
58C        and if CPOCO or CPOFA has set INFO .EQ. 0 .
59C
60C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
61C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
62C***ROUTINES CALLED  CAXPY, CSCAL
63C***REVISION HISTORY  (YYMMDD)
64C   780814  DATE WRITTEN
65C   890831  Modified array declarations.  (WRB)
66C   890831  REVISION DATE from Version 3.2
67C   891214  Prologue converted to Version 4.0 format.  (BAB)
68C   900326  Removed duplicate information from DESCRIPTION section.
69C           (WRB)
70C   920501  Reformatted the REFERENCES section.  (WRB)
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***FIRST EXECUTABLE STATEMENT  CPODI
80C
81C     COMPUTE DETERMINANT
82C
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)
89            IF (DET(1) .EQ. 0.0E0) GO TO 60
90   10       IF (DET(1) .GE. 1.0E0) GO TO 20
91               DET(1) = S*DET(1)
92               DET(2) = DET(2) - 1.0E0
93            GO TO 10
94   20       CONTINUE
95   30       IF (DET(1) .LT. S) GO TO 40
96               DET(1) = DET(1)/S
97               DET(2) = DET(2) + 1.0E0
98            GO TO 30
99   40       CONTINUE
100   50    CONTINUE
101   60    CONTINUE
102   70 CONTINUE
103C
104C     COMPUTE INVERSE(R)
105C
106      IF (MOD(JOB,10) .EQ. 0) GO TO 140
107         DO 100 K = 1, N
108            A(K,K) = (1.0E0,0.0E0)/A(K,K)
109            T = -A(K,K)
110            CALL CSCAL(K-1,T,A(1,K),1)
111            KP1 = K + 1
112            IF (N .LT. KP1) GO TO 90
113            DO 80 J = KP1, N
114               T = A(K,J)
115               A(K,J) = (0.0E0,0.0E0)
116               CALL CAXPY(K,T,A(1,K),1,A(1,J),1)
117   80       CONTINUE
118   90       CONTINUE
119  100    CONTINUE
120C
121C        FORM  INVERSE(R) * CTRANS(INVERSE(R))
122C
123         DO 130 J = 1, N
124            JM1 = J - 1
125            IF (JM1 .LT. 1) GO TO 120
126            DO 110 K = 1, JM1
127               T = CONJG(A(K,J))
128               CALL CAXPY(K,T,A(1,J),1,A(1,K),1)
129  110       CONTINUE
130  120       CONTINUE
131            T = CONJG(A(J,J))
132            CALL CSCAL(J,T,A(1,J),1)
133  130    CONTINUE
134  140 CONTINUE
135      RETURN
136      END
137