1      SUBROUTINE CHPCO(AP,N,KPVT,RCOND,Z)
2C***BEGIN PROLOGUE  CHPCO
3C***DATE WRITTEN   780814   (YYMMDD)
4C***REVISION DATE  820801   (YYMMDD)
5C***REVISION HISTORY  (YYMMDD)
6C   000330  Modified array declarations.  (JEC)
7C***CATEGORY NO.  D2D1A
8C***KEYWORDS  COMPLEX,CONDITION,FACTOR,HERMITIAN,LINEAR ALGEBRA,LINPACK,
9C             MATRIX,PACKED
10C***AUTHOR  MOLER, C. B., (U. OF NEW MEXICO)
11C***PURPOSE  Factors a COMPLEX HERMITIAN matrix (packed form) by
12C            elimination with symmetric pivoting and estimates the
13C            condition of the matrix.
14C***DESCRIPTION
15C
16C     CHPCO factors a complex Hermitian matrix stored in packed
17C     form by elimination with symmetric pivoting and estimates
18C     the condition of the matrix.
19C
20C     if  RCOND  is not needed, CHPFA is slightly faster.
21C     To solve  A*X = B , follow CHPCO by CHPSL.
22C     To compute  INVERSE(A)*C , follow CHPCO by CHPSL.
23C     To compute  INVERSE(A) , follow CHPCO by CHPDI.
24C     To compute  DETERMINANT(A) , follow CHPCO by CHPDI.
25C     To compute  INERTIA(A), follow CHPCO by CHPDI.
26C
27C     On Entry
28C
29C        AP      COMPLEX (N*(N+1)/2)
30C                the packed form of a Hermitian matrix  A .  The
31C                columns of the upper triangle are stored sequentially
32C                in a one-dimensional array of length  N*(N+1)/2 .
33C                See comments below for details.
34C
35C        N       INTEGER
36C                the order of the matrix  A .
37C
38C     Output
39C
40C        AP      a block diagonal matrix and the multipliers which
41C                were used to obtain it stored in packed form.
42C                The factorization can be written  A = U*D*CTRANS(U)
43C                where  U  is a product of permutation and unit
44C                upper triangular matrices , CTRANS(U) is the
45C                conjugate transpose of  U , and  D  is block diagonal
46C                with 1 by 1 and 2 by 2 blocks.
47C
48C        KVPT    INTEGER(N)
49C                an integer vector of pivot indices.
50C
51C        RCOND   REAL
52C                an estimate of the reciprocal condition of  A .
53C                For the system  A*X = B , relative perturbations
54C                in  A  and  B  of size  EPSILON  may cause
55C                relative perturbations in  X  of size  EPSILON/RCOND .
56C                If  RCOND  is so small that the logical expression
57C                           1.0 + RCOND .EQ. 1.0
58C                is true, then  A  may be singular to working
59C                precision.  In particular,  RCOND  is zero  if
60C                exact singularity is detected or the estimate
61C                underflows.
62C
63C        Z       COMPLEX(N)
64C                a work vector whose contents are usually unimportant.
65C                If  A  is close to a singular matrix, then  Z  is
66C                an approximate null vector in the sense that
67C                NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
68C
69C     Packed Storage
70C
71C          The following program segment will pack the upper
72C          triangle of a Hermitian matrix.
73C
74C                K = 0
75C                DO 20 J = 1, N
76C                   DO 10 I = 1, J
77C                      K = K + 1
78C                      AP(K) = A(I,J)
79C             10    CONTINUE
80C             20 CONTINUE
81C
82C     LINPACK.  This version dated 08/14/78 .
83C     Cleve Moler, University of New Mexico, Argonne National Lab.
84C
85C     Subroutines and Functions
86C
87C     LINPACK CHPFA
88C     BLAS CAXPY,CDOTC,CSSCAL,SCASUM
89C     Fortran ABS,AIMAG,AMAX1,CMPLX,CONJG,IABS,REAL
90C***REFERENCES  DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W.,
91C                 *LINPACK USERS  GUIDE*, SIAM, 1979.
92C***ROUTINES CALLED  CAXPY,CDOTC,CHPFA,CSSCAL,SCASUM
93C***END PROLOGUE  CHPCO
94      INTEGER N,KPVT(*)
95      COMPLEX AP(*),Z(*)
96      REAL RCOND
97C
98      COMPLEX AK,AKM1,BK,BKM1,CDOTC,DENOM,EK,T
99      REAL ANORM,S,SCASUM,YNORM
100      INTEGER I,IJ,IK,IKM1,IKP1,INFO,J,JM1,J1
101      INTEGER K,KK,KM1K,KM1KM1,KP,KPS,KS
102      COMPLEX ZDUM,ZDUM2,CSIGN1
103      REAL CABS1
104      CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM))
105      CSIGN1(ZDUM,ZDUM2) = CABS1(ZDUM)*(ZDUM2/CABS1(ZDUM2))
106C
107C     FIND NORM OF A USING ONLY UPPER HALF
108C
109C***FIRST EXECUTABLE STATEMENT  CHPCO
110      J1 = 1
111      DO 30 J = 1, N
112         Z(J) = CMPLX(SCASUM(J,AP(J1),1),0.0E0)
113         IJ = J1
114         J1 = J1 + J
115         JM1 = J - 1
116         IF (JM1 .LT. 1) GO TO 20
117         DO 10 I = 1, JM1
118            Z(I) = CMPLX(REAL(Z(I))+CABS1(AP(IJ)),0.0E0)
119            IJ = IJ + 1
120   10    CONTINUE
121   20    CONTINUE
122   30 CONTINUE
123      ANORM = 0.0E0
124      DO 40 J = 1, N
125         ANORM = AMAX1(ANORM,REAL(Z(J)))
126   40 CONTINUE
127C
128C     FACTOR
129C
130      CALL CHPFA(AP,N,KPVT,INFO)
131C
132C     RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
133C     ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  A*Y = E .
134C     THE COMPONENTS OF  E  ARE CHOSEN TO CAUSE MAXIMUM LOCAL
135C     GROWTH IN THE ELEMENTS OF W  WHERE  U*D*W = E .
136C     THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW.
137C
138C     SOLVE U*D*W = E
139C
140      EK = (1.0E0,0.0E0)
141      DO 50 J = 1, N
142         Z(J) = (0.0E0,0.0E0)
143   50 CONTINUE
144      K = N
145      IK = (N*(N - 1))/2
146   60 IF (K .EQ. 0) GO TO 120
147         KK = IK + K
148         IKM1 = IK - (K - 1)
149         KS = 1
150         IF (KPVT(K) .LT. 0) KS = 2
151         KP = IABS(KPVT(K))
152         KPS = K + 1 - KS
153         IF (KP .EQ. KPS) GO TO 70
154            T = Z(KPS)
155            Z(KPS) = Z(KP)
156            Z(KP) = T
157   70    CONTINUE
158         IF (CABS1(Z(K)) .NE. 0.0E0) EK = CSIGN1(EK,Z(K))
159         Z(K) = Z(K) + EK
160         CALL CAXPY(K-KS,Z(K),AP(IK+1),1,Z(1),1)
161         IF (KS .EQ. 1) GO TO 80
162            IF (CABS1(Z(K-1)) .NE. 0.0E0) EK = CSIGN1(EK,Z(K-1))
163            Z(K-1) = Z(K-1) + EK
164            CALL CAXPY(K-KS,Z(K-1),AP(IKM1+1),1,Z(1),1)
165   80    CONTINUE
166         IF (KS .EQ. 2) GO TO 100
167            IF (CABS1(Z(K)) .LE. CABS1(AP(KK))) GO TO 90
168               S = CABS1(AP(KK))/CABS1(Z(K))
169               CALL CSSCAL(N,S,Z,1)
170               EK = CMPLX(S,0.0E0)*EK
171   90       CONTINUE
172            IF (CABS1(AP(KK)) .NE. 0.0E0) Z(K) = Z(K)/AP(KK)
173            IF (CABS1(AP(KK)) .EQ. 0.0E0) Z(K) = (1.0E0,0.0E0)
174         GO TO 110
175  100    CONTINUE
176            KM1K = IK + K - 1
177            KM1KM1 = IKM1 + K - 1
178            AK = AP(KK)/CONJG(AP(KM1K))
179            AKM1 = AP(KM1KM1)/AP(KM1K)
180            BK = Z(K)/CONJG(AP(KM1K))
181            BKM1 = Z(K-1)/AP(KM1K)
182            DENOM = AK*AKM1 - 1.0E0
183            Z(K) = (AKM1*BK - BKM1)/DENOM
184            Z(K-1) = (AK*BKM1 - BK)/DENOM
185  110    CONTINUE
186         K = K - KS
187         IK = IK - K
188         IF (KS .EQ. 2) IK = IK - (K + 1)
189      GO TO 60
190  120 CONTINUE
191      S = 1.0E0/SCASUM(N,Z,1)
192      CALL CSSCAL(N,S,Z,1)
193C
194C     SOLVE CTRANS(U)*Y = W
195C
196      K = 1
197      IK = 0
198  130 IF (K .GT. N) GO TO 160
199         KS = 1
200         IF (KPVT(K) .LT. 0) KS = 2
201         IF (K .EQ. 1) GO TO 150
202            Z(K) = Z(K) + CDOTC(K-1,AP(IK+1),1,Z(1),1)
203            IKP1 = IK + K
204            IF (KS .EQ. 2)
205     1         Z(K+1) = Z(K+1) + CDOTC(K-1,AP(IKP1+1),1,Z(1),1)
206            KP = IABS(KPVT(K))
207            IF (KP .EQ. K) GO TO 140
208               T = Z(K)
209               Z(K) = Z(KP)
210               Z(KP) = T
211  140       CONTINUE
212  150    CONTINUE
213         IK = IK + K
214         IF (KS .EQ. 2) IK = IK + (K + 1)
215         K = K + KS
216      GO TO 130
217  160 CONTINUE
218      S = 1.0E0/SCASUM(N,Z,1)
219      CALL CSSCAL(N,S,Z,1)
220C
221      YNORM = 1.0E0
222C
223C     SOLVE U*D*V = Y
224C
225      K = N
226      IK = N*(N - 1)/2
227  170 IF (K .EQ. 0) GO TO 230
228         KK = IK + K
229         IKM1 = IK - (K - 1)
230         KS = 1
231         IF (KPVT(K) .LT. 0) KS = 2
232         IF (K .EQ. KS) GO TO 190
233            KP = IABS(KPVT(K))
234            KPS = K + 1 - KS
235            IF (KP .EQ. KPS) GO TO 180
236               T = Z(KPS)
237               Z(KPS) = Z(KP)
238               Z(KP) = T
239  180       CONTINUE
240            CALL CAXPY(K-KS,Z(K),AP(IK+1),1,Z(1),1)
241            IF (KS .EQ. 2) CALL CAXPY(K-KS,Z(K-1),AP(IKM1+1),1,Z(1),1)
242  190    CONTINUE
243         IF (KS .EQ. 2) GO TO 210
244            IF (CABS1(Z(K)) .LE. CABS1(AP(KK))) GO TO 200
245               S = CABS1(AP(KK))/CABS1(Z(K))
246               CALL CSSCAL(N,S,Z,1)
247               YNORM = S*YNORM
248  200       CONTINUE
249            IF (CABS1(AP(KK)) .NE. 0.0E0) Z(K) = Z(K)/AP(KK)
250            IF (CABS1(AP(KK)) .EQ. 0.0E0) Z(K) = (1.0E0,0.0E0)
251         GO TO 220
252  210    CONTINUE
253            KM1K = IK + K - 1
254            KM1KM1 = IKM1 + K - 1
255            AK = AP(KK)/CONJG(AP(KM1K))
256            AKM1 = AP(KM1KM1)/AP(KM1K)
257            BK = Z(K)/CONJG(AP(KM1K))
258            BKM1 = Z(K-1)/AP(KM1K)
259            DENOM = AK*AKM1 - 1.0E0
260            Z(K) = (AKM1*BK - BKM1)/DENOM
261            Z(K-1) = (AK*BKM1 - BK)/DENOM
262  220    CONTINUE
263         K = K - KS
264         IK = IK - K
265         IF (KS .EQ. 2) IK = IK - (K + 1)
266      GO TO 170
267  230 CONTINUE
268      S = 1.0E0/SCASUM(N,Z,1)
269      CALL CSSCAL(N,S,Z,1)
270      YNORM = S*YNORM
271C
272C     SOLVE CTRANS(U)*Z = V
273C
274      K = 1
275      IK = 0
276  240 IF (K .GT. N) GO TO 270
277         KS = 1
278         IF (KPVT(K) .LT. 0) KS = 2
279         IF (K .EQ. 1) GO TO 260
280            Z(K) = Z(K) + CDOTC(K-1,AP(IK+1),1,Z(1),1)
281            IKP1 = IK + K
282            IF (KS .EQ. 2)
283     1         Z(K+1) = Z(K+1) + CDOTC(K-1,AP(IKP1+1),1,Z(1),1)
284            KP = IABS(KPVT(K))
285            IF (KP .EQ. K) GO TO 250
286               T = Z(K)
287               Z(K) = Z(KP)
288               Z(KP) = T
289  250       CONTINUE
290  260    CONTINUE
291         IK = IK + K
292         IF (KS .EQ. 2) IK = IK + (K + 1)
293         K = K + KS
294      GO TO 240
295  270 CONTINUE
296C     MAKE ZNORM = 1.0
297      S = 1.0E0/SCASUM(N,Z,1)
298      CALL CSSCAL(N,S,Z,1)
299      YNORM = S*YNORM
300C
301      IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM
302      IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0
303      RETURN
304      END
305