LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 hptrs
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)
Definition cblat2.f:3285
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
subroutine cgeru(m, n, alpha, x, incx, y, incy, a, lda)
CGERU
Definition cgeru.f:130
subroutine chptrs(uplo, n, nrhs, ap, ipiv, b, ldb, info)
CHPTRS
Definition chptrs.f:115
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:74
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81