1      SUBROUTINE CSICO(A,LDA,N,KPVT,RCOND,Z)
2C***BEGIN PROLOGUE  CSICO
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             SYMMETRIC
10C***AUTHOR  MOLER, C. B., (U. OF NEW MEXICO)
11C***PURPOSE  Factors a COMPLEX SYMMETRIC matrix by elimination with
12C            symmetric pivoting and estimates the condition of the
13C            matrix.
14C***DESCRIPTION
15C
16C     CSICO factors a complex symmetric matrix by elimination with
17C     symmetric pivoting and estimates the condition of the matrix.
18C
19C     If  RCOND  is not needed, CSIFA is slightly faster.
20C     To solve  A*X = B , follow CSICO by CSISL.
21C     To compute  INVERSE(A)*C , follow CSICO by CSISL.
22C     To compute  INVERSE(A) , follow CSICO by CSIDI.
23C     To compute  DETERMINANT(A) , follow CSICO by CSIDI.
24C
25C     On Entry
26C
27C        A       COMPLEX(LDA, N)
28C                the symmetric matrix to be factored.
29C                Only the diagonal and upper triangle are used.
30C
31C        LDA     INTEGER
32C                the leading dimension of the array  A .
33C
34C        N       INTEGER
35C                the order of the matrix  A .
36C
37C     On Return
38C
39C        A       a block diagonal matrix and the multipliers which
40C                were used to obtain it.
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     LINPACK.  This version dated 08/14/78 .
69C     Cleve Moler, University of New Mexico, Argonne National Lab.
70C
71C     Subroutines and Functions
72C
73C     LINPACK CSIFA
74C     BLAS CAXPY,CDOTU,CSSCAL,SCASUM
75C     Fortran ABS,AIMAG,AMAX1,CMPLX,IABS,REAL
76C***REFERENCES  DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W.,
77C                 *LINPACK USERS  GUIDE*, SIAM, 1979.
78C***ROUTINES CALLED  CAXPY,CDOTU,CSIFA,CSSCAL,SCASUM
79C***END PROLOGUE  CSICO
80      INTEGER LDA,N,KPVT(*)
81      COMPLEX A(LDA,*),Z(*)
82      REAL RCOND
83C
84      COMPLEX AK,AKM1,BK,BKM1,CDOTU,DENOM,EK,T
85      REAL ANORM,S,SCASUM,YNORM
86      INTEGER I,INFO,J,JM1,K,KP,KPS,KS
87      COMPLEX ZDUM,ZDUM2,CSIGN1
88      REAL CABS1
89      CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM))
90      CSIGN1(ZDUM,ZDUM2) = CABS1(ZDUM)*(ZDUM2/CABS1(ZDUM2))
91C
92C     FIND NORM OF A USING ONLY UPPER HALF
93C
94C***FIRST EXECUTABLE STATEMENT  CSICO
95      DO 30 J = 1, N
96         Z(J) = CMPLX(SCASUM(J,A(1,J),1),0.0E0)
97         JM1 = J - 1
98         IF (JM1 .LT. 1) GO TO 20
99         DO 10 I = 1, JM1
100            Z(I) = CMPLX(REAL(Z(I))+CABS1(A(I,J)),0.0E0)
101   10    CONTINUE
102   20    CONTINUE
103   30 CONTINUE
104      ANORM = 0.0E0
105      DO 40 J = 1, N
106         ANORM = AMAX1(ANORM,REAL(Z(J)))
107   40 CONTINUE
108C
109C     FACTOR
110C
111      CALL CSIFA(A,LDA,N,KPVT,INFO)
112C
113C     RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
114C     ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  A*Y = E .
115C     THE COMPONENTS OF  E  ARE CHOSEN TO CAUSE MAXIMUM LOCAL
116C     GROWTH IN THE ELEMENTS OF W  WHERE  U*D*W = E .
117C     THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW.
118C
119C     SOLVE U*D*W = E
120C
121      EK = (1.0E0,0.0E0)
122      DO 50 J = 1, N
123         Z(J) = (0.0E0,0.0E0)
124   50 CONTINUE
125      K = N
126   60 IF (K .EQ. 0) GO TO 120
127         KS = 1
128         IF (KPVT(K) .LT. 0) KS = 2
129         KP = IABS(KPVT(K))
130         KPS = K + 1 - KS
131         IF (KP .EQ. KPS) GO TO 70
132            T = Z(KPS)
133            Z(KPS) = Z(KP)
134            Z(KP) = T
135   70    CONTINUE
136         IF (CABS1(Z(K)) .NE. 0.0E0) EK = CSIGN1(EK,Z(K))
137         Z(K) = Z(K) + EK
138         CALL CAXPY(K-KS,Z(K),A(1,K),1,Z(1),1)
139         IF (KS .EQ. 1) GO TO 80
140            IF (CABS1(Z(K-1)) .NE. 0.0E0) EK = CSIGN1(EK,Z(K-1))
141            Z(K-1) = Z(K-1) + EK
142            CALL CAXPY(K-KS,Z(K-1),A(1,K-1),1,Z(1),1)
143   80    CONTINUE
144         IF (KS .EQ. 2) GO TO 100
145            IF (CABS1(Z(K)) .LE. CABS1(A(K,K))) GO TO 90
146               S = CABS1(A(K,K))/CABS1(Z(K))
147               CALL CSSCAL(N,S,Z,1)
148               EK = CMPLX(S,0.0E0)*EK
149   90       CONTINUE
150            IF (CABS1(A(K,K)) .NE. 0.0E0) Z(K) = Z(K)/A(K,K)
151            IF (CABS1(A(K,K)) .EQ. 0.0E0) Z(K) = (1.0E0,0.0E0)
152         GO TO 110
153  100    CONTINUE
154            AK = A(K,K)/A(K-1,K)
155            AKM1 = A(K-1,K-1)/A(K-1,K)
156            BK = Z(K)/A(K-1,K)
157            BKM1 = Z(K-1)/A(K-1,K)
158            DENOM = AK*AKM1 - 1.0E0
159            Z(K) = (AKM1*BK - BKM1)/DENOM
160            Z(K-1) = (AK*BKM1 - BK)/DENOM
161  110    CONTINUE
162         K = K - KS
163      GO TO 60
164  120 CONTINUE
165      S = 1.0E0/SCASUM(N,Z,1)
166      CALL CSSCAL(N,S,Z,1)
167C
168C     SOLVE TRANS(U)*Y = W
169C
170      K = 1
171  130 IF (K .GT. N) GO TO 160
172         KS = 1
173         IF (KPVT(K) .LT. 0) KS = 2
174         IF (K .EQ. 1) GO TO 150
175            Z(K) = Z(K) + CDOTU(K-1,A(1,K),1,Z(1),1)
176            IF (KS .EQ. 2)
177     1         Z(K+1) = Z(K+1) + CDOTU(K-1,A(1,K+1),1,Z(1),1)
178            KP = IABS(KPVT(K))
179            IF (KP .EQ. K) GO TO 140
180               T = Z(K)
181               Z(K) = Z(KP)
182               Z(KP) = T
183  140       CONTINUE
184  150    CONTINUE
185         K = K + KS
186      GO TO 130
187  160 CONTINUE
188      S = 1.0E0/SCASUM(N,Z,1)
189      CALL CSSCAL(N,S,Z,1)
190C
191      YNORM = 1.0E0
192C
193C     SOLVE U*D*V = Y
194C
195      K = N
196  170 IF (K .EQ. 0) GO TO 230
197         KS = 1
198         IF (KPVT(K) .LT. 0) KS = 2
199         IF (K .EQ. KS) GO TO 190
200            KP = IABS(KPVT(K))
201            KPS = K + 1 - KS
202            IF (KP .EQ. KPS) GO TO 180
203               T = Z(KPS)
204               Z(KPS) = Z(KP)
205               Z(KP) = T
206  180       CONTINUE
207            CALL CAXPY(K-KS,Z(K),A(1,K),1,Z(1),1)
208            IF (KS .EQ. 2) CALL CAXPY(K-KS,Z(K-1),A(1,K-1),1,Z(1),1)
209  190    CONTINUE
210         IF (KS .EQ. 2) GO TO 210
211            IF (CABS1(Z(K)) .LE. CABS1(A(K,K))) GO TO 200
212               S = CABS1(A(K,K))/CABS1(Z(K))
213               CALL CSSCAL(N,S,Z,1)
214               YNORM = S*YNORM
215  200       CONTINUE
216            IF (CABS1(A(K,K)) .NE. 0.0E0) Z(K) = Z(K)/A(K,K)
217            IF (CABS1(A(K,K)) .EQ. 0.0E0) Z(K) = (1.0E0,0.0E0)
218         GO TO 220
219  210    CONTINUE
220            AK = A(K,K)/A(K-1,K)
221            AKM1 = A(K-1,K-1)/A(K-1,K)
222            BK = Z(K)/A(K-1,K)
223            BKM1 = Z(K-1)/A(K-1,K)
224            DENOM = AK*AKM1 - 1.0E0
225            Z(K) = (AKM1*BK - BKM1)/DENOM
226            Z(K-1) = (AK*BKM1 - BK)/DENOM
227  220    CONTINUE
228         K = K - KS
229      GO TO 170
230  230 CONTINUE
231      S = 1.0E0/SCASUM(N,Z,1)
232      CALL CSSCAL(N,S,Z,1)
233      YNORM = S*YNORM
234C
235C     SOLVE TRANS(U)*Z = V
236C
237      K = 1
238  240 IF (K .GT. N) GO TO 270
239         KS = 1
240         IF (KPVT(K) .LT. 0) KS = 2
241         IF (K .EQ. 1) GO TO 260
242            Z(K) = Z(K) + CDOTU(K-1,A(1,K),1,Z(1),1)
243            IF (KS .EQ. 2)
244     1         Z(K+1) = Z(K+1) + CDOTU(K-1,A(1,K+1),1,Z(1),1)
245            KP = IABS(KPVT(K))
246            IF (KP .EQ. K) GO TO 250
247               T = Z(K)
248               Z(K) = Z(KP)
249               Z(KP) = T
250  250       CONTINUE
251  260    CONTINUE
252         K = K + KS
253      GO TO 240
254  270 CONTINUE
255C     MAKE ZNORM = 1.0
256      S = 1.0E0/SCASUM(N,Z,1)
257      CALL CSSCAL(N,S,Z,1)
258      YNORM = S*YNORM
259C
260      IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM
261      IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0
262      RETURN
263      END
264