LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zhetrs_rook()

subroutine zhetrs_rook ( character uplo,
integer n,
integer nrhs,
complex*16, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
complex*16, dimension( ldb, * ) b,
integer ldb,
integer info )

ZHETRS_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)

Download ZHETRS_ROOK + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZHETRS_ROOK solves a system of linear equations A*X = B with a complex
!> Hermitian matrix A using the factorization A = U*D*U**H or
!> A = L*D*L**H computed by ZHETRF_ROOK.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          Specifies whether the details of the factorization are stored
!>          as an upper or lower triangular matrix.
!>          = 'U':  Upper triangular, form is A = U*D*U**H;
!>          = 'L':  Lower triangular, form is A = L*D*L**H.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!>
[in]A
!>          A is COMPLEX*16 array, dimension (LDA,N)
!>          The block diagonal matrix D and the multipliers used to
!>          obtain the factor U or L as computed by ZHETRF_ROOK.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          Details of the interchanges and the block structure of D
!>          as determined by ZHETRF_ROOK.
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          On entry, the right hand side matrix B.
!>          On exit, the solution matrix X.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
!>
!>  November 2013,  Igor Kozachenko,
!>                  Computer Science Division,
!>                  University of California, Berkeley
!>
!>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
!>                  School of Mathematics,
!>                  University of Manchester
!>
!>

Definition at line 132 of file zhetrs_rook.f.

