LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
chptrs.f
Go to the documentation of this file.
1 *> \brief \b CHPTRS
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CHPTRS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chptrs.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chptrs.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chptrs.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CHPTRS( UPLO, N, NRHS, AP, IPIV, B, LDB, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER UPLO
25 * INTEGER INFO, LDB, N, NRHS
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IPIV( * )
29 * COMPLEX AP( * ), B( LDB, * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> CHPTRS solves a system of linear equations A*X = B with a complex
39 *> Hermitian matrix A stored in packed format using the factorization
40 *> A = U*D*U**H or A = L*D*L**H computed by CHPTRF.
41 *> \endverbatim
42 *
43 * Arguments:
44 * ==========
45 *
46 *> \param[in] UPLO
47 *> \verbatim
48 *> UPLO is CHARACTER*1
49 *> Specifies whether the details of the factorization are stored
50 *> as an upper or lower triangular matrix.
51 *> = 'U': Upper triangular, form is A = U*D*U**H;
52 *> = 'L': Lower triangular, form is A = L*D*L**H.
53 *> \endverbatim
54 *>
55 *> \param[in] N
56 *> \verbatim
57 *> N is INTEGER
58 *> The order of the matrix A. N >= 0.
59 *> \endverbatim
60 *>
61 *> \param[in] NRHS
62 *> \verbatim
63 *> NRHS is INTEGER
64 *> The number of right hand sides, i.e., the number of columns
65 *> of the matrix B. NRHS >= 0.
66 *> \endverbatim
67 *>
68 *> \param[in] AP
69 *> \verbatim
70 *> AP is COMPLEX array, dimension (N*(N+1)/2)
71 *> The block diagonal matrix D and the multipliers used to
72 *> obtain the factor U or L as computed by CHPTRF, stored as a
73 *> packed triangular matrix.
74 *> \endverbatim
75 *>
76 *> \param[in] IPIV
77 *> \verbatim
78 *> IPIV is INTEGER array, dimension (N)
79 *> Details of the interchanges and the block structure of D
80 *> as determined by CHPTRF.
81 *> \endverbatim
82 *>
83 *> \param[in,out] B
84 *> \verbatim
85 *> B is COMPLEX array, dimension (LDB,NRHS)
86 *> On entry, the right hand side matrix B.
87 *> On exit, the solution matrix X.
88 *> \endverbatim
89 *>
90 *> \param[in] LDB
91 *> \verbatim
92 *> LDB is INTEGER
93 *> The leading dimension of the array B. LDB >= max(1,N).
94 *> \endverbatim
95 *>
96 *> \param[out] INFO
97 *> \verbatim
98 *> INFO is INTEGER
99 *> = 0: successful exit
100 *> < 0: if INFO = -i, the i-th argument had an illegal value
101 *> \endverbatim
102 *
103 * Authors:
104 * ========
105 *
106 *> \author Univ. of Tennessee
107 *> \author Univ. of California Berkeley
108 *> \author Univ. of Colorado Denver
109 *> \author NAG Ltd.
110 *
111 *> \ingroup complexOTHERcomputational
112 *
113 * =====================================================================
114  SUBROUTINE chptrs( UPLO, N, NRHS, AP, IPIV, B, LDB, INFO )
115 *
116 * -- LAPACK computational routine --
117 * -- LAPACK is a software package provided by Univ. of Tennessee, --
118 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
119 *
120 * .. Scalar Arguments ..
121  CHARACTER UPLO
122  INTEGER INFO, LDB, N, NRHS
123 * ..
124 * .. Array Arguments ..
125  INTEGER IPIV( * )
126  COMPLEX AP( * ), B( LDB, * )
127 * ..
128 *
129 * =====================================================================
130 *
131 * .. Parameters ..
132  COMPLEX ONE
133  parameter( one = ( 1.0e+0, 0.0e+0 ) )
134 * ..
135 * .. Local Scalars ..
136  LOGICAL UPPER
137  INTEGER J, K, KC, KP
138  REAL S
139  COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM
140 * ..
141 * .. External Functions ..
142  LOGICAL LSAME
143  EXTERNAL lsame
144 * ..
145 * .. External Subroutines ..
146  EXTERNAL cgemv, cgeru, clacgv, csscal, cswap, xerbla
147 * ..
148 * .. Intrinsic Functions ..
149  INTRINSIC conjg, max, real
150 * ..
151 * .. Executable Statements ..
152 *
153  info = 0
154  upper = lsame( uplo, 'U' )
155  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
156  info = -1
157  ELSE IF( n.LT.0 ) THEN
158  info = -2
159  ELSE IF( nrhs.LT.0 ) THEN
160  info = -3
161  ELSE IF( ldb.LT.max( 1, n ) ) THEN
162  info = -7
163  END IF
164  IF( info.NE.0 ) THEN
165  CALL xerbla( 'CHPTRS', -info )
166  RETURN
167  END IF
168 *
169 * Quick return if possible
170 *
171  IF( n.EQ.0 .OR. nrhs.EQ.0 )
172  $ RETURN
173 *
174  IF( upper ) THEN
175 *
176 * Solve A*X = B, where A = U*D*U**H.
177 *
178 * First solve U*D*X = B, overwriting B with X.
179 *
180 * K is the main loop index, decreasing from N to 1 in steps of
181 * 1 or 2, depending on the size of the diagonal blocks.
182 *
183  k = n
184  kc = n*( n+1 ) / 2 + 1
185  10 CONTINUE
186 *
187 * If K < 1, exit from loop.
188 *
189  IF( k.LT.1 )
190  $ GO TO 30
191 *
192  kc = kc - k
193  IF( ipiv( k ).GT.0 ) THEN
194 *
195 * 1 x 1 diagonal block
196 *
197 * Interchange rows K and IPIV(K).
198 *
199  kp = ipiv( k )
200  IF( kp.NE.k )
201  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
202 *
203 * Multiply by inv(U(K)), where U(K) is the transformation
204 * stored in column K of A.
205 *
206  CALL cgeru( k-1, nrhs, -one, ap( kc ), 1, b( k, 1 ), ldb,
207  $ b( 1, 1 ), ldb )
208 *
209 * Multiply by the inverse of the diagonal block.
210 *
211  s = real( one ) / real( ap( kc+k-1 ) )
212  CALL csscal( nrhs, s, b( k, 1 ), ldb )
213  k = k - 1
214  ELSE
215 *
216 * 2 x 2 diagonal block
217 *
218 * Interchange rows K-1 and -IPIV(K).
219 *
220  kp = -ipiv( k )
221  IF( kp.NE.k-1 )
222  $ CALL cswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
223 *
224 * Multiply by inv(U(K)), where U(K) is the transformation
225 * stored in columns K-1 and K of A.
226 *
227  CALL cgeru( k-2, nrhs, -one, ap( kc ), 1, b( k, 1 ), ldb,
228  $ b( 1, 1 ), ldb )
229  CALL cgeru( k-2, nrhs, -one, ap( kc-( k-1 ) ), 1,
230  $ b( k-1, 1 ), ldb, b( 1, 1 ), ldb )
231 *
232 * Multiply by the inverse of the diagonal block.
233 *
234  akm1k = ap( kc+k-2 )
235  akm1 = ap( kc-1 ) / akm1k
236  ak = ap( kc+k-1 ) / conjg( akm1k )
237  denom = akm1*ak - one
238  DO 20 j = 1, nrhs
239  bkm1 = b( k-1, j ) / akm1k
240  bk = b( k, j ) / conjg( akm1k )
241  b( k-1, j ) = ( ak*bkm1-bk ) / denom
242  b( k, j ) = ( akm1*bk-bkm1 ) / denom
243  20 CONTINUE
244  kc = kc - k + 1
245  k = k - 2
246  END IF
247 *
248  GO TO 10
249  30 CONTINUE
250 *
251 * Next solve U**H *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  40 CONTINUE
259 *
260 * If K > N, exit from loop.
261 *
262  IF( k.GT.n )
263  $ GO TO 50
264 *
265  IF( ipiv( k ).GT.0 ) THEN
266 *
267 * 1 x 1 diagonal block
268 *
269 * Multiply by inv(U**H(K)), where U(K) is the transformation
270 * stored in column K of A.
271 *
272  IF( k.GT.1 ) THEN
273  CALL clacgv( nrhs, b( k, 1 ), ldb )
274  CALL cgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
275  $ ldb, ap( kc ), 1, one, b( k, 1 ), ldb )
276  CALL clacgv( nrhs, b( k, 1 ), ldb )
277  END IF
278 *
279 * Interchange rows K and IPIV(K).
280 *
281  kp = ipiv( k )
282  IF( kp.NE.k )
283  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
284  kc = kc + k
285  k = k + 1
286  ELSE
287 *
288 * 2 x 2 diagonal block
289 *
290 * Multiply by inv(U**H(K+1)), where U(K+1) is the transformation
291 * stored in columns K and K+1 of A.
292 *
293  IF( k.GT.1 ) THEN
294  CALL clacgv( nrhs, b( k, 1 ), ldb )
295  CALL cgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
296  $ ldb, ap( kc ), 1, one, b( k, 1 ), ldb )
297  CALL clacgv( nrhs, b( k, 1 ), ldb )
298 *
299  CALL clacgv( nrhs, b( k+1, 1 ), ldb )
300  CALL cgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
301  $ ldb, ap( kc+k ), 1, one, b( k+1, 1 ), ldb )
302  CALL clacgv( nrhs, b( k+1, 1 ), ldb )
303  END IF
304 *
305 * Interchange rows K and -IPIV(K).
306 *
307  kp = -ipiv( k )
308  IF( kp.NE.k )
309  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
310  kc = kc + 2*k + 1
311  k = k + 2
312  END IF
313 *
314  GO TO 40
315  50 CONTINUE
316 *
317  ELSE
318 *
319 * Solve A*X = B, where A = L*D*L**H.
320 *
321 * First solve L*D*X = B, overwriting B with X.
322 *
323 * K is the main loop index, increasing from 1 to N in steps of
324 * 1 or 2, depending on the size of the diagonal blocks.
325 *
326  k = 1
327  kc = 1
328  60 CONTINUE
329 *
330 * If K > N, exit from loop.
331 *
332  IF( k.GT.n )
333  $ GO TO 80
334 *
335  IF( ipiv( k ).GT.0 ) THEN
336 *
337 * 1 x 1 diagonal block
338 *
339 * Interchange rows K and IPIV(K).
340 *
341  kp = ipiv( k )
342  IF( kp.NE.k )
343  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
344 *
345 * Multiply by inv(L(K)), where L(K) is the transformation
346 * stored in column K of A.
347 *
348  IF( k.LT.n )
349  $ CALL cgeru( n-k, nrhs, -one, ap( kc+1 ), 1, b( k, 1 ),
350  $ ldb, b( k+1, 1 ), ldb )
351 *
352 * Multiply by the inverse of the diagonal block.
353 *
354  s = real( one ) / real( ap( kc ) )
355  CALL csscal( nrhs, s, b( k, 1 ), ldb )
356  kc = kc + n - k + 1
357  k = k + 1
358  ELSE
359 *
360 * 2 x 2 diagonal block
361 *
362 * Interchange rows K+1 and -IPIV(K).
363 *
364  kp = -ipiv( k )
365  IF( kp.NE.k+1 )
366  $ CALL cswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
367 *
368 * Multiply by inv(L(K)), where L(K) is the transformation
369 * stored in columns K and K+1 of A.
370 *
371  IF( k.LT.n-1 ) THEN
372  CALL cgeru( n-k-1, nrhs, -one, ap( kc+2 ), 1, b( k, 1 ),
373  $ ldb, b( k+2, 1 ), ldb )
374  CALL cgeru( n-k-1, nrhs, -one, ap( kc+n-k+2 ), 1,
375  $ b( k+1, 1 ), ldb, b( k+2, 1 ), ldb )
376  END IF
377 *
378 * Multiply by the inverse of the diagonal block.
379 *
380  akm1k = ap( kc+1 )
381  akm1 = ap( kc ) / conjg( akm1k )
382  ak = ap( kc+n-k+1 ) / akm1k
383  denom = akm1*ak - one
384  DO 70 j = 1, nrhs
385  bkm1 = b( k, j ) / conjg( akm1k )
386  bk = b( k+1, j ) / akm1k
387  b( k, j ) = ( ak*bkm1-bk ) / denom
388  b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
389  70 CONTINUE
390  kc = kc + 2*( n-k ) + 1
391  k = k + 2
392  END IF
393 *
394  GO TO 60
395  80 CONTINUE
396 *
397 * Next solve L**H *X = B, overwriting B with X.
398 *
399 * K is the main loop index, decreasing from N to 1 in steps of
400 * 1 or 2, depending on the size of the diagonal blocks.
401 *
402  k = n
403  kc = n*( n+1 ) / 2 + 1
404  90 CONTINUE
405 *
406 * If K < 1, exit from loop.
407 *
408  IF( k.LT.1 )
409  $ GO TO 100
410 *
411  kc = kc - ( n-k+1 )
412  IF( ipiv( k ).GT.0 ) THEN
413 *
414 * 1 x 1 diagonal block
415 *
416 * Multiply by inv(L**H(K)), where L(K) is the transformation
417 * stored in column K of A.
418 *
419  IF( k.LT.n ) THEN
420  CALL clacgv( nrhs, b( k, 1 ), ldb )
421  CALL cgemv( 'Conjugate transpose', n-k, nrhs, -one,
422  $ b( k+1, 1 ), ldb, ap( kc+1 ), 1, one,
423  $ b( k, 1 ), ldb )
424  CALL clacgv( nrhs, b( k, 1 ), ldb )
425  END IF
426 *
427 * Interchange rows K and IPIV(K).
428 *
429  kp = ipiv( k )
430  IF( kp.NE.k )
431  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
432  k = k - 1
433  ELSE
434 *
435 * 2 x 2 diagonal block
436 *
437 * Multiply by inv(L**H(K-1)), where L(K-1) is the transformation
438 * stored in columns K-1 and K of A.
439 *
440  IF( k.LT.n ) THEN
441  CALL clacgv( nrhs, b( k, 1 ), ldb )
442  CALL cgemv( 'Conjugate transpose', n-k, nrhs, -one,
443  $ b( k+1, 1 ), ldb, ap( kc+1 ), 1, one,
444  $ b( k, 1 ), ldb )
445  CALL clacgv( nrhs, b( k, 1 ), ldb )
446 *
447  CALL clacgv( nrhs, b( k-1, 1 ), ldb )
448  CALL cgemv( 'Conjugate transpose', n-k, nrhs, -one,
449  $ b( k+1, 1 ), ldb, ap( kc-( n-k ) ), 1, one,
450  $ b( k-1, 1 ), ldb )
451  CALL clacgv( nrhs, b( k-1, 1 ), ldb )
452  END IF
453 *
454 * Interchange rows K and -IPIV(K).
455 *
456  kp = -ipiv( k )
457  IF( kp.NE.k )
458  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
459  kc = kc - ( n-k+2 )
460  k = k - 2
461  END IF
462 *
463  GO TO 90
464  100 CONTINUE
465  END IF
466 *
467  RETURN
468 *
469 * End of CHPTRS
470 *
471  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine cgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CGERU
Definition: cgeru.f:130
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:74
subroutine chptrs(UPLO, N, NRHS, AP, IPIV, B, LDB, INFO)
CHPTRS
Definition: chptrs.f:115