LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 *> \date November 2011
126 *
127 *> \ingroup double_lin
128 *
129 * =====================================================================
130  SUBROUTINE dlavsp( UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB,
131  $ info )
132 *
133 * -- LAPACK test routine (version 3.4.0) --
134 * -- LAPACK is a software package provided by Univ. of Tennessee, --
135 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136 * November 2011
137 *
138 * .. Scalar Arguments ..
139  CHARACTER DIAG, TRANS, UPLO
140  INTEGER INFO, LDB, N, NRHS
141 * ..
142 * .. Array Arguments ..
143  INTEGER IPIV( * )
144  DOUBLE PRECISION A( * ), B( ldb, * )
145 * ..
146 *
147 * =====================================================================
148 *
149 * .. Parameters ..
150  DOUBLE PRECISION ONE
151  parameter ( one = 1.0d+0 )
152 * ..
153 * .. Local Scalars ..
154  LOGICAL NOUNIT
155  INTEGER J, K, KC, KCNEXT, KP
156  DOUBLE PRECISION D11, D12, D21, D22, T1, T2
157 * ..
158 * .. External Functions ..
159  LOGICAL LSAME
160  EXTERNAL lsame
161 * ..
162 * .. External Subroutines ..
163  EXTERNAL dgemv, dger, dscal, dswap, xerbla
164 * ..
165 * .. Intrinsic Functions ..
166  INTRINSIC abs, max
167 * ..
168 * .. Executable Statements ..
169 *
170 * Test the input parameters.
171 *
172  info = 0
173  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
174  info = -1
175  ELSE IF( .NOT.lsame( trans, 'N' ) .AND. .NOT.
176  $ lsame( trans, 'T' ) .AND. .NOT.lsame( trans, 'C' ) ) THEN
177  info = -2
178  ELSE IF( .NOT.lsame( diag, 'U' ) .AND. .NOT.lsame( diag, 'N' ) )
179  $ THEN
180  info = -3
181  ELSE IF( n.LT.0 ) THEN
182  info = -4
183  ELSE IF( ldb.LT.max( 1, n ) ) THEN
184  info = -8
185  END IF
186  IF( info.NE.0 ) THEN
187  CALL xerbla( 'DLAVSP ', -info )
188  RETURN
189  END IF
190 *
191 * Quick return if possible.
192 *
193  IF( n.EQ.0 )
194  $ RETURN
195 *
196  nounit = lsame( diag, 'N' )
197 *------------------------------------------
198 *
199 * Compute B := A * B (No transpose)
200 *
201 *------------------------------------------
202  IF( lsame( trans, 'N' ) ) THEN
203 *
204 * Compute B := U*B
205 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
206 *
207  IF( lsame( uplo, 'U' ) ) THEN
208 *
209 * Loop forward applying the transformations.
210 *
211  k = 1
212  kc = 1
213  10 CONTINUE
214  IF( k.GT.n )
215  $ GO TO 30
216 *
217 * 1 x 1 pivot block
218 *
219  IF( ipiv( k ).GT.0 ) THEN
220 *
221 * Multiply by the diagonal element if forming U * D.
222 *
223  IF( nounit )
224  $ CALL dscal( nrhs, a( kc+k-1 ), b( k, 1 ), ldb )
225 *
226 * Multiply by P(K) * inv(U(K)) if K > 1.
227 *
228  IF( k.GT.1 ) THEN
229 *
230 * Apply the transformation.
231 *
232  CALL dger( k-1, nrhs, one, a( kc ), 1, b( k, 1 ), ldb,
233  $ b( 1, 1 ), ldb )
234 *
235 * Interchange if P(K) != I.
236 *
237  kp = ipiv( k )
238  IF( kp.NE.k )
239  $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
240  END IF
241  kc = kc + k
242  k = k + 1
243  ELSE
244 *
245 * 2 x 2 pivot block
246 *
247  kcnext = kc + k
248 *
249 * Multiply by the diagonal block if forming U * D.
250 *
251  IF( nounit ) THEN
252  d11 = a( kcnext-1 )
253  d22 = a( kcnext+k )
254  d12 = a( kcnext+k-1 )
255  d21 = d12
256  DO 20 j = 1, nrhs
257  t1 = b( k, j )
258  t2 = b( k+1, j )
259  b( k, j ) = d11*t1 + d12*t2
260  b( k+1, j ) = d21*t1 + d22*t2
261  20 CONTINUE
262  END IF
263 *
264 * Multiply by P(K) * inv(U(K)) if K > 1.
265 *
266  IF( k.GT.1 ) THEN
267 *
268 * Apply the transformations.
269 *
270  CALL dger( k-1, nrhs, one, a( kc ), 1, b( k, 1 ), ldb,
271  $ b( 1, 1 ), ldb )
272  CALL dger( k-1, nrhs, one, a( kcnext ), 1,
273  $ b( k+1, 1 ), ldb, b( 1, 1 ), ldb )
274 *
275 * Interchange if P(K) != I.
276 *
277  kp = abs( ipiv( k ) )
278  IF( kp.NE.k )
279  $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
280  END IF
281  kc = kcnext + k + 1
282  k = k + 2
283  END IF
284  GO TO 10
285  30 CONTINUE
286 *
287 * Compute B := L*B
288 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
289 *
290  ELSE
291 *
292 * Loop backward applying the transformations to B.
293 *
294  k = n
295  kc = n*( n+1 ) / 2 + 1
296  40 CONTINUE
297  IF( k.LT.1 )
298  $ GO TO 60
299  kc = kc - ( n-k+1 )
300 *
301 * Test the pivot index. If greater than zero, a 1 x 1
302 * pivot was used, otherwise a 2 x 2 pivot was used.
303 *
304  IF( ipiv( k ).GT.0 ) THEN
305 *
306 * 1 x 1 pivot block:
307 *
308 * Multiply by the diagonal element if forming L * D.
309 *
310  IF( nounit )
311  $ CALL dscal( nrhs, a( kc ), b( k, 1 ), ldb )
312 *
313 * Multiply by P(K) * inv(L(K)) if K < N.
314 *
315  IF( k.NE.n ) THEN
316  kp = ipiv( k )
317 *
318 * Apply the transformation.
319 *
320  CALL dger( n-k, nrhs, one, a( kc+1 ), 1, b( k, 1 ),
321  $ ldb, b( k+1, 1 ), ldb )
322 *
323 * Interchange if a permutation was applied at the
324 * K-th step of the factorization.
325 *
326  IF( kp.NE.k )
327  $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
328  END IF
329  k = k - 1
330 *
331  ELSE
332 *
333 * 2 x 2 pivot block:
334 *
335  kcnext = kc - ( n-k+2 )
336 *
337 * Multiply by the diagonal block if forming L * D.
338 *
339  IF( nounit ) THEN
340  d11 = a( kcnext )
341  d22 = a( kc )
342  d21 = a( kcnext+1 )
343  d12 = d21
344  DO 50 j = 1, nrhs
345  t1 = b( k-1, j )
346  t2 = b( k, j )
347  b( k-1, j ) = d11*t1 + d12*t2
348  b( k, j ) = d21*t1 + d22*t2
349  50 CONTINUE
350  END IF
351 *
352 * Multiply by P(K) * inv(L(K)) if K < N.
353 *
354  IF( k.NE.n ) THEN
355 *
356 * Apply the transformation.
357 *
358  CALL dger( n-k, nrhs, one, a( kc+1 ), 1, b( k, 1 ),
359  $ ldb, b( k+1, 1 ), ldb )
360  CALL dger( n-k, nrhs, one, a( kcnext+2 ), 1,
361  $ b( k-1, 1 ), ldb, b( k+1, 1 ), ldb )
362 *
363 * Interchange if a permutation was applied at the
364 * K-th step of the factorization.
365 *
366  kp = abs( ipiv( k ) )
367  IF( kp.NE.k )
368  $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
369  END IF
370  kc = kcnext
371  k = k - 2
372  END IF
373  GO TO 40
374  60 CONTINUE
375  END IF
376 *----------------------------------------
377 *
378 * Compute B := A' * B (transpose)
379 *
380 *----------------------------------------
381  ELSE
382 *
383 * Form B := U'*B
384 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
385 * and U' = inv(U'(1))*P(1)* ... *inv(U'(m))*P(m)
386 *
387  IF( lsame( uplo, 'U' ) ) THEN
388 *
389 * Loop backward applying the transformations.
390 *
391  k = n
392  kc = n*( n+1 ) / 2 + 1
393  70 CONTINUE
394  IF( k.LT.1 )
395  $ GO TO 90
396  kc = kc - k
397 *
398 * 1 x 1 pivot block.
399 *
400  IF( ipiv( k ).GT.0 ) THEN
401  IF( k.GT.1 ) THEN
402 *
403 * Interchange if P(K) != I.
404 *
405  kp = ipiv( k )
406  IF( kp.NE.k )
407  $ CALL dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
408 *
409 * Apply the transformation
410 *
411  CALL dgemv( 'Transpose', k-1, nrhs, one, b, ldb,
412  $ a( kc ), 1, one, b( k, 1 ), ldb )
413  END IF
414  IF( nounit )
415  $ CALL dscal( nrhs, a( kc+k-1 ), b( k, 1 ), ldb )
416  k = k - 1
417 *
418 * 2 x 2 pivot block.
419 *
420  ELSE
421  kcnext = kc - ( k-1 )
422  IF( k.GT.2 ) THEN
423 *
424 * Interchange if P(K) != I.
425 *
426  kp = abs( ipiv( k ) )
427  IF( kp.NE.k-1 )
428  $ CALL dswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
429  $ ldb )
430 *
431 * Apply the transformations
432 *
433  CALL dgemv( 'Transpose', k-2, nrhs, one, b, ldb,
434  $ a( kc ), 1, one, b( k, 1 ), ldb )
435  CALL dgemv( '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'*B
460 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
461 * and L' = 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 dswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
483 *
484 * Apply the transformation
485 *
486  CALL dgemv( '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 dscal( 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 dswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
505  $ ldb )
506 *
507 * Apply the transformation
508 *
509  CALL dgemv( 'Transpose', n-k-1, nrhs, one,
510  $ b( k+2, 1 ), ldb, a( kcnext+1 ), 1, one,
511  $ b( k+1, 1 ), ldb )
512  CALL dgemv( 'Transpose', n-k-1, nrhs, one,
513  $ b( k+2, 1 ), ldb, a( kc+2 ), 1, one,
514  $ b( k, 1 ), ldb )
515  END IF
516 *
517 * Multiply by the diagonal block if non-unit.
518 *
519  IF( nounit ) THEN
520  d11 = a( kc )
521  d22 = a( kcnext )
522  d21 = a( kc+1 )
523  d12 = d21
524  DO 110 j = 1, nrhs
525  t1 = b( k, j )
526  t2 = b( k+1, j )
527  b( k, j ) = d11*t1 + d12*t2
528  b( k+1, j ) = d21*t1 + d22*t2
529  110 CONTINUE
530  END IF
531  kc = kcnext + ( n-k )
532  k = k + 2
533  END IF
534  GO TO 100
535  120 CONTINUE
536  END IF
537 *
538  END IF
539  RETURN
540 *
541 * End of DLAVSP
542 *
543  END
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dswap
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:53
dger
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:132
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
dlavsp
subroutine dlavsp(UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB, INFO)
DLAVSP
Definition: dlavsp.f:132