134*
135* -- LAPACK computational routine --
136* -- LAPACK is a software package provided by Univ. of Tennessee, --
137* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138*
139* .. Scalar Arguments ..
140 CHARACTER UPLO
141 INTEGER INFO, LDA, LDB, N, NRHS
142* ..
143* .. Array Arguments ..
144 INTEGER IPIV( * )
145 COMPLEX*16 A( LDA, * ), B( LDB, * )
146* ..
147*
148* =====================================================================
149*
150* .. Parameters ..
151 COMPLEX*16 ONE
152 parameter( one = ( 1.0d+0, 0.0d+0 ) )
153* ..
154* .. Local Scalars ..
155 LOGICAL UPPER
156 INTEGER J, K, KP
157 DOUBLE PRECISION S
158 COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
159* ..
160* .. External Functions ..
161 LOGICAL LSAME
162 EXTERNAL lsame
163* ..
164* .. External Subroutines ..
165 EXTERNAL zgemv, zgeru, zlacgv, zdscal, zswap,
166 $ xerbla
167* ..
168* .. Intrinsic Functions ..
169 INTRINSIC dconjg, max, dble
170* ..
171* .. Executable Statements ..
172*
173 info = 0
174 upper = lsame( uplo, 'U' )
175 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
176 info = -1
177 ELSE IF( n.LT.0 ) THEN
178 info = -2
179 ELSE IF( nrhs.LT.0 ) THEN
180 info = -3
181 ELSE IF( lda.LT.max( 1, n ) ) THEN
182 info = -5
183 ELSE IF( ldb.LT.max( 1, n ) ) THEN
184 info = -8
185 END IF
186 IF( info.NE.0 ) THEN
187 CALL xerbla( 'ZHETRS_ROOK', -info )
188 RETURN
189 END IF
190*
191* Quick return if possible
192*
193 IF( n.EQ.0 .OR. nrhs.EQ.0 )
194 $ RETURN
195*
196 IF( upper ) THEN
197*
198* Solve A*X = B, where A = U*D*U**H.
199*
200* First solve U*D*X = B, overwriting B with X.
201*
202* K is the main loop index, decreasing from N to 1 in steps of
203* 1 or 2, depending on the size of the diagonal blocks.
204*
205 k = n
206 10 CONTINUE
207*
208* If K < 1, exit from loop.
209*
210 IF( k.LT.1 )
211 $ GO TO 30
212*
213 IF( ipiv( k ).GT.0 ) THEN
214*
215* 1 x 1 diagonal block
216*
217* Interchange rows K and IPIV(K).
218*
219 kp = ipiv( k )
220 IF( kp.NE.k )
221 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
222*
223* Multiply by inv(U(K)), where U(K) is the transformation
224* stored in column K of A.
225*
226 CALL zgeru( k-1, nrhs, -one, a( 1, k ), 1, b( k, 1 ),
227 $ ldb,
228 $ b( 1, 1 ), ldb )
229*
230* Multiply by the inverse of the diagonal block.
231*
232 s = dble( one ) / dble( a( k, k ) )
233 CALL zdscal( 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 zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
244*
245 kp = -ipiv( k-1)
246 IF( kp.NE.k-1 )
247 $ CALL zswap( 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 zgeru( k-2, nrhs, -one, a( 1, k ), 1, b( k, 1 ),
253 $ ldb,
254 $ b( 1, 1 ), ldb )
255 CALL zgeru( k-2, nrhs, -one, a( 1, k-1 ), 1, b( k-1, 1 ),
256 $ ldb, b( 1, 1 ), ldb )
257*
258* Multiply by the inverse of the diagonal block.
259*
260 akm1k = a( k-1, k )
261 akm1 = a( k-1, k-1 ) / akm1k
262 ak = a( k, k ) / dconjg( akm1k )
263 denom = akm1*ak - one
264 DO 20 j = 1, nrhs
265 bkm1 = b( k-1, j ) / akm1k
266 bk = b( k, j ) / dconjg( akm1k )
267 b( k-1, j ) = ( ak*bkm1-bk ) / denom
268 b( k, j ) = ( akm1*bk-bkm1 ) / denom
269 20 CONTINUE
270 k = k - 2
271 END IF
272*
273 GO TO 10
274 30 CONTINUE
275*
276* Next solve U**H *X = B, overwriting B with X.
277*
278* K is the main loop index, increasing from 1 to N in steps of
279* 1 or 2, depending on the size of the diagonal blocks.
280*
281 k = 1
282 40 CONTINUE
283*
284* If K > N, exit from loop.
285*
286 IF( k.GT.n )
287 $ GO TO 50
288*
289 IF( ipiv( k ).GT.0 ) THEN
290*
291* 1 x 1 diagonal block
292*
293* Multiply by inv(U**H(K)), where U(K) is the transformation
294* stored in column K of A.
295*
296 IF( k.GT.1 ) THEN
297 CALL zlacgv( nrhs, b( k, 1 ), ldb )
298 CALL zgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
299 $ ldb, a( 1, k ), 1, one, b( k, 1 ), ldb )
300 CALL zlacgv( nrhs, b( k, 1 ), ldb )
301 END IF
302*
303* Interchange rows K and IPIV(K).
304*
305 kp = ipiv( k )
306 IF( kp.NE.k )
307 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
308 k = k + 1
309 ELSE
310*
311* 2 x 2 diagonal block
312*
313* Multiply by inv(U**H(K+1)), where U(K+1) is the transformation
314* stored in columns K and K+1 of A.
315*
316 IF( k.GT.1 ) THEN
317 CALL zlacgv( nrhs, b( k, 1 ), ldb )
318 CALL zgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
319 $ ldb, a( 1, k ), 1, one, b( k, 1 ), ldb )
320 CALL zlacgv( nrhs, b( k, 1 ), ldb )
321*
322 CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
323 CALL zgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
324 $ ldb, a( 1, k+1 ), 1, one, b( k+1, 1 ), ldb )
325 CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
326 END IF
327*
328* Interchange rows K and -IPIV(K), then K+1 and -IPIV(K+1)
329*
330 kp = -ipiv( k )
331 IF( kp.NE.k )
332 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
333*
334 kp = -ipiv( k+1 )
335 IF( kp.NE.k+1 )
336 $ CALL zswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
337*
338 k = k + 2
339 END IF
340*
341 GO TO 40
342 50 CONTINUE
343*
344 ELSE
345*
346* Solve A*X = B, where A = L*D*L**H.
347*
348* First solve L*D*X = B, overwriting B with X.
349*
350* K is the main loop index, increasing from 1 to N in steps of
351* 1 or 2, depending on the size of the diagonal blocks.
352*
353 k = 1
354 60 CONTINUE
355*
356* If K > N, exit from loop.
357*
358 IF( k.GT.n )
359 $ GO TO 80
360*
361 IF( ipiv( k ).GT.0 ) THEN
362*
363* 1 x 1 diagonal block
364*
365* Interchange rows K and IPIV(K).
366*
367 kp = ipiv( k )
368 IF( kp.NE.k )
369 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
370*
371* Multiply by inv(L(K)), where L(K) is the transformation
372* stored in column K of A.
373*
374 IF( k.LT.n )
375 $ CALL zgeru( n-k, nrhs, -one, a( k+1, k ), 1, b( k,
376 $ 1 ),
377 $ ldb, b( k+1, 1 ), ldb )
378*
379* Multiply by the inverse of the diagonal block.
380*
381 s = dble( one ) / dble( a( k, k ) )
382 CALL zdscal( nrhs, s, b( k, 1 ), ldb )
383 k = k + 1
384 ELSE
385*
386* 2 x 2 diagonal block
387*
388* Interchange rows K and -IPIV(K), then K+1 and -IPIV(K+1)
389*
390 kp = -ipiv( k )
391 IF( kp.NE.k )
392 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
393*
394 kp = -ipiv( k+1 )
395 IF( kp.NE.k+1 )
396 $ CALL zswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
397*
398* Multiply by inv(L(K)), where L(K) is the transformation
399* stored in columns K and K+1 of A.
400*
401 IF( k.LT.n-1 ) THEN
402 CALL zgeru( n-k-1, nrhs, -one, a( k+2, k ), 1, b( k,
403 $ 1 ),
404 $ ldb, b( k+2, 1 ), ldb )
405 CALL zgeru( n-k-1, nrhs, -one, a( k+2, k+1 ), 1,
406 $ b( k+1, 1 ), ldb, b( k+2, 1 ), ldb )
407 END IF
408*
409* Multiply by the inverse of the diagonal block.
410*
411 akm1k = a( k+1, k )
412 akm1 = a( k, k ) / dconjg( akm1k )
413 ak = a( k+1, k+1 ) / akm1k
414 denom = akm1*ak - one
415 DO 70 j = 1, nrhs
416 bkm1 = b( k, j ) / dconjg( akm1k )
417 bk = b( k+1, j ) / akm1k
418 b( k, j ) = ( ak*bkm1-bk ) / denom
419 b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
420 70 CONTINUE
421 k = k + 2
422 END IF
423*
424 GO TO 60
425 80 CONTINUE
426*
427* Next solve L**H *X = B, overwriting B with X.
428*
429* K is the main loop index, decreasing from N to 1 in steps of
430* 1 or 2, depending on the size of the diagonal blocks.
431*
432 k = n
433 90 CONTINUE
434*
435* If K < 1, exit from loop.
436*
437 IF( k.LT.1 )
438 $ GO TO 100
439*
440 IF( ipiv( k ).GT.0 ) THEN
441*
442* 1 x 1 diagonal block
443*
444* Multiply by inv(L**H(K)), where L(K) is the transformation
445* stored in column K of A.
446*
447 IF( k.LT.n ) THEN
448 CALL zlacgv( nrhs, b( k, 1 ), ldb )
449 CALL zgemv( 'Conjugate transpose', n-k, nrhs, -one,
450 $ b( k+1, 1 ), ldb, a( k+1, k ), 1, one,
451 $ b( k, 1 ), ldb )
452 CALL zlacgv( nrhs, b( k, 1 ), ldb )
453 END IF
454*
455* Interchange rows K and IPIV(K).
456*
457 kp = ipiv( k )
458 IF( kp.NE.k )
459 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
460 k = k - 1
461 ELSE
462*
463* 2 x 2 diagonal block
464*
465* Multiply by inv(L**H(K-1)), where L(K-1) is the transformation
466* stored in columns K-1 and K of A.
467*
468 IF( k.LT.n ) THEN
469 CALL zlacgv( nrhs, b( k, 1 ), ldb )
470 CALL zgemv( 'Conjugate transpose', n-k, nrhs, -one,
471 $ b( k+1, 1 ), ldb, a( k+1, k ), 1, one,
472 $ b( k, 1 ), ldb )
473 CALL zlacgv( nrhs, b( k, 1 ), ldb )
474*
475 CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
476 CALL zgemv( 'Conjugate transpose', n-k, nrhs, -one,
477 $ b( k+1, 1 ), ldb, a( k+1, k-1 ), 1, one,
478 $ b( k-1, 1 ), ldb )
479 CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
480 END IF
481*
482* Interchange rows K and -IPIV(K), then K-1 and -IPIV(K-1)
483*
484 kp = -ipiv( k )
485 IF( kp.NE.k )
486 $ CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
487*
488 kp = -ipiv( k-1 )
489 IF( kp.NE.k-1 )
490 $ CALL zswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
491*
492 k = k - 2
493 END IF
494*
495 GO TO 90
496 100 CONTINUE
497 END IF
498*
499 RETURN
500*
501* End of ZHETRS_ROOK
502*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
zgeru
subroutine zgeru(m, n, alpha, x, incx, y, incy, a, lda)
ZGERU
Definition zgeru.f:130
zlacgv
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
Definition zlacgv.f:72
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
zdscal
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: