LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DORT01( ROWCOL, M, N, U, LDU, WORK, LWORK, RESID ) 00002 * 00003 * -- LAPACK test routine (version 3.1) -- 00004 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00005 * November 2006 00006 * 00007 * .. Scalar Arguments .. 00008 CHARACTER ROWCOL 00009 INTEGER LDU, LWORK, M, N 00010 DOUBLE PRECISION RESID 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION U( LDU, * ), WORK( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * DORT01 checks that the matrix U is orthogonal by computing the ratio 00020 * 00021 * RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R', 00022 * or 00023 * RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'. 00024 * 00025 * Alternatively, if there isn't sufficient workspace to form 00026 * I - U*U' or I - U'*U, the ratio is computed as 00027 * 00028 * RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R', 00029 * or 00030 * RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'. 00031 * 00032 * where EPS is the machine precision. ROWCOL is used only if m = n; 00033 * if m > n, ROWCOL is assumed to be 'C', and if m < n, ROWCOL is 00034 * assumed to be 'R'. 00035 * 00036 * Arguments 00037 * ========= 00038 * 00039 * ROWCOL (input) CHARACTER 00040 * Specifies whether the rows or columns of U should be checked 00041 * for orthogonality. Used only if M = N. 00042 * = 'R': Check for orthogonal rows of U 00043 * = 'C': Check for orthogonal columns of U 00044 * 00045 * M (input) INTEGER 00046 * The number of rows of the matrix U. 00047 * 00048 * N (input) INTEGER 00049 * The number of columns of the matrix U. 00050 * 00051 * U (input) DOUBLE PRECISION array, dimension (LDU,N) 00052 * The orthogonal matrix U. U is checked for orthogonal columns 00053 * if m > n or if m = n and ROWCOL = 'C'. U is checked for 00054 * orthogonal rows if m < n or if m = n and ROWCOL = 'R'. 00055 * 00056 * LDU (input) INTEGER 00057 * The leading dimension of the array U. LDU >= max(1,M). 00058 * 00059 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 00060 * 00061 * LWORK (input) INTEGER 00062 * The length of the array WORK. For best performance, LWORK 00063 * should be at least N*(N+1) if ROWCOL = 'C' or M*(M+1) if 00064 * ROWCOL = 'R', but the test will be done even if LWORK is 0. 00065 * 00066 * RESID (output) DOUBLE PRECISION 00067 * RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or 00068 * RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'. 00069 * 00070 * ===================================================================== 00071 * 00072 * .. Parameters .. 00073 DOUBLE PRECISION ZERO, ONE 00074 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00075 * .. 00076 * .. Local Scalars .. 00077 CHARACTER TRANSU 00078 INTEGER I, J, K, LDWORK, MNMIN 00079 DOUBLE PRECISION EPS, TMP 00080 * .. 00081 * .. External Functions .. 00082 LOGICAL LSAME 00083 DOUBLE PRECISION DDOT, DLAMCH, DLANSY 00084 EXTERNAL LSAME, DDOT, DLAMCH, DLANSY 00085 * .. 00086 * .. External Subroutines .. 00087 EXTERNAL DLASET, DSYRK 00088 * .. 00089 * .. Intrinsic Functions .. 00090 INTRINSIC ABS, DBLE, MAX, MIN 00091 * .. 00092 * .. Executable Statements .. 00093 * 00094 RESID = ZERO 00095 * 00096 * Quick return if possible 00097 * 00098 IF( M.LE.0 .OR. N.LE.0 ) 00099 $ RETURN 00100 * 00101 EPS = DLAMCH( 'Precision' ) 00102 IF( M.LT.N .OR. ( M.EQ.N .AND. LSAME( ROWCOL, 'R' ) ) ) THEN 00103 TRANSU = 'N' 00104 K = N 00105 ELSE 00106 TRANSU = 'T' 00107 K = M 00108 END IF 00109 MNMIN = MIN( M, N ) 00110 * 00111 IF( ( MNMIN+1 )*MNMIN.LE.LWORK ) THEN 00112 LDWORK = MNMIN 00113 ELSE 00114 LDWORK = 0 00115 END IF 00116 IF( LDWORK.GT.0 ) THEN 00117 * 00118 * Compute I - U*U' or I - U'*U. 00119 * 00120 CALL DLASET( 'Upper', MNMIN, MNMIN, ZERO, ONE, WORK, LDWORK ) 00121 CALL DSYRK( 'Upper', TRANSU, MNMIN, K, -ONE, U, LDU, ONE, WORK, 00122 $ LDWORK ) 00123 * 00124 * Compute norm( I - U*U' ) / ( K * EPS ) . 00125 * 00126 RESID = DLANSY( '1', 'Upper', MNMIN, WORK, LDWORK, 00127 $ WORK( LDWORK*MNMIN+1 ) ) 00128 RESID = ( RESID / DBLE( K ) ) / EPS 00129 ELSE IF( TRANSU.EQ.'T' ) THEN 00130 * 00131 * Find the maximum element in abs( I - U'*U ) / ( m * EPS ) 00132 * 00133 DO 20 J = 1, N 00134 DO 10 I = 1, J 00135 IF( I.NE.J ) THEN 00136 TMP = ZERO 00137 ELSE 00138 TMP = ONE 00139 END IF 00140 TMP = TMP - DDOT( M, U( 1, I ), 1, U( 1, J ), 1 ) 00141 RESID = MAX( RESID, ABS( TMP ) ) 00142 10 CONTINUE 00143 20 CONTINUE 00144 RESID = ( RESID / DBLE( M ) ) / EPS 00145 ELSE 00146 * 00147 * Find the maximum element in abs( I - U*U' ) / ( n * EPS ) 00148 * 00149 DO 40 J = 1, M 00150 DO 30 I = 1, J 00151 IF( I.NE.J ) THEN 00152 TMP = ZERO 00153 ELSE 00154 TMP = ONE 00155 END IF 00156 TMP = TMP - DDOT( N, U( J, 1 ), LDU, U( I, 1 ), LDU ) 00157 RESID = MAX( RESID, ABS( TMP ) ) 00158 30 CONTINUE 00159 40 CONTINUE 00160 RESID = ( RESID / DBLE( N ) ) / EPS 00161 END IF 00162 RETURN 00163 * 00164 * End of DORT01 00165 * 00166 END