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