001:       SUBROUTINE ZHETRS( 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:       COMPLEX*16         A( LDA, * ), B( LDB, * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  ZHETRS solves a system of linear equations A*X = B with a complex
020: *  Hermitian matrix A using the factorization A = U*D*U**H or
021: *  A = L*D*L**H computed by ZHETRF.
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**H;
030: *          = 'L':  Lower triangular, form is A = L*D*L**H.
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) COMPLEX*16 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 ZHETRF.
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 ZHETRF.
049: *
050: *  B       (input/output) COMPLEX*16 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:       COMPLEX*16         ONE
065:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
066: *     ..
067: *     .. Local Scalars ..
068:       LOGICAL            UPPER
069:       INTEGER            J, K, KP
070:       DOUBLE PRECISION   S
071:       COMPLEX*16         AK, AKM1, AKM1K, BK, BKM1, DENOM
072: *     ..
073: *     .. External Functions ..
074:       LOGICAL            LSAME
075:       EXTERNAL           LSAME
076: *     ..
077: *     .. External Subroutines ..
078:       EXTERNAL           XERBLA, ZDSCAL, ZGEMV, ZGERU, ZLACGV, ZSWAP
079: *     ..
080: *     .. Intrinsic Functions ..
081:       INTRINSIC          DBLE, DCONJG, 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( 'ZHETRS', -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 ZSWAP( 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 ZGERU( 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:             S = DBLE( ONE ) / DBLE( A( K, K ) )
144:             CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
145:             K = K - 1
146:          ELSE
147: *
148: *           2 x 2 diagonal block
149: *
150: *           Interchange rows K-1 and -IPIV(K).
151: *
152:             KP = -IPIV( K )
153:             IF( KP.NE.K-1 )
154:      $         CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
155: *
156: *           Multiply by inv(U(K)), where U(K) is the transformation
157: *           stored in columns K-1 and K of A.
158: *
159:             CALL ZGERU( K-2, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
160:      $                  B( 1, 1 ), LDB )
161:             CALL ZGERU( K-2, NRHS, -ONE, A( 1, K-1 ), 1, B( K-1, 1 ),
162:      $                  LDB, B( 1, 1 ), LDB )
163: *
164: *           Multiply by the inverse of the diagonal block.
165: *
166:             AKM1K = A( K-1, K )
167:             AKM1 = A( K-1, K-1 ) / AKM1K
168:             AK = A( K, K ) / DCONJG( AKM1K )
169:             DENOM = AKM1*AK - ONE
170:             DO 20 J = 1, NRHS
171:                BKM1 = B( K-1, J ) / AKM1K
172:                BK = B( K, J ) / DCONJG( AKM1K )
173:                B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
174:                B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
175:    20       CONTINUE
176:             K = K - 2
177:          END IF
178: *
179:          GO TO 10
180:    30    CONTINUE
181: *
182: *        Next solve U'*X = B, overwriting B with X.
183: *
184: *        K is the main loop index, increasing from 1 to N in steps of
185: *        1 or 2, depending on the size of the diagonal blocks.
186: *
187:          K = 1
188:    40    CONTINUE
189: *
190: *        If K > N, exit from loop.
191: *
192:          IF( K.GT.N )
193:      $      GO TO 50
194: *
195:          IF( IPIV( K ).GT.0 ) THEN
196: *
197: *           1 x 1 diagonal block
198: *
199: *           Multiply by inv(U'(K)), where U(K) is the transformation
200: *           stored in column K of A.
201: *
202:             IF( K.GT.1 ) THEN
203:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
204:                CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
205:      $                     LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
206:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
207:             END IF
208: *
209: *           Interchange rows K and IPIV(K).
210: *
211:             KP = IPIV( K )
212:             IF( KP.NE.K )
213:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
214:             K = K + 1
215:          ELSE
216: *
217: *           2 x 2 diagonal block
218: *
219: *           Multiply by inv(U'(K+1)), where U(K+1) is the transformation
220: *           stored in columns K and K+1 of A.
221: *
222:             IF( K.GT.1 ) THEN
223:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
224:                CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
225:      $                     LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
226:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
227: *
228:                CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
229:                CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
230:      $                     LDB, A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB )
231:                CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
232:             END IF
233: *
234: *           Interchange rows K and -IPIV(K).
235: *
236:             KP = -IPIV( K )
237:             IF( KP.NE.K )
238:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
239:             K = K + 2
240:          END IF
241: *
242:          GO TO 40
243:    50    CONTINUE
244: *
245:       ELSE
246: *
247: *        Solve A*X = B, where A = L*D*L'.
248: *
249: *        First solve L*D*X = B, overwriting B with X.
250: *
251: *        K is the main loop index, increasing from 1 to N in steps of
252: *        1 or 2, depending on the size of the diagonal blocks.
253: *
254:          K = 1
255:    60    CONTINUE
256: *
257: *        If K > N, exit from loop.
258: *
259:          IF( K.GT.N )
260:      $      GO TO 80
261: *
262:          IF( IPIV( K ).GT.0 ) THEN
263: *
264: *           1 x 1 diagonal block
265: *
266: *           Interchange rows K and IPIV(K).
267: *
268:             KP = IPIV( K )
269:             IF( KP.NE.K )
270:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
271: *
272: *           Multiply by inv(L(K)), where L(K) is the transformation
273: *           stored in column K of A.
274: *
275:             IF( K.LT.N )
276:      $         CALL ZGERU( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ),
277:      $                     LDB, B( K+1, 1 ), LDB )
278: *
279: *           Multiply by the inverse of the diagonal block.
280: *
281:             S = DBLE( ONE ) / DBLE( A( K, K ) )
282:             CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
283:             K = K + 1
284:          ELSE
285: *
286: *           2 x 2 diagonal block
287: *
288: *           Interchange rows K+1 and -IPIV(K).
289: *
290:             KP = -IPIV( K )
291:             IF( KP.NE.K+1 )
292:      $         CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
293: *
294: *           Multiply by inv(L(K)), where L(K) is the transformation
295: *           stored in columns K and K+1 of A.
296: *
297:             IF( K.LT.N-1 ) THEN
298:                CALL ZGERU( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ),
299:      $                     LDB, B( K+2, 1 ), LDB )
300:                CALL ZGERU( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1,
301:      $                     B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
302:             END IF
303: *
304: *           Multiply by the inverse of the diagonal block.
305: *
306:             AKM1K = A( K+1, K )
307:             AKM1 = A( K, K ) / DCONJG( AKM1K )
308:             AK = A( K+1, K+1 ) / AKM1K
309:             DENOM = AKM1*AK - ONE
310:             DO 70 J = 1, NRHS
311:                BKM1 = B( K, J ) / DCONJG( AKM1K )
312:                BK = B( K+1, J ) / AKM1K
313:                B( K, J ) = ( AK*BKM1-BK ) / DENOM
314:                B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
315:    70       CONTINUE
316:             K = K + 2
317:          END IF
318: *
319:          GO TO 60
320:    80    CONTINUE
321: *
322: *        Next solve L'*X = B, overwriting B with X.
323: *
324: *        K is the main loop index, decreasing from N to 1 in steps of
325: *        1 or 2, depending on the size of the diagonal blocks.
326: *
327:          K = N
328:    90    CONTINUE
329: *
330: *        If K < 1, exit from loop.
331: *
332:          IF( K.LT.1 )
333:      $      GO TO 100
334: *
335:          IF( IPIV( K ).GT.0 ) THEN
336: *
337: *           1 x 1 diagonal block
338: *
339: *           Multiply by inv(L'(K)), where L(K) is the transformation
340: *           stored in column K of A.
341: *
342:             IF( K.LT.N ) THEN
343:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
344:                CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
345:      $                     B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
346:      $                     B( K, 1 ), LDB )
347:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
348:             END IF
349: *
350: *           Interchange rows K and IPIV(K).
351: *
352:             KP = IPIV( K )
353:             IF( KP.NE.K )
354:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
355:             K = K - 1
356:          ELSE
357: *
358: *           2 x 2 diagonal block
359: *
360: *           Multiply by inv(L'(K-1)), where L(K-1) is the transformation
361: *           stored in columns K-1 and K of A.
362: *
363:             IF( K.LT.N ) THEN
364:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
365:                CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
366:      $                     B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
367:      $                     B( K, 1 ), LDB )
368:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
369: *
370:                CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
371:                CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
372:      $                     B( K+1, 1 ), LDB, A( K+1, K-1 ), 1, ONE,
373:      $                     B( K-1, 1 ), LDB )
374:                CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
375:             END IF
376: *
377: *           Interchange rows K and -IPIV(K).
378: *
379:             KP = -IPIV( K )
380:             IF( KP.NE.K )
381:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
382:             K = K - 2
383:          END IF
384: *
385:          GO TO 90
386:   100    CONTINUE
387:       END IF
388: *
389:       RETURN
390: *
391: *     End of ZHETRS
392: *
393:       END
394: