LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ dlavsp()

subroutine dlavsp ( character  UPLO,
character  TRANS,
character  DIAG,
integer  N,
integer  NRHS,
double precision, dimension( * )  A,
integer, dimension( * )  IPIV,
double precision, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

DLAVSP

Purpose:
 DLAVSP  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 DSPTRF.

 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 (N*(N+1)/2)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L, stored as a packed triangular
          matrix as computed by DSPTRF.
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          The pivot indices from DSPTRF.
[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.
Date
December 2016

Definition at line 132 of file dlavsp.f.

132 *
133 * -- LAPACK test routine (version 3.7.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 * December 2016
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 *
Here is the call graph for this function:
Here is the caller graph for this function:
dger
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:132
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
dswap
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:84
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81