```001:       SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       CHARACTER          UPLO
009:       INTEGER            INFO, N
010: *     ..
011: *     .. Array Arguments ..
012:       INTEGER            IPIV( * )
013:       DOUBLE PRECISION   AP( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  DSPTRF computes the factorization of a real symmetric matrix A stored
020: *  in packed format using the Bunch-Kaufman diagonal pivoting method:
021: *
022: *     A = U*D*U**T  or  A = L*D*L**T
023: *
024: *  where U (or L) is a product of permutation and unit upper (lower)
025: *  triangular matrices, and D is symmetric and block diagonal with
026: *  1-by-1 and 2-by-2 diagonal blocks.
027: *
028: *  Arguments
029: *  =========
030: *
031: *  UPLO    (input) CHARACTER*1
032: *          = 'U':  Upper triangle of A is stored;
033: *          = 'L':  Lower triangle of A is stored.
034: *
035: *  N       (input) INTEGER
036: *          The order of the matrix A.  N >= 0.
037: *
038: *  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
039: *          On entry, the upper or lower triangle of the symmetric matrix
040: *          A, packed columnwise in a linear array.  The j-th column of A
041: *          is stored in the array AP as follows:
042: *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
043: *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
044: *
045: *          On exit, the block diagonal matrix D and the multipliers used
046: *          to obtain the factor U or L, stored as a packed triangular
047: *          matrix overwriting A (see below for further details).
048: *
049: *  IPIV    (output) INTEGER array, dimension (N)
050: *          Details of the interchanges and the block structure of D.
051: *          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
052: *          interchanged and D(k,k) is a 1-by-1 diagonal block.
053: *          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
054: *          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
055: *          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
056: *          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
057: *          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
058: *
059: *  INFO    (output) INTEGER
060: *          = 0: successful exit
061: *          < 0: if INFO = -i, the i-th argument had an illegal value
062: *          > 0: if INFO = i, D(i,i) is exactly zero.  The factorization
063: *               has been completed, but the block diagonal matrix D is
064: *               exactly singular, and division by zero will occur if it
065: *               is used to solve a system of equations.
066: *
067: *  Further Details
068: *  ===============
069: *
070: *  5-96 - Based on modifications by J. Lewis, Boeing Computer Services
071: *         Company
072: *
073: *  If UPLO = 'U', then A = U*D*U', where
074: *     U = P(n)*U(n)* ... *P(k)U(k)* ...,
075: *  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
076: *  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
077: *  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
078: *  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
079: *  that if the diagonal block D(k) is of order s (s = 1 or 2), then
080: *
081: *             (   I    v    0   )   k-s
082: *     U(k) =  (   0    I    0   )   s
083: *             (   0    0    I   )   n-k
084: *                k-s   s   n-k
085: *
086: *  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
087: *  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
088: *  and A(k,k), and v overwrites A(1:k-2,k-1:k).
089: *
090: *  If UPLO = 'L', then A = L*D*L', where
091: *     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
092: *  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
093: *  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
094: *  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
095: *  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
096: *  that if the diagonal block D(k) is of order s (s = 1 or 2), then
097: *
098: *             (   I    0     0   )  k-1
099: *     L(k) =  (   0    I     0   )  s
100: *             (   0    v     I   )  n-k-s+1
101: *                k-1   s  n-k-s+1
102: *
103: *  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
104: *  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
105: *  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
106: *
107: *  =====================================================================
108: *
109: *     .. Parameters ..
110:       DOUBLE PRECISION   ZERO, ONE
111:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
112:       DOUBLE PRECISION   EIGHT, SEVTEN
113:       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
114: *     ..
115: *     .. Local Scalars ..
116:       LOGICAL            UPPER
117:       INTEGER            I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
118:      \$                   KSTEP, KX, NPP
119:       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
120:      \$                   ROWMAX, T, WK, WKM1, WKP1
121: *     ..
122: *     .. External Functions ..
123:       LOGICAL            LSAME
124:       INTEGER            IDAMAX
125:       EXTERNAL           LSAME, IDAMAX
126: *     ..
127: *     .. External Subroutines ..
128:       EXTERNAL           DSCAL, DSPR, DSWAP, XERBLA
129: *     ..
130: *     .. Intrinsic Functions ..
131:       INTRINSIC          ABS, MAX, SQRT
132: *     ..
133: *     .. Executable Statements ..
134: *
135: *     Test the input parameters.
136: *
137:       INFO = 0
138:       UPPER = LSAME( UPLO, 'U' )
139:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
140:          INFO = -1
141:       ELSE IF( N.LT.0 ) THEN
142:          INFO = -2
143:       END IF
144:       IF( INFO.NE.0 ) THEN
145:          CALL XERBLA( 'DSPTRF', -INFO )
146:          RETURN
147:       END IF
148: *
149: *     Initialize ALPHA for use in choosing pivot block size.
150: *
151:       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
152: *
153:       IF( UPPER ) THEN
154: *
155: *        Factorize A as U*D*U' using the upper triangle of A
156: *
157: *        K is the main loop index, decreasing from N to 1 in steps of
158: *        1 or 2
159: *
160:          K = N
161:          KC = ( N-1 )*N / 2 + 1
162:    10    CONTINUE
163:          KNC = KC
164: *
165: *        If K < 1, exit from loop
166: *
167:          IF( K.LT.1 )
168:      \$      GO TO 110
169:          KSTEP = 1
170: *
171: *        Determine rows and columns to be interchanged and whether
172: *        a 1-by-1 or 2-by-2 pivot block will be used
173: *
174:          ABSAKK = ABS( AP( KC+K-1 ) )
175: *
176: *        IMAX is the row-index of the largest off-diagonal element in
177: *        column K, and COLMAX is its absolute value
178: *
179:          IF( K.GT.1 ) THEN
180:             IMAX = IDAMAX( K-1, AP( KC ), 1 )
181:             COLMAX = ABS( AP( KC+IMAX-1 ) )
182:          ELSE
183:             COLMAX = ZERO
184:          END IF
185: *
186:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
187: *
188: *           Column K is zero: set INFO and continue
189: *
190:             IF( INFO.EQ.0 )
191:      \$         INFO = K
192:             KP = K
193:          ELSE
194:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
195: *
196: *              no interchange, use 1-by-1 pivot block
197: *
198:                KP = K
199:             ELSE
200: *
201: *              JMAX is the column-index of the largest off-diagonal
202: *              element in row IMAX, and ROWMAX is its absolute value
203: *
204:                ROWMAX = ZERO
205:                JMAX = IMAX
206:                KX = IMAX*( IMAX+1 ) / 2 + IMAX
207:                DO 20 J = IMAX + 1, K
208:                   IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
209:                      ROWMAX = ABS( AP( KX ) )
210:                      JMAX = J
211:                   END IF
212:                   KX = KX + J
213:    20          CONTINUE
214:                KPC = ( IMAX-1 )*IMAX / 2 + 1
215:                IF( IMAX.GT.1 ) THEN
216:                   JMAX = IDAMAX( IMAX-1, AP( KPC ), 1 )
217:                   ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-1 ) ) )
218:                END IF
219: *
220:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
221: *
222: *                 no interchange, use 1-by-1 pivot block
223: *
224:                   KP = K
225:                ELSE IF( ABS( AP( KPC+IMAX-1 ) ).GE.ALPHA*ROWMAX ) THEN
226: *
227: *                 interchange rows and columns K and IMAX, use 1-by-1
228: *                 pivot block
229: *
230:                   KP = IMAX
231:                ELSE
232: *
233: *                 interchange rows and columns K-1 and IMAX, use 2-by-2
234: *                 pivot block
235: *
236:                   KP = IMAX
237:                   KSTEP = 2
238:                END IF
239:             END IF
240: *
241:             KK = K - KSTEP + 1
242:             IF( KSTEP.EQ.2 )
243:      \$         KNC = KNC - K + 1
244:             IF( KP.NE.KK ) THEN
245: *
246: *              Interchange rows and columns KK and KP in the leading
247: *              submatrix A(1:k,1:k)
248: *
249:                CALL DSWAP( KP-1, AP( KNC ), 1, AP( KPC ), 1 )
250:                KX = KPC + KP - 1
251:                DO 30 J = KP + 1, KK - 1
252:                   KX = KX + J - 1
253:                   T = AP( KNC+J-1 )
254:                   AP( KNC+J-1 ) = AP( KX )
255:                   AP( KX ) = T
256:    30          CONTINUE
257:                T = AP( KNC+KK-1 )
258:                AP( KNC+KK-1 ) = AP( KPC+KP-1 )
259:                AP( KPC+KP-1 ) = T
260:                IF( KSTEP.EQ.2 ) THEN
261:                   T = AP( KC+K-2 )
262:                   AP( KC+K-2 ) = AP( KC+KP-1 )
263:                   AP( KC+KP-1 ) = T
264:                END IF
265:             END IF
266: *
267: *           Update the leading submatrix
268: *
269:             IF( KSTEP.EQ.1 ) THEN
270: *
271: *              1-by-1 pivot block D(k): column k now holds
272: *
273: *              W(k) = U(k)*D(k)
274: *
275: *              where U(k) is the k-th column of U
276: *
277: *              Perform a rank-1 update of A(1:k-1,1:k-1) as
278: *
279: *              A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)'
280: *
281:                R1 = ONE / AP( KC+K-1 )
282:                CALL DSPR( UPLO, K-1, -R1, AP( KC ), 1, AP )
283: *
284: *              Store U(k) in column k
285: *
286:                CALL DSCAL( K-1, R1, AP( KC ), 1 )
287:             ELSE
288: *
289: *              2-by-2 pivot block D(k): columns k and k-1 now hold
290: *
291: *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
292: *
293: *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
294: *              of U
295: *
296: *              Perform a rank-2 update of A(1:k-2,1:k-2) as
297: *
298: *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
299: *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )'
300: *
301:                IF( K.GT.2 ) THEN
302: *
303:                   D12 = AP( K-1+( K-1 )*K / 2 )
304:                   D22 = AP( K-1+( K-2 )*( K-1 ) / 2 ) / D12
305:                   D11 = AP( K+( K-1 )*K / 2 ) / D12
306:                   T = ONE / ( D11*D22-ONE )
307:                   D12 = T / D12
308: *
309:                   DO 50 J = K - 2, 1, -1
310:                      WKM1 = D12*( D11*AP( J+( K-2 )*( K-1 ) / 2 )-
311:      \$                      AP( J+( K-1 )*K / 2 ) )
312:                      WK = D12*( D22*AP( J+( K-1 )*K / 2 )-
313:      \$                    AP( J+( K-2 )*( K-1 ) / 2 ) )
314:                      DO 40 I = J, 1, -1
315:                         AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) -
316:      \$                     AP( I+( K-1 )*K / 2 )*WK -
317:      \$                     AP( I+( K-2 )*( K-1 ) / 2 )*WKM1
318:    40                CONTINUE
319:                      AP( J+( K-1 )*K / 2 ) = WK
320:                      AP( J+( K-2 )*( K-1 ) / 2 ) = WKM1
321:    50             CONTINUE
322: *
323:                END IF
324: *
325:             END IF
326:          END IF
327: *
328: *        Store details of the interchanges in IPIV
329: *
330:          IF( KSTEP.EQ.1 ) THEN
331:             IPIV( K ) = KP
332:          ELSE
333:             IPIV( K ) = -KP
334:             IPIV( K-1 ) = -KP
335:          END IF
336: *
337: *        Decrease K and return to the start of the main loop
338: *
339:          K = K - KSTEP
340:          KC = KNC - K
341:          GO TO 10
342: *
343:       ELSE
344: *
345: *        Factorize A as L*D*L' using the lower triangle of A
346: *
347: *        K is the main loop index, increasing from 1 to N in steps of
348: *        1 or 2
349: *
350:          K = 1
351:          KC = 1
352:          NPP = N*( N+1 ) / 2
353:    60    CONTINUE
354:          KNC = KC
355: *
356: *        If K > N, exit from loop
357: *
358:          IF( K.GT.N )
359:      \$      GO TO 110
360:          KSTEP = 1
361: *
362: *        Determine rows and columns to be interchanged and whether
363: *        a 1-by-1 or 2-by-2 pivot block will be used
364: *
365:          ABSAKK = ABS( AP( KC ) )
366: *
367: *        IMAX is the row-index of the largest off-diagonal element in
368: *        column K, and COLMAX is its absolute value
369: *
370:          IF( K.LT.N ) THEN
371:             IMAX = K + IDAMAX( N-K, AP( KC+1 ), 1 )
372:             COLMAX = ABS( AP( KC+IMAX-K ) )
373:          ELSE
374:             COLMAX = ZERO
375:          END IF
376: *
377:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
378: *
379: *           Column K is zero: set INFO and continue
380: *
381:             IF( INFO.EQ.0 )
382:      \$         INFO = K
383:             KP = K
384:          ELSE
385:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
386: *
387: *              no interchange, use 1-by-1 pivot block
388: *
389:                KP = K
390:             ELSE
391: *
392: *              JMAX is the column-index of the largest off-diagonal
393: *              element in row IMAX, and ROWMAX is its absolute value
394: *
395:                ROWMAX = ZERO
396:                KX = KC + IMAX - K
397:                DO 70 J = K, IMAX - 1
398:                   IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
399:                      ROWMAX = ABS( AP( KX ) )
400:                      JMAX = J
401:                   END IF
402:                   KX = KX + N - J
403:    70          CONTINUE
404:                KPC = NPP - ( N-IMAX+1 )*( N-IMAX+2 ) / 2 + 1
405:                IF( IMAX.LT.N ) THEN
406:                   JMAX = IMAX + IDAMAX( N-IMAX, AP( KPC+1 ), 1 )
407:                   ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-IMAX ) ) )
408:                END IF
409: *
410:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
411: *
412: *                 no interchange, use 1-by-1 pivot block
413: *
414:                   KP = K
415:                ELSE IF( ABS( AP( KPC ) ).GE.ALPHA*ROWMAX ) THEN
416: *
417: *                 interchange rows and columns K and IMAX, use 1-by-1
418: *                 pivot block
419: *
420:                   KP = IMAX
421:                ELSE
422: *
423: *                 interchange rows and columns K+1 and IMAX, use 2-by-2
424: *                 pivot block
425: *
426:                   KP = IMAX
427:                   KSTEP = 2
428:                END IF
429:             END IF
430: *
431:             KK = K + KSTEP - 1
432:             IF( KSTEP.EQ.2 )
433:      \$         KNC = KNC + N - K + 1
434:             IF( KP.NE.KK ) THEN
435: *
436: *              Interchange rows and columns KK and KP in the trailing
437: *              submatrix A(k:n,k:n)
438: *
439:                IF( KP.LT.N )
440:      \$            CALL DSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ),
441:      \$                        1 )
442:                KX = KNC + KP - KK
443:                DO 80 J = KK + 1, KP - 1
444:                   KX = KX + N - J + 1
445:                   T = AP( KNC+J-KK )
446:                   AP( KNC+J-KK ) = AP( KX )
447:                   AP( KX ) = T
448:    80          CONTINUE
449:                T = AP( KNC )
450:                AP( KNC ) = AP( KPC )
451:                AP( KPC ) = T
452:                IF( KSTEP.EQ.2 ) THEN
453:                   T = AP( KC+1 )
454:                   AP( KC+1 ) = AP( KC+KP-K )
455:                   AP( KC+KP-K ) = T
456:                END IF
457:             END IF
458: *
459: *           Update the trailing submatrix
460: *
461:             IF( KSTEP.EQ.1 ) THEN
462: *
463: *              1-by-1 pivot block D(k): column k now holds
464: *
465: *              W(k) = L(k)*D(k)
466: *
467: *              where L(k) is the k-th column of L
468: *
469:                IF( K.LT.N ) THEN
470: *
471: *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
472: *
473: *                 A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
474: *
475:                   R1 = ONE / AP( KC )
476:                   CALL DSPR( UPLO, N-K, -R1, AP( KC+1 ), 1,
477:      \$                       AP( KC+N-K+1 ) )
478: *
479: *                 Store L(k) in column K
480: *
481:                   CALL DSCAL( N-K, R1, AP( KC+1 ), 1 )
482:                END IF
483:             ELSE
484: *
485: *              2-by-2 pivot block D(k): columns K and K+1 now hold
486: *
487: *              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
488: *
489: *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
490: *              of L
491: *
492:                IF( K.LT.N-1 ) THEN
493: *
494: *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
495: *
496: *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )'
497: *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )'
498: *
499:                   D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 )
500:                   D11 = AP( K+1+K*( 2*N-K-1 ) / 2 ) / D21
501:                   D22 = AP( K+( K-1 )*( 2*N-K ) / 2 ) / D21
502:                   T = ONE / ( D11*D22-ONE )
503:                   D21 = T / D21
504: *
505:                   DO 100 J = K + 2, N
506:                      WK = D21*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-
507:      \$                    AP( J+K*( 2*N-K-1 ) / 2 ) )
508:                      WKP1 = D21*( D22*AP( J+K*( 2*N-K-1 ) / 2 )-
509:      \$                      AP( J+( K-1 )*( 2*N-K ) / 2 ) )
510: *
511:                      DO 90 I = J, N
512:                         AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )*
513:      \$                     ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) /
514:      \$                     2 )*WK - AP( I+K*( 2*N-K-1 ) / 2 )*WKP1
515:    90                CONTINUE
516: *
517:                      AP( J+( K-1 )*( 2*N-K ) / 2 ) = WK
518:                      AP( J+K*( 2*N-K-1 ) / 2 ) = WKP1
519: *
520:   100             CONTINUE
521:                END IF
522:             END IF
523:          END IF
524: *
525: *        Store details of the interchanges in IPIV
526: *
527:          IF( KSTEP.EQ.1 ) THEN
528:             IPIV( K ) = KP
529:          ELSE
530:             IPIV( K ) = -KP
531:             IPIV( K+1 ) = -KP
532:          END IF
533: *
534: *        Increase K and return to the start of the main loop
535: *
536:          K = K + KSTEP
537:          KC = KNC + N - K + 2
538:          GO TO 60
539: *
540:       END IF
541: *
542:   110 CONTINUE
543:       RETURN
544: *
545: *     End of DSPTRF
546: *
547:       END
548: ```