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