1*> \brief \b CLAVSY
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 CLAVSY( 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            A( LDA, * ), B( LDB, * )
21*       ..
22*
23*
24*> \par Purpose:
25*  =============
26*>
27*> \verbatim
28*>
29*> CLAVSY  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 CSYTRF.
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 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 CSYTRF.
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 CSYTRF or CHETRF.
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 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 complex_lin
150*
151*  =====================================================================
152      SUBROUTINE CLAVSY( 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            A( LDA, * ), B( LDB, * )
167*     ..
168*
169*  =====================================================================
170*
171*     .. Parameters ..
172      COMPLEX            CONE
173      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
174*     ..
175*     .. Local Scalars ..
176      LOGICAL            NOUNIT
177      INTEGER            J, K, KP
178      COMPLEX            D11, D12, D21, D22, T1, T2
179*     ..
180*     .. External Functions ..
181      LOGICAL            LSAME
182      EXTERNAL           LSAME
183*     ..
184*     .. External Subroutines ..
185      EXTERNAL           CGEMV, CGERU, CSCAL, CSWAP, 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.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( 'CLAVSY ', -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 CSCAL( 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 CGERU( 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 CSWAP( 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 CGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
290     $                        LDB, B( 1, 1 ), LDB )
291                  CALL CGERU( 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 CSWAP( 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 CSCAL( 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 CGERU( 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 CSWAP( 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 CGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
373     $                        B( K, 1 ), LDB, B( K+1, 1 ), LDB )
374                  CALL CGERU( 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 CSWAP( 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       IF( K.LT.1 )
406     $         GO TO 90
407*
408*           1 x 1 pivot block.
409*
410            IF( IPIV( K ).GT.0 ) THEN
411               IF( K.GT.1 ) THEN
412*
413*                 Interchange if P(K) != I.
414*
415                  KP = IPIV( K )
416                  IF( KP.NE.K )
417     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
418*
419*                 Apply the transformation
420*
421                  CALL CGEMV( 'Transpose', K-1, NRHS, CONE, B, LDB,
422     $                        A( 1, K ), 1, CONE, B( K, 1 ), LDB )
423               END IF
424               IF( NOUNIT )
425     $            CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
426               K = K - 1
427*
428*           2 x 2 pivot block.
429*
430            ELSE
431               IF( K.GT.2 ) THEN
432*
433*                 Interchange if P(K) != I.
434*
435                  KP = ABS( IPIV( K ) )
436                  IF( KP.NE.K-1 )
437     $               CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
438     $                           LDB )
439*
440*                 Apply the transformations
441*
442                  CALL CGEMV( 'Transpose', K-2, NRHS, CONE, B, LDB,
443     $                        A( 1, K ), 1, CONE, B( K, 1 ), LDB )
444                  CALL CGEMV( 'Transpose', K-2, NRHS, CONE, B, LDB,
445     $                        A( 1, K-1 ), 1, CONE, B( K-1, 1 ), LDB )
446               END IF
447*
448*              Multiply by the diagonal block if non-unit.
449*
450               IF( NOUNIT ) THEN
451                  D11 = A( K-1, K-1 )
452                  D22 = A( K, K )
453                  D12 = A( K-1, K )
454                  D21 = D12
455                  DO 80 J = 1, NRHS
456                     T1 = B( K-1, J )
457                     T2 = B( K, J )
458                     B( K-1, J ) = D11*T1 + D12*T2
459                     B( K, J ) = D21*T1 + D22*T2
460   80             CONTINUE
461               END IF
462               K = K - 2
463            END IF
464            GO TO 70
465   90       CONTINUE
466*
467*        Form  B := L'*B
468*        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
469*        and   L' = inv(L'(m))*P(m)* ... *inv(L'(1))*P(1)
470*
471         ELSE
472*
473*           Loop forward applying the L-transformations.
474*
475            K = 1
476  100       CONTINUE
477            IF( K.GT.N )
478     $         GO TO 120
479*
480*           1 x 1 pivot block
481*
482            IF( IPIV( K ).GT.0 ) THEN
483               IF( K.LT.N ) THEN
484*
485*                 Interchange if P(K) != I.
486*
487                  KP = IPIV( K )
488                  IF( KP.NE.K )
489     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
490*
491*                 Apply the transformation
492*
493                  CALL CGEMV( 'Transpose', N-K, NRHS, CONE, B( K+1, 1 ),
494     $                       LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB )
495               END IF
496               IF( NOUNIT )
497     $            CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
498               K = K + 1
499*
500*           2 x 2 pivot block.
501*
502            ELSE
503               IF( K.LT.N-1 ) THEN
504*
505*              Interchange if P(K) != I.
506*
507                  KP = ABS( IPIV( K ) )
508                  IF( KP.NE.K+1 )
509     $               CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
510     $                           LDB )
511*
512*                 Apply the transformation
513*
514                  CALL CGEMV( 'Transpose', N-K-1, NRHS, CONE,
515     $                        B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, CONE,
516     $                        B( K+1, 1 ), LDB )
517                  CALL CGEMV( 'Transpose', N-K-1, NRHS, CONE,
518     $                        B( K+2, 1 ), LDB, A( K+2, K ), 1, CONE,
519     $                        B( K, 1 ), LDB )
520               END IF
521*
522*              Multiply by the diagonal block if non-unit.
523*
524               IF( NOUNIT ) THEN
525                  D11 = A( K, K )
526                  D22 = A( K+1, K+1 )
527                  D21 = A( K+1, K )
528                  D12 = D21
529                  DO 110 J = 1, NRHS
530                     T1 = B( K, J )
531                     T2 = B( K+1, J )
532                     B( K, J ) = D11*T1 + D12*T2
533                     B( K+1, J ) = D21*T1 + D22*T2
534  110             CONTINUE
535               END IF
536               K = K + 2
537            END IF
538            GO TO 100
539  120       CONTINUE
540         END IF
541      END IF
542      RETURN
543*
544*     End of CLAVSY
545*
546      END
547