LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ chetrs_rook()

subroutine chetrs_rook ( character  UPLO,
integer  N,
integer  NRHS,
complex, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

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)

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

Purpose:
 CHETRS_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 CHETRF_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 array, dimension (LDA,N)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L as computed by CHETRF_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 CHETRF_ROOK.
[in,out]B
          B is COMPLEX 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 chetrs_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 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 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine cgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CGERU
Definition: cgeru.f:130
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:74
Here is the call graph for this function:
Here is the caller graph for this function: