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