LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlavsy.f
Go to the documentation of this file.
1*> \brief \b DLAVSY
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 DLAVSY( 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* DOUBLE PRECISION A( LDA, * ), B( LDB, * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> DLAVSY performs one of the matrix-vector operations
30*> x := A*x or x := A'*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 DSYTRF.
33*>
34*> If TRANS = 'N', multiplies by U or U * D (or L or L * D)
35*> If TRANS = 'T', multiplies by U' or D * U' (or L' or D * L')
36*> If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L')
37*> \endverbatim
38*
39* Arguments:
40* ==========
41*
42*> \param[in] UPLO
43*> \verbatim
44*> UPLO is CHARACTER*1
45*> Specifies whether the factor stored in A is upper or lower
46*> triangular.
47*> = 'U': Upper triangular
48*> = 'L': Lower triangular
49*> \endverbatim
50*>
51*> \param[in] TRANS
52*> \verbatim
53*> TRANS is CHARACTER*1
54*> Specifies the operation to be performed:
55*> = 'N': x := A*x
56*> = 'T': x := A'*x
57*> = 'C': x := A'*x
58*> \endverbatim
59*>
60*> \param[in] DIAG
61*> \verbatim
62*> DIAG is CHARACTER*1
63*> Specifies whether or not the diagonal blocks are unit
64*> matrices. If the diagonal blocks are assumed to be unit,
65*> then A = U or A = L, otherwise A = U*D or A = L*D.
66*> = 'U': Diagonal blocks are assumed to be unit matrices.
67*> = 'N': Diagonal blocks are assumed to be non-unit matrices.
68*> \endverbatim
69*>
70*> \param[in] N
71*> \verbatim
72*> N is INTEGER
73*> The number of rows and columns of the matrix A. N >= 0.
74*> \endverbatim
75*>
76*> \param[in] NRHS
77*> \verbatim
78*> NRHS is INTEGER
79*> The number of right hand sides, i.e., the number of vectors
80*> x to be multiplied by A. NRHS >= 0.
81*> \endverbatim
82*>
83*> \param[in] A
84*> \verbatim
85*> A is DOUBLE PRECISION array, dimension (LDA,N)
86*> The block diagonal matrix D and the multipliers used to
87*> obtain the factor U or L as computed by DSYTRF.
88*> Stored as a 2-D triangular matrix.
89*> \endverbatim
90*>
91*> \param[in] LDA
92*> \verbatim
93*> LDA is INTEGER
94*> The leading dimension of the array A. LDA >= max(1,N).
95*> \endverbatim
96*>
97*> \param[in] IPIV
98*> \verbatim
99*> IPIV is INTEGER array, dimension (N)
100*> Details of the interchanges and the block structure of D,
101*> as determined by DSYTRF.
102*>
103*> If UPLO = 'U':
104*> If IPIV(k) > 0, then rows and columns k and IPIV(k)
105*> were interchanged and D(k,k) is a 1-by-1 diagonal block.
106*> (If IPIV( k ) = k, no interchange was done).
107*>
108*> If IPIV(k) = IPIV(k-1) < 0, then rows and
109*> columns k-1 and -IPIV(k) were interchanged,
110*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
111*>
112*> If UPLO = 'L':
113*> If IPIV(k) > 0, then rows and columns k and IPIV(k)
114*> were interchanged and D(k,k) is a 1-by-1 diagonal block.
115*> (If IPIV( k ) = k, no interchange was done).
116*>
117*> If IPIV(k) = IPIV(k+1) < 0, then rows and
118*> columns k+1 and -IPIV(k) were interchanged,
119*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
120*> \endverbatim
121*>
122*> \param[in,out] B
123*> \verbatim
124*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
125*> On entry, B contains NRHS vectors of length N.
126*> On exit, B is overwritten with the product A * B.
127*> \endverbatim
128*>
129*> \param[in] LDB
130*> \verbatim
131*> LDB is INTEGER
132*> The leading dimension of the array B. LDB >= max(1,N).
133*> \endverbatim
134*>
135*> \param[out] INFO
136*> \verbatim
137*> INFO is INTEGER
138*> = 0: successful exit
139*> < 0: if INFO = -k, the k-th argument had an illegal value
140*> \endverbatim
141*
142* Authors:
143* ========
144*
145*> \author Univ. of Tennessee
146*> \author Univ. of California Berkeley
147*> \author Univ. of Colorado Denver
148*> \author NAG Ltd.
149*
150*> \ingroup double_lin
151*
152* =====================================================================
153 SUBROUTINE dlavsy( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B,
154 $ LDB, INFO )
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*
548 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:130
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:156
subroutine dlavsy(UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
DLAVSY
Definition: dlavsy.f:155