1       SUBROUTINE SLAVSY( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B,
  2      $                   LDB, INFO )
  3 *
  4 *  -- LAPACK auxiliary routine (version 3.1) --
  5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  6 *     November 2006
  7 *
  8 *     .. Scalar Arguments ..
  9       CHARACTER          DIAG, TRANS, UPLO
 10       INTEGER            INFO, LDA, LDB, N, NRHS
 11 *     ..
 12 *     .. Array Arguments ..
 13       INTEGER            IPIV( * )
 14       REAL               A( LDA, * ), B( LDB, * )
 15 *     ..
 16 *
 17 *  Purpose
 18 *  =======
 19 *
 20 *  SLAVSY  performs one of the matrix-vector operations
 21 *     x := A*x  or  x := A'*x,
 22 *  where x is an N element vector and A is one of the factors
 23 *  from the block U*D*U' or L*D*L' factorization computed by SSYTRF.
 24 *
 25 *  If TRANS = 'N', multiplies by U  or U * D  (or L  or L * D)
 26 *  If TRANS = 'T', multiplies by U' or D * U' (or L' or D * L')
 27 *  If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L')
 28 *
 29 *  Arguments
 30 *  =========
 31 *
 32 *  UPLO    (input) CHARACTER*1
 33 *          Specifies whether the factor stored in A is upper or lower
 34 *          triangular.
 35 *          = 'U':  Upper triangular
 36 *          = 'L':  Lower triangular
 37 *
 38 *  TRANS   (input) CHARACTER*1
 39 *          Specifies the operation to be performed:
 40 *          = 'N':  x := A*x
 41 *          = 'T':  x := A'*x
 42 *          = 'C':  x := A'*x
 43 *
 44 *  DIAG    (input) CHARACTER*1
 45 *          Specifies whether or not the diagonal blocks are unit
 46 *          matrices.  If the diagonal blocks are assumed to be unit,
 47 *          then A = U or A = L, otherwise A = U*D or A = L*D.
 48 *          = 'U':  Diagonal blocks are assumed to be unit matrices.
 49 *          = 'N':  Diagonal blocks are assumed to be non-unit matrices.
 50 *
 51 *  N       (input) INTEGER
 52 *          The number of rows and columns of the matrix A.  N >= 0.
 53 *
 54 *  NRHS    (input) INTEGER
 55 *          The number of right hand sides, i.e., the number of vectors
 56 *          x to be multiplied by A.  NRHS >= 0.
 57 *
 58 *  A       (input) REAL array, dimension (LDA,N)
 59 *          The block diagonal matrix D and the multipliers used to
 60 *          obtain the factor U or L as computed by SSYTRF.
 61 *
 62 *  LDA     (input) INTEGER
 63 *          The leading dimension of the array A.  LDA >= max(1,N).
 64 *
 65 *  IPIV    (input) INTEGER array, dimension (N)
 66 *          The pivot indices from SSYTRF.
 67 *
 68 *  B       (input/output) REAL array, dimension (LDB,NRHS)
 69 *          On entry, B contains NRHS vectors of length N.
 70 *          On exit, B is overwritten with the product A * B.
 71 *
 72 *  LDB     (input) INTEGER
 73 *          The leading dimension of the array B.  LDB >= max(1,N).
 74 *
 75 *  INFO    (output) INTEGER
 76 *          = 0: successful exit
 77 *          < 0: if INFO = -k, the k-th argument had an illegal value
 78 *
 79 *  =====================================================================
 80 *
 81 *     .. Parameters ..
 82       REAL               ONE
 83       PARAMETER          ( ONE = 1.0E+0 )
 84 *     ..
 85 *     .. Local Scalars ..
 86       LOGICAL            NOUNIT
 87       INTEGER            J, K, KP
 88       REAL               D11, D12, D21, D22, T1, T2
 89 *     ..
 90 *     .. External Functions ..
 91       LOGICAL            LSAME
 92       EXTERNAL           LSAME
 93 *     ..
 94 *     .. External Subroutines ..
 95       EXTERNAL           SGEMV, SGER, SSCAL, SSWAP, XERBLA
 96 *     ..
 97 *     .. Intrinsic Functions ..
 98       INTRINSIC          ABSMAX
 99 *     ..
