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