1*> \brief \b ZHET01_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 ZHET01_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*       DOUBLE PRECISION   RESID
18*       ..
19*       .. Array Arguments ..
20*       INTEGER            IPIV( * )
21*       DOUBLE PRECISION   RWORK( * )
22*       COMPLEX*16         A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
23*                          E( * )
24*       ..
25*
26*
27*> \par Purpose:
28*  =============
29*>
30*> \verbatim
31*>
32*> ZHET01_3 reconstructs a Hermitian indefinite matrix A from its
33*> block L*D*L' or U*D*U' factorization computed by ZHETRF_RK
34*> (or ZHETRF_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*16 array, dimension (LDAFAC,N)
72*>          Diagonal of the block diagonal matrix D and factors U or L
73*>          as computed by ZHETRF_RK and ZHETRF_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*16 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 ZHETRF_RK (or ZHETRF_BK).
103*> \endverbatim
104*>
105*> \param[out] C
106*> \verbatim
107*>          C is COMPLEX*16 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 DOUBLE PRECISION array, dimension (N)
119*> \endverbatim
120*>
121*> \param[out] RESID
122*> \verbatim
123*>          RESID is DOUBLE PRECISION
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*> \ingroup complex16_lin
137*
138*  =====================================================================
139      SUBROUTINE ZHET01_3( UPLO, N, A, LDA, AFAC, LDAFAC, E, IPIV, C,
140     $                     LDC, RWORK, RESID )
141*
142*  -- LAPACK test routine --
143*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
144*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
145*
146*     .. Scalar Arguments ..
147      CHARACTER          UPLO
148      INTEGER            LDA, LDAFAC, LDC, N
149      DOUBLE PRECISION   RESID
150*     ..
151*     .. Array Arguments ..
152      INTEGER            IPIV( * )
153      DOUBLE PRECISION   RWORK( * )
154      COMPLEX*16         A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
155     $                   E( * )
156*     ..
157*
158*  =====================================================================
159*
160*     .. Parameters ..
161      DOUBLE PRECISION   ZERO, ONE
162      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
163      COMPLEX*16         CZERO, CONE
164      PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
165     $                   CONE = ( 1.0D+0, 0.0D+0 ) )
166*     ..
167*     .. Local Scalars ..
168      INTEGER            I, INFO, J
169      DOUBLE PRECISION   ANORM, EPS
170*     ..
171*     .. External Functions ..
172      LOGICAL            LSAME
173      DOUBLE PRECISION   ZLANHE, DLAMCH
174      EXTERNAL           LSAME, ZLANHE, DLAMCH
175*     ..
176*     .. External Subroutines ..
177      EXTERNAL           ZLASET, ZLAVHE_ROOK, ZSYCONVF_ROOK
178*     ..
179*     .. Intrinsic Functions ..
180      INTRINSIC          DIMAG, DBLE
181*     ..
182*     .. Executable Statements ..
183*
184*     Quick exit if N = 0.
185*
186      IF( N.LE.0 ) THEN
187         RESID = ZERO
188         RETURN
189      END IF
190*
191*     a) Revert to multiplyers of L
192*
193      CALL ZSYCONVF_ROOK( UPLO, 'R', N, AFAC, LDAFAC, E, IPIV, INFO )
194*
195*     1) Determine EPS and the norm of A.
196*
197      EPS = DLAMCH( 'Epsilon' )
198      ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
199*
200*     Check the imaginary parts of the diagonal elements and return with
201*     an error code if any are nonzero.
202*
203      DO J = 1, N
204         IF( DIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
205            RESID = ONE / EPS
206            RETURN
207         END IF
208      END DO
209*
210*     2) Initialize C to the identity matrix.
211*
212      CALL ZLASET( 'Full', N, N, CZERO, CONE, C, LDC )
213*
214*     3) Call ZLAVHE_ROOK to form the product D * U' (or D * L' ).
215*
216      CALL ZLAVHE_ROOK( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC,
217     $                  LDAFAC, IPIV, C, LDC, INFO )
218*
219*     4) Call ZLAVHE_RK again to multiply by U (or L ).
220*
221      CALL ZLAVHE_ROOK( UPLO, 'No transpose', 'Unit', N, N, AFAC,
222     $                  LDAFAC, IPIV, C, LDC, INFO )
223*
224*     5) Compute the difference  C - A .
225*
226      IF( LSAME( UPLO, 'U' ) ) THEN
227         DO J = 1, N
228            DO I = 1, J - 1
229               C( I, J ) = C( I, J ) - A( I, J )
230            END DO
231            C( J, J ) = C( J, J ) - DBLE( A( J, J ) )
232         END DO
233      ELSE
234         DO J = 1, N
235            C( J, J ) = C( J, J ) - DBLE( A( J, J ) )
236            DO I = J + 1, N
237               C( I, J ) = C( I, J ) - A( I, J )
238            END DO
239         END DO
240      END IF
241*
242*     6) Compute norm( C - A ) / ( N * norm(A) * EPS )
243*
244      RESID = ZLANHE( '1', UPLO, N, C, LDC, RWORK )
245*
246      IF( ANORM.LE.ZERO ) THEN
247         IF( RESID.NE.ZERO )
248     $      RESID = ONE / EPS
249      ELSE
250         RESID = ( ( RESID/DBLE( N ) )/ANORM ) / EPS
251      END IF
252*
253*     b) Convert to factor of L (or U)
254*
255      CALL ZSYCONVF_ROOK( UPLO, 'C', N, AFAC, LDAFAC, E, IPIV, INFO )
256*
257      RETURN
258*
259*     End of ZHET01_3
260*
261      END
262