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

◆ dlavsy()

subroutine dlavsy ( character  uplo,
character  trans,
character  diag,
integer  n,
integer  nrhs,
double precision, dimension( lda, * )  a,
integer  lda,
integer, dimension( * )  ipiv,
double precision, dimension( ldb, * )  b,
integer  ldb,
integer  info 
)

DLAVSY

Purpose:
 DLAVSY  performs one of the matrix-vector operations
    x := A*x  or  x := A'*x,
 where x is an N element vector and A is one of the factors
 from the block U*D*U' or L*D*L' factorization computed by DSYTRF.

 If TRANS = 'N', multiplies by U  or U * D  (or L  or L * D)
 If TRANS = 'T', multiplies by U' or D * U' (or L' or D * L')
 If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L')
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the factor stored in A is upper or lower
          triangular.
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]TRANS
          TRANS is CHARACTER*1
          Specifies the operation to be performed:
          = 'N':  x := A*x
          = 'T':  x := A'*x
          = 'C':  x := A'*x
[in]DIAG
          DIAG is CHARACTER*1
          Specifies whether or not the diagonal blocks are unit
          matrices.  If the diagonal blocks are assumed to be unit,
          then A = U or A = L, otherwise A = U*D or A = L*D.
          = 'U':  Diagonal blocks are assumed to be unit matrices.
          = 'N':  Diagonal blocks are assumed to be non-unit matrices.
[in]N
          N is INTEGER
          The number of rows and columns of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of vectors
          x to be multiplied by A.  NRHS >= 0.
[in]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L as computed by DSYTRF.
          Stored as a 2-D triangular matrix.
