LAPACK  3.9.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.
Date
November 2013
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 138 of file chetrs_rook.f.

138 *
139 * -- LAPACK computational routine (version 3.5.0) --
140 * -- LAPACK is a software package provided by Univ. of Tennessee, --
141 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142 * November 2013
143 *
144 * .. Scalar Arguments ..
145  CHARACTER UPLO
146  INTEGER INFO, LDA, LDB, N, NRHS
147 * ..
148 * .. Array Arguments ..
149  INTEGER IPIV( * )
150  COMPLEX A( LDA, * ), B( LDB, * )
151 * ..
152 *
153 * =====================================================================
154 *
155 * .. Parameters ..
156  COMPLEX ONE
157  parameter( one = ( 1.0e+0, 0.0e+0 ) )
158 * ..
159 * .. Local Scalars ..
160  LOGICAL UPPER
161  INTEGER J, K, KP
162  REAL S
163  COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM
164 * ..
165 * .. External Functions ..
166  LOGICAL LSAME
167  EXTERNAL lsame
168 * ..
169 * .. External Subroutines ..
170  EXTERNAL cgemv, cgeru, clacgv, csscal, cswap, xerbla
171 * ..
172 * .. Intrinsic Functions ..
173  INTRINSIC conjg, max, real
174 * ..
175 * .. Executable Statements ..
176 *
177  info = 0
178  upper = lsame( uplo, 'U' )
179  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
180  info = -1
181  ELSE IF( n.LT.0 ) THEN
182  info = -2
183  ELSE IF( nrhs.LT.0 ) THEN
184  info = -3
185  ELSE IF( lda.LT.max( 1, n ) ) THEN
186  info = -5
187  ELSE IF( ldb.LT.max( 1, n ) ) THEN
188  info = -8
189  END IF
190  IF( info.NE.0 ) THEN
191  CALL xerbla( 'CHETRS_ROOK', -info )
192  RETURN
193  END IF
194 *
195 * Quick return if possible
196 *
197  IF( n.EQ.0 .OR. nrhs.EQ.0 )
198  $ RETURN
199 *
200  IF( upper ) THEN
201 *
202 * Solve A*X = B, where A = U*D*U**H.
203 *
204 * First solve U*D*X = B, overwriting B with X.
205 *
206 * K is the main loop index, decreasing from N to 1 in steps of
207 * 1 or 2, depending on the size of the diagonal blocks.
208 *
209  k = n
210  10 CONTINUE
211 *
212 * If K < 1, exit from loop.
213 *
214  IF( k.LT.1 )
215  $ GO TO 30
216 *
217  IF( ipiv( k ).GT.0 ) THEN
218 *
219 * 1 x 1 diagonal block
220 *
221 * Interchange rows K and IPIV(K).
222 *
223  kp = ipiv( k )
224  IF( kp.NE.k )
225  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
226 *
227 * Multiply by inv(U(K)), where U(K) is the transformation
228 * stored in column K of A.
229 *
230  CALL cgeru( k-1, nrhs, -one, a( 1, k ), 1, b( k, 1 ), ldb,
231  $ b( 1, 1 ), ldb )
232 *
233 * Multiply by the inverse of the diagonal block.
234 *
235  s = real( one ) / real( a( k, k ) )
236  CALL csscal( nrhs, s, b( k, 1 ), ldb )
237  k = k - 1
238  ELSE
239 *
240 * 2 x 2 diagonal block
241 *
242 * Interchange rows K and -IPIV(K), then K-1 and -IPIV(K-1)
243 *
244  kp = -ipiv( k )
245  IF( kp.NE.k )
246  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
247 *
248  kp = -ipiv( k-1)
249  IF( kp.NE.k-1 )
250  $ CALL cswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
251 *
252 * Multiply by inv(U(K)), where U(K) is the transformation
253 * stored in columns K-1 and K of A.
254 *
255  CALL cgeru( k-2, nrhs, -one, a( 1, k ), 1, b( k, 1 ), ldb,
256  $ b( 1, 1 ), ldb )
257  CALL cgeru( k-2, nrhs, -one, a( 1, k-1 ), 1, b( k-1, 1 ),
258  $ ldb, b( 1, 1 ), ldb )
259 *
260 * Multiply by the inverse of the diagonal block.
261 *
262  akm1k = a( k-1, k )
263  akm1 = a( k-1, k-1 ) / akm1k
264  ak = a( k, k ) / conjg( akm1k )
265  denom = akm1*ak - one
266  DO 20 j = 1, nrhs
267  bkm1 = b( k-1, j ) / akm1k
268  bk = b( k, j ) / conjg( akm1k )
269  b( k-1, j ) = ( ak*bkm1-bk ) / denom
270  b( k, j ) = ( akm1*bk-bkm1 ) / denom
271  20 CONTINUE
272  k = k - 2
273  END IF
274 *
275  GO TO 10
276  30 CONTINUE
277 *
278 * Next solve U**H *X = B, overwriting B with X.
279 *
280 * K is the main loop index, increasing from 1 to N in steps of
281 * 1 or 2, depending on the size of the diagonal blocks.
282 *
283  k = 1
284  40 CONTINUE
285 *
286 * If K > N, exit from loop.
287 *
288  IF( k.GT.n )
289  $ GO TO 50
290 *
291  IF( ipiv( k ).GT.0 ) THEN
292 *
293 * 1 x 1 diagonal block
294 *
295 * Multiply by inv(U**H(K)), where U(K) is the transformation
296 * stored in column K of A.
297 *
298  IF( k.GT.1 ) THEN
299  CALL clacgv( nrhs, b( k, 1 ), ldb )
300  CALL cgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
301  $ ldb, a( 1, k ), 1, one, b( k, 1 ), ldb )
302  CALL clacgv( nrhs, b( k, 1 ), ldb )
303  END IF
304 *
305 * Interchange rows K and IPIV(K).
306 *
307  kp = ipiv( k )
308  IF( kp.NE.k )
309  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
310  k = k + 1
311  ELSE
312 *
313 * 2 x 2 diagonal block
314 *
315 * Multiply by inv(U**H(K+1)), where U(K+1) is the transformation
316 * stored in columns K and K+1 of A.
317 *
318  IF( k.GT.1 ) THEN
319  CALL clacgv( nrhs, b( k, 1 ), ldb )
320  CALL cgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
321  $ ldb, a( 1, k ), 1, one, b( k, 1 ), ldb )
322  CALL clacgv( nrhs, b( k, 1 ), ldb )
323 *
324  CALL clacgv( nrhs, b( k+1, 1 ), ldb )
325  CALL cgemv( 'Conjugate transpose', k-1, nrhs, -one, b,
326  $ ldb, a( 1, k+1 ), 1, one, b( k+1, 1 ), ldb )
327  CALL clacgv( nrhs, b( k+1, 1 ), ldb )
328  END IF
329 *
330 * Interchange rows K and -IPIV(K), then K+1 and -IPIV(K+1)
331 *
332  kp = -ipiv( k )
333  IF( kp.NE.k )
334  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
335 *
336  kp = -ipiv( k+1 )
337  IF( kp.NE.k+1 )
338  $ CALL cswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
339 *
340  k = k + 2
341  END IF
342 *
343  GO TO 40
344  50 CONTINUE
345 *
346  ELSE
347 *
348 * Solve A*X = B, where A = L*D*L**H.
349 *
350 * First solve L*D*X = B, overwriting B with X.
351 *
352 * K is the main loop index, increasing from 1 to N in steps of
353 * 1 or 2, depending on the size of the diagonal blocks.
354 *
355  k = 1
356  60 CONTINUE
357 *
358 * If K > N, exit from loop.
359 *
360  IF( k.GT.n )
361  $ GO TO 80
362 *
363  IF( ipiv( k ).GT.0 ) THEN
364 *
365 * 1 x 1 diagonal block
366 *
367 * Interchange rows K and IPIV(K).
368 *
369  kp = ipiv( k )
370  IF( kp.NE.k )
371  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
372 *
373 * Multiply by inv(L(K)), where L(K) is the transformation
374 * stored in column K of A.
375 *
376  IF( k.LT.n )
377  $ CALL cgeru( n-k, nrhs, -one, a( k+1, k ), 1, b( k, 1 ),
378  $ ldb, b( k+1, 1 ), ldb )
379 *
380 * Multiply by the inverse of the diagonal block.
381 *
382  s = real( one ) / real( a( k, k ) )
383  CALL csscal( nrhs, s, b( k, 1 ), ldb )
384  k = k + 1
385  ELSE
386 *
387 * 2 x 2 diagonal block
388 *
389 * Interchange rows K and -IPIV(K), then K+1 and -IPIV(K+1)
390 *
391  kp = -ipiv( k )
392  IF( kp.NE.k )
393  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
394 *
395  kp = -ipiv( k+1 )
396  IF( kp.NE.k+1 )
397  $ CALL cswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
398 *
399 * Multiply by inv(L(K)), where L(K) is the transformation
400 * stored in columns K and K+1 of A.
401 *
402  IF( k.LT.n-1 ) THEN
403  CALL cgeru( n-k-1, nrhs, -one, a( k+2, k ), 1, b( k, 1 ),
404  $ ldb, b( k+2, 1 ), ldb )
405  CALL cgeru( 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 ) / conjg( 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 ) / conjg( 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 clacgv( nrhs, b( k, 1 ), ldb )
449  CALL cgemv( '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 clacgv( 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 cswap( 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 clacgv( nrhs, b( k, 1 ), ldb )
470  CALL cgemv( '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 clacgv( nrhs, b( k, 1 ), ldb )
474 *
475  CALL clacgv( nrhs, b( k-1, 1 ), ldb )
476  CALL cgemv( '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 clacgv( 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 cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
487 *
488  kp = -ipiv( k-1 )
489  IF( kp.NE.k-1 )
490  $ CALL cswap( 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 CHETRS_ROOK
502 *
Here is the call graph for this function:
Here is the caller graph for this function:
csscal
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:80
cgemv
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
clacgv
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:76
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:83
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
cgeru
subroutine cgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CGERU
Definition: cgeru.f:132