001:       SUBROUTINE CHETRS( 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: *  CHETRS solves a system of linear equations A*X = B with a complex
021: *  Hermitian matrix A using the factorization A = U*D*U**H or
022: *  A = L*D*L**H computed by CHETRF.
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**H;
031: *          = 'L':  Lower triangular, form is A = L*D*L**H.
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 CHETRF.
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 CHETRF.
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:       REAL               S
072:       COMPLEX            AK, AKM1, AKM1K, BK, BKM1, DENOM
073: *     ..
074: *     .. External Functions ..
075:       LOGICAL            LSAME
076:       EXTERNAL           LSAME
077: *     ..
078: *     .. External Subroutines ..
079:       EXTERNAL           CGEMV, CGERU, CLACGV, CSSCAL, CSWAP, XERBLA
080: *     ..
081: *     .. Intrinsic Functions ..
082:       INTRINSIC          CONJG, MAX, REAL
083: *     ..
084: *     .. Executable Statements ..
085: *
086:       INFO = 0
087:       UPPER = LSAME( UPLO, 'U' )
088:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
089:          INFO = -1
090:       ELSE IF( N.LT.0 ) THEN
091:          INFO = -2
092:       ELSE IF( NRHS.LT.0 ) THEN
093:          INFO = -3
094:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
095:          INFO = -5
096:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
097:          INFO = -8
098:       END IF
099:       IF( INFO.NE.0 ) THEN
100:          CALL XERBLA( 'CHETRS', -INFO )
101:          RETURN
102:       END IF
103: *
104: *     Quick return if possible
105: *
106:       IF( N.EQ.0 .OR. NRHS.EQ.0 )
107:      $   RETURN
108: *
109:       IF( UPPER ) THEN
110: *
111: *        Solve A*X = B, where A = U*D*U'.
112: *
113: *        First solve U*D*X = B, overwriting B with X.
114: *
115: *        K is the main loop index, decreasing from N to 1 in steps of
116: *        1 or 2, depending on the size of the diagonal blocks.
117: *
118:          K = N
119:    10    CONTINUE
120: *
121: *        If K < 1, exit from loop.
122: *
123:          IF( K.LT.1 )
124:      $      GO TO 30
125: *
126:          IF( IPIV( K ).GT.0 ) THEN
127: *
128: *           1 x 1 diagonal block
129: *
130: *           Interchange rows K and IPIV(K).
131: *
132:             KP = IPIV( K )
133:             IF( KP.NE.K )
134:      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
135: *
136: *           Multiply by inv(U(K)), where U(K) is the transformation
137: *           stored in column K of A.
138: *
139:             CALL CGERU( K-1, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
140:      $                  B( 1, 1 ), LDB )
141: *
142: *           Multiply by the inverse of the diagonal block.
143: *
144:             S = REAL( ONE ) / REAL( A( K, K ) )
145:             CALL CSSCAL( NRHS, S, B( K, 1 ), LDB )
146:             K = K - 1
147:          ELSE
148: *
149: *           2 x 2 diagonal block
150: *
151: *           Interchange rows K-1 and -IPIV(K).
152: *
153:             KP = -IPIV( K )
154:             IF( KP.NE.K-1 )
155:      $         CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
156: *
157: *           Multiply by inv(U(K)), where U(K) is the transformation
158: *           stored in columns K-1 and K of A.
159: *
160:             CALL CGERU( K-2, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
161:      $                  B( 1, 1 ), LDB )
162:             CALL CGERU( K-2, NRHS, -ONE, A( 1, K-1 ), 1, B( K-1, 1 ),
163:      $                  LDB, B( 1, 1 ), LDB )
164: *
165: *           Multiply by the inverse of the diagonal block.
166: *
167:             AKM1K = A( K-1, K )
168:             AKM1 = A( K-1, K-1 ) / AKM1K
169:             AK = A( K, K ) / CONJG( AKM1K )
170:             DENOM = AKM1*AK - ONE
171:             DO 20 J = 1, NRHS
172:                BKM1 = B( K-1, J ) / AKM1K
173:                BK = B( K, J ) / CONJG( AKM1K )
174:                B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
175:                B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
176:    20       CONTINUE
177:             K = K - 2
178:          END IF
179: *
180:          GO TO 10
181:    30    CONTINUE
182: *
183: *        Next solve U'*X = B, overwriting B with X.
184: *
185: *        K is the main loop index, increasing from 1 to N in steps of
186: *        1 or 2, depending on the size of the diagonal blocks.
187: *
188:          K = 1
189:    40    CONTINUE
190: *
191: *        If K > N, exit from loop.
192: *
193:          IF( K.GT.N )
194:      $      GO TO 50
195: *
196:          IF( IPIV( K ).GT.0 ) THEN
197: *
198: *           1 x 1 diagonal block
199: *
200: *           Multiply by inv(U'(K)), where U(K) is the transformation
201: *           stored in column K of A.
202: *
203:             IF( K.GT.1 ) THEN
204:                CALL CLACGV( NRHS, B( K, 1 ), LDB )
205:                CALL CGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
206:      $                     LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
207:                CALL CLACGV( NRHS, B( K, 1 ), LDB )
208:             END IF
209: *
210: *           Interchange rows K and IPIV(K).
211: *
212:             KP = IPIV( K )
213:             IF( KP.NE.K )
214:      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
215:             K = K + 1
216:          ELSE
217: *
218: *           2 x 2 diagonal block
219: *
220: *           Multiply by inv(U'(K+1)), where U(K+1) is the transformation
221: *           stored in columns K and K+1 of A.
222: *
223:             IF( K.GT.1 ) THEN
224:                CALL CLACGV( NRHS, B( K, 1 ), LDB )
225:                CALL CGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
226:      $                     LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
227:                CALL CLACGV( NRHS, B( K, 1 ), LDB )
228: *
229:                CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
230:                CALL CGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
231:      $                     LDB, A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB )
232:                CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
233:             END IF
234: *
235: *           Interchange rows K and -IPIV(K).
236: *
237:             KP = -IPIV( K )
238:             IF( KP.NE.K )
239:      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
240:             K = K + 2
241:          END IF
242: *
243:          GO TO 40
244:    50    CONTINUE
245: *
246:       ELSE
247: *
248: *        Solve A*X = B, where A = L*D*L'.
249: *
250: *        First solve L*D*X = B, overwriting B with X.
251: *
252: *        K is the main loop index, increasing from 1 to N in steps of
253: *        1 or 2, depending on the size of the diagonal blocks.
254: *
255:          K = 1
256:    60    CONTINUE
257: *
258: *        If K > N, exit from loop.
259: *
260:          IF( K.GT.N )
261:      $      GO TO 80
262: *
263:          IF( IPIV( K ).GT.0 ) THEN
264: *
265: *           1 x 1 diagonal block
266: *
267: *           Interchange rows K and IPIV(K).
268: *
269:             KP = IPIV( K )
270:             IF( KP.NE.K )
271:      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
272: *
273: *           Multiply by inv(L(K)), where L(K) is the transformation
274: *           stored in column K of A.
275: *
276:             IF( K.LT.N )
277:      $         CALL CGERU( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ),
278:      $                     LDB, B( K+1, 1 ), LDB )
279: *
280: *           Multiply by the inverse of the diagonal block.
281: *
282:             S = REAL( ONE ) / REAL( A( K, K ) )
283:             CALL CSSCAL( NRHS, S, B( K, 1 ), LDB )
284:             K = K + 1
285:          ELSE
286: *
287: *           2 x 2 diagonal block
288: *
289: *           Interchange rows K+1 and -IPIV(K).
290: *
291:             KP = -IPIV( K )
292:             IF( KP.NE.K+1 )
293:      $         CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
294: *
295: *           Multiply by inv(L(K)), where L(K) is the transformation
296: *           stored in columns K and K+1 of A.
297: *
298:             IF( K.LT.N-1 ) THEN
299:                CALL CGERU( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ),
300:      $                     LDB, B( K+2, 1 ), LDB )
301:                CALL CGERU( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1,
302:      $                     B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
303:             END IF
304: *
305: *           Multiply by the inverse of the diagonal block.
306: *
307:             AKM1K = A( K+1, K )
308:             AKM1 = A( K, K ) / CONJG( AKM1K )
309:             AK = A( K+1, K+1 ) / AKM1K
310:             DENOM = AKM1*AK - ONE
311:             DO 70 J = 1, NRHS
312:                BKM1 = B( K, J ) / CONJG( AKM1K )
313:                BK = B( K+1, J ) / AKM1K
314:                B( K, J ) = ( AK*BKM1-BK ) / DENOM
315:                B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
316:    70       CONTINUE
317:             K = K + 2
318:          END IF
319: *
320:          GO TO 60
321:    80    CONTINUE
322: *
323: *        Next solve L'*X = B, overwriting B with X.
324: *
325: *        K is the main loop index, decreasing from N to 1 in steps of
326: *        1 or 2, depending on the size of the diagonal blocks.
327: *
328:          K = N
329:    90    CONTINUE
330: *
331: *        If K < 1, exit from loop.
332: *
333:          IF( K.LT.1 )
334:      $      GO TO 100
335: *
336:          IF( IPIV( K ).GT.0 ) THEN
337: *
338: *           1 x 1 diagonal block
339: *
340: *           Multiply by inv(L'(K)), where L(K) is the transformation
341: *           stored in column K of A.
342: *
343:             IF( K.LT.N ) THEN
344:                CALL CLACGV( NRHS, B( K, 1 ), LDB )
345:                CALL CGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
346:      $                     B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
347:      $                     B( K, 1 ), LDB )
348:                CALL CLACGV( NRHS, B( K, 1 ), LDB )
349:             END IF
350: *
351: *           Interchange rows K and IPIV(K).
352: *
353:             KP = IPIV( K )
354:             IF( KP.NE.K )
355:      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
356:             K = K - 1
357:          ELSE
358: *
359: *           2 x 2 diagonal block
360: *
361: *           Multiply by inv(L'(K-1)), where L(K-1) is the transformation
362: *           stored in columns K-1 and K of A.
363: *
364:             IF( K.LT.N ) THEN
365:                CALL CLACGV( NRHS, B( K, 1 ), LDB )
366:                CALL CGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
367:      $                     B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
368:      $                     B( K, 1 ), LDB )
369:                CALL CLACGV( NRHS, B( K, 1 ), LDB )
370: *
371:                CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
372:                CALL CGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
373:      $                     B( K+1, 1 ), LDB, A( K+1, K-1 ), 1, ONE,
374:      $                     B( K-1, 1 ), LDB )
375:                CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
376:             END IF
377: *
378: *           Interchange rows K and -IPIV(K).
379: *
380:             KP = -IPIV( K )
381:             IF( KP.NE.K )
382:      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
383:             K = K - 2
384:          END IF
385: *
386:          GO TO 90
387:   100    CONTINUE
388:       END IF
389: *
390:       RETURN
391: *
392: *     End of CHETRS
393: *
394:       END
395: