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*>          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 CSYTRF.
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 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*> \date November 2013
149*
150*> \ingroup complex_lin
151*
152*  =====================================================================
153      SUBROUTINE CLAVSY( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B,
154     $                   LDB, INFO )
155*
156*  -- LAPACK test routine (version 3.5.0) --
157*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
158*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
159*     November 2013
160*
161*     .. Scalar Arguments ..
162      CHARACTER          DIAG, TRANS, UPLO
163      INTEGER            INFO, LDA, LDB, N, NRHS
164*     ..
165*     .. Array Arguments ..
166      INTEGER            IPIV( * )
167      COMPLEX            A( LDA, * ), B( LDB, * )
168*     ..
169*
170*  =====================================================================
171*
172*     .. Parameters ..
173      COMPLEX            CONE
174      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
175*     ..
176*     .. Local Scalars ..
177      LOGICAL            NOUNIT
178      INTEGER            J, K, KP
179      COMPLEX            D11, D12, D21, D22, T1, T2
180*     ..
181*     .. External Functions ..
182      LOGICAL            LSAME
183      EXTERNAL           LSAME
184*     ..
185*     .. External Subroutines ..
186      EXTERNAL           CGEMV, CGERU, CSCAL, CSWAP, XERBLA
187*     ..
188*     .. Intrinsic Functions ..
189      INTRINSIC          ABS, MAX
190*     ..
191*     .. Executable Statements ..
192*
193*     Test the input parameters.
194*
195      INFO = 0
196      IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
197         INFO = -1
198      ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'T' ) )
199     $          THEN
200         INFO = -2
201      ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) )
202     $          THEN
203         INFO = -3
204      ELSE IF( N.LT.0 ) THEN
205         INFO = -4
206      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
207         INFO = -6
208      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
209         INFO = -9
210      END IF
211      IF( INFO.NE.0 ) THEN
212         CALL XERBLA( 'CLAVSY ', -INFO )
213         RETURN
214      END IF
215*
216*     Quick return if possible.
217*
218      IF( N.EQ.0 )
219     $   RETURN
220*
221      NOUNIT = LSAME( DIAG, 'N' )
222*------------------------------------------
223*
224*     Compute  B := A * B  (No transpose)
225*
226*------------------------------------------
227      IF( LSAME( TRANS, 'N' ) ) THEN
228*
229*        Compute  B := U*B
230*        where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
231*
232         IF( LSAME( UPLO, 'U' ) ) THEN
233*
234*        Loop forward applying the transformations.
235*
236            K = 1
237   10       CONTINUE
238            IF( K.GT.N )
239     $         GO TO 30
240            IF( IPIV( K ).GT.0 ) THEN
241*
242*              1 x 1 pivot block
243*
244*              Multiply by the diagonal element if forming U * D.
245*
246               IF( NOUNIT )
247     $            CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
248*
249*              Multiply by  P(K) * inv(U(K))  if K > 1.
250*
251               IF( K.GT.1 ) THEN
252*
253*                 Apply the transformation.
254*
255                  CALL CGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
256     $                        LDB, B( 1, 1 ), LDB )
257*
258*                 Interchange if P(K) != I.
259*
260                  KP = IPIV( K )
261                  IF( KP.NE.K )
262     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
263               END IF
264               K = K + 1
265            ELSE
266*
267*              2 x 2 pivot block
268*
269*              Multiply by the diagonal block if forming U * D.
270*
271               IF( NOUNIT ) THEN
272                  D11 = A( K, K )
273                  D22 = A( K+1, K+1 )
274                  D12 = A( K, K+1 )
275                  D21 = D12
276                  DO 20 J = 1, NRHS
277                     T1 = B( K, J )
278                     T2 = B( K+1, J )
279                     B( K, J ) = D11*T1 + D12*T2
280                     B( K+1, J ) = D21*T1 + D22*T2
281   20             CONTINUE
282               END IF
283*
284*              Multiply by  P(K) * inv(U(K))  if K > 1.
285*
286               IF( K.GT.1 ) THEN
287*
288*                 Apply the transformations.
289*
290                  CALL CGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
291     $                        LDB, B( 1, 1 ), LDB )
292                  CALL CGERU( K-1, NRHS, CONE, A( 1, K+1 ), 1,
293     $                        B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
294*
295*                 Interchange if P(K) != I.
296*
297                  KP = ABS( IPIV( K ) )
298                  IF( KP.NE.K )
299     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
300               END IF
301               K = K + 2
302            END IF
303            GO TO 10
304   30       CONTINUE
305*
306*        Compute  B := L*B
307*        where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
308*
309         ELSE
310*
311*           Loop backward applying the transformations to B.
312*
313            K = N
314   40       CONTINUE
315            IF( K.LT.1 )
316     $         GO TO 60
317*
318*           Test the pivot index.  If greater than zero, a 1 x 1
319*           pivot was used, otherwise a 2 x 2 pivot was used.
320*
321            IF( IPIV( K ).GT.0 ) THEN
322*
323*              1 x 1 pivot block:
324*
325*              Multiply by the diagonal element if forming L * D.
326*
327               IF( NOUNIT )
328     $            CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
329*
330*              Multiply by  P(K) * inv(L(K))  if K < N.
331*
332               IF( K.NE.N ) THEN
333                  KP = IPIV( K )
334*
335*                 Apply the transformation.
336*
337                  CALL CGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
338     $                        B( K, 1 ), LDB, B( K+1, 1 ), LDB )
339*
340*                 Interchange if a permutation was applied at the
341*                 K-th step of the factorization.
342*
343                  IF( KP.NE.K )
344     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
345               END IF
346               K = K - 1
347*
348            ELSE
349*
350*              2 x 2 pivot block:
351*
352*              Multiply by the diagonal block if forming L * D.
353*
354               IF( NOUNIT ) THEN
355                  D11 = A( K-1, K-1 )
356                  D22 = A( K, K )
357                  D21 = A( K, K-1 )
358                  D12 = D21
359                  DO 50 J = 1, NRHS
360                     T1 = B( K-1, J )
361                     T2 = B( K, J )
362                     B( K-1, J ) = D11*T1 + D12*T2
363                     B( K, J ) = D21*T1 + D22*T2
364   50             CONTINUE
365               END IF
366*
367*              Multiply by  P(K) * inv(L(K))  if K < N.
368*
369               IF( K.NE.N ) THEN
370*
371*                 Apply the transformation.
372*
373                  CALL CGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
374     $                        B( K, 1 ), LDB, B( K+1, 1 ), LDB )
375                  CALL CGERU( N-K, NRHS, CONE, A( K+1, K-1 ), 1,
376     $                        B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
377*
378*                 Interchange if a permutation was applied at the
379*                 K-th step of the factorization.
380*
381                  KP = ABS( IPIV( K ) )
382                  IF( KP.NE.K )
383     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
384               END IF
385               K = K - 2
386            END IF
387            GO TO 40
388   60       CONTINUE
389         END IF
390*----------------------------------------
391*
392*     Compute  B := A' * B  (transpose)
393*
394*----------------------------------------
395      ELSE IF( LSAME( TRANS, 'T' ) ) THEN
396*
397*        Form  B := U'*B
398*        where U  = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
399*        and   U' = inv(U'(1))*P(1)* ... *inv(U'(m))*P(m)
400*
401         IF( LSAME( UPLO, 'U' ) ) THEN
402*
403*           Loop backward applying the transformations.
404*
405            K = N
406   70       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 CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
419*
420*                 Apply the transformation
421*
422                  CALL CGEMV( '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 CSCAL( 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 CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
439     $                           LDB )
440*
441*                 Apply the transformations
442*
443                  CALL CGEMV( 'Transpose', K-2, NRHS, CONE, B, LDB,
444     $                        A( 1, K ), 1, CONE, B( K, 1 ), LDB )
445                  CALL CGEMV( '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 CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
491*
492*                 Apply the transformation
493*
494                  CALL CGEMV( '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 CSCAL( 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 CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
511     $                           LDB )
512*
513*                 Apply the transformation
514*
515                  CALL CGEMV( '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 CGEMV( '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 CLAVSY
546*
547      END
548