00001 SUBROUTINE SORT03( RC, MU, MV, N, K, U, LDU, V, LDV, WORK, LWORK, 00002 $ RESULT, INFO ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER*( * ) RC 00010 INTEGER INFO, K, LDU, LDV, LWORK, MU, MV, N 00011 REAL RESULT 00012 * .. 00013 * .. Array Arguments .. 00014 REAL U( LDU, * ), V( LDV, * ), WORK( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * SORT03 compares two orthogonal matrices U and V to see if their 00021 * corresponding rows or columns span the same spaces. The rows are 00022 * checked if RC = 'R', and the columns are checked if RC = 'C'. 00023 * 00024 * RESULT is the maximum of 00025 * 00026 * | V*V' - I | / ( MV ulp ), if RC = 'R', or 00027 * 00028 * | V'*V - I | / ( MV ulp ), if RC = 'C', 00029 * 00030 * and the maximum over rows (or columns) 1 to K of 00031 * 00032 * | U(i) - S*V(i) |/ ( N ulp ) 00033 * 00034 * where S is +-1 (chosen to minimize the expression), U(i) is the i-th 00035 * row (column) of U, and V(i) is the i-th row (column) of V. 00036 * 00037 * Arguments 00038 * ========== 00039 * 00040 * RC (input) CHARACTER*1 00041 * If RC = 'R' the rows of U and V are to be compared. 00042 * If RC = 'C' the columns of U and V are to be compared. 00043 * 00044 * MU (input) INTEGER 00045 * The number of rows of U if RC = 'R', and the number of 00046 * columns if RC = 'C'. If MU = 0 SORT03 does nothing. 00047 * MU must be at least zero. 00048 * 00049 * MV (input) INTEGER 00050 * The number of rows of V if RC = 'R', and the number of 00051 * columns if RC = 'C'. If MV = 0 SORT03 does nothing. 00052 * MV must be at least zero. 00053 * 00054 * N (input) INTEGER 00055 * If RC = 'R', the number of columns in the matrices U and V, 00056 * and if RC = 'C', the number of rows in U and V. If N = 0 00057 * SORT03 does nothing. N must be at least zero. 00058 * 00059 * K (input) INTEGER 00060 * The number of rows or columns of U and V to compare. 00061 * 0 <= K <= max(MU,MV). 00062 * 00063 * U (input) REAL array, dimension (LDU,N) 00064 * The first matrix to compare. If RC = 'R', U is MU by N, and 00065 * if RC = 'C', U is N by MU. 00066 * 00067 * LDU (input) INTEGER 00068 * The leading dimension of U. If RC = 'R', LDU >= max(1,MU), 00069 * and if RC = 'C', LDU >= max(1,N). 00070 * 00071 * V (input) REAL array, dimension (LDV,N) 00072 * The second matrix to compare. If RC = 'R', V is MV by N, and 00073 * if RC = 'C', V is N by MV. 00074 * 00075 * LDV (input) INTEGER 00076 * The leading dimension of V. If RC = 'R', LDV >= max(1,MV), 00077 * and if RC = 'C', LDV >= max(1,N). 00078 * 00079 * WORK (workspace) REAL array, dimension (LWORK) 00080 * 00081 * LWORK (input) INTEGER 00082 * The length of the array WORK. For best performance, LWORK 00083 * should be at least N*N if RC = 'C' or M*M if RC = 'R', but 00084 * the tests will be done even if LWORK is 0. 00085 * 00086 * RESULT (output) REAL 00087 * The value computed by the test described above. RESULT is 00088 * limited to 1/ulp to avoid overflow. 00089 * 00090 * INFO (output) INTEGER 00091 * 0 indicates a successful exit 00092 * -k indicates the k-th parameter had an illegal value 00093 * 00094 * ===================================================================== 00095 * 00096 * .. Parameters .. 00097 REAL ZERO, ONE 00098 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00099 * .. 00100 * .. Local Scalars .. 00101 INTEGER I, IRC, J, LMX 00102 REAL RES1, RES2, S, ULP 00103 * .. 00104 * .. External Functions .. 00105 LOGICAL LSAME 00106 INTEGER ISAMAX 00107 REAL SLAMCH 00108 EXTERNAL LSAME, ISAMAX, SLAMCH 00109 * .. 00110 * .. Intrinsic Functions .. 00111 INTRINSIC ABS, MAX, MIN, REAL, SIGN 00112 * .. 00113 * .. External Subroutines .. 00114 EXTERNAL SORT01, XERBLA 00115 * .. 00116 * .. Executable Statements .. 00117 * 00118 * Check inputs 00119 * 00120 INFO = 0 00121 IF( LSAME( RC, 'R' ) ) THEN 00122 IRC = 0 00123 ELSE IF( LSAME( RC, 'C' ) ) THEN 00124 IRC = 1 00125 ELSE 00126 IRC = -1 00127 END IF 00128 IF( IRC.EQ.-1 ) THEN 00129 INFO = -1 00130 ELSE IF( MU.LT.0 ) THEN 00131 INFO = -2 00132 ELSE IF( MV.LT.0 ) THEN 00133 INFO = -3 00134 ELSE IF( N.LT.0 ) THEN 00135 INFO = -4 00136 ELSE IF( K.LT.0 .OR. K.GT.MAX( MU, MV ) ) THEN 00137 INFO = -5 00138 ELSE IF( ( IRC.EQ.0 .AND. LDU.LT.MAX( 1, MU ) ) .OR. 00139 $ ( IRC.EQ.1 .AND. LDU.LT.MAX( 1, N ) ) ) THEN 00140 INFO = -7 00141 ELSE IF( ( IRC.EQ.0 .AND. LDV.LT.MAX( 1, MV ) ) .OR. 00142 $ ( IRC.EQ.1 .AND. LDV.LT.MAX( 1, N ) ) ) THEN 00143 INFO = -9 00144 END IF 00145 IF( INFO.NE.0 ) THEN 00146 CALL XERBLA( 'SORT03', -INFO ) 00147 RETURN 00148 END IF 00149 * 00150 * Initialize result 00151 * 00152 RESULT = ZERO 00153 IF( MU.EQ.0 .OR. MV.EQ.0 .OR. N.EQ.0 ) 00154 $ RETURN 00155 * 00156 * Machine constants 00157 * 00158 ULP = SLAMCH( 'Precision' ) 00159 * 00160 IF( IRC.EQ.0 ) THEN 00161 * 00162 * Compare rows 00163 * 00164 RES1 = ZERO 00165 DO 20 I = 1, K 00166 LMX = ISAMAX( N, U( I, 1 ), LDU ) 00167 S = SIGN( ONE, U( I, LMX ) )*SIGN( ONE, V( I, LMX ) ) 00168 DO 10 J = 1, N 00169 RES1 = MAX( RES1, ABS( U( I, J )-S*V( I, J ) ) ) 00170 10 CONTINUE 00171 20 CONTINUE 00172 RES1 = RES1 / ( REAL( N )*ULP ) 00173 * 00174 * Compute orthogonality of rows of V. 00175 * 00176 CALL SORT01( 'Rows', MV, N, V, LDV, WORK, LWORK, RES2 ) 00177 * 00178 ELSE 00179 * 00180 * Compare columns 00181 * 00182 RES1 = ZERO 00183 DO 40 I = 1, K 00184 LMX = ISAMAX( N, U( 1, I ), 1 ) 00185 S = SIGN( ONE, U( LMX, I ) )*SIGN( ONE, V( LMX, I ) ) 00186 DO 30 J = 1, N 00187 RES1 = MAX( RES1, ABS( U( J, I )-S*V( J, I ) ) ) 00188 30 CONTINUE 00189 40 CONTINUE 00190 RES1 = RES1 / ( REAL( N )*ULP ) 00191 * 00192 * Compute orthogonality of columns of V. 00193 * 00194 CALL SORT01( 'Columns', N, MV, V, LDV, WORK, LWORK, RES2 ) 00195 END IF 00196 * 00197 RESULT = MIN( MAX( RES1, RES2 ), ONE / ULP ) 00198 RETURN 00199 * 00200 * End of SORT03 00201 * 00202 END
1.7.2