LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
clavhe_rook.f
Go to the documentation of this file.
1 *> \brief \b CLAVHE_ROOK
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE CLAVHE_ROOK( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B,
12 * LDB, INFO )
13 *
14 * .. Scalar Arguments ..
15 * CHARACTER DIAG, TRANS, UPLO
16 * INTEGER INFO, LDA, LDB, N, NRHS
17 * ..
18 * .. Array Arguments ..
19 * INTEGER IPIV( * )
20 * COMPLEX A( LDA, * ), B( LDB, * )
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> CLAVHE_ROOK performs one of the matrix-vector operations
30 *> x := A*x or x := A^H*x,
31 *> where x is an N element vector and A is one of the factors
32 *> from the block U*D*U' or L*D*L' factorization computed by CHETRF_ROOK.
33 *>
34 *> If TRANS = 'N', multiplies by U or U * D (or L or L * D)
35 *> If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L')
36 *> \endverbatim
37 *
38 * Arguments:
39 * ==========
40 *
41 *> \param[in] UPLO
42 *> \verbatim
43 *> UPLO is CHARACTER*1
44 *> Specifies whether the factor stored in A is upper or lower
45 *> triangular.
46 *> = 'U': Upper triangular
47 *> = 'L': Lower triangular
48 *> \endverbatim
49 *>
50 *> \param[in] TRANS
51 *> \verbatim
52 *> TRANS is CHARACTER*1
53 *> Specifies the operation to be performed:
54 *> = 'N': x := A*x
55 *> = 'C': x := A^H*x
56 *> \endverbatim
57 *>
58 *> \param[in] DIAG
59 *> \verbatim
60 *> DIAG is CHARACTER*1
61 *> Specifies whether or not the diagonal blocks are unit
62 *> matrices. If the diagonal blocks are assumed to be unit,
63 *> then A = U or A = L, otherwise A = U*D or A = L*D.
64 *> = 'U': Diagonal blocks are assumed to be unit matrices.
65 *> = 'N': Diagonal blocks are assumed to be non-unit matrices.
66 *> \endverbatim
67 *>
68 *> \param[in] N
69 *> \verbatim
70 *> N is INTEGER
71 *> The number of rows and columns of the matrix A. N >= 0.
72 *> \endverbatim
73 *>
74 *> \param[in] NRHS
75 *> \verbatim
76 *> NRHS is INTEGER
77 *> The number of right hand sides, i.e., the number of vectors
78 *> x to be multiplied by A. NRHS >= 0.
79 *> \endverbatim
80 *>
81 *> \param[in] A
82 *> \verbatim
83 *> A is COMPLEX array, dimension (LDA,N)
84 *> The block diagonal matrix D and the multipliers used to
85 *> obtain the factor U or L as computed by CHETRF_ROOK.
86 *> Stored as a 2-D triangular matrix.
87 *> \endverbatim
88 *>
89 *> \param[in] LDA
90 *> \verbatim
91 *> LDA is INTEGER
92 *> The leading dimension of the array A. LDA >= max(1,N).
93 *> \endverbatim
94 *>
95 *> \param[out] IPIV
96 *> \verbatim
97 *> IPIV is INTEGER array, dimension (N)
98 *> Details of the interchanges and the block structure of D,
99 *> as determined by CHETRF_ROOK.
100 *> If UPLO = 'U':
101 *> Only the last KB elements of IPIV are set.
102 *>
103 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
104 *> interchanged and D(k,k) is a 1-by-1 diagonal block.
105 *>
106 *> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
107 *> columns k and -IPIV(k) were interchanged and rows and
108 *> columns k-1 and -IPIV(k-1) were inerchaged,
109 *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
110 *>
111 *> If UPLO = 'L':
112 *> Only the first KB elements of IPIV are set.
113 *>
114 *> If IPIV(k) > 0, then rows and columns k and IPIV(k)
115 *> were interchanged and D(k,k) is a 1-by-1 diagonal block.
116 *>
117 *> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
118 *> columns k and -IPIV(k) were interchanged and rows and
119 *> columns k+1 and -IPIV(k+1) were inerchaged,
120 *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
121 *> \endverbatim
122 *>
123 *> \param[in,out] B
124 *> \verbatim
125 *> B is COMPLEX array, dimension (LDB,NRHS)
126 *> On entry, B contains NRHS vectors of length N.
127 *> On exit, B is overwritten with the product A * B.
128 *> \endverbatim
129 *>
130 *> \param[in] LDB
131 *> \verbatim
132 *> LDB is INTEGER
133 *> The leading dimension of the array B. LDB >= max(1,N).
134 *> \endverbatim
135 *>
136 *> \param[out] INFO
137 *> \verbatim
138 *> INFO is INTEGER
139 *> = 0: successful exit
140 *> < 0: if INFO = -k, the k-th argument had an illegal value
141 *> \endverbatim
142 *
143 * Authors:
144 * ========
145 *
146 *> \author Univ. of Tennessee
147 *> \author Univ. of California Berkeley
148 *> \author Univ. of Colorado Denver
149 *> \author NAG Ltd.
150 *
151 *> \date November 2013
152 *
153 *> \ingroup complex_lin
154 *
155 * =====================================================================
156  SUBROUTINE clavhe_rook( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV,
157  $ b, ldb, info )
158 *
159 * -- LAPACK test routine (version 3.5.0) --
160 * -- LAPACK is a software package provided by Univ. of Tennessee, --
161 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
162 * November 2013
163 *
164 * .. Scalar Arguments ..
165  CHARACTER DIAG, TRANS, UPLO
166  INTEGER INFO, LDA, LDB, N, NRHS
167 * ..
168 * .. Array Arguments ..
169  INTEGER IPIV( * )
170  COMPLEX A( lda, * ), B( ldb, * )
171 * ..
172 *
173 * =====================================================================
174 *
175 * .. Parameters ..
176  COMPLEX CONE
177  parameter ( cone = ( 1.0e+0, 0.0e+0 ) )
178 * ..
179 * .. Local Scalars ..
180  LOGICAL NOUNIT
181  INTEGER J, K, KP
182  COMPLEX D11, D12, D21, D22, T1, T2
183 * ..
184 * .. External Functions ..
185  LOGICAL LSAME
186  EXTERNAL lsame
187 * ..
188 * .. External Subroutines ..
189  EXTERNAL cgemv, cgeru, clacgv, cscal, cswap, xerbla
190 * ..
191 * .. Intrinsic Functions ..
192  INTRINSIC abs, conjg, max
193 * ..
194 * .. Executable Statements ..
195 *
196 * Test the input parameters.
197 *
198  info = 0
199  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
200  info = -1
201  ELSE IF( .NOT.lsame( trans, 'N' ) .AND. .NOT.lsame( trans, 'C' ) )
202  $ THEN
203  info = -2
204  ELSE IF( .NOT.lsame( diag, 'U' ) .AND. .NOT.lsame( diag, 'N' ) )
205  $ THEN
206  info = -3
207  ELSE IF( n.LT.0 ) THEN
208  info = -4
209  ELSE IF( lda.LT.max( 1, n ) ) THEN
210  info = -6
211  ELSE IF( ldb.LT.max( 1, n ) ) THEN
212  info = -9
213  END IF
214  IF( info.NE.0 ) THEN
215  CALL xerbla( 'CLAVHE_ROOK ', -info )
216  RETURN
217  END IF
218 *
219 * Quick return if possible.
220 *
221  IF( n.EQ.0 )
222  $ RETURN
223 *
224  nounit = lsame( diag, 'N' )
225 *------------------------------------------
226 *
227 * Compute B := A * B (No transpose)
228 *
229 *------------------------------------------
230  IF( lsame( trans, 'N' ) ) THEN
231 *
232 * Compute B := U*B
233 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
234 *
235  IF( lsame( uplo, 'U' ) ) THEN
236 *
237 * Loop forward applying the transformations.
238 *
239  k = 1
240  10 CONTINUE
241  IF( k.GT.n )
242  $ GO TO 30
243  IF( ipiv( k ).GT.0 ) THEN
244 *
245 * 1 x 1 pivot block
246 *
247 * Multiply by the diagonal element if forming U * D.
248 *
249  IF( nounit )
250  $ CALL cscal( nrhs, a( k, k ), b( k, 1 ), ldb )
251 *
252 * Multiply by P(K) * inv(U(K)) if K > 1.
253 *
254  IF( k.GT.1 ) THEN
255 *
256 * Apply the transformation.
257 *
258  CALL cgeru( k-1, nrhs, cone, a( 1, k ), 1, b( k, 1 ),
259  $ ldb, b( 1, 1 ), ldb )
260 *
261 * Interchange if P(K) != I.
262 *
263  kp = ipiv( k )
264  IF( kp.NE.k )
265  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
266  END IF
267  k = k + 1
268  ELSE
269 *
270 * 2 x 2 pivot block
271 *
272 * Multiply by the diagonal block if forming U * D.
273 *
274  IF( nounit ) THEN
275  d11 = a( k, k )
276  d22 = a( k+1, k+1 )
277  d12 = a( k, k+1 )
278  d21 = conjg( d12 )
279  DO 20 j = 1, nrhs
280  t1 = b( k, j )
281  t2 = b( k+1, j )
282  b( k, j ) = d11*t1 + d12*t2
283  b( k+1, j ) = d21*t1 + d22*t2
284  20 CONTINUE
285  END IF
286 *
287 * Multiply by P(K) * inv(U(K)) if K > 1.
288 *
289  IF( k.GT.1 ) THEN
290 *
291 * Apply the transformations.
292 *
293  CALL cgeru( k-1, nrhs, cone, a( 1, k ), 1, b( k, 1 ),
294  $ ldb, b( 1, 1 ), ldb )
295  CALL cgeru( k-1, nrhs, cone, a( 1, k+1 ), 1,
296  $ b( k+1, 1 ), ldb, b( 1, 1 ), ldb )
297 *
298 * Interchange if a permutation was applied at the
299 * K-th step of the factorization.
300 *
301 * Swap the first of pair with IMAXth
302 *
303  kp = abs( ipiv( k ) )
304  IF( kp.NE.k )
305  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
306 *
307 * NOW swap the first of pair with Pth
308 *
309  kp = abs( ipiv( k+1 ) )
310  IF( kp.NE.k+1 )
311  $ CALL cswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
312  $ ldb )
313  END IF
314  k = k + 2
315  END IF
316  GO TO 10
317  30 CONTINUE
318 *
319 * Compute B := L*B
320 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
321 *
322  ELSE
323 *
324 * Loop backward applying the transformations to B.
325 *
326  k = n
327  40 CONTINUE
328  IF( k.LT.1 )
329  $ GO TO 60
330 *
331 * Test the pivot index. If greater than zero, a 1 x 1
332 * pivot was used, otherwise a 2 x 2 pivot was used.
333 *
334  IF( ipiv( k ).GT.0 ) THEN
335 *
336 * 1 x 1 pivot block:
337 *
338 * Multiply by the diagonal element if forming L * D.
339 *
340  IF( nounit )
341  $ CALL cscal( nrhs, a( k, k ), b( k, 1 ), ldb )
342 *
343 * Multiply by P(K) * inv(L(K)) if K < N.
344 *
345  IF( k.NE.n ) THEN
346  kp = ipiv( k )
347 *
348 * Apply the transformation.
349 *
350  CALL cgeru( n-k, nrhs, cone, a( k+1, k ), 1,
351  $ b( k, 1 ), ldb, b( k+1, 1 ), ldb )
352 *
353 * Interchange if a permutation was applied at the
354 * K-th step of the factorization.
355 *
356  IF( kp.NE.k )
357  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
358  END IF
359  k = k - 1
360 *
361  ELSE
362 *
363 * 2 x 2 pivot block:
364 *
365 * Multiply by the diagonal block if forming L * D.
366 *
367  IF( nounit ) THEN
368  d11 = a( k-1, k-1 )
369  d22 = a( k, k )
370  d21 = a( k, k-1 )
371  d12 = conjg( d21 )
372  DO 50 j = 1, nrhs
373  t1 = b( k-1, j )
374  t2 = b( k, j )
375  b( k-1, j ) = d11*t1 + d12*t2
376  b( k, j ) = d21*t1 + d22*t2
377  50 CONTINUE
378  END IF
379 *
380 * Multiply by P(K) * inv(L(K)) if K < N.
381 *
382  IF( k.NE.n ) THEN
383 *
384 * Apply the transformation.
385 *
386  CALL cgeru( n-k, nrhs, cone, a( k+1, k ), 1,
387  $ b( k, 1 ), ldb, b( k+1, 1 ), ldb )
388  CALL cgeru( n-k, nrhs, cone, a( k+1, k-1 ), 1,
389  $ b( k-1, 1 ), ldb, b( k+1, 1 ), ldb )
390 *
391 * Interchange if a permutation was applied at the
392 * K-th step of the factorization.
393 *
394 *
395 * Swap the second of pair with IMAXth
396 *
397  kp = abs( ipiv( k ) )
398  IF( kp.NE.k )
399  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
400 *
401 * NOW swap the first of pair with Pth
402 *
403  kp = abs( ipiv( k-1 ) )
404  IF( kp.NE.k-1 )
405  $ CALL cswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
406  $ ldb )
407 *
408  END IF
409  k = k - 2
410  END IF
411  GO TO 40
412  60 CONTINUE
413  END IF
414 *--------------------------------------------------
415 *
416 * Compute B := A^H * B (conjugate transpose)
417 *
418 *--------------------------------------------------
419  ELSE
420 *
421 * Form B := U^H*B
422 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
423 * and U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m)
424 *
425  IF( lsame( uplo, 'U' ) ) THEN
426 *
427 * Loop backward applying the transformations.
428 *
429  k = n
430  70 IF( k.LT.1 )
431  $ GO TO 90
432 *
433 * 1 x 1 pivot block.
434 *
435  IF( ipiv( k ).GT.0 ) THEN
436  IF( k.GT.1 ) THEN
437 *
438 * Interchange if P(K) != I.
439 *
440  kp = ipiv( k )
441  IF( kp.NE.k )
442  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
443 *
444 * Apply the transformation
445 * y = y - B' conjg(x),
446 * where x is a column of A and y is a row of B.
447 *
448  CALL clacgv( nrhs, b( k, 1 ), ldb )
449  CALL cgemv( 'Conjugate', k-1, nrhs, cone, b, ldb,
450  $ a( 1, k ), 1, cone, b( k, 1 ), ldb )
451  CALL clacgv( nrhs, b( k, 1 ), ldb )
452  END IF
453  IF( nounit )
454  $ CALL cscal( nrhs, a( k, k ), b( k, 1 ), ldb )
455  k = k - 1
456 *
457 * 2 x 2 pivot block.
458 *
459  ELSE
460  IF( k.GT.2 ) THEN
461 *
462 * Swap the second of pair with Pth
463 *
464  kp = abs( ipiv( k ) )
465  IF( kp.NE.k )
466  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
467 *
468 * Now swap the first of pair with IMAX(r)th
469 *
470  kp = abs( ipiv( k-1 ) )
471  IF( kp.NE.k-1 )
472  $ CALL cswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
473  $ ldb )
474 *
475 * Apply the transformations
476 * y = y - B' conjg(x),
477 * where x is a block column of A and y is a block
478 * row of B.
479 *
480  CALL clacgv( nrhs, b( k, 1 ), ldb )
481  CALL cgemv( 'Conjugate', k-2, nrhs, cone, b, ldb,
482  $ a( 1, k ), 1, cone, b( k, 1 ), ldb )
483  CALL clacgv( nrhs, b( k, 1 ), ldb )
484 *
485  CALL clacgv( nrhs, b( k-1, 1 ), ldb )
486  CALL cgemv( 'Conjugate', k-2, nrhs, cone, b, ldb,
487  $ a( 1, k-1 ), 1, cone, b( k-1, 1 ), ldb )
488  CALL clacgv( nrhs, b( k-1, 1 ), ldb )
489  END IF
490 *
491 * Multiply by the diagonal block if non-unit.
492 *
493  IF( nounit ) THEN
494  d11 = a( k-1, k-1 )
495  d22 = a( k, k )
496  d12 = a( k-1, k )
497  d21 = conjg( d12 )
498  DO 80 j = 1, nrhs
499  t1 = b( k-1, j )
500  t2 = b( k, j )
501  b( k-1, j ) = d11*t1 + d12*t2
502  b( k, j ) = d21*t1 + d22*t2
503  80 CONTINUE
504  END IF
505  k = k - 2
506  END IF
507  GO TO 70
508  90 CONTINUE
509 *
510 * Form B := L^H*B
511 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
512 * and L^H = inv(L^H(m))*P(m)* ... *inv(L^H(1))*P(1)
513 *
514  ELSE
515 *
516 * Loop forward applying the L-transformations.
517 *
518  k = 1
519  100 CONTINUE
520  IF( k.GT.n )
521  $ GO TO 120
522 *
523 * 1 x 1 pivot block
524 *
525  IF( ipiv( k ).GT.0 ) THEN
526  IF( k.LT.n ) THEN
527 *
528 * Interchange if P(K) != I.
529 *
530  kp = ipiv( k )
531  IF( kp.NE.k )
532  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
533 *
534 * Apply the transformation
535 *
536  CALL clacgv( nrhs, b( k, 1 ), ldb )
537  CALL cgemv( 'Conjugate', n-k, nrhs, cone, b( k+1, 1 ),
538  $ ldb, a( k+1, k ), 1, cone, b( k, 1 ), ldb )
539  CALL clacgv( nrhs, b( k, 1 ), ldb )
540  END IF
541  IF( nounit )
542  $ CALL cscal( nrhs, a( k, k ), b( k, 1 ), ldb )
543  k = k + 1
544 *
545 * 2 x 2 pivot block.
546 *
547  ELSE
548  IF( k.LT.n-1 ) THEN
549 *
550 * Swap the first of pair with Pth
551 *
552  kp = abs( ipiv( k ) )
553  IF( kp.NE.k )
554  $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
555 *
556 * Now swap the second of pair with IMAX(r)th
557 *
558  kp = abs( ipiv( k+1 ) )
559  IF( kp.NE.k+1 )
560  $ CALL cswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
561  $ ldb )
562 *
563 * Apply the transformation
564 *
565  CALL clacgv( nrhs, b( k+1, 1 ), ldb )
566  CALL cgemv( 'Conjugate', n-k-1, nrhs, cone,
567  $ b( k+2, 1 ), ldb, a( k+2, k+1 ), 1, cone,
568  $ b( k+1, 1 ), ldb )
569  CALL clacgv( nrhs, b( k+1, 1 ), ldb )
570 *
571  CALL clacgv( nrhs, b( k, 1 ), ldb )
572  CALL cgemv( 'Conjugate', n-k-1, nrhs, cone,
573  $ b( k+2, 1 ), ldb, a( k+2, k ), 1, cone,
574  $ b( k, 1 ), ldb )
575  CALL clacgv( nrhs, b( k, 1 ), ldb )
576  END IF
577 *
578 * Multiply by the diagonal block if non-unit.
579 *
580  IF( nounit ) THEN
581  d11 = a( k, k )
582  d22 = a( k+1, k+1 )
583  d21 = a( k+1, k )
584  d12 = conjg( d21 )
585  DO 110 j = 1, nrhs
586  t1 = b( k, j )
587  t2 = b( k+1, j )
588  b( k, j ) = d11*t1 + d12*t2
589  b( k+1, j ) = d21*t1 + d22*t2
590  110 CONTINUE
591  END IF
592  k = k + 2
593  END IF
594  GO TO 100
595  120 CONTINUE
596  END IF
597 *
598  END IF
599  RETURN
600 *
601 * End of CLAVHE_ROOK
602 *
603  END
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cscal
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54
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
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
cgeru
subroutine cgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CGERU
Definition: cgeru.f:132
clavhe_rook
subroutine clavhe_rook(UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
CLAVHE_ROOK
Definition: clavhe_rook.f:158