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