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