1*> \brief \b ZHETRS2
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZHETRS2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrs2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrs2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrs2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZHETRS2( 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*> ZHETRS2 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 ZHETRF 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**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*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 ZHETRF.
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 ZHETRF.
87*> \endverbatim
88*>
89*> \param[in,out] B
90*> \verbatim
91*>          B is COMPLEX*16 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*16 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*> \date June 2016
123*
124*> \ingroup complex16HEcomputational
125*
126*  =====================================================================
127      SUBROUTINE ZHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
128     $                    WORK, INFO )
129*
130*  -- LAPACK computational routine (version 3.7.0) --
131*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
132*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
133*     June 2016
134*
135*     .. Scalar Arguments ..
136      CHARACTER          UPLO
137      INTEGER            INFO, LDA, LDB, N, NRHS
138*     ..
139*     .. Array Arguments ..
140      INTEGER            IPIV( * )
141      COMPLEX*16       A( LDA, * ), B( LDB, * ), WORK( * )
142*     ..
143*
144*  =====================================================================
145*
146*     .. Parameters ..
147      COMPLEX*16         ONE
148      PARAMETER          ( ONE = (1.0D+0,0.0D+0) )
149*     ..
150*     .. Local Scalars ..
151      LOGICAL            UPPER
152      INTEGER            I, IINFO, J, K, KP
153      DOUBLE PRECISION   S
154      COMPLEX*16         AK, AKM1, AKM1K, BK, BKM1, DENOM
155*     ..
156*     .. External Functions ..
157      LOGICAL            LSAME
158      EXTERNAL           LSAME
159*     ..
160*     .. External Subroutines ..
161      EXTERNAL           ZDSCAL, ZSYCONV, ZSWAP, ZTRSM, XERBLA
162*     ..
163*     .. Intrinsic Functions ..
164      INTRINSIC          DBLE, DCONJG, MAX
165*     ..
166*     .. Executable Statements ..
167*
168      INFO = 0
169      UPPER = LSAME( UPLO, 'U' )
170      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
171         INFO = -1
172      ELSE IF( N.LT.0 ) THEN
173         INFO = -2
174      ELSE IF( NRHS.LT.0 ) THEN
175         INFO = -3
176      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
177         INFO = -5
178      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
179         INFO = -8
180      END IF
181      IF( INFO.NE.0 ) THEN
182         CALL XERBLA( 'ZHETRS2', -INFO )
183         RETURN
184      END IF
185*
186*     Quick return if possible
187*
188      IF( N.EQ.0 .OR. NRHS.EQ.0 )
189     $   RETURN
190*
191*     Convert A
192*
193      CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
194*
195      IF( UPPER ) THEN
196*
197*        Solve A*X = B, where A = U*D*U**H.
198*
199*       P**T * B
200        K=N
201        DO WHILE ( K .GE. 1 )
202         IF( IPIV( K ).GT.0 ) THEN
203*           1 x 1 diagonal block
204*           Interchange rows K and IPIV(K).
205            KP = IPIV( K )
206            IF( KP.NE.K )
207     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
208            K=K-1
209         ELSE
210*           2 x 2 diagonal block
211*           Interchange rows K-1 and -IPIV(K).
212            KP = -IPIV( K )
213            IF( KP.EQ.-IPIV( K-1 ) )
214     $         CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
215            K=K-2
216         END IF
217        END DO
218*
219*  Compute (U \P**T * B) -> B    [ (U \P**T * B) ]
220*
221        CALL ZTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
222*
223*  Compute D \ B -> B   [ D \ (U \P**T * B) ]
224*
225         I=N
226         DO WHILE ( I .GE. 1 )
227            IF( IPIV(I) .GT. 0 ) THEN
228              S = DBLE( ONE ) / DBLE( A( I, I ) )
229              CALL ZDSCAL( NRHS, S, 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 ) / DCONJG( AKM1K )
235                  DENOM = AKM1*AK - ONE
236                  DO 15 J = 1, NRHS
237                     BKM1 = B( I-1, J ) / AKM1K
238                     BK = B( I, J ) / DCONJG( 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**H \ B) -> B   [ U**H \ (D \ (U \P**T * B) ) ]
249*
250         CALL ZTRSM('L','U','C','U',N,NRHS,ONE,A,LDA,B,LDB)
251*
252*       P * B  [ P * (U**H \ (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**H.
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              S = DBLE( ONE ) / DBLE( A( I, I ) )
307              CALL ZDSCAL( NRHS, S, B( I, 1 ), LDB )
308            ELSE
309                  AKM1K = WORK(I)
310                  AKM1 = A( I, I ) / DCONJG( AKM1K )
311                  AK = A( I+1, I+1 ) / AKM1K
312                  DENOM = AKM1*AK - ONE
313                  DO 25 J = 1, NRHS
314                     BKM1 = B( I, J ) / DCONJG( AKM1K )
315                     BK = B( I+1, J ) / AKM1K
316                     B( I, J ) = ( AK*BKM1-BK ) / DENOM
317                     B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
318 25              CONTINUE
319                  I = I + 1
320            ENDIF
321            I = I + 1
322         END DO
323*
324*  Compute (L**H \ B) -> B   [ L**H \ (D \ (L \P**T * B) ) ]
325*
326        CALL ZTRSM('L','L','C','U',N,NRHS,ONE,A,LDA,B,LDB)
327*
328*       P * B  [ P * (L**H \ (D \ (L \P**T * B) )) ]
329*
330        K=N
331        DO WHILE ( K .GE. 1 )
332         IF( IPIV( K ).GT.0 ) THEN
333*           1 x 1 diagonal block
334*           Interchange rows K and IPIV(K).
335            KP = IPIV( K )
336            IF( KP.NE.K )
337     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
338            K=K-1
339         ELSE
340*           2 x 2 diagonal block
341*           Interchange rows K-1 and -IPIV(K).
342            KP = -IPIV( K )
343            IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
344     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
345            K=K-2
346         ENDIF
347        END DO
348*
349      END IF
350*
351*     Revert A
352*
353      CALL ZSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
354*
355      RETURN
356*
357*     End of ZHETRS2
358*
359      END
360