LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
dlavsp.f
Go to the documentation of this file.
1*> \brief \b DLAVSP
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 DLAVSP( 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* DOUBLE PRECISION A( * ), B( LDB, * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> DLAVSP 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 DSPTRF.
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 (N*(N+1)/2)
86*> The block diagonal matrix D and the multipliers used to
87*> obtain the factor U or L, stored as a packed triangular
88*> matrix as computed by DSPTRF.
89*> \endverbatim
90*>
91*> \param[in] IPIV
92*> \verbatim
93*> IPIV is INTEGER array, dimension (N)
94*> The pivot indices from DSPTRF.
95*> \endverbatim
96*>
97*> \param[in,out] B
98*> \verbatim
99*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
100*> On entry, B contains NRHS vectors of length N.
101*> On exit, B is overwritten with the product A * B.
102*> \endverbatim
103*>
104*> \param[in] LDB
105*> \verbatim
106*> LDB is INTEGER
107*> The leading dimension of the array B. LDB >= max(1,N).
108*> \endverbatim
109*>
110*> \param[out] INFO
111*> \verbatim
112*> INFO is INTEGER
113*> = 0: successful exit
114*> < 0: if INFO = -k, the k-th argument had an illegal value
115*> \endverbatim
116*
117* Authors:
118* ========
119*
120*> \author Univ. of Tennessee
121*> \author Univ. of California Berkeley
122*> \author Univ. of Colorado Denver
123*> \author NAG Ltd.
124*
125*> \ingroup double_lin
126*
127* =====================================================================
128 SUBROUTINE dlavsp( UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB,
129 $ INFO )
130*
131* -- LAPACK test routine --
132* -- LAPACK is a software package provided by Univ. of Tennessee, --
133* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
134*
135* .. Scalar Arguments ..
136 CHARACTER DIAG, TRANS, UPLO
137 INTEGER INFO, LDB, N, NRHS
138* ..
139* .. Array Arguments ..
140 INTEGER IPIV( * )
141 DOUBLE PRECISION A( * ), B( LDB, * )
142* ..
143*
144* =====================================================================
145*
146* .. Parameters ..
147 DOUBLE PRECISION ONE
148 parameter( one = 1.0d+0 )
149* ..
150* .. Local Scalars ..
151 LOGICAL NOUNIT
152 INTEGER J, K, KC, KCNEXT, KP
153 DOUBLE PRECISION D11, D12, D21, D22, T1, T2
154* ..
155* .. External Functions ..
156 LOGICAL LSAME
157 EXTERNAL lsame
158* ..
159* .. External Subroutines ..
160 EXTERNAL dgemv, dger, dscal, dswap, xerbla
161* ..
162* .. Intrinsic Functions ..
163 INTRINSIC abs, max
164* ..
165* .. Executable Statements ..
166*
167* Test the input parameters.
168*
169 info = 0
170 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
171 info = -1
172 ELSE IF( .NOT.lsame( trans, 'N' ) .AND. .NOT.
173 $ lsame( trans, 'T' ) .AND. .NOT.lsame( trans, 'C' ) ) THEN
174 info = -2
175 ELSE IF( .NOT.lsame( diag, 'U' ) .AND. .NOT.lsame( diag, 'N' ) )
176 $ THEN
177 info = -3
178 ELSE IF( n.LT.0 ) THEN
179 info = -4
180 ELSE IF( ldb.LT.max( 1, n ) ) THEN
181 info = -8
182 END IF
183 IF( info.NE.0 ) THEN
184 CALL xerbla( 'DLAVSP ', -info )
185 RETURN
186 END IF
187*
188* Quick return if possible.
189*
190 IF( n.EQ.0 )
191 $ RETURN
192*
193 nounit = lsame( diag, 'N' )
194*------------------------------------------
195*
196* Compute B := A * B (No transpose)
197*
198*------------------------------------------
199 IF( lsame( trans, 'N' ) ) THEN
200*
201* Compute B := U*B
202* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
203*
204 IF( lsame( uplo, 'U' ) ) THEN
205*
206* Loop forward applying the transformations.
207*
208 k = 1
209 kc = 1
210 10 CONTINUE
211 IF( k.GT.n )
212 $ GO TO 30
213*
214* 1 x 1 pivot block
215*
216 IF( ipiv( k ).GT.0 ) THEN
217*
218* Multiply by the diagonal element if forming U * D.
219*
220 IF( nounit )
221 $ CALL dscal( nrhs, a( kc+k-1 ), b( k, 1 ), ldb )
222*
223* Multiply by P(K) * inv(U(K)) if K > 1.
224*
225 IF( k.GT.1 ) THEN
226*
227* Apply the transformation.
228*
229 CALL dger( k-1, nrhs, one, a( kc ), 1, b( k, 1 ), ldb,
230 $ b( 1, 1 ), ldb )
231*
232* Interchange if P(K) != I.
233*
234 kp = ipiv( k )
235 IF( kp.NE.k )
236 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
237 END IF
238 kc = kc + k
239 k = k + 1
240 ELSE
241*
242* 2 x 2 pivot block
243*
244 kcnext = kc + k
245*
246* Multiply by the diagonal block if forming U * D.
247*
248 IF( nounit ) THEN
249 d11 = a( kcnext-1 )
250 d22 = a( kcnext+k )
251 d12 = a( kcnext+k-1 )
252 d21 = d12
253 DO 20 j = 1, nrhs
254 t1 = b( k, j )
255 t2 = b( k+1, j )
256 b( k, j ) = d11*t1 + d12*t2
257 b( k+1, j ) = d21*t1 + d22*t2
258 20 CONTINUE
259 END IF
260*
261* Multiply by P(K) * inv(U(K)) if K > 1.
262*
263 IF( k.GT.1 ) THEN
264*
265* Apply the transformations.
266*
267 CALL dger( k-1, nrhs, one, a( kc ), 1, b( k, 1 ), ldb,
268 $ b( 1, 1 ), ldb )
269 CALL dger( k-1, nrhs, one, a( kcnext ), 1,
270 $ b( k+1, 1 ), ldb, b( 1, 1 ), ldb )
271*
272* Interchange if P(K) != I.
273*
274 kp = abs( ipiv( k ) )
275 IF( kp.NE.k )
276 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
277 END IF
278 kc = kcnext + k + 1
279 k = k + 2
280 END IF
281 GO TO 10
282 30 CONTINUE
283*
284* Compute B := L*B
285* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
286*
287 ELSE
288*
289* Loop backward applying the transformations to B.
290*
291 k = n
292 kc = n*( n+1 ) / 2 + 1
293 40 CONTINUE
294 IF( k.LT.1 )
295 $ GO TO 60
296 kc = kc - ( n-k+1 )
297*
298* Test the pivot index. If greater than zero, a 1 x 1
299* pivot was used, otherwise a 2 x 2 pivot was used.
300*
301 IF( ipiv( k ).GT.0 ) THEN
302*
303* 1 x 1 pivot block:
304*
305* Multiply by the diagonal element if forming L * D.
306*
307 IF( nounit )
308 $ CALL dscal( nrhs, a( kc ), b( k, 1 ), ldb )
309*
310* Multiply by P(K) * inv(L(K)) if K < N.
311*
312 IF( k.NE.n ) THEN
313 kp = ipiv( k )
314*
315* Apply the transformation.
316*
317 CALL dger( n-k, nrhs, one, a( kc+1 ), 1, b( k, 1 ),
318 $ ldb, b( k+1, 1 ), ldb )
319*
320* Interchange if a permutation was applied at the
321* K-th step of the factorization.
322*
323 IF( kp.NE.k )
324 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
325 END IF
326 k = k - 1
327*
328 ELSE
329*
330* 2 x 2 pivot block:
331*
332 kcnext = kc - ( n-k+2 )
333*
334* Multiply by the diagonal block if forming L * D.
335*
336 IF( nounit ) THEN
337 d11 = a( kcnext )
338 d22 = a( kc )
339 d21 = a( kcnext+1 )
340 d12 = d21
341 DO 50 j = 1, nrhs
342 t1 = b( k-1, j )
343 t2 = b( k, j )
344 b( k-1, j ) = d11*t1 + d12*t2
345 b( k, j ) = d21*t1 + d22*t2
346 50 CONTINUE
347 END IF
348*
349* Multiply by P(K) * inv(L(K)) if K < N.
350*
351 IF( k.NE.n ) THEN
352*
353* Apply the transformation.
354*
355 CALL dger( n-k, nrhs, one, a( kc+1 ), 1, b( k, 1 ),
356 $ ldb, b( k+1, 1 ), ldb )
357 CALL dger( n-k, nrhs, one, a( kcnext+2 ), 1,
358 $ b( k-1, 1 ), ldb, b( k+1, 1 ), ldb )
359*
360* Interchange if a permutation was applied at the
361* K-th step of the factorization.
362*
363 kp = abs( ipiv( k ) )
364 IF( kp.NE.k )
365 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
366 END IF
367 kc = kcnext
368 k = k - 2
369 END IF
370 GO TO 40
371 60 CONTINUE
372 END IF
373*----------------------------------------
374*
375* Compute B := A' * B (transpose)
376*
377*----------------------------------------
378 ELSE
379*
380* Form B := U'*B
381* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
382* and U' = inv(U'(1))*P(1)* ... *inv(U'(m))*P(m)
383*
384 IF( lsame( uplo, 'U' ) ) THEN
385*
386* Loop backward applying the transformations.
387*
388 k = n
389 kc = n*( n+1 ) / 2 + 1
390 70 CONTINUE
391 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 dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
405*
406* Apply the transformation
407*
408 CALL dgemv( 'Transpose', k-1, nrhs, one, b, ldb,
409 $ a( kc ), 1, one, b( k, 1 ), ldb )
410 END IF
411 IF( nounit )
412 $ CALL dscal( nrhs, a( kc+k-1 ), b( k, 1 ), ldb )
413 k = k - 1
414*
415* 2 x 2 pivot block.
416*
417 ELSE
418 kcnext = kc - ( k-1 )
419 IF( k.GT.2 ) THEN
420*
421* Interchange if P(K) != I.
422*
423 kp = abs( ipiv( k ) )
424 IF( kp.NE.k-1 )
425 $ CALL dswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
426 $ ldb )
427*
428* Apply the transformations
429*
430 CALL dgemv( 'Transpose', k-2, nrhs, one, b, ldb,
431 $ a( kc ), 1, one, b( k, 1 ), ldb )
432 CALL dgemv( 'Transpose', k-2, nrhs, one, b, ldb,
433 $ a( kcnext ), 1, one, b( k-1, 1 ), ldb )
434 END IF
435*
436* Multiply by the diagonal block if non-unit.
437*
438 IF( nounit ) THEN
439 d11 = a( kc-1 )
440 d22 = a( kc+k-1 )
441 d12 = a( kc+k-2 )
442 d21 = d12
443 DO 80 j = 1, nrhs
444 t1 = b( k-1, j )
445 t2 = b( k, j )
446 b( k-1, j ) = d11*t1 + d12*t2
447 b( k, j ) = d21*t1 + d22*t2
448 80 CONTINUE
449 END IF
450 kc = kcnext
451 k = k - 2
452 END IF
453 GO TO 70
454 90 CONTINUE
455*
456* Form B := L'*B
457* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
458* and L' = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
459*
460 ELSE
461*
462* Loop forward applying the L-transformations.
463*
464 k = 1
465 kc = 1
466 100 CONTINUE
467 IF( k.GT.n )
468 $ GO TO 120
469*
470* 1 x 1 pivot block
471*
472 IF( ipiv( k ).GT.0 ) THEN
473 IF( k.LT.n ) THEN
474*
475* Interchange if P(K) != I.
476*
477 kp = ipiv( k )
478 IF( kp.NE.k )
479 $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
480*
481* Apply the transformation
482*
483 CALL dgemv( 'Transpose', n-k, nrhs, one, b( k+1, 1 ),
484 $ ldb, a( kc+1 ), 1, one, b( k, 1 ), ldb )
485 END IF
486 IF( nounit )
487 $ CALL dscal( nrhs, a( kc ), b( k, 1 ), ldb )
488 kc = kc + n - k + 1
489 k = k + 1
490*
491* 2 x 2 pivot block.
492*
493 ELSE
494 kcnext = kc + n - k + 1
495 IF( k.LT.n-1 ) THEN
496*
497* Interchange if P(K) != I.
498*
499 kp = abs( ipiv( k ) )
500 IF( kp.NE.k+1 )
501 $ CALL dswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
502 $ ldb )
503*
504* Apply the transformation
505*
506 CALL dgemv( 'Transpose', n-k-1, nrhs, one,
507 $ b( k+2, 1 ), ldb, a( kcnext+1 ), 1, one,
508 $ b( k+1, 1 ), ldb )
509 CALL dgemv( 'Transpose', n-k-1, nrhs, one,
510 $ b( k+2, 1 ), ldb, a( kc+2 ), 1, one,
511 $ b( k, 1 ), ldb )
512 END IF
513*
514* Multiply by the diagonal block if non-unit.
515*
516 IF( nounit ) THEN
517 d11 = a( kc )
518 d22 = a( kcnext )
519 d21 = a( kc+1 )
520 d12 = d21
521 DO 110 j = 1, nrhs
522 t1 = b( k, j )
523 t2 = b( k+1, j )
524 b( k, j ) = d11*t1 + d12*t2
525 b( k+1, j ) = d21*t1 + d22*t2
526 110 CONTINUE
527 END IF
528 kc = kcnext + ( n-k )
529 k = k + 2
530 END IF
531 GO TO 100
532 120 CONTINUE
533 END IF
534*
535 END IF
536 RETURN
537*
538* End of DLAVSP
539*
540 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dlavsp
subroutine dlavsp(uplo, trans, diag, n, nrhs, a, ipiv, b, ldb, info)
DLAVSP
Definition dlavsp.f:130
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
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