LAPACK 3.12.0
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 134 of file zhetrs_rook.f.

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*16 A( LDA, * ), B( LDB, * )
148* ..
149*
150* =====================================================================
151*
152* .. Parameters ..
153 COMPLEX*16 ONE
154 parameter( one = ( 1.0d+0, 0.0d+0 ) )
155* ..
156* .. Local Scalars ..
157 LOGICAL UPPER
158 INTEGER J, K, KP
159 DOUBLE PRECISION S
160 COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
161* ..
162* .. External Functions ..
163 LOGICAL LSAME
164 EXTERNAL lsame
165* ..
166* .. External Subroutines ..
167 EXTERNAL zgemv, zgeru, zlacgv, zdscal, zswap, xerbla
168* ..
169* .. Intrinsic Functions ..
170 INTRINSIC dconjg, max, dble
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( 'ZHETRS_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 zswap( 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 zgeru( 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 = 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 ), ldb,
253 $ b( 1, 1 ), ldb )
254 CALL zgeru( 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 ) / dconjg( akm1k )
262 denom = akm1*ak - one
263 DO 20 j = 1, nrhs
264 bkm1 = b( k-1, j ) / akm1k
265 bk = b( k, j ) / dconjg( 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 zlacgv( nrhs, b( k, 1 ), ldb )
297 CALL zgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
298 $ ldb, a( 1, k ), 1, one, b( k, 1 ), ldb )
299 CALL zlacgv( 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 zswap( 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 zlacgv( nrhs, b( k, 1 ), ldb )
317 CALL zgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
318 $ ldb, a( 1, k ), 1, one, b( k, 1 ), ldb )
319 CALL zlacgv( nrhs, b( k, 1 ), ldb )
320*
321 CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
322 CALL zgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
323 $ ldb, a( 1, k+1 ), 1, one, b( k+1, 1 ), ldb )
324 CALL zlacgv( 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 zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
332*
333 kp = -ipiv( k+1 )
334 IF( kp.NE.k+1 )
335 $ CALL zswap( 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 zswap( 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 zgeru( 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 = dble( one ) / dble( a( k, k ) )
380 CALL zdscal( 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 zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
391*
392 kp = -ipiv( k+1 )
393 IF( kp.NE.k+1 )
394 $ CALL zswap( 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 zgeru( n-k-1, nrhs, -one, a( k+2, k ), 1, b( k, 1 ),
401 $ ldb, b( k+2, 1 ), ldb )
402 CALL zgeru( 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 ) / dconjg( 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 ) / dconjg( 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 zlacgv( nrhs, b( k, 1 ), ldb )
446 CALL zgemv( '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 zlacgv( 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 zswap( 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 zlacgv( nrhs, b( k, 1 ), ldb )
467 CALL zgemv( '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 zlacgv( nrhs, b( k, 1 ), ldb )
471*
472 CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
473 CALL zgemv( '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 zlacgv( 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 zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
484*
485 kp = -ipiv( k-1 )
486 IF( kp.NE.k-1 )
487 $ CALL zswap( 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 ZHETRS_ROOK
499*
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:74
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: