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