LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
chetrs_rook.f
Go to the documentation of this file.
1*> \brief \b CHETRS_ROOK computes the solution to a system of linear equations A * X = B for HE matrices using factorization obtained with one of the bounded diagonal pivoting methods (max 2 interchanges)
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CHETRS_ROOK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetrs_rook.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetrs_rook.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetrs_rook.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CHETRS_ROOK( 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_ROOK 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_ROOK.
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_ROOK.
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_ROOK.
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_rook
117*
118*> \par Contributors:
119* ==================
120*>
121*> \verbatim
122*>
123*> November 2013, Igor Kozachenko,
124*> Computer Science Division,
125*> University of California, Berkeley
126*>
127*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
128*> School of Mathematics,
129*> University of Manchester
130*>
131*> \endverbatim
132*
133* =====================================================================
134 SUBROUTINE chetrs_rook( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
135 $ INFO )
136*
137* -- LAPACK computational routine --
138* -- LAPACK is a software package provided by Univ. of Tennessee, --
139* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140*
141* .. Scalar Arguments ..
142 CHARACTER UPLO
143 INTEGER INFO, LDA, LDB, N, NRHS
144* ..
145* .. Array Arguments ..
146 INTEGER IPIV( * )
147 COMPLEX A( LDA, * ), B( LDB, * )
148* ..
149*
150* =====================================================================
151*
152* .. Parameters ..
153 COMPLEX ONE
154 parameter( one = ( 1.0e+0, 0.0e+0 ) )
155* ..
156* .. Local Scalars ..
157 LOGICAL UPPER
158 INTEGER J, K, KP
159 REAL S
160 COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM
161* ..
162* .. External Functions ..
163 LOGICAL LSAME
164 EXTERNAL lsame
165* ..
166* .. External Subroutines ..
167 EXTERNAL cgemv, cgeru, clacgv, csscal, cswap, xerbla
168* ..
169* .. Intrinsic Functions ..
170 INTRINSIC conjg, max, real
171* ..
172* .. Executable Statements ..
173*
174 info = 0
175 upper = lsame( uplo, 'U' )
176 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
177 info = -1
178 ELSE IF( n.LT.0 ) THEN
179 info = -2
180 ELSE IF( nrhs.LT.0 ) THEN
181 info = -3
182 ELSE IF( lda.LT.max( 1, n ) ) THEN
183 info = -5
184 ELSE IF( ldb.LT.max( 1, n ) ) THEN
185 info = -8
186 END IF
187 IF( info.NE.0 ) THEN
188 CALL xerbla( 'CHETRS_ROOK', -info )
189 RETURN
190 END IF
191*
192* Quick return if possible
193*
194 IF( n.EQ.0 .OR. nrhs.EQ.0 )
195 $ RETURN
196*
197 IF( upper ) THEN
198*
199* Solve A*X = B, where A = U*D*U**H.
200*
201* First solve U*D*X = B, overwriting B with X.
202*
203* K is the main loop index, decreasing from N to 1 in steps of
204* 1 or 2, depending on the size of the diagonal blocks.
205*
206 k = n
207 10 CONTINUE
208*
209* If K < 1, exit from loop.
210*
211 IF( k.LT.1 )
212 $ GO TO 30
213*
214 IF( ipiv( k ).GT.0 ) THEN
215*
216* 1 x 1 diagonal block
217*
218* Interchange rows K and IPIV(K).
219*
220 kp = ipiv( k )
221 IF( kp.NE.k )
222 $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
223*
224* Multiply by inv(U(K)), where U(K) is the transformation
225* stored in column K of A.
226*
227 CALL cgeru( k-1, nrhs, -one, a( 1, k ), 1, b( k, 1 ), ldb,
228 $ b( 1, 1 ), ldb )
229*
230* Multiply by the inverse of the diagonal block.
231*
232 s = real( one ) / real( a( k, k ) )
233 CALL csscal( nrhs, s, b( k, 1 ), ldb )
234 k = k - 1
235 ELSE
236*
237* 2 x 2 diagonal block
238*
239* Interchange rows K and -IPIV(K), then K-1 and -IPIV(K-1)
240*
241 kp = -ipiv( k )
242 IF( kp.NE.k )
243 $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
244*
245 kp = -ipiv( k-1)
246 IF( kp.NE.k-1 )
247 $ CALL cswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
248*
249* Multiply by inv(U(K)), where U(K) is the transformation
250* stored in columns K-1 and K of A.
251*
252 CALL cgeru( k-2, nrhs, -one, a( 1, k ), 1, b( k, 1 ), ldb,
253 $ b( 1, 1 ), ldb )
254 CALL cgeru( k-2, nrhs, -one, a( 1, k-1 ), 1, b( k-1, 1 ),
255 $ ldb, b( 1, 1 ), ldb )
256*
257* Multiply by the inverse of the diagonal block.
258*
259 akm1k = a( k-1, k )
260 akm1 = a( k-1, k-1 ) / akm1k
261 ak = a( k, k ) / conjg( akm1k )
262 denom = akm1*ak - one
263 DO 20 j = 1, nrhs
264 bkm1 = b( k-1, j ) / akm1k
265 bk = b( k, j ) / conjg( akm1k )
266 b( k-1, j ) = ( ak*bkm1-bk ) / denom
267 b( k, j ) = ( akm1*bk-bkm1 ) / denom
268 20 CONTINUE
269 k = k - 2
270 END IF
271*
272 GO TO 10
273 30 CONTINUE
274*
275* Next solve U**H *X = B, overwriting B with X.
276*
277* K is the main loop index, increasing from 1 to N in steps of
278* 1 or 2, depending on the size of the diagonal blocks.
279*
280 k = 1
281 40 CONTINUE
282*
283* If K > N, exit from loop.
284*
285 IF( k.GT.n )
286 $ GO TO 50
287*
288 IF( ipiv( k ).GT.0 ) THEN
289*
290* 1 x 1 diagonal block
291*
292* Multiply by inv(U**H(K)), where U(K) is the transformation
293* stored in column K 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 END IF
301*
302* Interchange rows K and IPIV(K).
303*
304 kp = ipiv( k )
305 IF( kp.NE.k )
306 $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
307 k = k + 1
308 ELSE
309*
310* 2 x 2 diagonal block
311*
312* Multiply by inv(U**H(K+1)), where U(K+1) is the transformation
313* stored in columns K and K+1 of A.
314*
315 IF( k.GT.1 ) THEN
316 CALL clacgv( nrhs, b( k, 1 ), ldb )
317 CALL cgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
318 $ ldb, a( 1, k ), 1, one, b( k, 1 ), ldb )
319 CALL clacgv( nrhs, b( k, 1 ), ldb )
320*
321 CALL clacgv( nrhs, b( k+1, 1 ), ldb )
322 CALL cgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
323 $ ldb, a( 1, k+1 ), 1, one, b( k+1, 1 ), ldb )
324 CALL clacgv( nrhs, b( k+1, 1 ), ldb )
325 END IF
326*
327* Interchange rows K and -IPIV(K), then K+1 and -IPIV(K+1)
328*
329 kp = -ipiv( k )
330 IF( kp.NE.k )
331 $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
332*
333 kp = -ipiv( k+1 )
334 IF( kp.NE.k+1 )
335 $ CALL cswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
336*
337 k = k + 2
338 END IF
339*
340 GO TO 40
341 50 CONTINUE
342*
343 ELSE
344*
345* Solve A*X = B, where A = L*D*L**H.
346*
347* First solve L*D*X = B, overwriting B with X.
348*
349* K is the main loop index, increasing from 1 to N in steps of
350* 1 or 2, depending on the size of the diagonal blocks.
351*
352 k = 1
353 60 CONTINUE
354*
355* If K > N, exit from loop.
356*
357 IF( k.GT.n )
358 $ GO TO 80
359*
360 IF( ipiv( k ).GT.0 ) THEN
361*
362* 1 x 1 diagonal block
363*
364* Interchange rows K and IPIV(K).
365*
366 kp = ipiv( k )
367 IF( kp.NE.k )
368 $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
369*
370* Multiply by inv(L(K)), where L(K) is the transformation
371* stored in column K of A.
372*
373 IF( k.LT.n )
374 $ CALL cgeru( n-k, nrhs, -one, a( k+1, k ), 1, b( k, 1 ),
375 $ ldb, b( k+1, 1 ), ldb )
376*
377* Multiply by the inverse of the diagonal block.
378*
379 s = real( one ) / real( a( k, k ) )
380 CALL csscal( nrhs, s, b( k, 1 ), ldb )
381 k = k + 1
382 ELSE
383*
384* 2 x 2 diagonal block
385*
386* Interchange rows K and -IPIV(K), then K+1 and -IPIV(K+1)
387*
388 kp = -ipiv( k )
389 IF( kp.NE.k )
390 $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
391*
392 kp = -ipiv( k+1 )
393 IF( kp.NE.k+1 )
394 $ CALL cswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
395*
396* Multiply by inv(L(K)), where L(K) is the transformation
397* stored in columns K and K+1 of A.
398*
399 IF( k.LT.n-1 ) THEN
400 CALL cgeru( n-k-1, nrhs, -one, a( k+2, k ), 1, b( k, 1 ),
401 $ ldb, b( k+2, 1 ), ldb )
402 CALL cgeru( n-k-1, nrhs, -one, a( k+2, k+1 ), 1,
403 $ b( k+1, 1 ), ldb, b( k+2, 1 ), ldb )
404 END IF
405*
406* Multiply by the inverse of the diagonal block.
407*
408 akm1k = a( k+1, k )
409 akm1 = a( k, k ) / conjg( akm1k )
410 ak = a( k+1, k+1 ) / akm1k
411 denom = akm1*ak - one
412 DO 70 j = 1, nrhs
413 bkm1 = b( k, j ) / conjg( akm1k )
414 bk = b( k+1, j ) / akm1k
415 b( k, j ) = ( ak*bkm1-bk ) / denom
416 b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
417 70 CONTINUE
418 k = k + 2
419 END IF
420*
421 GO TO 60
422 80 CONTINUE
423*
424* Next solve L**H *X = B, overwriting B with X.
425*
426* K is the main loop index, decreasing from N to 1 in steps of
427* 1 or 2, depending on the size of the diagonal blocks.
428*
429 k = n
430 90 CONTINUE
431*
432* If K < 1, exit from loop.
433*
434 IF( k.LT.1 )
435 $ GO TO 100
436*
437 IF( ipiv( k ).GT.0 ) THEN
438*
439* 1 x 1 diagonal block
440*
441* Multiply by inv(L**H(K)), where L(K) is the transformation
442* stored in column K of A.
443*
444 IF( k.LT.n ) THEN
445 CALL clacgv( nrhs, b( k, 1 ), ldb )
446 CALL cgemv( 'Conjugate transpose', n-k, nrhs, -one,
447 $ b( k+1, 1 ), ldb, a( k+1, k ), 1, one,
448 $ b( k, 1 ), ldb )
449 CALL clacgv( nrhs, b( k, 1 ), ldb )
450 END IF
451*
452* Interchange rows K and IPIV(K).
453*
454 kp = ipiv( k )
455 IF( kp.NE.k )
456 $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
457 k = k - 1
458 ELSE
459*
460* 2 x 2 diagonal block
461*
462* Multiply by inv(L**H(K-1)), where L(K-1) is the transformation
463* stored in columns K-1 and K of A.
464*
465 IF( k.LT.n ) THEN
466 CALL clacgv( nrhs, b( k, 1 ), ldb )
467 CALL cgemv( 'Conjugate transpose', n-k, nrhs, -one,
468 $ b( k+1, 1 ), ldb, a( k+1, k ), 1, one,
469 $ b( k, 1 ), ldb )
470 CALL clacgv( nrhs, b( k, 1 ), ldb )
471*
472 CALL clacgv( nrhs, b( k-1, 1 ), ldb )
473 CALL cgemv( 'Conjugate transpose', n-k, nrhs, -one,
474 $ b( k+1, 1 ), ldb, a( k+1, k-1 ), 1, one,
475 $ b( k-1, 1 ), ldb )
476 CALL clacgv( nrhs, b( k-1, 1 ), ldb )
477 END IF
478*
479* Interchange rows K and -IPIV(K), then K-1 and -IPIV(K-1)
480*
481 kp = -ipiv( k )
482 IF( kp.NE.k )
483 $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
484*
485 kp = -ipiv( k-1 )
486 IF( kp.NE.k-1 )
487 $ CALL cswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
488*
489 k = k - 2
490 END IF
491*
492 GO TO 90
493 100 CONTINUE
494 END IF
495*
496 RETURN
497*
498* End of CHETRS_ROOK
499*
500 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_rook(uplo, n, nrhs, a, lda, ipiv, b, ldb, info)
CHETRS_ROOK computes the solution to a system of linear equations A * X = B for HE matrices using fac...
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