1*> \brief \b ZSYT01_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 ZSYT01_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*> ZSYT01_3 reconstructs a symmetric indefinite matrix A from its
33*> block L*D*L' or U*D*U' factorization computed by ZSYTRF_RK
34*> (or ZSYTRF_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*>          symmetric 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 symmetric 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 ZSYTRF_RK and ZSYTRF_BK:
74*>            a) ONLY diagonal elements of the symmetric 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 symmetric 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 ZSYTRF_RK (or ZSYTRF_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*> \date June 2017
137*
138*> \ingroup complex16_lin
139*
140*  =====================================================================
141      SUBROUTINE ZSYT01_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      DOUBLE PRECISION   RESID
153*     ..
154*     .. Array Arguments ..
155      INTEGER            IPIV( * )
156      DOUBLE PRECISION   RWORK( * )
157      COMPLEX*16         A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
158     $                   E( * )
159*     ..
160*
161*  =====================================================================
162*
163*     .. Parameters ..
164      DOUBLE PRECISION   ZERO, ONE
165      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
166      COMPLEX*16         CZERO, CONE
167      PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
168     $                   CONE = ( 1.0D+0, 0.0D+0 ) )
169*     ..
170*     .. Local Scalars ..
171      INTEGER            I, INFO, J
172      DOUBLE PRECISION   ANORM, EPS
173*     ..
174*     .. External Functions ..
175      LOGICAL            LSAME
176      DOUBLE PRECISION   DLAMCH, ZLANSY
177      EXTERNAL           LSAME, DLAMCH, ZLANSY
178*     ..
179*     .. External Subroutines ..
180      EXTERNAL           ZLASET, ZLAVSY_ROOK, ZSYCONVF_ROOK
181*     ..
182*     .. Intrinsic Functions ..
183      INTRINSIC          DBLE
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 ZSYCONVF_ROOK( UPLO, 'R', N, AFAC, LDAFAC, E, IPIV, INFO )
197*
198*     1) Determine EPS and the norm of A.
199*
200      EPS = DLAMCH( 'Epsilon' )
201      ANORM = ZLANSY( '1', UPLO, N, A, LDA, RWORK )
202*
203*     2) Initialize C to the identity matrix.
204*
205      CALL ZLASET( 'Full', N, N, CZERO, CONE, C, LDC )
206*
207*     3) Call ZLAVSY_ROOK to form the product D * U' (or D * L' ).
208*
209      CALL ZLAVSY_ROOK( UPLO, 'Transpose', 'Non-unit', N, N, AFAC,
210     $                  LDAFAC, IPIV, C, LDC, INFO )
211*
212*     4) Call ZLAVSY_ROOK again to multiply by U (or L ).
213*
214      CALL ZLAVSY_ROOK( UPLO, 'No transpose', 'Unit', N, N, AFAC,
215     $                  LDAFAC, IPIV, C, LDC, INFO )
216*
217*     5) Compute the difference  C - A .
218*
219      IF( LSAME( UPLO, 'U' ) ) THEN
220         DO J = 1, N
221            DO I = 1, J
222               C( I, J ) = C( I, J ) - A( I, J )
223            END DO
224         END DO
225      ELSE
226         DO J = 1, N
227            DO I = J, N
228               C( I, J ) = C( I, J ) - A( I, J )
229            END DO
230         END DO
231      END IF
232*
233*     6) Compute norm( C - A ) / ( N * norm(A) * EPS )
234*
235      RESID = ZLANSY( '1', UPLO, N, C, LDC, RWORK )
236*
237      IF( ANORM.LE.ZERO ) THEN
238         IF( RESID.NE.ZERO )
239     $      RESID = ONE / EPS
240      ELSE
241         RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
242      END IF
243
244*
245*     b) Convert to factor of L (or U)
246*
247      CALL ZSYCONVF_ROOK( UPLO, 'C', N, AFAC, LDAFAC, E, IPIV, INFO )
248*
249      RETURN
250*
251*     End of ZSYT01_3
252*
253      END
254