1*> \brief \b CHET01_3
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 CHET01_3( UPLO, N, A, LDA, AFAC, LDAFAC, E, IPIV, C,
12*                            LDC, RWORK, RESID )
13*
14*       .. Scalar Arguments ..
15*       CHARACTER          UPLO
16*       INTEGER            LDA, LDAFAC, LDC, N
17*       REAL               RESID
18*       ..
19*       .. Array Arguments ..
20*       INTEGER            IPIV( * )
21*       REAL               RWORK( * )
22*       COMPLEX            A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
23*                          E( * )
24*       ..
25*
26*
27*> \par Purpose:
28*  =============
29*>
30*> \verbatim
31*>
32*> CHET01_3 reconstructs a Hermitian indefinite matrix A from its
33*> block L*D*L' or U*D*U' factorization computed by CHETRF_RK
34*> (or CHETRF_BK) and computes the residual
35*>    norm( C - A ) / ( N * norm(A) * EPS ),
36*> where C is the reconstructed matrix and EPS is the machine epsilon.
37*> \endverbatim
38*
39*  Arguments:
40*  ==========
41*
42*> \param[in] UPLO
43*> \verbatim
44*>          UPLO is CHARACTER*1
45*>          Specifies whether the upper or lower triangular part of the
46*>          Hermitian matrix A is stored:
47*>          = 'U':  Upper triangular
48*>          = 'L':  Lower triangular
49*> \endverbatim
50*>
51*> \param[in] N
52*> \verbatim
53*>          N is INTEGER
54*>          The number of rows and columns of the matrix A.  N >= 0.
55*> \endverbatim
56*>
57*> \param[in] A
58*> \verbatim
59*>          A is COMPLEX*16 array, dimension (LDA,N)
60*>          The original Hermitian matrix A.
61*> \endverbatim
62*>
63*> \param[in] LDA
64*> \verbatim
65*>          LDA is INTEGER
66*>          The leading dimension of the array A.  LDA >= max(1,N)
67*> \endverbatim
68*>
69*> \param[in] AFAC
70*> \verbatim
71*>          AFAC is COMPLEX array, dimension (LDAFAC,N)
72*>          Diagonal of the block diagonal matrix D and factors U or L
73*>          as computed by CHETRF_RK and CHETRF_BK:
74*>            a) ONLY diagonal elements of the Hermitian block diagonal
75*>               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
76*>               (superdiagonal (or subdiagonal) elements of D
77*>                should be provided on entry in array E), and
78*>            b) If UPLO = 'U': factor U in the superdiagonal part of A.
79*>               If UPLO = 'L': factor L in the subdiagonal part of A.
80*> \endverbatim
81*>
82*> \param[in] LDAFAC
83*> \verbatim
84*>          LDAFAC is INTEGER
85*>          The leading dimension of the array AFAC.
86*>          LDAFAC >= max(1,N).
87*> \endverbatim
88*>
89*> \param[in] E
90*> \verbatim
91*>          E is COMPLEX array, dimension (N)
92*>          On entry, contains the superdiagonal (or subdiagonal)
93*>          elements of the Hermitian block diagonal matrix D
94*>          with 1-by-1 or 2-by-2 diagonal blocks, where
95*>          If UPLO = 'U': E(i) = D(i-1,i),i=2:N, E(1) not referenced;
96*>          If UPLO = 'L': E(i) = D(i+1,i),i=1:N-1, E(N) not referenced.
97*> \endverbatim
98*>
99*> \param[in] IPIV
100*> \verbatim
101*>          IPIV is INTEGER array, dimension (N)
102*>          The pivot indices from CHETRF_RK (or CHETRF_BK).
103*> \endverbatim
104*>
105*> \param[out] C
106*> \verbatim
107*>          C is COMPLEX array, dimension (LDC,N)
108*> \endverbatim
109*>
110*> \param[in] LDC
111*> \verbatim
112*>          LDC is INTEGER
113*>          The leading dimension of the array C.  LDC >= max(1,N).
114*> \endverbatim
115*>
116*> \param[out] RWORK
117*> \verbatim
118*>          RWORK is REAL array, dimension (N)
119*> \endverbatim
120*>
121*> \param[out] RESID
122*> \verbatim
123*>          RESID is REAL
124*>          If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
125*>          If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
126*> \endverbatim
127*
128*  Authors:
129*  ========
130*
131*> \author Univ. of Tennessee
132*> \author Univ. of California Berkeley
133*> \author Univ. of Colorado Denver
134*> \author NAG Ltd.
135*
136*> \date June 2017
137*
138*> \ingroup complex_lin
139*
140*  =====================================================================
141      SUBROUTINE CHET01_3( UPLO, N, A, LDA, AFAC, LDAFAC, E, IPIV, C,
142     $                     LDC, RWORK, RESID )
143*
144*  -- LAPACK test routine (version 3.7.1) --
145*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
146*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
147*     June 2017
148*
149*     .. Scalar Arguments ..
150      CHARACTER          UPLO
151      INTEGER            LDA, LDAFAC, LDC, N
152      REAL               RESID
153*     ..
154*     .. Array Arguments ..
155      INTEGER            IPIV( * )
156      REAL               RWORK( * )
157      COMPLEX            A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
158     $                   E( * )
159*     ..
160*
161*  =====================================================================
162*
163*     .. Parameters ..
164      REAL               ZERO, ONE
165      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
166      COMPLEX            CZERO, CONE
167      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
168     $                   CONE = ( 1.0E+0, 0.0E+0 ) )
169*     ..
170*     .. Local Scalars ..
171      INTEGER            I, INFO, J
172      REAL               ANORM, EPS
173*     ..
174*     .. External Functions ..
175      LOGICAL            LSAME
176      REAL               CLANHE, SLAMCH
177      EXTERNAL           LSAME, CLANHE, SLAMCH
178*     ..
179*     .. External Subroutines ..
180      EXTERNAL           CLASET, CLAVHE_ROOK, CSYCONVF_ROOK
181*     ..
182*     .. Intrinsic Functions ..
183      INTRINSIC          AIMAG, REAL
184*     ..
185*     .. Executable Statements ..
186*
187*     Quick exit if N = 0.
188*
189      IF( N.LE.0 ) THEN
190         RESID = ZERO
191         RETURN
192      END IF
193*
194*     a) Revert to multiplyers of L
195*
196      CALL CSYCONVF_ROOK( UPLO, 'R', N, AFAC, LDAFAC, E, IPIV, INFO )
197*
198*     1) Determine EPS and the norm of A.
199*
200      EPS = SLAMCH( 'Epsilon' )
201      ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK )
202*
203*     Check the imaginary parts of the diagonal elements and return with
204*     an error code if any are nonzero.
205*
206      DO J = 1, N
207         IF( AIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
208            RESID = ONE / EPS
209            RETURN
210         END IF
211      END DO
212*
213*     2) Initialize C to the identity matrix.
214*
215      CALL CLASET( 'Full', N, N, CZERO, CONE, C, LDC )
216*
217*     3) Call CLAVHE_ROOK to form the product D * U' (or D * L' ).
218*
219      CALL CLAVHE_ROOK( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC,
220     $                  LDAFAC, IPIV, C, LDC, INFO )
221*
222*     4) Call ZLAVHE_RK again to multiply by U (or L ).
223*
224      CALL CLAVHE_ROOK( UPLO, 'No transpose', 'Unit', N, N, AFAC,
225     $                  LDAFAC, IPIV, C, LDC, INFO )
226*
227*     5) Compute the difference  C - A .
228*
229      IF( LSAME( UPLO, 'U' ) ) THEN
230         DO J = 1, N
231            DO I = 1, J - 1
232               C( I, J ) = C( I, J ) - A( I, J )
233            END DO
234            C( J, J ) = C( J, J ) - REAL( A( J, J ) )
235         END DO
236      ELSE
237         DO J = 1, N
238            C( J, J ) = C( J, J ) - REAL( A( J, J ) )
239            DO I = J + 1, N
240               C( I, J ) = C( I, J ) - A( I, J )
241            END DO
242         END DO
243      END IF
244*
245*     6) Compute norm( C - A ) / ( N * norm(A) * EPS )
246*
247      RESID = CLANHE( '1', UPLO, N, C, LDC, RWORK )
248*
249      IF( ANORM.LE.ZERO ) THEN
250         IF( RESID.NE.ZERO )
251     $      RESID = ONE / EPS
252      ELSE
253         RESID = ( ( RESID/REAL( N ) )/ANORM ) / EPS
254      END IF
255*
256*     b) Convert to factor of L (or U)
257*
258      CALL CSYCONVF_ROOK( UPLO, 'C', N, AFAC, LDAFAC, E, IPIV, INFO )
259*
260      RETURN
261*
262*     End of CHET01_3
263*
264      END
265