![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CUNT01 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 CUNT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK, 00012 * RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER ROWCOL 00016 * INTEGER LDU, LWORK, M, N 00017 * REAL RESID 00018 * .. 00019 * .. Array Arguments .. 00020 * REAL RWORK( * ) 00021 * COMPLEX U( LDU, * ), WORK( * ) 00022 * .. 00023 * 00024 * 00025 *> \par Purpose: 00026 * ============= 00027 *> 00028 *> \verbatim 00029 *> 00030 *> CUNT01 checks that the matrix U is unitary by computing the ratio 00031 *> 00032 *> RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R', 00033 *> or 00034 *> RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'. 00035 *> 00036 *> Alternatively, if there isn't sufficient workspace to form 00037 *> I - U*U' or I - U'*U, the ratio is computed as 00038 *> 00039 *> RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R', 00040 *> or 00041 *> RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'. 00042 *> 00043 *> where EPS is the machine precision. ROWCOL is used only if m = n; 00044 *> if m > n, ROWCOL is assumed to be 'C', and if m < n, ROWCOL is 00045 *> assumed to be 'R'. 00046 *> \endverbatim 00047 * 00048 * Arguments: 00049 * ========== 00050 * 00051 *> \param[in] ROWCOL 00052 *> \verbatim 00053 *> ROWCOL is CHARACTER 00054 *> Specifies whether the rows or columns of U should be checked 00055 *> for orthogonality. Used only if M = N. 00056 *> = 'R': Check for orthogonal rows of U 00057 *> = 'C': Check for orthogonal columns of U 00058 *> \endverbatim 00059 *> 00060 *> \param[in] M 00061 *> \verbatim 00062 *> M is INTEGER 00063 *> The number of rows of the matrix U. 00064 *> \endverbatim 00065 *> 00066 *> \param[in] N 00067 *> \verbatim 00068 *> N is INTEGER 00069 *> The number of columns of the matrix U. 00070 *> \endverbatim 00071 *> 00072 *> \param[in] U 00073 *> \verbatim 00074 *> U is COMPLEX array, dimension (LDU,N) 00075 *> The unitary matrix U. U is checked for orthogonal columns 00076 *> if m > n or if m = n and ROWCOL = 'C'. U is checked for 00077 *> orthogonal rows if m < n or if m = n and ROWCOL = 'R'. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] LDU 00081 *> \verbatim 00082 *> LDU is INTEGER 00083 *> The leading dimension of the array U. LDU >= max(1,M). 00084 *> \endverbatim 00085 *> 00086 *> \param[out] WORK 00087 *> \verbatim 00088 *> WORK is COMPLEX array, dimension (LWORK) 00089 *> \endverbatim 00090 *> 00091 *> \param[in] LWORK 00092 *> \verbatim 00093 *> LWORK is INTEGER 00094 *> The length of the array WORK. For best performance, LWORK 00095 *> should be at least N*N if ROWCOL = 'C' or M*M if 00096 *> ROWCOL = 'R', but the test will be done even if LWORK is 0. 00097 *> \endverbatim 00098 *> 00099 *> \param[out] RWORK 00100 *> \verbatim 00101 *> RWORK is REAL array, dimension (min(M,N)) 00102 *> Used only if LWORK is large enough to use the Level 3 BLAS 00103 *> code. 00104 *> \endverbatim 00105 *> 00106 *> \param[out] RESID 00107 *> \verbatim 00108 *> RESID is REAL 00109 *> RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or 00110 *> RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'. 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 complex_eig 00124 * 00125 * ===================================================================== 00126 SUBROUTINE CUNT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK, 00127 $ 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 ROWCOL 00136 INTEGER LDU, LWORK, M, N 00137 REAL RESID 00138 * .. 00139 * .. Array Arguments .. 00140 REAL RWORK( * ) 00141 COMPLEX U( LDU, * ), WORK( * ) 00142 * .. 00143 * 00144 * ===================================================================== 00145 * 00146 * .. Parameters .. 00147 REAL ZERO, ONE 00148 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00149 * .. 00150 * .. Local Scalars .. 00151 CHARACTER TRANSU 00152 INTEGER I, J, K, LDWORK, MNMIN 00153 REAL EPS 00154 COMPLEX TMP, ZDUM 00155 * .. 00156 * .. External Functions .. 00157 LOGICAL LSAME 00158 REAL CLANSY, SLAMCH 00159 COMPLEX CDOTC 00160 EXTERNAL LSAME, CLANSY, SLAMCH, CDOTC 00161 * .. 00162 * .. External Subroutines .. 00163 EXTERNAL CHERK, CLASET 00164 * .. 00165 * .. Intrinsic Functions .. 00166 INTRINSIC ABS, AIMAG, CMPLX, MAX, MIN, REAL 00167 * .. 00168 * .. Statement Functions .. 00169 REAL CABS1 00170 * .. 00171 * .. Statement Function definitions .. 00172 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) 00173 * .. 00174 * .. Executable Statements .. 00175 * 00176 RESID = ZERO 00177 * 00178 * Quick return if possible 00179 * 00180 IF( M.LE.0 .OR. N.LE.0 ) 00181 $ RETURN 00182 * 00183 EPS = SLAMCH( 'Precision' ) 00184 IF( M.LT.N .OR. ( M.EQ.N .AND. LSAME( ROWCOL, 'R' ) ) ) THEN 00185 TRANSU = 'N' 00186 K = N 00187 ELSE 00188 TRANSU = 'C' 00189 K = M 00190 END IF 00191 MNMIN = MIN( M, N ) 00192 * 00193 IF( ( MNMIN+1 )*MNMIN.LE.LWORK ) THEN 00194 LDWORK = MNMIN 00195 ELSE 00196 LDWORK = 0 00197 END IF 00198 IF( LDWORK.GT.0 ) THEN 00199 * 00200 * Compute I - U*U' or I - U'*U. 00201 * 00202 CALL CLASET( 'Upper', MNMIN, MNMIN, CMPLX( ZERO ), 00203 $ CMPLX( ONE ), WORK, LDWORK ) 00204 CALL CHERK( 'Upper', TRANSU, MNMIN, K, -ONE, U, LDU, ONE, WORK, 00205 $ LDWORK ) 00206 * 00207 * Compute norm( I - U*U' ) / ( K * EPS ) . 00208 * 00209 RESID = CLANSY( '1', 'Upper', MNMIN, WORK, LDWORK, RWORK ) 00210 RESID = ( RESID / REAL( K ) ) / EPS 00211 ELSE IF( TRANSU.EQ.'C' ) THEN 00212 * 00213 * Find the maximum element in abs( I - U'*U ) / ( m * EPS ) 00214 * 00215 DO 20 J = 1, N 00216 DO 10 I = 1, J 00217 IF( I.NE.J ) THEN 00218 TMP = ZERO 00219 ELSE 00220 TMP = ONE 00221 END IF 00222 TMP = TMP - CDOTC( M, U( 1, I ), 1, U( 1, J ), 1 ) 00223 RESID = MAX( RESID, CABS1( TMP ) ) 00224 10 CONTINUE 00225 20 CONTINUE 00226 RESID = ( RESID / REAL( M ) ) / EPS 00227 ELSE 00228 * 00229 * Find the maximum element in abs( I - U*U' ) / ( n * EPS ) 00230 * 00231 DO 40 J = 1, M 00232 DO 30 I = 1, J 00233 IF( I.NE.J ) THEN 00234 TMP = ZERO 00235 ELSE 00236 TMP = ONE 00237 END IF 00238 TMP = TMP - CDOTC( N, U( J, 1 ), LDU, U( I, 1 ), LDU ) 00239 RESID = MAX( RESID, CABS1( TMP ) ) 00240 30 CONTINUE 00241 40 CONTINUE 00242 RESID = ( RESID / REAL( N ) ) / EPS 00243 END IF 00244 RETURN 00245 * 00246 * End of CUNT01 00247 * 00248 END