![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZSYT03 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE ZSYT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK, 00012 * RWORK, RCOND, RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER UPLO 00016 * INTEGER LDA, LDAINV, LDWORK, N 00017 * DOUBLE PRECISION RCOND, RESID 00018 * .. 00019 * .. Array Arguments .. 00020 * DOUBLE PRECISION RWORK( * ) 00021 * COMPLEX*16 A( LDA, * ), AINV( LDAINV, * ), 00022 * $ WORK( LDWORK, * ) 00023 * .. 00024 * 00025 * 00026 *> \par Purpose: 00027 * ============= 00028 *> 00029 *> \verbatim 00030 *> 00031 *> ZSYT03 computes the residual for a complex symmetric matrix times 00032 *> its inverse: 00033 *> norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ) 00034 *> where EPS is the machine epsilon. 00035 *> \endverbatim 00036 * 00037 * Arguments: 00038 * ========== 00039 * 00040 *> \param[in] UPLO 00041 *> \verbatim 00042 *> UPLO is CHARACTER*1 00043 *> Specifies whether the upper or lower triangular part of the 00044 *> complex symmetric matrix A is stored: 00045 *> = 'U': Upper triangular 00046 *> = 'L': Lower triangular 00047 *> \endverbatim 00048 *> 00049 *> \param[in] N 00050 *> \verbatim 00051 *> N is INTEGER 00052 *> The number of rows and columns of the matrix A. N >= 0. 00053 *> \endverbatim 00054 *> 00055 *> \param[in] A 00056 *> \verbatim 00057 *> A is COMPLEX*16 array, dimension (LDA,N) 00058 *> The original complex symmetric matrix A. 00059 *> \endverbatim 00060 *> 00061 *> \param[in] LDA 00062 *> \verbatim 00063 *> LDA is INTEGER 00064 *> The leading dimension of the array A. LDA >= max(1,N) 00065 *> \endverbatim 00066 *> 00067 *> \param[in,out] AINV 00068 *> \verbatim 00069 *> AINV is COMPLEX*16 array, dimension (LDAINV,N) 00070 *> On entry, the inverse of the matrix A, stored as a symmetric 00071 *> matrix in the same format as A. 00072 *> In this version, AINV is expanded into a full matrix and 00073 *> multiplied by A, so the opposing triangle of AINV will be 00074 *> changed; i.e., if the upper triangular part of AINV is 00075 *> stored, the lower triangular part will be used as work space. 00076 *> \endverbatim 00077 *> 00078 *> \param[in] LDAINV 00079 *> \verbatim 00080 *> LDAINV is INTEGER 00081 *> The leading dimension of the array AINV. LDAINV >= max(1,N). 00082 *> \endverbatim 00083 *> 00084 *> \param[out] WORK 00085 *> \verbatim 00086 *> WORK is COMPLEX*16 array, dimension (LDWORK,N) 00087 *> \endverbatim 00088 *> 00089 *> \param[in] LDWORK 00090 *> \verbatim 00091 *> LDWORK is INTEGER 00092 *> The leading dimension of the array WORK. LDWORK >= max(1,N). 00093 *> \endverbatim 00094 *> 00095 *> \param[out] RWORK 00096 *> \verbatim 00097 *> RWORK is DOUBLE PRECISION array, dimension (N) 00098 *> \endverbatim 00099 *> 00100 *> \param[out] RCOND 00101 *> \verbatim 00102 *> RCOND is DOUBLE PRECISION 00103 *> The reciprocal of the condition number of A, computed as 00104 *> RCOND = 1/ (norm(A) * norm(AINV)). 00105 *> \endverbatim 00106 *> 00107 *> \param[out] RESID 00108 *> \verbatim 00109 *> RESID is DOUBLE PRECISION 00110 *> norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS ) 00111 *> \endverbatim 00112 * 00113 * Authors: 00114 * ======== 00115 * 00116 *> \author Univ. of Tennessee 00117 *> \author Univ. of California Berkeley 00118 *> \author Univ. of Colorado Denver 00119 *> \author NAG Ltd. 00120 * 00121 *> \date November 2011 00122 * 00123 *> \ingroup complex16_lin 00124 * 00125 * ===================================================================== 00126 SUBROUTINE ZSYT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK, 00127 $ RWORK, RCOND, RESID ) 00128 * 00129 * -- LAPACK test routine (version 3.4.0) -- 00130 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00131 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00132 * November 2011 00133 * 00134 * .. Scalar Arguments .. 00135 CHARACTER UPLO 00136 INTEGER LDA, LDAINV, LDWORK, N 00137 DOUBLE PRECISION RCOND, RESID 00138 * .. 00139 * .. Array Arguments .. 00140 DOUBLE PRECISION RWORK( * ) 00141 COMPLEX*16 A( LDA, * ), AINV( LDAINV, * ), 00142 $ WORK( LDWORK, * ) 00143 * .. 00144 * 00145 * ===================================================================== 00146 * 00147 * 00148 * .. Parameters .. 00149 DOUBLE PRECISION ZERO, ONE 00150 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00151 COMPLEX*16 CZERO, CONE 00152 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00153 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00154 * .. 00155 * .. Local Scalars .. 00156 INTEGER I, J 00157 DOUBLE PRECISION AINVNM, ANORM, EPS 00158 * .. 00159 * .. External Functions .. 00160 LOGICAL LSAME 00161 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY 00162 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANSY 00163 * .. 00164 * .. External Subroutines .. 00165 EXTERNAL ZSYMM 00166 * .. 00167 * .. Intrinsic Functions .. 00168 INTRINSIC DBLE 00169 * .. 00170 * .. Executable Statements .. 00171 * 00172 * Quick exit if N = 0 00173 * 00174 IF( N.LE.0 ) THEN 00175 RCOND = ONE 00176 RESID = ZERO 00177 RETURN 00178 END IF 00179 * 00180 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. 00181 * 00182 EPS = DLAMCH( 'Epsilon' ) 00183 ANORM = ZLANSY( '1', UPLO, N, A, LDA, RWORK ) 00184 AINVNM = ZLANSY( '1', UPLO, N, AINV, LDAINV, RWORK ) 00185 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00186 RCOND = ZERO 00187 RESID = ONE / EPS 00188 RETURN 00189 END IF 00190 RCOND = ( ONE / ANORM ) / AINVNM 00191 * 00192 * Expand AINV into a full matrix and call ZSYMM to multiply 00193 * AINV on the left by A (store the result in WORK). 00194 * 00195 IF( LSAME( UPLO, 'U' ) ) THEN 00196 DO 20 J = 1, N 00197 DO 10 I = 1, J - 1 00198 AINV( J, I ) = AINV( I, J ) 00199 10 CONTINUE 00200 20 CONTINUE 00201 ELSE 00202 DO 40 J = 1, N 00203 DO 30 I = J + 1, N 00204 AINV( J, I ) = AINV( I, J ) 00205 30 CONTINUE 00206 40 CONTINUE 00207 END IF 00208 CALL ZSYMM( 'Left', UPLO, N, N, -CONE, A, LDA, AINV, LDAINV, 00209 $ CZERO, WORK, LDWORK ) 00210 * 00211 * Add the identity matrix to WORK . 00212 * 00213 DO 50 I = 1, N 00214 WORK( I, I ) = WORK( I, I ) + CONE 00215 50 CONTINUE 00216 * 00217 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS) 00218 * 00219 RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK ) 00220 * 00221 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N ) 00222 * 00223 RETURN 00224 * 00225 * End of ZSYT03 00226 * 00227 END