LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slavsp()

subroutine slavsp ( character  UPLO,
character  TRANS,
character  DIAG,
integer  N,
integer  NRHS,
real, dimension( * )  A,
integer, dimension( * )  IPIV,
real, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

SLAVSP

Purpose:
 SLAVSP  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 SSPTRF.

 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 REAL 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 SSPTRF.
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          The pivot indices from SSPTRF.
[in,out]B
          B is REAL 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 128 of file slavsp.f.

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  REAL A( * ), B( LDB, * )
142 * ..
143 *
144 * =====================================================================
145 *
146 * .. Parameters ..
147  REAL ONE
148  parameter( one = 1.0e+0 )
149 * ..
150 * .. Local Scalars ..
151  LOGICAL NOUNIT
152  INTEGER J, K, KC, KCNEXT, KP
153  REAL D11, D12, D21, D22, T1, T2
154 * ..
155 * .. External Functions ..
156  LOGICAL LSAME
157  EXTERNAL lsame
158 * ..
159 * .. External Subroutines ..
160  EXTERNAL sgemv, sger, sscal, sswap, 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( 'SLAVSP ', -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 sscal( 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 sger( 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 sswap( 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 sger( k-1, nrhs, one, a( kc ), 1, b( k, 1 ), ldb,
268  $ b( 1, 1 ), ldb )
269  CALL sger( 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 sswap( 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 sscal( 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 sger( 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 sswap( 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 sger( n-k, nrhs, one, a( kc+1 ), 1, b( k, 1 ),
356  $ ldb, b( k+1, 1 ), ldb )
357  CALL sger( 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 sswap( 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 sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
405 *
406 * Apply the transformation
407 *
408  CALL sgemv( 'Transpose', k-1, nrhs, one, b, ldb,
409  $ a( kc ), 1, one, b( k, 1 ), ldb )
410  END IF
411  IF( nounit )
412  $ CALL sscal( 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 sswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ),
426  $ ldb )
427 *
428 * Apply the transformations
429 *
430  CALL sgemv( 'Transpose', k-2, nrhs, one, b, ldb,
431  $ a( kc ), 1, one, b( k, 1 ), ldb )
432  CALL sgemv( '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 sswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
480 *
481 * Apply the transformation
482 *
483  CALL sgemv( '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 sscal( 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 sswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ),
502  $ ldb )
503 *
504 * Apply the transformation
505 *
506  CALL sgemv( '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 sgemv( '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 SLAVSP
539 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine sger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SGER
Definition: sger.f:130
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
Here is the call graph for this function:
Here is the caller graph for this function: