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*> \endverbatim
87*>
88*> \param[in] LDA
89*> \verbatim
90*>          LDA is INTEGER
91*>          The leading dimension of the array A.  LDA >= max(1,N).
92*> \endverbatim
93*>
94*> \param[in] IPIV
95*> \verbatim
96*>          IPIV is INTEGER array, dimension (N)
97*>          Details of the interchanges and the block structure of D,
98*>          as determined by ZSYTRF or ZHETRF.
99*>
100*>          If UPLO = 'U':
101*>               If IPIV(k) > 0, then rows and columns k and IPIV(k)
102*>               were interchanged and D(k,k) is a 1-by-1 diagonal block.
103*>               (If IPIV( k ) = k, no interchange was done).
104*>
105*>               If IPIV(k) = IPIV(k-1) < 0, then rows and
106*>               columns k-1 and -IPIV(k) were interchanged,
107*>               D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
108*>
109*>          If UPLO = 'L':
110*>               If IPIV(k) > 0, then rows and columns k and IPIV(k)
111*>               were interchanged and D(k,k) is a 1-by-1 diagonal block.
112*>               (If IPIV( k ) = k, no interchange was done).
113*>
114*>               If IPIV(k) = IPIV(k+1) < 0, then rows and
115*>               columns k+1 and -IPIV(k) were interchanged,
116*>               D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
117*> \endverbatim
118*>
119*> \param[in,out] B
120*> \verbatim
121*>          B is COMPLEX*16 array, dimension (LDB,NRHS)
122*>          On entry, B contains NRHS vectors of length N.
123*>          On exit, B is overwritten with the product A * B.
124*> \endverbatim
125*>
126*> \param[in] LDB
127*> \verbatim
128*>          LDB is INTEGER
129*>          The leading dimension of the array B.  LDB >= max(1,N).
130*> \endverbatim
131*>
132*> \param[out] INFO
133*> \verbatim
134*>          INFO is INTEGER
135*>          = 0: successful exit
136*>          < 0: if INFO = -k, the k-th argument had an illegal value
137*> \endverbatim
138*
139*  Authors:
140*  ========
141*
142*> \author Univ. of Tennessee
143*> \author Univ. of California Berkeley
144*> \author Univ. of Colorado Denver
145*> \author NAG Ltd.
146*
147*> \date April 2012
148*
149*> \ingroup complex16_lin
150*
151*  =====================================================================
152      SUBROUTINE ZLAVSY( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B,
153     $                   LDB, INFO )
154*
155*  -- LAPACK test routine (version 3.4.1) --
156*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
157*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
158*     April 2012
159*
160*     .. Scalar Arguments ..
161      CHARACTER          DIAG, TRANS, UPLO
162      INTEGER            INFO, LDA, LDB, N, NRHS
163*     ..
164*     .. Array Arguments ..
165      INTEGER            IPIV( * )
166      COMPLEX*16         A( LDA, * ), B( LDB, * )
167*     ..
168*
169*  =====================================================================
170*
171*     .. Parameters ..
172      COMPLEX*16         CONE
173      PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
174*     ..
175*     .. Local Scalars ..
176      LOGICAL            NOUNIT
177      INTEGER            J, K, KP
178      COMPLEX*16         D11, D12, D21, D22, T1, T2
179*     ..
180*     .. External Functions ..
181      LOGICAL            LSAME
182      EXTERNAL           LSAME
183*     ..
184*     .. External Subroutines ..
185      EXTERNAL           XERBLA, ZGEMV, ZGERU, ZSCAL, ZSWAP
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.LSAME( TRANS, 'T' ) )
198     $          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( 'ZLAVSY ', -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 ZSCAL( 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 ZGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
255     $                        LDB, B( 1, 1 ), LDB )
256*
257*                 Interchange if P(K) != I.
258*
259                  KP = IPIV( K )
260                  IF( KP.NE.K )
261     $               CALL ZSWAP( 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 ZGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
290     $                        LDB, B( 1, 1 ), LDB )
291                  CALL ZGERU( K-1, NRHS, CONE, A( 1, K+1 ), 1,
292     $                        B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
293*
294*                 Interchange if P(K) != I.
295*
296                  KP = ABS( IPIV( K ) )
297                  IF( KP.NE.K )
298     $               CALL ZSWAP( 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 ZSCAL( 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 ZGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
337     $                         B( K, 1 ), 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 ZSWAP( 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 ZGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
373     $                        B( K, 1 ), LDB, B( K+1, 1 ), LDB )
374                  CALL ZGERU( N-K, NRHS, CONE, 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 ZSWAP( 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 IF( LSAME( TRANS, 'T' ) ) THEN
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) != I.
415*
416                  KP = IPIV( K )
417                  IF( KP.NE.K )
418     $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
419*
420*                 Apply the transformation
421*
422                  CALL ZGEMV( 'Transpose', K-1, NRHS, CONE, B, LDB,
423     $                        A( 1, K ), 1, CONE, B( K, 1 ), LDB )
424               END IF
425               IF( NOUNIT )
426     $            CALL ZSCAL( 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) != I.
435*
436                  KP = ABS( IPIV( K ) )
437                  IF( KP.NE.K-1 )
438     $               CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
439     $                           LDB )
440*
441*                 Apply the transformations
442*
443                  CALL ZGEMV( 'Transpose', K-2, NRHS, CONE, B, LDB,
444     $                        A( 1, K ), 1, CONE, B( K, 1 ), LDB )
445                  CALL ZGEMV( 'Transpose', K-2, NRHS, CONE, B, LDB,
446     $                        A( 1, K-1 ), 1, CONE, 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) != I.
487*
488                  KP = IPIV( K )
489                  IF( KP.NE.K )
490     $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
491*
492*                 Apply the transformation
493*
494                  CALL ZGEMV( 'Transpose', N-K, NRHS, CONE, B( K+1, 1 ),
495     $                       LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB )
496               END IF
497               IF( NOUNIT )
498     $            CALL ZSCAL( 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) != I.
507*
508                  KP = ABS( IPIV( K ) )
509                  IF( KP.NE.K+1 )
510     $               CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
511     $                           LDB )
512*
513*                 Apply the transformation
514*
515                  CALL ZGEMV( 'Transpose', N-K-1, NRHS, CONE,
516     $                        B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, CONE,
517     $                        B( K+1, 1 ), LDB )
518                  CALL ZGEMV( 'Transpose', N-K-1, NRHS, CONE,
519     $                        B( K+2, 1 ), LDB, A( K+2, K ), 1, CONE,
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      END IF
543      RETURN
544*
545*     End of ZLAVSY
546*
547      END
548