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