1*DECK DSPFA
2      SUBROUTINE DSPFA (AP, N, KPVT, INFO)
3C***BEGIN PROLOGUE  DSPFA
4C***PURPOSE  Factor a real symmetric matrix stored in packed form by
5C            elimination with symmetric pivoting.
6C***LIBRARY   SLATEC (LINPACK)
7C***CATEGORY  D2B1A
8C***TYPE      DOUBLE PRECISION (SSPFA-S, DSPFA-D, CHPFA-C, CSPFA-C)
9C***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION, PACKED,
10C             SYMMETRIC
11C***AUTHOR  Bunch, J., (UCSD)
12C***DESCRIPTION
13C
14C     DSPFA factors a double precision symmetric matrix stored in
15C     packed form by elimination with symmetric pivoting.
16C
17C     To solve  A*X = B , follow DSPFA by DSPSL.
18C     To compute  INVERSE(A)*C , follow DSPFA by DSPSL.
19C     To compute  DETERMINANT(A) , follow DSPFA by DSPDI.
20C     To compute  INERTIA(A) , follow DSPFA by DSPDI.
21C     To compute  INVERSE(A) , follow DSPFA by DSPDI.
22C
23C     On Entry
24C
25C        AP      DOUBLE PRECISION (N*(N+1)/2)
26C                the packed form of a symmetric matrix  A .  The
27C                columns of the upper triangle are stored sequentially
28C                in a one-dimensional array of length  N*(N+1)/2 .
29C                See comments below for details.
30C
31C        N       INTEGER
32C                the order of the matrix  A .
33C
34C     Output
35C
36C        AP      a block diagonal matrix and the multipliers which
37C                were used to obtain it stored in packed form.
38C                The factorization can be written  A = U*D*TRANS(U)
39C                where  U  is a product of permutation and unit
40C                upper triangular matrices, TRANS(U) is the
41C                transpose of  U , and  D  is block diagonal
42C                with 1 by 1 and 2 by 2 blocks.
43C
44C        KPVT    INTEGER(N)
45C                an integer vector of pivot indices.
46C
47C        INFO    INTEGER
48C                = 0  normal value.
49C                = K  if the K-th pivot block is singular.  This is
50C                     not an error condition for this subroutine,
51C                     but it does indicate that DSPSL or DSPDI may
52C                     divide by zero if called.
53C
54C     Packed Storage
55C
56C          The following program segment will pack the upper
57C          triangle of a symmetric matrix.
58C
59C                K = 0
60C                DO 20 J = 1, N
61C                   DO 10 I = 1, J
62C                      K = K + 1
63C                      AP(K)  = A(I,J)
64C             10    CONTINUE
65C             20 CONTINUE
66C
67C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
68C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
69C***ROUTINES CALLED  DAXPY, DSWAP, IDAMAX
70C***REVISION HISTORY  (YYMMDD)
71C   780814  DATE WRITTEN
72C   890531  Changed all specific intrinsics to generic.  (WRB)
73C   890831  Modified array declarations.  (WRB)
74C   891107  Modified routine equivalence list.  (WRB)
75C   891107  REVISION DATE from Version 3.2
76C   891214  Prologue converted to Version 4.0 format.  (BAB)
77C   900326  Removed duplicate information from DESCRIPTION section.
78C           (WRB)
79C   920501  Reformatted the REFERENCES section.  (WRB)
80C***END PROLOGUE  DSPFA
81      INTEGER N,KPVT(*),INFO
82      DOUBLE PRECISION AP(*)
83C
84      DOUBLE PRECISION AK,AKM1,BK,BKM1,DENOM,MULK,MULKM1,T
85      DOUBLE PRECISION ABSAKK,ALPHA,COLMAX,ROWMAX
86      INTEGER IDAMAX,IJ,IK,IKM1,IM,IMAX,IMAXP1,IMIM,IMJ,IMK
87      INTEGER J,JJ,JK,JKM1,JMAX,JMIM,K,KK,KM1,KM1K,KM1KM1,KM2,KSTEP
88      LOGICAL SWAP
89C***FIRST EXECUTABLE STATEMENT  DSPFA
90C
91C     INITIALIZE
92C
93C     ALPHA IS USED IN CHOOSING PIVOT BLOCK SIZE.
94C
95      ALPHA = (1.0D0 + SQRT(17.0D0))/8.0D0
96C
97      INFO = 0
98C
99C     MAIN LOOP ON K, WHICH GOES FROM N TO 1.
100C
101      K = N
102      IK = (N*(N - 1))/2
103   10 CONTINUE
104C
105C        LEAVE THE LOOP IF K=0 OR K=1.
106C
107         IF (K .EQ. 0) GO TO 200
108         IF (K .GT. 1) GO TO 20
109            KPVT(1) = 1
110            IF (AP(1) .EQ. 0.0D0) INFO = 1
111            GO TO 200
112   20    CONTINUE
113C
114C        THIS SECTION OF CODE DETERMINES THE KIND OF
115C        ELIMINATION TO BE PERFORMED.  WHEN IT IS COMPLETED,
116C        KSTEP WILL BE SET TO THE SIZE OF THE PIVOT BLOCK, AND
117C        SWAP WILL BE SET TO .TRUE. IF AN INTERCHANGE IS
118C        REQUIRED.
119C
120         KM1 = K - 1
121         KK = IK + K
122         ABSAKK = ABS(AP(KK))
123C
124C        DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN
125C        COLUMN K.
126C
127         IMAX = IDAMAX(K-1,AP(IK+1),1)
128         IMK = IK + IMAX
129         COLMAX = ABS(AP(IMK))
130         IF (ABSAKK .LT. ALPHA*COLMAX) GO TO 30
131            KSTEP = 1
132            SWAP = .FALSE.
133         GO TO 90
134   30    CONTINUE
135C
136C           DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN
137C           ROW IMAX.
138C
139            ROWMAX = 0.0D0
140            IMAXP1 = IMAX + 1
141            IM = IMAX*(IMAX - 1)/2
142            IMJ = IM + 2*IMAX
143            DO 40 J = IMAXP1, K
144               ROWMAX = MAX(ROWMAX,ABS(AP(IMJ)))
145               IMJ = IMJ + J
146   40       CONTINUE
147            IF (IMAX .EQ. 1) GO TO 50
148               JMAX = IDAMAX(IMAX-1,AP(IM+1),1)
149               JMIM = JMAX + IM
150               ROWMAX = MAX(ROWMAX,ABS(AP(JMIM)))
151   50       CONTINUE
152            IMIM = IMAX + IM
153            IF (ABS(AP(IMIM)) .LT. ALPHA*ROWMAX) GO TO 60
154               KSTEP = 1
155               SWAP = .TRUE.
156            GO TO 80
157   60       CONTINUE
158            IF (ABSAKK .LT. ALPHA*COLMAX*(COLMAX/ROWMAX)) GO TO 70
159               KSTEP = 1
160               SWAP = .FALSE.
161            GO TO 80
162   70       CONTINUE
163               KSTEP = 2
164               SWAP = IMAX .NE. KM1
165   80       CONTINUE
166   90    CONTINUE
167         IF (MAX(ABSAKK,COLMAX) .NE. 0.0D0) GO TO 100
168C
169C           COLUMN K IS ZERO.  SET INFO AND ITERATE THE LOOP.
170C
171            KPVT(K) = K
172            INFO = K
173         GO TO 190
174  100    CONTINUE
175         IF (KSTEP .EQ. 2) GO TO 140
176C
177C           1 X 1 PIVOT BLOCK.
178C
179            IF (.NOT.SWAP) GO TO 120
180C
181C              PERFORM AN INTERCHANGE.
182C
183               CALL DSWAP(IMAX,AP(IM+1),1,AP(IK+1),1)
184               IMJ = IK + IMAX
185               DO 110 JJ = IMAX, K
186                  J = K + IMAX - JJ
187                  JK = IK + J
188                  T = AP(JK)
189                  AP(JK) = AP(IMJ)
190                  AP(IMJ) = T
191                  IMJ = IMJ - (J - 1)
192  110          CONTINUE
193  120       CONTINUE
194C
195C           PERFORM THE ELIMINATION.
196C
197            IJ = IK - (K - 1)
198            DO 130 JJ = 1, KM1
199               J = K - JJ
200               JK = IK + J
201               MULK = -AP(JK)/AP(KK)
202               T = MULK
203               CALL DAXPY(J,T,AP(IK+1),1,AP(IJ+1),1)
204               AP(JK) = MULK
205               IJ = IJ - (J - 1)
206  130       CONTINUE
207C
208C           SET THE PIVOT ARRAY.
209C
210            KPVT(K) = K
211            IF (SWAP) KPVT(K) = IMAX
212         GO TO 190
213  140    CONTINUE
214C
215C           2 X 2 PIVOT BLOCK.
216C
217            KM1K = IK + K - 1
218            IKM1 = IK - (K - 1)
219            IF (.NOT.SWAP) GO TO 160
220C
221C              PERFORM AN INTERCHANGE.
222C
223               CALL DSWAP(IMAX,AP(IM+1),1,AP(IKM1+1),1)
224               IMJ = IKM1 + IMAX
225               DO 150 JJ = IMAX, KM1
226                  J = KM1 + IMAX - JJ
227                  JKM1 = IKM1 + J
228                  T = AP(JKM1)
229                  AP(JKM1) = AP(IMJ)
230                  AP(IMJ) = T
231                  IMJ = IMJ - (J - 1)
232  150          CONTINUE
233               T = AP(KM1K)
234               AP(KM1K) = AP(IMK)
235               AP(IMK) = T
236  160       CONTINUE
237C
238C           PERFORM THE ELIMINATION.
239C
240            KM2 = K - 2
241            IF (KM2 .EQ. 0) GO TO 180
242               AK = AP(KK)/AP(KM1K)
243               KM1KM1 = IKM1 + K - 1
244               AKM1 = AP(KM1KM1)/AP(KM1K)
245               DENOM = 1.0D0 - AK*AKM1
246               IJ = IK - (K - 1) - (K - 2)
247               DO 170 JJ = 1, KM2
248                  J = KM1 - JJ
249                  JK = IK + J
250                  BK = AP(JK)/AP(KM1K)
251                  JKM1 = IKM1 + J
252                  BKM1 = AP(JKM1)/AP(KM1K)
253                  MULK = (AKM1*BK - BKM1)/DENOM
254                  MULKM1 = (AK*BKM1 - BK)/DENOM
255                  T = MULK
256                  CALL DAXPY(J,T,AP(IK+1),1,AP(IJ+1),1)
257                  T = MULKM1
258                  CALL DAXPY(J,T,AP(IKM1+1),1,AP(IJ+1),1)
259                  AP(JK) = MULK
260                  AP(JKM1) = MULKM1
261                  IJ = IJ - (J - 1)
262  170          CONTINUE
263  180       CONTINUE
264C
265C           SET THE PIVOT ARRAY.
266C
267            KPVT(K) = 1 - K
268            IF (SWAP) KPVT(K) = -IMAX
269            KPVT(K-1) = KPVT(K)
270  190    CONTINUE
271         IK = IK - (K - 1)
272         IF (KSTEP .EQ. 2) IK = IK - (K - 2)
273         K = K - KSTEP
274      GO TO 10
275  200 CONTINUE
276      RETURN
277      END
278