1*> \brief \b CHETRS_3
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CHETRS_3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetrs_3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetrs_3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetrs_3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CHETRS_3( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB,
22*                            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, * ), E( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*> CHETRS_3 solves a system of linear equations A * X = B with a complex
39*> Hermitian matrix A using the factorization computed
40*> by CHETRF_RK or CHETRF_BK:
41*>
42*>    A = P*U*D*(U**H)*(P**T) or A = P*L*D*(L**H)*(P**T),
43*>
44*> where U (or L) is unit upper (or lower) triangular matrix,
45*> U**H (or L**H) is the conjugate of U (or L), P is a permutation
46*> matrix, P**T is the transpose of P, and D is Hermitian and block
47*> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
48*>
49*> This algorithm is using Level 3 BLAS.
50*> \endverbatim
51*
52*  Arguments:
53*  ==========
54*
55*> \param[in] UPLO
56*> \verbatim
57*>          UPLO is CHARACTER*1
58*>          Specifies whether the details of the factorization are
59*>          stored as an upper or lower triangular matrix:
60*>          = 'U':  Upper triangular, form is A = P*U*D*(U**H)*(P**T);
61*>          = 'L':  Lower triangular, form is A = P*L*D*(L**H)*(P**T).
62*> \endverbatim
63*>
64*> \param[in] N
65*> \verbatim
66*>          N is INTEGER
67*>          The order of the matrix A.  N >= 0.
68*> \endverbatim
69*>
70*> \param[in] NRHS
71*> \verbatim
72*>          NRHS is INTEGER
73*>          The number of right hand sides, i.e., the number of columns
74*>          of the matrix B.  NRHS >= 0.
75*> \endverbatim
76*>
77*> \param[in] A
78*> \verbatim
79*>          A is COMPLEX array, dimension (LDA,N)
80*>          Diagonal of the block diagonal matrix D and factors U or L
81*>          as computed by CHETRF_RK and CHETRF_BK:
82*>            a) ONLY diagonal elements of the Hermitian block diagonal
83*>               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
84*>               (superdiagonal (or subdiagonal) elements of D
85*>                should be provided on entry in array E), and
86*>            b) If UPLO = 'U': factor U in the superdiagonal part of A.
87*>               If UPLO = 'L': factor L in the subdiagonal part of A.
88*> \endverbatim
89*>
90*> \param[in] LDA
91*> \verbatim
92*>          LDA is INTEGER
93*>          The leading dimension of the array A.  LDA >= max(1,N).
94*> \endverbatim
95*>
96*> \param[in] E
97*> \verbatim
98*>          E is COMPLEX array, dimension (N)
99*>          On entry, contains the superdiagonal (or subdiagonal)
100*>          elements of the Hermitian block diagonal matrix D
101*>          with 1-by-1 or 2-by-2 diagonal blocks, where
102*>          If UPLO = 'U': E(i) = D(i-1,i),i=2:N, E(1) not referenced;
103*>          If UPLO = 'L': E(i) = D(i+1,i),i=1:N-1, E(N) not referenced.
104*>
105*>          NOTE: For 1-by-1 diagonal block D(k), where
106*>          1 <= k <= N, the element E(k) is not referenced in both
107*>          UPLO = 'U' or UPLO = 'L' cases.
108*> \endverbatim
109*>
110*> \param[in] IPIV
111*> \verbatim
112*>          IPIV is INTEGER array, dimension (N)
113*>          Details of the interchanges and the block structure of D
114*>          as determined by CHETRF_RK or CHETRF_BK.
115*> \endverbatim
116*>
117*> \param[in,out] B
118*> \verbatim
119*>          B is COMPLEX array, dimension (LDB,NRHS)
120*>          On entry, the right hand side matrix B.
121*>          On exit, the solution matrix X.
122*> \endverbatim
123*>
124*> \param[in] LDB
125*> \verbatim
126*>          LDB is INTEGER
127*>          The leading dimension of the array B.  LDB >= max(1,N).
128*> \endverbatim
129*>
130*> \param[out] INFO
131*> \verbatim
132*>          INFO is INTEGER
133*>          = 0:  successful exit
134*>          < 0:  if INFO = -i, the i-th argument had an illegal value
135*> \endverbatim
136*
137*  Authors:
138*  ========
139*
140*> \author Univ. of Tennessee
141*> \author Univ. of California Berkeley
142*> \author Univ. of Colorado Denver
143*> \author NAG Ltd.
144*
145*> \ingroup complexHEcomputational
146*
147*> \par Contributors:
148*  ==================
149*>
150*> \verbatim
151*>
152*>  June 2017,  Igor Kozachenko,
153*>                  Computer Science Division,
154*>                  University of California, Berkeley
155*>
156*>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
157*>                  School of Mathematics,
158*>                  University of Manchester
159*>
160*> \endverbatim
161*
162*  =====================================================================
163      SUBROUTINE CHETRS_3( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB,
164     $                     INFO )
165*
166*  -- LAPACK computational routine --
167*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
168*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169*
170*     .. Scalar Arguments ..
171      CHARACTER          UPLO
172      INTEGER            INFO, LDA, LDB, N, NRHS
173*     ..
174*     .. Array Arguments ..
175      INTEGER            IPIV( * )
176      COMPLEX            A( LDA, * ), B( LDB, * ), E( * )
177*     ..
178*
179*  =====================================================================
180*
181*     .. Parameters ..
182      COMPLEX            ONE
183      PARAMETER          ( ONE = ( 1.0E+0,0.0E+0 ) )
184*     ..
185*     .. Local Scalars ..
186      LOGICAL            UPPER
187      INTEGER            I, J, K, KP
188      REAL               S
189      COMPLEX            AK, AKM1, AKM1K, BK, BKM1, DENOM
190*     ..
191*     .. External Functions ..
192      LOGICAL            LSAME
193      EXTERNAL           LSAME
194*     ..
195*     .. External Subroutines ..
196      EXTERNAL           CSSCAL, CSWAP, CTRSM, XERBLA
197*     ..
198*     .. Intrinsic Functions ..
199      INTRINSIC          ABS, CONJG, MAX, REAL
200*     ..
201*     .. Executable Statements ..
202*
203      INFO = 0
204      UPPER = LSAME( UPLO, 'U' )
205      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
206         INFO = -1
207      ELSE IF( N.LT.0 ) THEN
208         INFO = -2
209      ELSE IF( NRHS.LT.0 ) THEN
210         INFO = -3
211      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
212         INFO = -5
213      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
214         INFO = -9
215      END IF
216      IF( INFO.NE.0 ) THEN
217         CALL XERBLA( 'CHETRS_3', -INFO )
218         RETURN
219      END IF
220*
221*     Quick return if possible
222*
223      IF( N.EQ.0 .OR. NRHS.EQ.0 )
224     $   RETURN
225*
226      IF( UPPER ) THEN
227*
228*        Begin Upper
229*
230*        Solve A*X = B, where A = U*D*U**H.
231*
232*        P**T * B
233*
234*        Interchange rows K and IPIV(K) of matrix B in the same order
235*        that the formation order of IPIV(I) vector for Upper case.
236*
237*        (We can do the simple loop over IPIV with decrement -1,
238*        since the ABS value of IPIV(I) represents the row index
239*        of the interchange with row i in both 1x1 and 2x2 pivot cases)
240*
241         DO K = N, 1, -1
242            KP = ABS( IPIV( K ) )
243            IF( KP.NE.K ) THEN
244               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
245            END IF
246         END DO
247*
248*        Compute (U \P**T * B) -> B    [ (U \P**T * B) ]
249*
250         CALL CTRSM( 'L', 'U', 'N', 'U', N, NRHS, ONE, A, LDA, B, LDB )
251*
252*        Compute D \ B -> B   [ D \ (U \P**T * B) ]
253*
254         I = N
255         DO WHILE ( I.GE.1 )
256            IF( IPIV( I ).GT.0 ) THEN
257               S = REAL( ONE ) / REAL( A( I, I ) )
258               CALL CSSCAL( NRHS, S, B( I, 1 ), LDB )
259            ELSE IF ( I.GT.1 ) THEN
260               AKM1K = E( I )
261               AKM1 = A( I-1, I-1 ) / AKM1K
262               AK = A( I, I ) / CONJG( AKM1K )
263               DENOM = AKM1*AK - ONE
264               DO J = 1, NRHS
265                  BKM1 = B( I-1, J ) / AKM1K
266                  BK = B( I, J ) / CONJG( AKM1K )
267                  B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
268                  B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
269               END DO
270               I = I - 1
271            END IF
272            I = I - 1
273         END DO
274*
275*        Compute (U**H \ B) -> B   [ U**H \ (D \ (U \P**T * B) ) ]
276*
277         CALL CTRSM( 'L', 'U', 'C', 'U', N, NRHS, ONE, A, LDA, B, LDB )
278*
279*        P * B  [ P * (U**H \ (D \ (U \P**T * B) )) ]
280*
281*        Interchange rows K and IPIV(K) of matrix B in reverse order
282*        from the formation order of IPIV(I) vector for Upper case.
283*
284*        (We can do the simple loop over IPIV with increment 1,
285*        since the ABS value of IPIV(I) represents the row index
286*        of the interchange with row i in both 1x1 and 2x2 pivot cases)
287*
288         DO K = 1, N, 1
289            KP = ABS( IPIV( K ) )
290            IF( KP.NE.K ) THEN
291               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
292            END IF
293         END DO
294*
295      ELSE
296*
297*        Begin Lower
298*
299*        Solve A*X = B, where A = L*D*L**H.
300*
301*        P**T * B
302*        Interchange rows K and IPIV(K) of matrix B in the same order
303*        that the formation order of IPIV(I) vector for Lower case.
304*
305*        (We can do the simple loop over IPIV with increment 1,
306*        since the ABS value of IPIV(I) represents the row index
307*        of the interchange with row i in both 1x1 and 2x2 pivot cases)
308*
309         DO K = 1, N, 1
310            KP = ABS( IPIV( K ) )
311            IF( KP.NE.K ) THEN
312               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
313            END IF
314         END DO
315*
316*        Compute (L \P**T * B) -> B    [ (L \P**T * B) ]
317*
318         CALL CTRSM( 'L', 'L', 'N', 'U', N, NRHS, ONE, A, LDA, B, LDB )
319*
320*        Compute D \ B -> B   [ D \ (L \P**T * B) ]
321*
322         I = 1
323         DO WHILE ( I.LE.N )
324            IF( IPIV( I ).GT.0 ) THEN
325               S = REAL( ONE ) / REAL( A( I, I ) )
326               CALL CSSCAL( NRHS, S, B( I, 1 ), LDB )
327            ELSE IF( I.LT.N ) THEN
328               AKM1K = E( I )
329               AKM1 = A( I, I ) / CONJG( AKM1K )
330               AK = A( I+1, I+1 ) / AKM1K
331               DENOM = AKM1*AK - ONE
332               DO  J = 1, NRHS
333                  BKM1 = B( I, J ) / CONJG( AKM1K )
334                  BK = B( I+1, J ) / AKM1K
335                  B( I, J ) = ( AK*BKM1-BK ) / DENOM
336                  B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
337               END DO
338               I = I + 1
339            END IF
340            I = I + 1
341         END DO
342*
343*        Compute (L**H \ B) -> B   [ L**H \ (D \ (L \P**T * B) ) ]
344*
345         CALL CTRSM('L', 'L', 'C', 'U', N, NRHS, ONE, A, LDA, B, LDB )
346*
347*        P * B  [ P * (L**H \ (D \ (L \P**T * B) )) ]
348*
349*        Interchange rows K and IPIV(K) of matrix B in reverse order
350*        from the formation order of IPIV(I) vector for Lower case.
351*
352*        (We can do the simple loop over IPIV with decrement -1,
353*        since the ABS value of IPIV(I) represents the row index
354*        of the interchange with row i in both 1x1 and 2x2 pivot cases)
355*
356         DO K = N, 1, -1
357            KP = ABS( IPIV( K ) )
358            IF( KP.NE.K ) THEN
359               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
360            END IF
361         END DO
362*
363*        END Lower
364*
365      END IF
366*
367      RETURN
368*
369*     End of CHETRS_3
370*
371      END
372