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