LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ clavsp()

subroutine clavsp ( character uplo,
character trans,
character diag,
integer n,
integer nrhs,
complex, dimension( * ) a,
integer, dimension( * ) ipiv,
complex, dimension( ldb, * ) b,
integer ldb,
integer info )

CLAVSP

Purpose:
!>
!>    CLAVSP  performs one of the matrix-vector operations
!>       x := A*x  or  x := A^T*x,
!>    where x is an N element vector and  A is one of the factors
!>    from the symmetric factorization computed by CSPTRF.
!>    CSPTRF produces a factorization of the form
!>         U * D * U^T     or     L * D * L^T,
!>    where U (or L) is a product of permutation and unit upper (lower)
!>    triangular matrices, U^T (or L^T) is the transpose of
!>    U (or L), and D is symmetric and block diagonal with 1 x 1 and
!>    2 x 2 diagonal blocks.  The multipliers for the transformations
!>    and the upper or lower triangular parts of the diagonal blocks
!>    are stored columnwise in packed format in the linear array A.
!>
!>    If TRANS = 'N' or 'n', CLAVSP multiplies either by U or U * D
!>    (or L or L * D).
!>    If TRANS = 'C' or 'c', CLAVSP multiplies either by U^T or D * U^T
!>    (or L^T or D * L^T ).
!> 
!>  UPLO   - CHARACTER*1
!>           On entry, UPLO specifies whether the triangular matrix
!>           stored in A is upper or lower triangular.
!>              UPLO = 'U' or 'u'   The matrix is upper triangular.
!>              UPLO = 'L' or 'l'   The matrix is lower triangular.
!>           Unchanged on exit.
!>
!>  TRANS  - CHARACTER*1
!>           On entry, TRANS specifies the operation to be performed as
!>           follows:
!>              TRANS = 'N' or 'n'   x := A*x.
!>              TRANS = 'T' or 't'   x := A^T*x.
!>           Unchanged on exit.
!>
!>  DIAG   - CHARACTER*1
!>           On entry, DIAG specifies whether the diagonal blocks are
!>           assumed to be unit matrices, as follows:
!>              DIAG = 'U' or 'u'   Diagonal blocks are unit matrices.
!>              DIAG = 'N' or 'n'   Diagonal blocks are non-unit.
!>           Unchanged on exit.
!>
!>  N      - INTEGER
!>           On entry, N specifies the order of the matrix A.
!>           N must be at least zero.
!>           Unchanged on exit.
!>
!>  NRHS   - INTEGER
!>           On entry, NRHS specifies the number of right hand sides,
!>           i.e., the number of vectors x to be multiplied by A.
!>           NRHS must be at least zero.
!>           Unchanged on exit.
!>
!>  A      - COMPLEX array, dimension( N*(N+1)/2 )
!>           On entry, A contains a block diagonal matrix and the
!>           multipliers of the transformations used to obtain it,
!>           stored as a packed triangular matrix.
!>           Unchanged on exit.
!>
!>  IPIV   - INTEGER array, dimension( N )
!>           On entry, IPIV contains the vector of pivot indices as
!>           determined by CSPTRF.
!>           If IPIV( K ) = K, no interchange was done.
!>           If IPIV( K ) <> K but IPIV( K ) > 0, then row K was inter-
!>           changed with row IPIV( K ) and a 1 x 1 pivot block was used.
!>           If IPIV( K ) < 0 and UPLO = 'U', then row K-1 was exchanged
!>           with row | IPIV( K ) | and a 2 x 2 pivot block was used.
!>           If IPIV( K ) < 0 and UPLO = 'L', then row K+1 was exchanged
!>           with row | IPIV( K ) | and a 2 x 2 pivot block was used.
!>
!>  B      - COMPLEX array, dimension( LDB, NRHS )
!>           On entry, B contains NRHS vectors of length N.
!>           On exit, B is overwritten with the product A * B.
!>
!>  LDB    - INTEGER
!>           On entry, LDB contains the leading dimension of B as
!>           declared in the calling program.  LDB must be at least
!>           max( 1, N ).
!>           Unchanged on exit.
!>
!>  INFO   - INTEGER
!>           INFO is the error flag.
!>           On exit, a value of 0 indicates a successful exit.
!>           A negative value, say -K, indicates that the K-th argument
!>           has an illegal value.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 129 of file clavsp.f.

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, cscal, cswap, xerbla
162* ..
163* .. Intrinsic Functions ..
164 INTRINSIC abs, 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 ELSE IF( .NOT.lsame( trans, 'N' ) .AND. .NOT.lsame( trans, 'T' ) )
174 $ THEN
175 info = -2
176 ELSE IF( .NOT.lsame( diag, 'U' ) .AND. .NOT.lsame( diag, 'N' ) )
177 $ THEN
178 info = -3
179 ELSE IF( n.LT.0 ) THEN
180 info = -4
181 ELSE IF( ldb.LT.max( 1, n ) ) THEN
182 info = -8
183 END IF
184 IF( info.NE.0 ) THEN
185 CALL xerbla( 'CLAVSP ', -info )
186 RETURN
187 END IF
188*
189* Quick return if possible.
190*
191 IF( n.EQ.0 )
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 IF( k.GT.n )
213 $ GO TO 30
214*
215* 1 x 1 pivot block
216*
217 IF( ipiv( k ).GT.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 IF( k.GT.1 ) 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 IF( kp.NE.k )
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 = 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 IF( k.GT.1 ) 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 IF( kp.NE.k )
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 IF( k.LT.1 )
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 IF( ipiv( k ).GT.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 IF( k.NE.n ) 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 IF( kp.NE.k )
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 = 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 IF( k.NE.n ) 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 IF( kp.NE.k )
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^T * B (transpose)
377*
378*-------------------------------------------------
379 ELSE
380*
381* Form B := U^T*B
382* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
383* and U^T = inv(U^T(1))*P(1)* ... *inv(U^T(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 70 IF( k.LT.1 )
392 $ GO TO 90
393 kc = kc - k
394*
395* 1 x 1 pivot block.
396*
397 IF( ipiv( k ).GT.0 ) THEN
398 IF( k.GT.1 ) THEN
399*
400* Interchange if P(K) != I.
401*
402 kp = ipiv( k )
403 IF( kp.NE.k )
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 cgemv( 'Transpose', k-1, nrhs, one, b, ldb,
411 $ a( kc ), 1, one, b( k, 1 ), ldb )
412 END IF
413 IF( nounit )
414 $ CALL cscal( nrhs, a( kc+k-1 ), b( k, 1 ), ldb )
415 k = k - 1
416*
417* 2 x 2 pivot block.
418*
419 ELSE
420 kcnext = kc - ( k-1 )
421 IF( k.GT.2 ) THEN
422*
423* Interchange if P(K) != I.
424*
425 kp = abs( ipiv( k ) )
426 IF( kp.NE.k-1 )
427 $ CALL cswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
428 $ ldb )
429*
430* Apply the transformations.
431*
432 CALL cgemv( 'Transpose', k-2, nrhs, one, b, ldb,
433 $ a( kc ), 1, one, b( k, 1 ), ldb )
434*
435 CALL cgemv( 'Transpose', k-2, nrhs, one, b, ldb,
436 $ a( kcnext ), 1, one, b( k-1, 1 ), ldb )
437 END IF
438*
439* Multiply by the diagonal block if non-unit.
440*
441 IF( nounit ) THEN
442 d11 = a( kc-1 )
443 d22 = a( kc+k-1 )
444 d12 = a( kc+k-2 )
445 d21 = d12
446 DO 80 j = 1, nrhs
447 t1 = b( k-1, j )
448 t2 = b( k, j )
449 b( k-1, j ) = d11*t1 + d12*t2
450 b( k, j ) = d21*t1 + d22*t2
451 80 CONTINUE
452 END IF
453 kc = kcnext
454 k = k - 2
455 END IF
456 GO TO 70
457 90 CONTINUE
458*
459* Form B := L^T*B
460* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
461* and L^T = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
462*
463 ELSE
464*
465* Loop forward applying the L-transformations.
466*
467 k = 1
468 kc = 1
469 100 CONTINUE
470 IF( k.GT.n )
471 $ GO TO 120
472*
473* 1 x 1 pivot block
474*
475 IF( ipiv( k ).GT.0 ) THEN
476 IF( k.LT.n ) THEN
477*
478* Interchange if P(K) != I.
479*
480 kp = ipiv( k )
481 IF( kp.NE.k )
482 $ CALL cswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
483*
484* Apply the transformation
485*
486 CALL cgemv( 'Transpose', n-k, nrhs, one, b( k+1, 1 ),
487 $ ldb, a( kc+1 ), 1, one, b( k, 1 ), ldb )
488 END IF
489 IF( nounit )
490 $ CALL cscal( nrhs, a( kc ), b( k, 1 ), ldb )
491 kc = kc + n - k + 1
492 k = k + 1
493*
494* 2 x 2 pivot block.
495*
496 ELSE
497 kcnext = kc + n - k + 1
498 IF( k.LT.n-1 ) THEN
499*
500* Interchange if P(K) != I.
501*
502 kp = abs( ipiv( k ) )
503 IF( kp.NE.k+1 )
504 $ CALL cswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
505 $ ldb )
506*
507* Apply the transformation
508*
509 CALL cgemv( 'Transpose', n-k-1, nrhs, one,
510 $ b( k+2, 1 ), ldb, a( kcnext+1 ), 1, one,
511 $ b( k+1, 1 ), ldb )
512*
513 CALL cgemv( 'Transpose', n-k-1, nrhs, one,
514 $ b( k+2, 1 ), ldb, a( kc+2 ), 1, one,
515 $ b( k, 1 ), ldb )
516 END IF
517*
518* Multiply by the diagonal block if non-unit.
519*
520 IF( nounit ) THEN
521 d11 = a( kc )
522 d22 = a( kcnext )
523 d21 = a( kc+1 )
524 d12 = d21
525 DO 110 j = 1, nrhs
526 t1 = b( k, j )
527 t2 = b( k+1, j )
528 b( k, j ) = d11*t1 + d12*t2
529 b( k+1, j ) = d21*t1 + d22*t2
530 110 CONTINUE
531 END IF
532 kc = kcnext + ( n-k )
533 k = k + 2
534 END IF
535 GO TO 100
536 120 CONTINUE
537 END IF
538*
539 END IF
540 RETURN
541*
542* End of CLAVSP
543*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
cgeru
subroutine cgeru(m, n, alpha, x, incx, y, incy, a, lda)
CGERU
Definition cgeru.f:130
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: