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