[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 DSYTRF.

          If UPLO = 'U':
               If IPIV(k) > 0, then rows and columns k and IPIV(k)
               were interchanged and D(k,k) is a 1-by-1 diagonal block.
               (If IPIV( k ) = k, no interchange was done).

               If IPIV(k) = IPIV(k-1) < 0, then rows and
               columns k-1 and -IPIV(k) were interchanged,
               D(k-1:k,k-1:k) is a 2-by-2 diagonal block.

          If UPLO = 'L':
               If IPIV(k) > 0, then rows and columns k and IPIV(k)
               were interchanged and D(k,k) is a 1-by-1 diagonal block.
               (If IPIV( k ) = k, no interchange was done).

               If IPIV(k) = IPIV(k+1) < 0, then rows and
               columns k+1 and -IPIV(k) were interchanged,
               D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[in,out]B
          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
          On entry, B contains NRHS vectors of length N.
          On exit, B is overwritten with the product A * B.
[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 = -k, the k-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 153 of file dlavsy.f.

155*
156* -- LAPACK test routine --
157* -- LAPACK is a software package provided by Univ. of Tennessee, --
158* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
159*
160* .. Scalar Arguments ..
161 CHARACTER DIAG, TRANS, UPLO
162 INTEGER INFO, LDA, LDB, N, NRHS
163* ..
164* .. Array Arguments ..
165 INTEGER IPIV( * )
166 DOUBLE PRECISION A( LDA, * ), B( LDB, * )
167* ..
168*
169* =====================================================================
170*
171* .. Parameters ..
172 DOUBLE PRECISION ONE
173 parameter( one = 1.0d+0 )
174* ..
175* .. Local Scalars ..
176 LOGICAL NOUNIT
177 INTEGER J, K, KP
178 DOUBLE PRECISION D11, D12, D21, D22, T1, T2
179* ..
180* .. External Functions ..
181 LOGICAL LSAME
182 EXTERNAL lsame
183* ..
184* .. External Subroutines ..
185 EXTERNAL dgemv, dger, dscal, dswap, xerbla
186* ..
187* .. Intrinsic Functions ..
188 INTRINSIC abs, max
189* ..
190* .. Executable Statements ..
191*
192* Test the input parameters.
193*
194 info = 0
195 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
196 info = -1
197 ELSE IF( .NOT.lsame( trans, 'N' ) .AND. .NOT.
198 $ lsame( trans, 'T' ) .AND. .NOT.lsame( trans, 'C' ) ) THEN
199 info = -2
200 ELSE IF( .NOT.lsame( diag, 'U' ) .AND. .NOT.lsame( diag, 'N' ) )
201 $ THEN
202 info = -3
203 ELSE IF( n.LT.0 ) THEN
204 info = -4
205 ELSE IF( lda.LT.max( 1, n ) ) THEN
206 info = -6
207 ELSE IF( ldb.LT.max( 1, n ) ) THEN
208 info = -9
209 END IF
210 IF( info.NE.0 ) THEN
211 CALL xerbla( 'DLAVSY ', -info )
212 RETURN
213 END IF
214*
215* Quick return if possible.
216*
217 IF( n.EQ.0 )
218 $ RETURN
219*
220 nounit = lsame( diag, 'N' )
221*------------------------------------------
222*
223* Compute B := A * B (No transpose)
224*
225*------------------------------------------
226 IF( lsame( trans, 'N' ) ) THEN
227*
228* Compute B := U*B
229* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
230*
231 IF( lsame( uplo, 'U' ) ) THEN
232*
233* Loop forward applying the transformations.
234*
235 k = 1
236 10 CONTINUE
237 IF( k.GT.n )
238 $ GO TO 30
239 IF( ipiv( k ).GT.0 ) THEN
240*
241* 1 x 1 pivot block
242*
243* Multiply by the diagonal element if forming U * D.
244*
245 IF( nounit )
246 $ CALL dscal( nrhs, a( k, k ), b( k, 1 ), ldb )
247*
248* Multiply by P(K) * inv(U(K)) if K > 1.
249*
250 IF( k.GT.1 ) THEN
251*
252* Apply the transformation.
253*
254 CALL dger( k-1, nrhs, one, a( 1, k ), 1, b( k, 1 ),
255 $ ldb, b( 1, 1 ), ldb )
256*
257* Interchange if P(K) .ne. I.
258*
259 kp = ipiv( k )
260 IF( kp.NE.k )
261 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
262 END IF
263 k = k + 1
264 ELSE
265*
266* 2 x 2 pivot block
267*
268* Multiply by the diagonal block if forming U * D.
269*
270 IF( nounit ) THEN
271 d11 = a( k, k )
272 d22 = a( k+1, k+1 )
273 d12 = a( k, k+1 )
274 d21 = d12
275 DO 20 j = 1, nrhs
276 t1 = b( k, j )
277 t2 = b( k+1, j )
278 b( k, j ) = d11*t1 + d12*t2
279 b( k+1, j ) = d21*t1 + d22*t2
280 20 CONTINUE
281 END IF
282*
283* Multiply by P(K) * inv(U(K)) if K > 1.
284*
285 IF( k.GT.1 ) THEN
286*
287* Apply the transformations.
288*
289 CALL dger( k-1, nrhs, one, a( 1, k ), 1, b( k, 1 ),
290 $ ldb, b( 1, 1 ), ldb )
291 CALL dger( k-1, nrhs, one, a( 1, k+1 ), 1,
292 $ b( k+1, 1 ), ldb, b( 1, 1 ), ldb )
293*
294* Interchange if P(K) .ne. I.
295*
296 kp = abs( ipiv( k ) )
297 IF( kp.NE.k )
298 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
299 END IF
300 k = k + 2
301 END IF
302 GO TO 10
303 30 CONTINUE
304*
305* Compute B := L*B
306* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
307*
308 ELSE
309*
310* Loop backward applying the transformations to B.
311*
312 k = n
313 40 CONTINUE
314 IF( k.LT.1 )
315 $ GO TO 60
316*
317* Test the pivot index. If greater than zero, a 1 x 1
318* pivot was used, otherwise a 2 x 2 pivot was used.
319*
320 IF( ipiv( k ).GT.0 ) THEN
321*
322* 1 x 1 pivot block:
323*
324* Multiply by the diagonal element if forming L * D.
325*
326 IF( nounit )
327 $ CALL dscal( nrhs, a( k, k ), b( k, 1 ), ldb )
328*
329* Multiply by P(K) * inv(L(K)) if K < N.
330*
331 IF( k.NE.n ) THEN
332 kp = ipiv( k )
333*
334* Apply the transformation.
335*
336 CALL dger( n-k, nrhs, one, a( k+1, k ), 1, b( k, 1 ),
337 $ ldb, b( k+1, 1 ), ldb )
338*
339* Interchange if a permutation was applied at the
340* K-th step of the factorization.
341*
342 IF( kp.NE.k )
343 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
344 END IF
345 k = k - 1
346*
347 ELSE
348*
349* 2 x 2 pivot block:
350*
351* Multiply by the diagonal block if forming L * D.
352*
353 IF( nounit ) THEN
354 d11 = a( k-1, k-1 )
355 d22 = a( k, k )
356 d21 = a( k, k-1 )
357 d12 = d21
358 DO 50 j = 1, nrhs
359 t1 = b( k-1, j )
360 t2 = b( k, j )
361 b( k-1, j ) = d11*t1 + d12*t2
362 b( k, j ) = d21*t1 + d22*t2
363 50 CONTINUE
364 END IF
365*
366* Multiply by P(K) * inv(L(K)) if K < N.
367*
368 IF( k.NE.n ) THEN
369*
370* Apply the transformation.
371*
372 CALL dger( n-k, nrhs, one, a( k+1, k ), 1, b( k, 1 ),
373 $ ldb, b( k+1, 1 ), ldb )
374 CALL dger( n-k, nrhs, one, a( k+1, k-1 ), 1,
375 $ b( k-1, 1 ), ldb, b( k+1, 1 ), ldb )
376*
377* Interchange if a permutation was applied at the
378* K-th step of the factorization.
379*
380 kp = abs( ipiv( k ) )
381 IF( kp.NE.k )
382 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
383 END IF
384 k = k - 2
385 END IF
386 GO TO 40
387 60 CONTINUE
388 END IF
389*----------------------------------------
390*
391* Compute B := A' * B (transpose)
392*
393*----------------------------------------
394 ELSE
395*
396* Form B := U'*B
397* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
398* and U' = inv(U'(1))*P(1)* ... *inv(U'(m))*P(m)
399*
400 IF( lsame( uplo, 'U' ) ) THEN
401*
402* Loop backward applying the transformations.
403*
404 k = n
405 70 CONTINUE
406 IF( k.LT.1 )
407 $ GO TO 90
408*
409* 1 x 1 pivot block.
410*
411 IF( ipiv( k ).GT.0 ) THEN
412 IF( k.GT.1 ) THEN
413*
414* Interchange if P(K) .ne. I.
415*
416 kp = ipiv( k )
417 IF( kp.NE.k )
418 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
419*
420* Apply the transformation
421*
422 CALL dgemv( 'Transpose', k-1, nrhs, one, b, ldb,
423 $ a( 1, k ), 1, one, b( k, 1 ), ldb )
424 END IF
425 IF( nounit )
426 $ CALL dscal( nrhs, a( k, k ), b( k, 1 ), ldb )
427 k = k - 1
428*
429* 2 x 2 pivot block.
430*
431 ELSE
432 IF( k.GT.2 ) THEN
433*
434* Interchange if P(K) .ne. I.
435*
436 kp = abs( ipiv( k ) )
437 IF( kp.NE.k-1 )
438 $ CALL dswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
439 $ ldb )
440*
441* Apply the transformations
442*
443 CALL dgemv( 'Transpose', k-2, nrhs, one, b, ldb,
444 $ a( 1, k ), 1, one, b( k, 1 ), ldb )
445 CALL dgemv( 'Transpose', k-2, nrhs, one, b, ldb,
446 $ a( 1, k-1 ), 1, one, b( k-1, 1 ), ldb )
447 END IF
448*
449* Multiply by the diagonal block if non-unit.
450*
451 IF( nounit ) THEN
452 d11 = a( k-1, k-1 )
453 d22 = a( k, k )
454 d12 = a( k-1, k )
455 d21 = d12
456 DO 80 j = 1, nrhs
457 t1 = b( k-1, j )
458 t2 = b( k, j )
459 b( k-1, j ) = d11*t1 + d12*t2
460 b( k, j ) = d21*t1 + d22*t2
461 80 CONTINUE
462 END IF
463 k = k - 2
464 END IF
465 GO TO 70
466 90 CONTINUE
467*
468* Form B := L'*B
469* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
470* and L' = inv(L'(m))*P(m)* ... *inv(L'(1))*P(1)
471*
472 ELSE
473*
474* Loop forward applying the L-transformations.
475*
476 k = 1
477 100 CONTINUE
478 IF( k.GT.n )
479 $ GO TO 120
480*
481* 1 x 1 pivot block
482*
483 IF( ipiv( k ).GT.0 ) THEN
484 IF( k.LT.n ) THEN
485*
486* Interchange if P(K) .ne. I.
487*
488 kp = ipiv( k )
489 IF( kp.NE.k )
490 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
491*
492* Apply the transformation
493*
494 CALL dgemv( 'Transpose', n-k, nrhs, one, b( k+1, 1 ),
495 $ ldb, a( k+1, k ), 1, one, b( k, 1 ), ldb )
496 END IF
497 IF( nounit )
498 $ CALL dscal( nrhs, a( k, k ), b( k, 1 ), ldb )
499 k = k + 1
500*
501* 2 x 2 pivot block.
502*
503 ELSE
504 IF( k.LT.n-1 ) THEN
505*
506* Interchange if P(K) .ne. I.
507*
508 kp = abs( ipiv( k ) )
509 IF( kp.NE.k+1 )
510 $ CALL dswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
511 $ ldb )
512*
513* Apply the transformation
514*
515 CALL dgemv( 'Transpose', n-k-1, nrhs, one,
516 $ b( k+2, 1 ), ldb, a( k+2, k+1 ), 1, one,
517 $ b( k+1, 1 ), ldb )
518 CALL dgemv( 'Transpose', n-k-1, nrhs, one,
519 $ b( k+2, 1 ), ldb, a( k+2, k ), 1, one,
520 $ b( k, 1 ), ldb )
521 END IF
522*
523* Multiply by the diagonal block if non-unit.
524*
525 IF( nounit ) THEN
526 d11 = a( k, k )
527 d22 = a( k+1, k+1 )
528 d21 = a( k+1, k )
529 d12 = d21
530 DO 110 j = 1, nrhs
531 t1 = b( k, j )
532 t2 = b( k+1, j )
533 b( k, j ) = d11*t1 + d12*t2
534 b( k+1, j ) = d21*t1 + d22*t2
535 110 CONTINUE
536 END IF
537 k = k + 2
538 END IF
539 GO TO 100
540 120 CONTINUE
541 END IF
542*
543 END IF
544 RETURN
545*
546* End of DLAVSY
547*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
dger
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: