OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
clavhp.f
Go to the documentation of this file.
1*> \brief \b CLAVHP
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 CLAVHP( UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB,
12* INFO )
13*
14* .. Scalar Arguments ..
15* CHARACTER DIAG, TRANS, UPLO
16* INTEGER INFO, LDB, N, NRHS
17* ..
18* .. Array Arguments ..
19* INTEGER IPIV( * )
20* COMPLEX A( * ), B( LDB, * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> CLAVHP 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 symmetric factorization computed by CHPTRF.
33*> CHPTRF produces a factorization of the form
34*> U * D * U^H or L * D * L^H,
35*> where U (or L) is a product of permutation and unit upper (lower)
36*> triangular matrices, U^H (or L^H) is the conjugate transpose of
37*> U (or L), and D is Hermitian and block diagonal with 1 x 1 and
38*> 2 x 2 diagonal blocks. The multipliers for the transformations
39*> and the upper or lower triangular parts of the diagonal blocks
40*> are stored columnwise in packed format in the linear array A.
41*>
42*> If TRANS = 'N' or 'n', CLAVHP multiplies either by U or U * D
43*> (or L or L * D).
44*> If TRANS = 'C' or 'c', CLAVHP multiplies either by U^H or D * U^H
45*> (or L^H or D * L^H ).
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \verbatim
52*> UPLO - CHARACTER*1
53*> On entry, UPLO specifies whether the triangular matrix
54*> stored in A is upper or lower triangular.
55*> UPLO = 'U' or 'u' The matrix is upper triangular.
56*> UPLO = 'L' or 'l' The matrix is lower triangular.
57*> Unchanged on exit.
58*>
59*> TRANS - CHARACTER*1
60*> On entry, TRANS specifies the operation to be performed as
61*> follows:
62*> TRANS = 'N' or 'n' x := A*x.
63*> TRANS = 'C' or 'c' x := A^H*x.
64*> Unchanged on exit.
65*>
66*> DIAG - CHARACTER*1
67*> On entry, DIAG specifies whether the diagonal blocks are
68*> assumed to be unit matrices, as follows:
69*> DIAG = 'U' or 'u' Diagonal blocks are unit matrices.
70*> DIAG = 'N' or 'n' Diagonal blocks are non-unit.
71*> Unchanged on exit.
72*>
73*> N - INTEGER
74*> On entry, N specifies the order of the matrix A.
75*> N must be at least zero.
76*> Unchanged on exit.
77*>
78*> NRHS - INTEGER
79*> On entry, NRHS specifies the number of right hand sides,
80*> i.e., the number of vectors x to be multiplied by A.
81*> NRHS must be at least zero.
82*> Unchanged on exit.
83*>
84*> A - COMPLEX array, dimension( N*(N+1)/2 )
85*> On entry, A contains a block diagonal matrix and the
86*> multipliers of the transformations used to obtain it,
87*> stored as a packed triangular matrix.
88*> Unchanged on exit.
89*>
90*> IPIV - INTEGER array, dimension( N )
91*> On entry, IPIV contains the vector of pivot indices as
92*> determined by CSPTRF or CHPTRF.
93*> If IPIV( K ) = K, no interchange was done.
94*> If IPIV( K ) <> K but IPIV( K ) > 0, then row K was inter-
95*> changed with row IPIV( K ) and a 1 x 1 pivot block was used.
96*> If IPIV( K ) < 0 and UPLO = 'U', then row K-1 was exchanged
97*> with row | IPIV( K ) | and a 2 x 2 pivot block was used.
98*> If IPIV( K ) < 0 and UPLO = 'L', then row K+1 was exchanged
99*> with row | IPIV( K ) | and a 2 x 2 pivot block was used.
100*>
101*> B - COMPLEX array, dimension( LDB, NRHS )
102*> On entry, B contains NRHS vectors of length N.
103*> On exit, B is overwritten with the product A * B.
104*>
105*> LDB - INTEGER
106*> On entry, LDB contains the leading dimension of B as
107*> declared in the calling program. LDB must be at least
108*> max( 1, N ).
109*> Unchanged on exit.
110*>
111*> INFO - INTEGER
112*> INFO is the error flag.
113*> On exit, a value of 0 indicates a successful exit.
114*> A negative value, say -K, indicates that the K-th argument
115*> has an illegal value.
116*> \endverbatim
117*
118* Authors:
119* ========
120*
121*> \author Univ. of Tennessee
122*> \author Univ. of California Berkeley
123*> \author Univ. of Colorado Denver
124*> \author NAG Ltd.
125*
126*> \ingroup complex_lin
127*
128* =====================================================================
129 SUBROUTINE clavhp( UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB,
130 $ INFO )
131*
132* -- LAPACK test routine --
133* -- LAPACK is a software package provided by Univ. of Tennessee, --
134* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
135*
136* .. Scalar Arguments ..
137 CHARACTER DIAG, TRANS, UPLO
138 INTEGER INFO, LDB, N, NRHS
139* ..
140* .. Array Arguments ..
141 INTEGER IPIV( * )
142 COMPLEX A( * ), B( LDB, * )
143* ..
144*
145* =====================================================================
146*
147* .. Parameters ..
148 COMPLEX ONE
149 parameter( one = ( 1.0e+0, 0.0e+0 ) )
150* ..
151* .. Local Scalars ..
152 LOGICAL NOUNIT
153 INTEGER J, K, KC, KCNEXT, KP
154 COMPLEX D11, D12, D21, D22, T1, T2
155* ..
156* .. External Functions ..
157 LOGICAL LSAME
158 EXTERNAL lsame
159* ..
160* .. External Subroutines ..
161 EXTERNAL cgemv, cgeru, clacgv, cscal, cswap, xerbla
162* ..
163* .. Intrinsic Functions ..
164 INTRINSIC abs, conjg, max
165* ..
166* .. Executable Statements ..
167*
168* Test the input parameters.
169*
170 info = 0
171 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'l' ) ) THEN
172 INFO = -1
173.NOT. ELSE IF( LSAME( TRANS, 'n.AND..NOT.' ) LSAME( TRANS, 'c' ) )
174 $ THEN
175 INFO = -2
176.NOT. ELSE IF( LSAME( DIAG, 'u.AND..NOT.' ) LSAME( DIAG, 'n' ) )
177 $ THEN
178 INFO = -3
179.LT. ELSE IF( N0 ) THEN
180 INFO = -4
181.LT. ELSE IF( LDBMAX( 1, N ) ) THEN
182 INFO = -8
183 END IF
184.NE. IF( INFO0 ) THEN
185 CALL XERBLA( 'clavhp ', -INFO )
186 RETURN
187 END IF
188*
189* Quick return if possible.
190*
191.EQ. IF( N0 )
192 $ RETURN
193*
194 NOUNIT = LSAME( DIAG, 'n' )
195*------------------------------------------
196*
197* Compute B := A * B (No transpose)
198*
199*------------------------------------------
200 IF( LSAME( TRANS, 'n' ) ) THEN
201*
202* Compute B := U*B
203* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
204*
205 IF( LSAME( UPLO, 'u' ) ) THEN
206*
207* Loop forward applying the transformations.
208*
209 K = 1
210 KC = 1
211 10 CONTINUE
212.GT. IF( KN )
213 $ GO TO 30
214*
215* 1 x 1 pivot block
216*
217.GT. IF( IPIV( K )0 ) THEN
218*
219* Multiply by the diagonal element if forming U * D.
220*
221 IF( NOUNIT )
222 $ CALL CSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
223*
224* Multiply by P(K) * inv(U(K)) if K > 1.
225*
226.GT. IF( K1 ) THEN
227*
228* Apply the transformation.
229*
230 CALL CGERU( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ),
231 $ LDB, B( 1, 1 ), LDB )
232*
233* Interchange if P(K) != I.
234*
235 KP = IPIV( K )
236.NE. IF( KPK )
237 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
238 END IF
239 KC = KC + K
240 K = K + 1
241 ELSE
242*
243* 2 x 2 pivot block
244*
245 KCNEXT = KC + K
246*
247* Multiply by the diagonal block if forming U * D.
248*
249 IF( NOUNIT ) THEN
250 D11 = A( KCNEXT-1 )
251 D22 = A( KCNEXT+K )
252 D12 = A( KCNEXT+K-1 )
253 D21 = CONJG( D12 )
254 DO 20 J = 1, NRHS
255 T1 = B( K, J )
256 T2 = B( K+1, J )
257 B( K, J ) = D11*T1 + D12*T2
258 B( K+1, J ) = D21*T1 + D22*T2
259 20 CONTINUE
260 END IF
261*
262* Multiply by P(K) * inv(U(K)) if K > 1.
263*
264.GT. IF( K1 ) THEN
265*
266* Apply the transformations.
267*
268 CALL CGERU( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ),
269 $ LDB, B( 1, 1 ), LDB )
270 CALL CGERU( K-1, NRHS, ONE, A( KCNEXT ), 1,
271 $ B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
272*
273* Interchange if P(K) != I.
274*
275 KP = ABS( IPIV( K ) )
276.NE. IF( KPK )
277 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
278 END IF
279 KC = KCNEXT + K + 1
280 K = K + 2
281 END IF
282 GO TO 10
283 30 CONTINUE
284*
285* Compute B := L*B
286* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
287*
288 ELSE
289*
290* Loop backward applying the transformations to B.
291*
292 K = N
293 KC = N*( N+1 ) / 2 + 1
294 40 CONTINUE
295.LT. IF( K1 )
296 $ GO TO 60
297 KC = KC - ( N-K+1 )
298*
299* Test the pivot index. If greater than zero, a 1 x 1
300* pivot was used, otherwise a 2 x 2 pivot was used.
301*
302.GT. IF( IPIV( K )0 ) THEN
303*
304* 1 x 1 pivot block:
305*
306* Multiply by the diagonal element if forming L * D.
307*
308 IF( NOUNIT )
309 $ CALL CSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
310*
311* Multiply by P(K) * inv(L(K)) if K < N.
312*
313.NE. IF( KN ) THEN
314 KP = IPIV( K )
315*
316* Apply the transformation.
317*
318 CALL CGERU( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
319 $ LDB, B( K+1, 1 ), LDB )
320*
321* Interchange if a permutation was applied at the
322* K-th step of the factorization.
323*
324.NE. IF( KPK )
325 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
326 END IF
327 K = K - 1
328*
329 ELSE
330*
331* 2 x 2 pivot block:
332*
333 KCNEXT = KC - ( N-K+2 )
334*
335* Multiply by the diagonal block if forming L * D.
336*
337 IF( NOUNIT ) THEN
338 D11 = A( KCNEXT )
339 D22 = A( KC )
340 D21 = A( KCNEXT+1 )
341 D12 = CONJG( D21 )
342 DO 50 J = 1, NRHS
343 T1 = B( K-1, J )
344 T2 = B( K, J )
345 B( K-1, J ) = D11*T1 + D12*T2
346 B( K, J ) = D21*T1 + D22*T2
347 50 CONTINUE
348 END IF
349*
350* Multiply by P(K) * inv(L(K)) if K < N.
351*
352.NE. IF( KN ) THEN
353*
354* Apply the transformation.
355*
356 CALL CGERU( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
357 $ LDB, B( K+1, 1 ), LDB )
358 CALL CGERU( N-K, NRHS, ONE, A( KCNEXT+2 ), 1,
359 $ B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
360*
361* Interchange if a permutation was applied at the
362* K-th step of the factorization.
363*
364 KP = ABS( IPIV( K ) )
365.NE. IF( KPK )
366 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
367 END IF
368 KC = KCNEXT
369 K = K - 2
370 END IF
371 GO TO 40
372 60 CONTINUE
373 END IF
374*-------------------------------------------------
375*
376* Compute B := A^H * B (conjugate transpose)
377*
378*-------------------------------------------------
379 ELSE
380*
381* Form B := U^H*B
382* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
383* and U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m)
384*
385 IF( LSAME( UPLO, 'u' ) ) THEN
386*
387* Loop backward applying the transformations.
388*
389 K = N
390 KC = N*( N+1 ) / 2 + 1
391.LT. 70 IF( K1 )
392 $ GO TO 90
393 KC = KC - K
394*
395* 1 x 1 pivot block.
396*
397.GT. IF( IPIV( K )0 ) THEN
398.GT. IF( K1 ) THEN
399*
400* Interchange if P(K) != I.
401*
402 KP = IPIV( K )
403.NE. IF( KPK )
404 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
405*
406* Apply the transformation:
407* y := y - B' * conjg(x)
408* where x is a column of A and y is a row of B.
409*
410 CALL CLACGV( NRHS, B( K, 1 ), LDB )
411 CALL CGEMV( 'conjugate', K-1, NRHS, ONE, B, LDB,
412 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
413 CALL CLACGV( NRHS, B( K, 1 ), LDB )
414 END IF
415 IF( NOUNIT )
416 $ CALL CSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
417 K = K - 1
418*
419* 2 x 2 pivot block.
420*
421 ELSE
422 KCNEXT = KC - ( K-1 )
423.GT. IF( K2 ) THEN
424*
425* Interchange if P(K) != I.
426*
427 KP = ABS( IPIV( K ) )
428.NE. IF( KPK-1 )
429 $ CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
430 $ LDB )
431*
432* Apply the transformations.
433*
434 CALL CLACGV( NRHS, B( K, 1 ), LDB )
435 CALL CGEMV( 'conjugate', K-2, NRHS, ONE, B, LDB,
436 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
437 CALL CLACGV( NRHS, B( K, 1 ), LDB )
438*
439 CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
440 CALL CGEMV( 'conjugate', K-2, NRHS, ONE, B, LDB,
441 $ A( KCNEXT ), 1, ONE, B( K-1, 1 ), LDB )
442 CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
443 END IF
444*
445* Multiply by the diagonal block if non-unit.
446*
447 IF( NOUNIT ) THEN
448 D11 = A( KC-1 )
449 D22 = A( KC+K-1 )
450 D12 = A( KC+K-2 )
451 D21 = CONJG( D12 )
452 DO 80 J = 1, NRHS
453 T1 = B( K-1, J )
454 T2 = B( K, J )
455 B( K-1, J ) = D11*T1 + D12*T2
456 B( K, J ) = D21*T1 + D22*T2
457 80 CONTINUE
458 END IF
459 KC = KCNEXT
460 K = K - 2
461 END IF
462 GO TO 70
463 90 CONTINUE
464*
465* Form B := L^H*B
466* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
467* and L^H = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
468*
469 ELSE
470*
471* Loop forward applying the L-transformations.
472*
473 K = 1
474 KC = 1
475 100 CONTINUE
476.GT. IF( KN )
477 $ GO TO 120
478*
479* 1 x 1 pivot block
480*
481.GT. IF( IPIV( K )0 ) THEN
482.LT. IF( KN ) THEN
483*
484* Interchange if P(K) != I.
485*
486 KP = IPIV( K )
487.NE. IF( KPK )
488 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
489*
490* Apply the transformation
491*
492 CALL CLACGV( NRHS, B( K, 1 ), LDB )
493 CALL CGEMV( 'conjugate', N-K, NRHS, ONE, B( K+1, 1 ),
494 $ LDB, A( KC+1 ), 1, ONE, B( K, 1 ), LDB )
495 CALL CLACGV( NRHS, B( K, 1 ), LDB )
496 END IF
497 IF( NOUNIT )
498 $ CALL CSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
499 KC = KC + N - K + 1
500 K = K + 1
501*
502* 2 x 2 pivot block.
503*
504 ELSE
505 KCNEXT = KC + N - K + 1
506.LT. IF( KN-1 ) THEN
507*
508* Interchange if P(K) != I.
509*
510 KP = ABS( IPIV( K ) )
511.NE. IF( KPK+1 )
512 $ CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
513 $ LDB )
514*
515* Apply the transformation
516*
517 CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
518 CALL CGEMV( 'conjugate', N-K-1, NRHS, ONE,
519 $ B( K+2, 1 ), LDB, A( KCNEXT+1 ), 1, ONE,
520 $ B( K+1, 1 ), LDB )
521 CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
522*
523 CALL CLACGV( NRHS, B( K, 1 ), LDB )
524 CALL CGEMV( 'conjugate', N-K-1, NRHS, ONE,
525 $ B( K+2, 1 ), LDB, A( KC+2 ), 1, ONE,
526 $ B( K, 1 ), LDB )
527 CALL CLACGV( NRHS, B( K, 1 ), LDB )
528 END IF
529*
530* Multiply by the diagonal block if non-unit.
531*
532 IF( NOUNIT ) THEN
533 D11 = A( KC )
534 D22 = A( KCNEXT )
535 D21 = A( KC+1 )
536 D12 = CONJG( D21 )
537 DO 110 J = 1, NRHS
538 T1 = B( K, J )
539 T2 = B( K+1, J )
540 B( K, J ) = D11*T1 + D12*T2
541 B( K+1, J ) = D21*T1 + D22*T2
542 110 CONTINUE
543 END IF
544 KC = KCNEXT + ( N-K )
545 K = K + 2
546 END IF
547 GO TO 100
548 120 CONTINUE
549 END IF
550*
551 END IF
552 RETURN
553*
554* End of CLAVHP
555*
556 END
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
clacgv
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:74
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:158
cgeru
subroutine cgeru(m, n, alpha, x, incx, y, incy, a, lda)
CGERU
Definition cgeru.f:130
clavhp
subroutine clavhp(uplo, trans, diag, n, nrhs, a, ipiv, b, ldb, info)
CLAVHP
Definition clavhp.f:131
max
#define max(a, b)
Definition macros.h:21