001:       SUBROUTINE ZHPTRS( 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:       COMPLEX*16         AP( * ), B( LDB, * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  ZHPTRS solves a system of linear equations A*X = B with a complex
021: *  Hermitian matrix A stored in packed format using the factorization
022: *  A = U*D*U**H or A = L*D*L**H computed by ZHPTRF.
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: *  AP      (input) COMPLEX*16 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 ZHPTRF, 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 ZHPTRF.
048: *
049: *  B       (input/output) COMPLEX*16 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:       COMPLEX*16         ONE
064:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
065: *     ..
066: *     .. Local Scalars ..
067:       LOGICAL            UPPER
068:       INTEGER            J, K, KC, KP
069:       DOUBLE PRECISION   S
070:       COMPLEX*16         AK, AKM1, AKM1K, BK, BKM1, DENOM
071: *     ..
072: *     .. External Functions ..
073:       LOGICAL            LSAME
074:       EXTERNAL           LSAME
075: *     ..
076: *     .. External Subroutines ..
077:       EXTERNAL           XERBLA, ZDSCAL, ZGEMV, ZGERU, ZLACGV, ZSWAP
078: *     ..
079: *     .. Intrinsic Functions ..
080:       INTRINSIC          DBLE, DCONJG, 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( LDB.LT.MAX( 1, N ) ) THEN
093:          INFO = -7
094:       END IF
095:       IF( INFO.NE.0 ) THEN
096:          CALL XERBLA( 'ZHPTRS', -INFO )
097:          RETURN
098:       END IF
099: *
100: *     Quick return if possible
101: *
102:       IF( N.EQ.0 .OR. NRHS.EQ.0 )
103:      $   RETURN
104: *
105:       IF( UPPER ) THEN
106: *
107: *        Solve A*X = B, where A = U*D*U'.
108: *
109: *        First solve U*D*X = B, overwriting B with X.
110: *
111: *        K is the main loop index, decreasing from N to 1 in steps of
112: *        1 or 2, depending on the size of the diagonal blocks.
113: *
114:          K = N
115:          KC = N*( N+1 ) / 2 + 1
116:    10    CONTINUE
117: *
118: *        If K < 1, exit from loop.
119: *
120:          IF( K.LT.1 )
121:      $      GO TO 30
122: *
123:          KC = KC - K
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 ZSWAP( 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 ZGERU( K-1, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
138:      $                  B( 1, 1 ), LDB )
139: *
140: *           Multiply by the inverse of the diagonal block.
141: *
142:             S = DBLE( ONE ) / DBLE( AP( KC+K-1 ) )
143:             CALL ZDSCAL( NRHS, S, 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 ZSWAP( 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 ZGERU( K-2, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
159:      $                  B( 1, 1 ), LDB )
160:             CALL ZGERU( K-2, NRHS, -ONE, AP( KC-( K-1 ) ), 1,
161:      $                  B( K-1, 1 ), LDB, B( 1, 1 ), LDB )
162: *
163: *           Multiply by the inverse of the diagonal block.
164: *
165:             AKM1K = AP( KC+K-2 )
166:             AKM1 = AP( KC-1 ) / AKM1K
167:             AK = AP( KC+K-1 ) / DCONJG( AKM1K )
168:             DENOM = AKM1*AK - ONE
169:             DO 20 J = 1, NRHS
170:                BKM1 = B( K-1, J ) / AKM1K
171:                BK = B( K, J ) / DCONJG( AKM1K )
172:                B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
173:                B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
174:    20       CONTINUE
175:             KC = KC - K + 1
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:          KC = 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 ZLACGV( NRHS, B( K, 1 ), LDB )
205:                CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
206:      $                     LDB, AP( KC ), 1, ONE, B( K, 1 ), LDB )
207:                CALL ZLACGV( 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 ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
215:             KC = KC + K
216:             K = K + 1
217:          ELSE
218: *
219: *           2 x 2 diagonal block
220: *
221: *           Multiply by inv(U'(K+1)), where U(K+1) is the transformation
222: *           stored in columns K and K+1 of A.
223: *
224:             IF( K.GT.1 ) THEN
225:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
226:                CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
227:      $                     LDB, AP( KC ), 1, ONE, B( K, 1 ), LDB )
228:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
229: *
230:                CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
231:                CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
232:      $                     LDB, AP( KC+K ), 1, ONE, B( K+1, 1 ), LDB )
233:                CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
234:             END IF
235: *
236: *           Interchange rows K and -IPIV(K).
237: *
238:             KP = -IPIV( K )
239:             IF( KP.NE.K )
240:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
241:             KC = KC + 2*K + 1
242:             K = K + 2
243:          END IF
244: *
245:          GO TO 40
246:    50    CONTINUE
247: *
248:       ELSE
249: *
250: *        Solve A*X = B, where A = L*D*L'.
251: *
252: *        First solve L*D*X = B, overwriting B with X.
253: *
254: *        K is the main loop index, increasing from 1 to N in steps of
255: *        1 or 2, depending on the size of the diagonal blocks.
256: *
257:          K = 1
258:          KC = 1
259:    60    CONTINUE
260: *
261: *        If K > N, exit from loop.
262: *
263:          IF( K.GT.N )
264:      $      GO TO 80
265: *
266:          IF( IPIV( K ).GT.0 ) THEN
267: *
268: *           1 x 1 diagonal block
269: *
270: *           Interchange rows K and IPIV(K).
271: *
272:             KP = IPIV( K )
273:             IF( KP.NE.K )
274:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
275: *
276: *           Multiply by inv(L(K)), where L(K) is the transformation
277: *           stored in column K of A.
278: *
279:             IF( K.LT.N )
280:      $         CALL ZGERU( N-K, NRHS, -ONE, AP( KC+1 ), 1, B( K, 1 ),
281:      $                     LDB, B( K+1, 1 ), LDB )
282: *
283: *           Multiply by the inverse of the diagonal block.
284: *
285:             S = DBLE( ONE ) / DBLE( AP( KC ) )
286:             CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
287:             KC = KC + N - K + 1
288:             K = K + 1
289:          ELSE
290: *
291: *           2 x 2 diagonal block
292: *
293: *           Interchange rows K+1 and -IPIV(K).
294: *
295:             KP = -IPIV( K )
296:             IF( KP.NE.K+1 )
297:      $         CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
298: *
299: *           Multiply by inv(L(K)), where L(K) is the transformation
300: *           stored in columns K and K+1 of A.
301: *
302:             IF( K.LT.N-1 ) THEN
303:                CALL ZGERU( N-K-1, NRHS, -ONE, AP( KC+2 ), 1, B( K, 1 ),
304:      $                     LDB, B( K+2, 1 ), LDB )
305:                CALL ZGERU( N-K-1, NRHS, -ONE, AP( KC+N-K+2 ), 1,
306:      $                     B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
307:             END IF
308: *
309: *           Multiply by the inverse of the diagonal block.
310: *
311:             AKM1K = AP( KC+1 )
312:             AKM1 = AP( KC ) / DCONJG( AKM1K )
313:             AK = AP( KC+N-K+1 ) / AKM1K
314:             DENOM = AKM1*AK - ONE
315:             DO 70 J = 1, NRHS
316:                BKM1 = B( K, J ) / DCONJG( AKM1K )
317:                BK = B( K+1, J ) / AKM1K
318:                B( K, J ) = ( AK*BKM1-BK ) / DENOM
319:                B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
320:    70       CONTINUE
321:             KC = KC + 2*( N-K ) + 1
322:             K = K + 2
323:          END IF
324: *
325:          GO TO 60
326:    80    CONTINUE
327: *
328: *        Next solve L'*X = B, overwriting B with X.
329: *
330: *        K is the main loop index, decreasing from N to 1 in steps of
331: *        1 or 2, depending on the size of the diagonal blocks.
332: *
333:          K = N
334:          KC = N*( N+1 ) / 2 + 1
335:    90    CONTINUE
336: *
337: *        If K < 1, exit from loop.
338: *
339:          IF( K.LT.1 )
340:      $      GO TO 100
341: *
342:          KC = KC - ( N-K+1 )
343:          IF( IPIV( K ).GT.0 ) THEN
344: *
345: *           1 x 1 diagonal block
346: *
347: *           Multiply by inv(L'(K)), where L(K) is the transformation
348: *           stored in column K of A.
349: *
350:             IF( K.LT.N ) THEN
351:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
352:                CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
353:      $                     B( K+1, 1 ), LDB, AP( KC+1 ), 1, ONE,
354:      $                     B( K, 1 ), LDB )
355:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
356:             END IF
357: *
358: *           Interchange rows K and IPIV(K).
359: *
360:             KP = IPIV( K )
361:             IF( KP.NE.K )
362:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
363:             K = K - 1
364:          ELSE
365: *
366: *           2 x 2 diagonal block
367: *
368: *           Multiply by inv(L'(K-1)), where L(K-1) is the transformation
369: *           stored in columns K-1 and K of A.
370: *
371:             IF( K.LT.N ) THEN
372:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
373:                CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
374:      $                     B( K+1, 1 ), LDB, AP( KC+1 ), 1, ONE,
375:      $                     B( K, 1 ), LDB )
376:                CALL ZLACGV( NRHS, B( K, 1 ), LDB )
377: *
378:                CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
379:                CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
380:      $                     B( K+1, 1 ), LDB, AP( KC-( N-K ) ), 1, ONE,
381:      $                     B( K-1, 1 ), LDB )
382:                CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
383:             END IF
384: *
385: *           Interchange rows K and -IPIV(K).
386: *
387:             KP = -IPIV( K )
388:             IF( KP.NE.K )
389:      $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
390:             KC = KC - ( N-K+2 )
391:             K = K - 2
392:          END IF
393: *
394:          GO TO 90
395:   100    CONTINUE
396:       END IF
397: *
398:       RETURN
399: *
400: *     End of ZHPTRS
401: *
402:       END
403: