3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
11 * SUBROUTINE ZSYT01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC,
14 * .. Scalar Arguments ..
16 * INTEGER LDA, LDAFAC, LDC, N
17 * DOUBLE PRECISION RESID
19 * .. Array Arguments ..
21 * DOUBLE PRECISION RWORK( * )
22 * COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
31 *> ZSYT01 reconstructs a hermitian indefinite matrix A from its
32 *> block L*D*L' or U*D*U' factorization and computes the residual
33 *> norm( C - A ) / ( N * norm(A) * EPS ),
34 *> where C is the reconstructed matrix and EPS is the machine epsilon.
42 *> UPLO is CHARACTER*1
43 *> Specifies whether the upper or lower triangular part of the
44 *> hermitian matrix A is stored:
45 *> = 'U': Upper triangular
46 *> = 'L': Lower triangular
52 *> The number of rows and columns of the matrix A. N >= 0.
57 *> A is COMPLEX*16 array, dimension (LDA,N)
58 *> The original hermitian matrix A.
64 *> The leading dimension of the array A. LDA >= max(1,N)
69 *> AFAC is COMPLEX*16 array, dimension (LDAFAC,N)
70 *> The factored form of the matrix A. AFAC contains the block
71 *> diagonal matrix D and the multipliers used to obtain the
72 *> factor L or U from the block L*D*L' or U*D*U' factorization
73 *> as computed by ZSYTRF.
79 *> The leading dimension of the array AFAC. LDAFAC >= max(1,N).
84 *> IPIV is INTEGER array, dimension (N)
85 *> The pivot indices from ZSYTRF.
90 *> C is COMPLEX*16 array, dimension (LDC,N)
96 *> The leading dimension of the array C. LDC >= max(1,N).
101 *> RWORK is COMPLEX*16 array, dimension (N)
106 *> RESID is COMPLEX*16
107 *> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
108 *> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
114 *> \author Univ. of Tennessee
115 *> \author Univ. of California Berkeley
116 *> \author Univ. of Colorado Denver
119 *> \date November 2016
121 * @generated from LIN/dsyt01_aa.f, fortran d -> z, Thu Nov 17 13:01:50 2016
123 *> \ingroup complex16_lin
125 * =====================================================================
126 SUBROUTINE ZSYT01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C,
127 $ LDC, RWORK, RESID )
129 * -- LAPACK test routine (version 3.5.0) --
130 * -- LAPACK is a software package provided by Univ. of Tennessee, --
131 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
134 * .. Scalar Arguments ..
136 INTEGER LDA, LDAFAC, LDC, N
137 DOUBLE PRECISION RESID
139 * .. Array Arguments ..
141 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
142 DOUBLE PRECISION RWORK( * )
145 * =====================================================================
148 DOUBLE PRECISION ZERO, ONE
149 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
150 COMPLEX*16 CZERO, CONE
151 PARAMETER ( CZERO = 0.0E+0, CONE = 1.0E+0 )
153 * .. Local Scalars ..
155 DOUBLE PRECISION ANORM, EPS
157 * .. External Functions ..
159 DOUBLE PRECISION DLAMCH, ZLANSY
160 EXTERNAL LSAME, DLAMCH, ZLANSY
162 * .. External Subroutines ..
163 EXTERNAL ZLASET, ZLAVSY
165 * .. Intrinsic Functions ..
168 * .. Executable Statements ..
170 * Quick exit if N = 0.
177 * Determine EPS and the norm of A.
179 EPS = DLAMCH( 'Epsilon' )
180 ANORM = ZLANSY( '1', UPLO, N, A, LDA, RWORK )
182 * Initialize C to the tridiagonal matrix T.
184 CALL ZLASET( 'Full', N, N, CZERO, CZERO, C, LDC )
185 CALL ZLACPY( 'F', 1, N, AFAC( 1, 1 ), LDAFAC+1, C( 1, 1 ), LDC+1 )
187 IF( LSAME( UPLO, 'U' ) ) THEN
188 CALL ZLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 1, 2 ),
190 CALL ZLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 2, 1 ),
193 CALL ZLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 1, 2 ),
195 CALL ZLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ),
199 * Call ZTRMM to form the product U' * D (or L * D ).
201 IF( LSAME( UPLO, 'U' ) ) THEN
202 CALL ZTRMM( 'Left', UPLO, 'Transpose', 'Unit', N-1, N,
203 $ CONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC )
205 CALL ZTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N,
206 $ CONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC )
209 * Call ZTRMM again to multiply by U (or L ).
211 IF( LSAME( UPLO, 'U' ) ) THEN
212 CALL ZTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1,
213 $ CONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC )
215 CALL ZTRMM( 'Right', UPLO, 'Transpose', 'Unit', N, N-1,
216 $ CONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC )
220 * Apply symmetric pivots
225 $ CALL ZSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC )
230 $ CALL ZSWAP( N, C( 1, J ), 1, C( 1, I ), 1 )
234 * Compute the difference C - A .
236 IF( LSAME( UPLO, 'U' ) ) THEN
239 C( I, J ) = C( I, J ) - A( I, J )
245 C( I, J ) = C( I, J ) - A( I, J )
250 * Compute norm( C - A ) / ( N * norm(A) * EPS )
252 RESID = ZLANSY( '1', UPLO, N, C, LDC, RWORK )
254 IF( ANORM.LE.ZERO ) THEN
258 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS