1*> \brief \b ZSYTRS2
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZSYTRS2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytrs2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytrs2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytrs2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
22*                           WORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          UPLO
26*       INTEGER            INFO, LDA, LDB, N, NRHS
27*       ..
28*       .. Array Arguments ..
29*       INTEGER            IPIV( * )
30*       COMPLEX*16       A( LDA, * ), B( LDB, * ), WORK( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> ZSYTRS2 solves a system of linear equations A*X = B with a complex
40*> symmetric matrix A using the factorization A = U*D*U**T or
41*> A = L*D*L**T computed by ZSYTRF and converted by ZSYCONV.
42*> \endverbatim
43*
44*  Arguments:
45*  ==========
46*
47*> \param[in] UPLO
48*> \verbatim
49*>          UPLO is CHARACTER*1
50*>          Specifies whether the details of the factorization are stored
51*>          as an upper or lower triangular matrix.
52*>          = 'U':  Upper triangular, form is A = U*D*U**T;
53*>          = 'L':  Lower triangular, form is A = L*D*L**T.
54*> \endverbatim
55*>
56*> \param[in] N
57*> \verbatim
58*>          N is INTEGER
59*>          The order of the matrix A.  N >= 0.
60*> \endverbatim
61*>
62*> \param[in] NRHS
63*> \verbatim
64*>          NRHS is INTEGER
65*>          The number of right hand sides, i.e., the number of columns
66*>          of the matrix B.  NRHS >= 0.
67*> \endverbatim
68*>
69*> \param[in,out] A
70*> \verbatim
71*>          A is COMPLEX*16 array, dimension (LDA,N)
72*>          The block diagonal matrix D and the multipliers used to
73*>          obtain the factor U or L as computed by ZSYTRF.
74*>          Note that A is input / output. This might be counter-intuitive,
75*>          and one may think that A is input only. A is input / output. This
76*>          is because, at the start of the subroutine, we permute A in a
77*>          "better" form and then we permute A back to its original form at
78*>          the end.
79*> \endverbatim
80*>
81*> \param[in] LDA
82*> \verbatim
83*>          LDA is INTEGER
84*>          The leading dimension of the array A.  LDA >= max(1,N).
85*> \endverbatim
86*>
87*> \param[in] IPIV
88*> \verbatim
89*>          IPIV is INTEGER array, dimension (N)
90*>          Details of the interchanges and the block structure of D
91*>          as determined by ZSYTRF.
92*> \endverbatim
93*>
94*> \param[in,out] B
95*> \verbatim
96*>          B is COMPLEX*16 array, dimension (LDB,NRHS)
97*>          On entry, the right hand side matrix B.
98*>          On exit, the solution matrix X.
99*> \endverbatim
100*>
101*> \param[in] LDB
102*> \verbatim
103*>          LDB is INTEGER
104*>          The leading dimension of the array B.  LDB >= max(1,N).
105*> \endverbatim
106*>
107*> \param[out] WORK
108*> \verbatim
109*>          WORK is COMPLEX*16 array, dimension (N)
110*> \endverbatim
111*>
112*> \param[out] INFO
113*> \verbatim
114*>          INFO is INTEGER
115*>          = 0:  successful exit
116*>          < 0:  if INFO = -i, the i-th argument had an illegal value
117*> \endverbatim
118*
119*  Authors:
120*  ========
121*
122*> \author Univ. of Tennessee
123*> \author Univ. of California Berkeley
124*> \author Univ. of Colorado Denver
125*> \author NAG Ltd.
126*
127*> \ingroup complex16SYcomputational
128*
129*  =====================================================================
130      SUBROUTINE ZSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
131     $                    WORK, INFO )
132*
133*  -- LAPACK computational routine --
134*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
135*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136*
137*     .. Scalar Arguments ..
138      CHARACTER          UPLO
139      INTEGER            INFO, LDA, LDB, N, NRHS
140*     ..
141*     .. Array Arguments ..
142      INTEGER            IPIV( * )
143      COMPLEX*16       A( LDA, * ), B( LDB, * ), WORK( * )
144*     ..
145*
146*  =====================================================================
147*
148*     .. Parameters ..
149      COMPLEX*16         ONE
150      PARAMETER          ( ONE = (1.0D+0,0.0D+0) )
151*     ..
152*     .. Local Scalars ..
153      LOGICAL            UPPER
154      INTEGER            I, IINFO, J, K, KP
155      COMPLEX*16         AK, AKM1, AKM1K, BK, BKM1, DENOM
156*     ..
157*     .. External Functions ..
158      LOGICAL            LSAME
159      EXTERNAL           LSAME
160*     ..
161*     .. External Subroutines ..
162      EXTERNAL           ZSCAL, ZSYCONV, ZSWAP, ZTRSM, XERBLA
163*     ..
164*     .. Intrinsic Functions ..
165      INTRINSIC          MAX
166*     ..
167*     .. Executable Statements ..
168*
169      INFO = 0
170      UPPER = LSAME( UPLO, 'U' )
171      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
172         INFO = -1
173      ELSE IF( N.LT.0 ) THEN
174         INFO = -2
175      ELSE IF( NRHS.LT.0 ) THEN
176         INFO = -3
177      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
178         INFO = -5
179      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
180         INFO = -8
181      END IF
182      IF( INFO.NE.0 ) THEN
183         CALL XERBLA( 'ZSYTRS2', -INFO )
184         RETURN
185      END IF
186*
187*     Quick return if possible
188*
189      IF( N.EQ.0 .OR. NRHS.EQ.0 )
190     $   RETURN
191*
192*     Convert A
193*
194      CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
195*
196      IF( UPPER ) THEN
197*
198*        Solve A*X = B, where A = U*D*U**T.
199*
200*       P**T * B
201        K=N
202        DO WHILE ( K .GE. 1 )
203         IF( IPIV( K ).GT.0 ) THEN
204*           1 x 1 diagonal block
205*           Interchange rows K and IPIV(K).
206            KP = IPIV( K )
207            IF( KP.NE.K )
208     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
209            K=K-1
210         ELSE
211*           2 x 2 diagonal block
212*           Interchange rows K-1 and -IPIV(K).
213            KP = -IPIV( K )
214            IF( KP.EQ.-IPIV( K-1 ) )
215     $         CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
216            K=K-2
217         END IF
218        END DO
219*
220*  Compute (U \P**T * B) -> B    [ (U \P**T * B) ]
221*
222        CALL ZTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
223*
224*  Compute D \ B -> B   [ D \ (U \P**T * B) ]
225*
226         I=N
227         DO WHILE ( I .GE. 1 )
228            IF( IPIV(I) .GT. 0 ) THEN
229              CALL ZSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
230            ELSEIF ( I .GT. 1) THEN
231               IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
232                  AKM1K = WORK(I)
233                  AKM1 = A( I-1, I-1 ) / AKM1K
234                  AK = A( I, I ) / AKM1K
235                  DENOM = AKM1*AK - ONE
236                  DO 15 J = 1, NRHS
237                     BKM1 = B( I-1, J ) / AKM1K
238                     BK = B( I, J ) / AKM1K
239                     B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
240                     B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
241 15              CONTINUE
242               I = I - 1
243               ENDIF
244            ENDIF
245            I = I - 1
246         END DO
247*
248*      Compute (U**T \ B) -> B   [ U**T \ (D \ (U \P**T * B) ) ]
249*
250         CALL ZTRSM('L','U','T','U',N,NRHS,ONE,A,LDA,B,LDB)
251*
252*       P * B  [ P * (U**T \ (D \ (U \P**T * B) )) ]
253*
254        K=1
255        DO WHILE ( K .LE. N )
256         IF( IPIV( K ).GT.0 ) THEN
257*           1 x 1 diagonal block
258*           Interchange rows K and IPIV(K).
259            KP = IPIV( K )
260            IF( KP.NE.K )
261     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
262            K=K+1
263         ELSE
264*           2 x 2 diagonal block
265*           Interchange rows K-1 and -IPIV(K).
266            KP = -IPIV( K )
267            IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
268     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
269            K=K+2
270         ENDIF
271        END DO
272*
273      ELSE
274*
275*        Solve A*X = B, where A = L*D*L**T.
276*
277*       P**T * B
278        K=1
279        DO WHILE ( K .LE. N )
280         IF( IPIV( K ).GT.0 ) THEN
281*           1 x 1 diagonal block
282*           Interchange rows K and IPIV(K).
283            KP = IPIV( K )
284            IF( KP.NE.K )
285     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
286            K=K+1
287         ELSE
288*           2 x 2 diagonal block
289*           Interchange rows K and -IPIV(K+1).
290            KP = -IPIV( K+1 )
291            IF( KP.EQ.-IPIV( K ) )
292     $         CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
293            K=K+2
294         ENDIF
295        END DO
296*
297*  Compute (L \P**T * B) -> B    [ (L \P**T * B) ]
298*
299        CALL ZTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
300*
301*  Compute D \ B -> B   [ D \ (L \P**T * B) ]
302*
303         I=1
304         DO WHILE ( I .LE. N )
305            IF( IPIV(I) .GT. 0 ) THEN
306              CALL ZSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
307            ELSE
308                  AKM1K = WORK(I)
309                  AKM1 = A( I, I ) / AKM1K
310                  AK = A( I+1, I+1 ) / AKM1K
311                  DENOM = AKM1*AK - ONE
312                  DO 25 J = 1, NRHS
313                     BKM1 = B( I, J ) / AKM1K
314                     BK = B( I+1, J ) / AKM1K
315                     B( I, J ) = ( AK*BKM1-BK ) / DENOM
316                     B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
317 25              CONTINUE
318                  I = I + 1
319            ENDIF
320            I = I + 1
321         END DO
322*
323*  Compute (L**T \ B) -> B   [ L**T \ (D \ (L \P**T * B) ) ]
324*
325        CALL ZTRSM('L','L','T','U',N,NRHS,ONE,A,LDA,B,LDB)
326*
327*       P * B  [ P * (L**T \ (D \ (L \P**T * B) )) ]
328*
329        K=N
330        DO WHILE ( K .GE. 1 )
331         IF( IPIV( K ).GT.0 ) THEN
332*           1 x 1 diagonal block
333*           Interchange rows K and IPIV(K).
334            KP = IPIV( K )
335            IF( KP.NE.K )
336     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
337            K=K-1
338         ELSE
339*           2 x 2 diagonal block
340*           Interchange rows K-1 and -IPIV(K).
341            KP = -IPIV( K )
342            IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
343     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
344            K=K-2
345         ENDIF
346        END DO
347*
348      END IF
349*
350*     Revert A
351*
352      CALL ZSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
353*
354      RETURN
355*
356*     End of ZSYTRS2
357*
358      END
359