100 *     .. Executable Statements ..
101 *
102 *     Test the input parameters.
103 *
104       INFO = 0
105       IF.NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
106          INFO = -1
107       ELSE IF.NOT.LSAME( TRANS, 'N' ) .AND. .NOT.
108      $         LSAME( TRANS, 'T' ) .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
109          INFO = -2
110       ELSE IF.NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) )
111      $          THEN
112          INFO = -3
113       ELSE IF( N.LT.0 ) THEN
114          INFO = -4
115       ELSE IF( LDA.LT.MAX1, N ) ) THEN
116          INFO = -6
117       ELSE IF( LDB.LT.MAX1, N ) ) THEN
118          INFO = -9
119       END IF
120       IF( INFO.NE.0 ) THEN
121          CALL XERBLA( 'SLAVSY '-INFO )
122          RETURN
123       END IF
124 *
125 *     Quick return if possible.
126 *
127       IF( N.EQ.0 )
128      $   RETURN
129 *
130       NOUNIT = LSAME( DIAG, 'N' )
131 *------------------------------------------
132 *
133 *     Compute  B := A * B  (No transpose)
134 *
135 *------------------------------------------
136       IF( LSAME( TRANS, 'N' ) ) THEN
137 *
138 *        Compute  B := U*B
139 *        where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
140 *
141          IF( LSAME( UPLO, 'U' ) ) THEN
142 *
143 *        Loop forward applying the transformations.
144 *
145             K = 1
146    10       CONTINUE
147             IF( K.GT.N )
148      $         GO TO 30
149             IF( IPIV( K ).GT.0 ) THEN
150 *
151 *              1 x 1 pivot block
152 *
153 *              Multiply by the diagonal element if forming U * D.
154 *
155                IF( NOUNIT )
156      $            CALL SSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
157 *
158 *              Multiply by  P(K) * inv(U(K))  if K > 1.
159 *
160                IF( K.GT.1 ) THEN
161 *
162 *                 Apply the transformation.
163 *
164                   CALL SGER( K-1, NRHS, ONE, A( 1, K ), 1, B( K, 1 ),
165      $                       LDB, B( 11 ), LDB )
166 *
167 *                 Interchange if P(K) .ne. I.
168 *
169                   KP = IPIV( K )
170                   IF( KP.NE.K )
171      $               CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
172                END IF
173                K = K + 1
174             ELSE
175 *
176 *              2 x 2 pivot block
177 *
178 *              Multiply by the diagonal block if forming U * D.
179 *
180                IF( NOUNIT ) THEN
181                   D11 = A( K, K )
182                   D22 = A( K+1, K+1 )
183                   D12 = A( K, K+1 )
184                   D21 = D12
185                   DO 20 J = 1, NRHS
186                      T1 = B( K, J )
187                      T2 = B( K+1, J )
188                      B( K, J ) = D11*T1 + D12*T2
189                      B( K+1, J ) = D21*T1 + D22*T2
190    20             CONTINUE
191                END IF
192 *
193 *              Multiply by  P(K) * inv(U(K))  if K > 1.
194 *
195                IF( K.GT.1 ) THEN
196 *
197 *                 Apply the transformations.
198 *
199                   CALL SGER( K-1, NRHS, ONE, A( 1, K ), 1, B( K, 1 ),
200      $                       LDB, B( 11 ), LDB )
201                   CALL SGER( K-1, NRHS, ONE, A( 1, K+1 ), 1,
202      $                       B( K+11 ), LDB, B( 11 ), LDB )
203 *
204 *                 Interchange if P(K) .ne. I.
205 *
206                   KP = ABS( IPIV( K ) )
207                   IF( KP.NE.K )
208      $               CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
209                END IF
210                K = K + 2
211             END IF
212             GO TO 10
213    30       CONTINUE
214 *
215 *        Compute  B := L*B
216 *        where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
217 *
218          ELSE
219 *
220 *           Loop backward applying the transformations to B.
221 *
222             K = N
223    40       CONTINUE
224             IF( K.LT.1 )
225      $         GO TO 60
226 *
227 *           Test the pivot index.  If greater than zero, a 1 x 1
228 *           pivot was used, otherwise a 2 x 2 pivot was used.
229 *
230             IF( IPIV( K ).GT.0 ) THEN
231 *
232 *              1 x 1 pivot block:
233 *
234 *              Multiply by the diagonal element if forming L * D.
235 *
236                IF( NOUNIT )
237      $            CALL SSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
238 *
239 *              Multiply by  P(K) * inv(L(K))  if K < N.
240 *
241                IF( K.NE.N ) THEN
242                   KP = IPIV( K )
243 *
244 *                 Apply the transformation.
245 *
246                   CALL SGER( N-K, NRHS, ONE, A( K+1, K ), 1, B( K, 1 ),
247      $                       LDB, B( K+11 ), LDB )
248 *
249 *                 Interchange if a permutation was applied at the
250 *                 K-th step of the factorization.
251 *
252                   IF( KP.NE.K )
253      $               CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
254                END IF
255                K = K - 1
256 *
257             ELSE
258 *
259 *              2 x 2 pivot block:
260 *
261 *              Multiply by the diagonal block if forming L * D.
262 *
263                IF( NOUNIT ) THEN
264                   D11 = A( K-1, K-1 )
265                   D22 = A( K, K )
266                   D21 = A( K, K-1 )
267                   D12 = D21
268                   DO 50 J = 1, NRHS
269                      T1 = B( K-1, J )
270                      T2 = B( K, J )
271                      B( K-1, J ) = D11*T1 + D12*T2
272                      B( K, J ) = D21*T1 + D22*T2
273    50             CONTINUE
274                END IF
275 *
276 *              Multiply by  P(K) * inv(L(K))  if K < N.
277 *
278                IF( K.NE.N ) THEN
279 *
280 *                 Apply the transformation.
281 *
282                   CALL SGER( N-K, NRHS, ONE, A( K+1, K ), 1, B( K, 1 ),
283      $                       LDB, B( K+11 ), LDB )
284                   CALL SGER( N-K, NRHS, ONE, A( K+1, K-1 ), 1,
285      $                       B( K-11 ), LDB, B( K+11 ), LDB )
286 *
287 *                 Interchange if a permutation was applied at the
288 *                 K-th step of the factorization.
289 *
290                   KP = ABS( IPIV( K ) )
291                   IF( KP.NE.K )
292      $               CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
293                END IF
294                K = K - 2
295             END IF
296             GO TO 40
297    60       CONTINUE
298          END IF
299 *----------------------------------------
300 *
301 *     Compute  B := A' * B  (transpose)
302 *
303 *----------------------------------------
304       ELSE
305 *
306 *        Form  B := U'*B
307 *        where U  = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
308 *        and   U' = inv(U'(1))*P(1)* ... *inv(U'(m))*P(m)
309 *
310          IF( LSAME( UPLO, 'U' ) ) THEN
311 *
312 *           Loop backward applying the transformations.
313 *
314             K = N
315    70       CONTINUE
316             IF( K.LT.1 )
317      $         GO TO 90
318 *
319 *           1 x 1 pivot block.
320 *
321             IF( IPIV( K ).GT.0 ) THEN
322                IF( K.GT.1 ) THEN
323 *
324 *                 Interchange if P(K) .ne. I.
325 *
326                   KP = IPIV( K )
327                   IF( KP.NE.K )
328      $               CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
329 *
330 *                 Apply the transformation
331 *
332                   CALL SGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB,
333      $                        A( 1, K ), 1, ONE, B( K, 1 ), LDB )
334                END IF
335                IF( NOUNIT )
336      $            CALL SSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
337                K = K - 1
338 *
339 *           2 x 2 pivot block.
340 *
341             ELSE
342                IF( K.GT.2 ) THEN
343 *
344 *                 Interchange if P(K) .ne. I.
345 *
346                   KP = ABS( IPIV( K ) )
347                   IF( KP.NE.K-1 )
348      $               CALL SSWAP( NRHS, B( K-11 ), LDB, B( KP, 1 ),
349      $                           LDB )
350 *
351 *                 Apply the transformations
352 *
353                   CALL SGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
354      $                        A( 1, K ), 1, ONE, B( K, 1 ), LDB )
355                   CALL SGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
356      $                        A( 1, K-1 ), 1, ONE, B( K-11 ), LDB )
357                END IF
358 *
359 *              Multiply by the diagonal block if non-unit.
360 *
361                IF( NOUNIT ) THEN
362                   D11 = A( K-1, K-1 )
363                   D22 = A( K, K )
364                   D12 = A( K-1, K )
365                   D21 = D12
366                   DO 80 J = 1, NRHS
367                      T1 = B( K-1, J )
368                      T2 = B( K, J )
369                      B( K-1, J ) = D11*T1 + D12*T2
370                      B( K, J ) = D21*T1 + D22*T2
371    80             CONTINUE
372                END IF
373                K = K - 2
374             END IF
375             GO TO 70
376    90       CONTINUE
377 *
378 *        Form  B := L'*B
379 *        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
380 *        and   L' = inv(L'(m))*P(m)* ... *inv(L'(1))*P(1)
381 *
382          ELSE
383 *
384 *           Loop forward applying the L-transformations.
385 *
386             K = 1
387   100       CONTINUE
388             IF( K.GT.N )
389      $         GO TO 120
390 *
391 *           1 x 1 pivot block
392 *
393             IF( IPIV( K ).GT.0 ) THEN
394                IF( K.LT.N ) THEN
395 *
396 *                 Interchange if P(K) .ne. I.
397 *
398                   KP = IPIV( K )
399                   IF( KP.NE.K )
400      $               CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
401 *
402 *                 Apply the transformation
403 *
404                   CALL SGEMV( 'Transpose', N-K, NRHS, ONE, B( K+11 ),
405      $                        LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB )
406                END IF
407                IF( NOUNIT )
408      $            CALL SSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
409                K = K + 1
410 *
411 *           2 x 2 pivot block.
412 *
413             ELSE
414                IF( K.LT.N-1 ) THEN
415 *
416 *              Interchange if P(K) .ne. I.
417 *
418                   KP = ABS( IPIV( K ) )
419                   IF( KP.NE.K+1 )
420      $               CALL SSWAP( NRHS, B( K+11 ), LDB, B( KP, 1 ),
421      $                           LDB )
422 *
423 *                 Apply the transformation
424 *
425                   CALL SGEMV( 'Transpose', N-K-1, NRHS, ONE,
426      $                        B( K+21 ), LDB, A( K+2, K+1 ), 1, ONE,
427      $                        B( K+11 ), LDB )
428                   CALL SGEMV( 'Transpose', N-K-1, NRHS, ONE,
429      $                        B( K+21 ), LDB, A( K+2, K ), 1, ONE,
430      $                        B( K, 1 ), LDB )
431                END IF
432 *
433 *              Multiply by the diagonal block if non-unit.
434 *
435                IF( NOUNIT ) THEN
436                   D11 = A( K, K )
437                   D22 = A( K+1, K+1 )
438                   D21 = A( K+1, K )
439                   D12 = D21
440                   DO 110 J = 1, NRHS
441                      T1 = B( K, J )
442                      T2 = B( K+1, J )
443                      B( K, J ) = D11*T1 + D12*T2
444                      B( K+1, J ) = D21*T1 + D22*T2
445   110             CONTINUE
446                END IF
447                K = K + 2
448             END IF
449             GO TO 100
450   120       CONTINUE
451          END IF
452 *
453       END IF
454       RETURN
455 *
456 *     End of SLAVSY
457 *
458       END