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