LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
subroutine zgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERU
Definition: zgeru.f:130
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:74
Here is the call graph for this function:
Here is the caller graph for this function: