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