```001:       SUBROUTINE DSPTRS( UPLO, N, NRHS, AP, IPIV, B, LDB, 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, LDB, N, NRHS
010: *     ..
011: *     .. Array Arguments ..
012:       INTEGER            IPIV( * )
013:       DOUBLE PRECISION   AP( * ), B( LDB, * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  DSPTRS solves a system of linear equations A*X = B with a real
020: *  symmetric matrix A stored in packed format using the factorization
021: *  A = U*D*U**T or A = L*D*L**T computed by DSPTRF.
022: *
023: *  Arguments
024: *  =========
025: *
026: *  UPLO    (input) CHARACTER*1
027: *          Specifies whether the details of the factorization are stored
028: *          as an upper or lower triangular matrix.
029: *          = 'U':  Upper triangular, form is A = U*D*U**T;
030: *          = 'L':  Lower triangular, form is A = L*D*L**T.
031: *
032: *  N       (input) INTEGER
033: *          The order of the matrix A.  N >= 0.
034: *
035: *  NRHS    (input) INTEGER
036: *          The number of right hand sides, i.e., the number of columns
037: *          of the matrix B.  NRHS >= 0.
038: *
039: *  AP      (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
040: *          The block diagonal matrix D and the multipliers used to
041: *          obtain the factor U or L as computed by DSPTRF, stored as a
042: *          packed triangular matrix.
043: *
044: *  IPIV    (input) INTEGER array, dimension (N)
045: *          Details of the interchanges and the block structure of D
046: *          as determined by DSPTRF.
047: *
048: *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
049: *          On entry, the right hand side matrix B.
050: *          On exit, the solution matrix X.
051: *
052: *  LDB     (input) INTEGER
053: *          The leading dimension of the array B.  LDB >= max(1,N).
054: *
055: *  INFO    (output) INTEGER
056: *          = 0:  successful exit
057: *          < 0: if INFO = -i, the i-th argument had an illegal value
058: *
059: *  =====================================================================
060: *
061: *     .. Parameters ..
062:       DOUBLE PRECISION   ONE
063:       PARAMETER          ( ONE = 1.0D+0 )
064: *     ..
065: *     .. Local Scalars ..
066:       LOGICAL            UPPER
067:       INTEGER            J, K, KC, KP
068:       DOUBLE PRECISION   AK, AKM1, AKM1K, BK, BKM1, DENOM
069: *     ..
070: *     .. External Functions ..
071:       LOGICAL            LSAME
072:       EXTERNAL           LSAME
073: *     ..
074: *     .. External Subroutines ..
075:       EXTERNAL           DGEMV, DGER, DSCAL, DSWAP, XERBLA
076: *     ..
077: *     .. Intrinsic Functions ..
078:       INTRINSIC          MAX
079: *     ..
080: *     .. Executable Statements ..
081: *
082:       INFO = 0
083:       UPPER = LSAME( UPLO, 'U' )
084:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
085:          INFO = -1
086:       ELSE IF( N.LT.0 ) THEN
087:          INFO = -2
088:       ELSE IF( NRHS.LT.0 ) THEN
089:          INFO = -3
090:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
091:          INFO = -7
092:       END IF
093:       IF( INFO.NE.0 ) THEN
094:          CALL XERBLA( 'DSPTRS', -INFO )
095:          RETURN
096:       END IF
097: *
098: *     Quick return if possible
099: *
100:       IF( N.EQ.0 .OR. NRHS.EQ.0 )
101:      \$   RETURN
102: *
103:       IF( UPPER ) THEN
104: *
105: *        Solve A*X = B, where A = U*D*U'.
106: *
107: *        First solve U*D*X = B, overwriting B with X.
108: *
109: *        K is the main loop index, decreasing from N to 1 in steps of
110: *        1 or 2, depending on the size of the diagonal blocks.
111: *
112:          K = N
113:          KC = N*( N+1 ) / 2 + 1
114:    10    CONTINUE
115: *
116: *        If K < 1, exit from loop.
117: *
118:          IF( K.LT.1 )
119:      \$      GO TO 30
120: *
121:          KC = KC - K
122:          IF( IPIV( K ).GT.0 ) THEN
123: *
124: *           1 x 1 diagonal block
125: *
126: *           Interchange rows K and IPIV(K).
127: *
128:             KP = IPIV( K )
129:             IF( KP.NE.K )
130:      \$         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
131: *
132: *           Multiply by inv(U(K)), where U(K) is the transformation
133: *           stored in column K of A.
134: *
135:             CALL DGER( K-1, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
136:      \$                 B( 1, 1 ), LDB )
137: *
138: *           Multiply by the inverse of the diagonal block.
139: *
140:             CALL DSCAL( NRHS, ONE / AP( KC+K-1 ), B( K, 1 ), LDB )
141:             K = K - 1
142:          ELSE
143: *
144: *           2 x 2 diagonal block
145: *
146: *           Interchange rows K-1 and -IPIV(K).
147: *
148:             KP = -IPIV( K )
149:             IF( KP.NE.K-1 )
150:      \$         CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
151: *
152: *           Multiply by inv(U(K)), where U(K) is the transformation
153: *           stored in columns K-1 and K of A.
154: *
155:             CALL DGER( K-2, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
156:      \$                 B( 1, 1 ), LDB )
157:             CALL DGER( K-2, NRHS, -ONE, AP( KC-( K-1 ) ), 1,
158:      \$                 B( K-1, 1 ), LDB, B( 1, 1 ), LDB )
159: *
160: *           Multiply by the inverse of the diagonal block.
161: *
162:             AKM1K = AP( KC+K-2 )
163:             AKM1 = AP( KC-1 ) / AKM1K
164:             AK = AP( KC+K-1 ) / AKM1K
165:             DENOM = AKM1*AK - ONE
166:             DO 20 J = 1, NRHS
167:                BKM1 = B( K-1, J ) / AKM1K
168:                BK = B( K, J ) / AKM1K
169:                B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
170:                B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
171:    20       CONTINUE
172:             KC = KC - K + 1
173:             K = K - 2
174:          END IF
175: *
176:          GO TO 10
177:    30    CONTINUE
178: *
179: *        Next solve U'*X = B, overwriting B with X.
180: *
181: *        K is the main loop index, increasing from 1 to N in steps of
182: *        1 or 2, depending on the size of the diagonal blocks.
183: *
184:          K = 1
185:          KC = 1
186:    40    CONTINUE
187: *
188: *        If K > N, exit from loop.
189: *
190:          IF( K.GT.N )
191:      \$      GO TO 50
192: *
193:          IF( IPIV( K ).GT.0 ) THEN
194: *
195: *           1 x 1 diagonal block
196: *
197: *           Multiply by inv(U'(K)), where U(K) is the transformation
198: *           stored in column K of A.
199: *
200:             CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, AP( KC ),
201:      \$                  1, ONE, B( K, 1 ), LDB )
202: *
203: *           Interchange rows K and IPIV(K).
204: *
205:             KP = IPIV( K )
206:             IF( KP.NE.K )
207:      \$         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
208:             KC = KC + K
209:             K = K + 1
210:          ELSE
211: *
212: *           2 x 2 diagonal block
213: *
214: *           Multiply by inv(U'(K+1)), where U(K+1) is the transformation
215: *           stored in columns K and K+1 of A.
216: *
217:             CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, AP( KC ),
218:      \$                  1, ONE, B( K, 1 ), LDB )
219:             CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB,
220:      \$                  AP( KC+K ), 1, ONE, B( K+1, 1 ), LDB )
221: *
222: *           Interchange rows K and -IPIV(K).
223: *
224:             KP = -IPIV( K )
225:             IF( KP.NE.K )
226:      \$         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
227:             KC = KC + 2*K + 1
228:             K = K + 2
229:          END IF
230: *
231:          GO TO 40
232:    50    CONTINUE
233: *
234:       ELSE
235: *
236: *        Solve A*X = B, where A = L*D*L'.
237: *
238: *        First solve L*D*X = B, overwriting B with X.
239: *
240: *        K is the main loop index, increasing from 1 to N in steps of
241: *        1 or 2, depending on the size of the diagonal blocks.
242: *
243:          K = 1
244:          KC = 1
245:    60    CONTINUE
246: *
247: *        If K > N, exit from loop.
248: *
249:          IF( K.GT.N )
250:      \$      GO TO 80
251: *
252:          IF( IPIV( K ).GT.0 ) THEN
253: *
254: *           1 x 1 diagonal block
255: *
256: *           Interchange rows K and IPIV(K).
257: *
258:             KP = IPIV( K )
259:             IF( KP.NE.K )
260:      \$         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
261: *
262: *           Multiply by inv(L(K)), where L(K) is the transformation
263: *           stored in column K of A.
264: *
265:             IF( K.LT.N )
266:      \$         CALL DGER( N-K, NRHS, -ONE, AP( KC+1 ), 1, B( K, 1 ),
267:      \$                    LDB, B( K+1, 1 ), LDB )
268: *
269: *           Multiply by the inverse of the diagonal block.
270: *
271:             CALL DSCAL( NRHS, ONE / AP( KC ), B( K, 1 ), LDB )
272:             KC = KC + N - K + 1
273:             K = K + 1
274:          ELSE
275: *
276: *           2 x 2 diagonal block
277: *
278: *           Interchange rows K+1 and -IPIV(K).
279: *
280:             KP = -IPIV( K )
281:             IF( KP.NE.K+1 )
282:      \$         CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
283: *
284: *           Multiply by inv(L(K)), where L(K) is the transformation
285: *           stored in columns K and K+1 of A.
286: *
287:             IF( K.LT.N-1 ) THEN
288:                CALL DGER( N-K-1, NRHS, -ONE, AP( KC+2 ), 1, B( K, 1 ),
289:      \$                    LDB, B( K+2, 1 ), LDB )
290:                CALL DGER( N-K-1, NRHS, -ONE, AP( KC+N-K+2 ), 1,
291:      \$                    B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
292:             END IF
293: *
294: *           Multiply by the inverse of the diagonal block.
295: *
296:             AKM1K = AP( KC+1 )
297:             AKM1 = AP( KC ) / AKM1K
298:             AK = AP( KC+N-K+1 ) / AKM1K
299:             DENOM = AKM1*AK - ONE
300:             DO 70 J = 1, NRHS
301:                BKM1 = B( K, J ) / AKM1K
302:                BK = B( K+1, J ) / AKM1K
303:                B( K, J ) = ( AK*BKM1-BK ) / DENOM
304:                B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
305:    70       CONTINUE
306:             KC = KC + 2*( N-K ) + 1
307:             K = K + 2
308:          END IF
309: *
310:          GO TO 60
311:    80    CONTINUE
312: *
313: *        Next solve L'*X = B, overwriting B with X.
314: *
315: *        K is the main loop index, decreasing from N to 1 in steps of
316: *        1 or 2, depending on the size of the diagonal blocks.
317: *
318:          K = N
319:          KC = N*( N+1 ) / 2 + 1
320:    90    CONTINUE
321: *
322: *        If K < 1, exit from loop.
323: *
324:          IF( K.LT.1 )
325:      \$      GO TO 100
326: *
327:          KC = KC - ( N-K+1 )
328:          IF( IPIV( K ).GT.0 ) THEN
329: *
330: *           1 x 1 diagonal block
331: *
332: *           Multiply by inv(L'(K)), where L(K) is the transformation
333: *           stored in column K of A.
334: *
335:             IF( K.LT.N )
336:      \$         CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
337:      \$                     LDB, AP( KC+1 ), 1, ONE, B( K, 1 ), LDB )
338: *
339: *           Interchange rows K and IPIV(K).
340: *
341:             KP = IPIV( K )
342:             IF( KP.NE.K )
343:      \$         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
344:             K = K - 1
345:          ELSE
346: *
347: *           2 x 2 diagonal block
348: *
349: *           Multiply by inv(L'(K-1)), where L(K-1) is the transformation
350: *           stored in columns K-1 and K of A.
351: *
352:             IF( K.LT.N ) THEN
353:                CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
354:      \$                     LDB, AP( KC+1 ), 1, ONE, B( K, 1 ), LDB )
355:                CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
356:      \$                     LDB, AP( KC-( N-K ) ), 1, ONE, B( K-1, 1 ),
357:      \$                     LDB )
358:             END IF
359: *
360: *           Interchange rows K and -IPIV(K).
361: *
362:             KP = -IPIV( K )
363:             IF( KP.NE.K )
364:      \$         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
365:             KC = KC - ( N-K+2 )
366:             K = K - 2
367:          END IF
368: *
369:          GO TO 90
370:   100    CONTINUE
371:       END IF
372: *
373:       RETURN
374: *
375: *     End of DSPTRS
376: *
377:       END
378: ```