 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DCKCSD( NM, MVAL, PVAL, QVAL, NMATS, ISEED, THRESH, 00002 $ MMAX, X, XF, U1, U2, V1T, V2T, THETA, IWORK, 00003 $ WORK, RWORK, NIN, NOUT, INFO ) 00004 IMPLICIT NONE 00005 * 00006 * Originally DCKGSV 00007 * -- LAPACK test routine (version 3.1) -- 00008 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00009 * November 2006 00010 * 00011 * Adapted to DCKCSD 00012 * July 2010 00013 * 00014 * .. Scalar Arguments .. 00015 INTEGER INFO, NIN, NM, NMATS, MMAX, NOUT 00016 DOUBLE PRECISION THRESH 00017 * .. 00018 * .. Array Arguments .. 00019 INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), PVAL( * ), 00020 $ QVAL( * ) 00021 DOUBLE PRECISION RWORK( * ), THETA( * ) 00022 DOUBLE PRECISION U1( * ), U2( * ), V1T( * ), V2T( * ), 00023 $ WORK( * ), X( * ), XF( * ) 00024 * .. 00025 * 00026 * Purpose 00027 * ======= 00028 * 00029 * DCKCSD tests DORCSD: 00030 * the CSD for an M-by-M orthogonal matrix X partitioned as 00031 * [ X11 X12; X21 X22 ]. X11 is P-by-Q. 00032 * 00033 * Arguments 00034 * ========= 00035 * 00036 * NM (input) INTEGER 00037 * The number of values of M contained in the vector MVAL. 00038 * 00039 * MVAL (input) INTEGER array, dimension (NM) 00040 * The values of the matrix row dimension M. 00041 * 00042 * PVAL (input) INTEGER array, dimension (NM) 00043 * The values of the matrix row dimension P. 00044 * 00045 * QVAL (input) INTEGER array, dimension (NM) 00046 * The values of the matrix column dimension Q. 00047 * 00048 * NMATS (input) INTEGER 00049 * The number of matrix types to be tested for each combination 00050 * of matrix dimensions. If NMATS >= NTYPES (the maximum 00051 * number of matrix types), then all the different types are 00052 * generated for testing. If NMATS < NTYPES, another input line 00053 * is read to get the numbers of the matrix types to be used. 00054 * 00055 * ISEED (input/output) INTEGER array, dimension (4) 00056 * On entry, the seed of the random number generator. The array 00057 * elements should be between 0 and 4095, otherwise they will be 00058 * reduced mod 4096, and ISEED(4) must be odd. 00059 * On exit, the next seed in the random number sequence after 00060 * all the test matrices have been generated. 00061 * 00062 * THRESH (input) DOUBLE PRECISION 00063 * The threshold value for the test ratios. A result is 00064 * included in the output file if RESULT >= THRESH. To have 00065 * every test ratio printed, use THRESH = 0. 00066 * 00067 * MMAX (input) INTEGER 00068 * The maximum value permitted for M, used in dimensioning the 00069 * work arrays. 00070 * 00071 * X (workspace) DOUBLE PRECISION array, dimension (MMAX*MMAX) 00072 * 00073 * XF (workspace) DOUBLE PRECISION array, dimension (MMAX*MMAX) 00074 * 00075 * U1 (workspace) DOUBLE PRECISION array, dimension (MMAX*MMAX) 00076 * 00077 * U2 (workspace) DOUBLE PRECISION array, dimension (MMAX*MMAX) 00078 * 00079 * V1T (workspace) DOUBLE PRECISION array, dimension (MMAX*MMAX) 00080 * 00081 * V2T (workspace) DOUBLE PRECISION array, dimension (MMAX*MMAX) 00082 * 00083 * THETA (workspace) DOUBLE PRECISION array, dimension (MMAX) 00084 * 00085 * IWORK (workspace) INTEGER array, dimension (MMAX) 00086 * 00087 * WORK (workspace) DOUBLE PRECISION array 00088 * 00089 * RWORK (workspace) DOUBLE PRECISION array 00090 * 00091 * NIN (input) INTEGER 00092 * The unit number for input. 00093 * 00094 * NOUT (input) INTEGER 00095 * The unit number for output. 00096 * 00097 * INFO (output) INTEGER 00098 * = 0 : successful exit 00099 * > 0 : If DLAROR returns an error code, the absolute value 00100 * of it is returned. 00101 * 00102 * ===================================================================== 00103 * 00104 * .. Parameters .. 00105 INTEGER NTESTS 00106 PARAMETER ( NTESTS = 9 ) 00107 INTEGER NTYPES 00108 PARAMETER ( NTYPES = 3 ) 00109 DOUBLE PRECISION GAPDIGIT, ORTH, PIOVER2, TEN 00110 PARAMETER ( GAPDIGIT = 18.0D0, ORTH = 1.0D-12, 00111 $ PIOVER2 = 1.57079632679489662D0, 00112 $ TEN = 10.0D0 ) 00113 DOUBLE PRECISION ZERO, ONE 00114 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00115 * .. 00116 * .. Local Scalars .. 00117 LOGICAL FIRSTT 00118 CHARACTER*3 PATH 00119 INTEGER I, IINFO, IM, IMAT, J, LDU1, LDU2, LDV1T, 00120 $ LDV2T, LDX, LWORK, M, NFAIL, NRUN, NT, P, Q, R 00121 * .. 00122 * .. Local Arrays .. 00123 LOGICAL DOTYPE( NTYPES ) 00124 DOUBLE PRECISION RESULT( NTESTS ) 00125 * .. 00126 * .. External Subroutines .. 00127 EXTERNAL ALAHDG, ALAREQ, ALASUM, DCSDTS, DLACSG, DLAROR, 00128 $ DLASET 00129 * .. 00130 * .. Intrinsic Functions .. 00131 INTRINSIC ABS, COS, MIN, SIN 00132 * .. 00133 * .. External Functions .. 00134 DOUBLE PRECISION DLANGE, DLARND 00135 EXTERNAL DLANGE, DLARND 00136 * .. 00137 * .. Executable Statements .. 00138 * 00139 * Initialize constants and the random number seed. 00140 * 00141 PATH( 1: 3 ) = 'CSD' 00142 INFO = 0 00143 NRUN = 0 00144 NFAIL = 0 00145 FIRSTT = .TRUE. 00146 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) 00147 LDX = MMAX 00148 LDU1 = MMAX 00149 LDU2 = MMAX 00150 LDV1T = MMAX 00151 LDV2T = MMAX 00152 LWORK = MMAX*MMAX 00153 * 00154 * Do for each value of M in MVAL. 00155 * 00156 DO 30 IM = 1, NM 00157 M = MVAL( IM ) 00158 P = PVAL( IM ) 00159 Q = QVAL( IM ) 00160 * 00161 DO 20 IMAT = 1, NTYPES 00162 * 00163 * Do the tests only if DOTYPE( IMAT ) is true. 00164 * 00165 IF( .NOT.DOTYPE( IMAT ) ) 00166 $ GO TO 20 00167 * 00168 * Generate X 00169 * 00170 IF( IMAT.EQ.1 ) THEN 00171 CALL DLAROR( 'L', 'I', M, M, X, LDX, ISEED, WORK, IINFO ) 00172 IF( M .NE. 0 .AND. IINFO .NE. 0 ) THEN 00173 WRITE( NOUT, FMT = 9999 ) M, IINFO 00174 INFO = ABS( IINFO ) 00175 GO TO 20 00176 END IF 00177 ELSE IF( IMAT.EQ.2 ) THEN 00178 R = MIN( P, M-P, Q, M-Q ) 00179 DO I = 1, R 00180 THETA(I) = PIOVER2 * DLARND( 1, ISEED ) 00181 END DO 00182 CALL DLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK ) 00183 DO I = 1, M 00184 DO J = 1, M 00185 X(I+(J-1)*LDX) = X(I+(J-1)*LDX) + 00186 $ ORTH*DLARND(2,ISEED) 00187 END DO 00188 END DO 00189 ELSE 00190 R = MIN( P, M-P, Q, M-Q ) 00191 DO I = 1, R+1 00192 THETA(I) = TEN**(-DLARND(1,ISEED)*GAPDIGIT) 00193 END DO 00194 DO I = 2, R+1 00195 THETA(I) = THETA(I-1) + THETA(I) 00196 END DO 00197 DO I = 1, R 00198 THETA(I) = PIOVER2 * THETA(I) / THETA(R+1) 00199 END DO 00200 CALL DLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK ) 00201 END IF 00202 * 00203 NT = 9 00204 * 00205 CALL DCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T, 00206 $ LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK, 00207 $ RWORK, RESULT ) 00208 * 00209 * Print information about the tests that did not 00210 * pass the threshold. 00211 * 00212 DO 10 I = 1, NT 00213 IF( RESULT( I ).GE.THRESH ) THEN 00214 IF( NFAIL.EQ.0 .AND. FIRSTT ) THEN 00215 FIRSTT = .FALSE. 00216 CALL ALAHDG( NOUT, PATH ) 00217 END IF 00218 WRITE( NOUT, FMT = 9998 )M, P, Q, IMAT, I, 00219 $ RESULT( I ) 00220 NFAIL = NFAIL + 1 00221 END IF 00222 10 CONTINUE 00223 NRUN = NRUN + NT 00224 20 CONTINUE 00225 30 CONTINUE 00226 * 00227 * Print a summary of the results. 00228 * 00229 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, 0 ) 00230 * 00231 9999 FORMAT( ' DLAROR in DCKCSD: M = ', I5, ', INFO = ', I15 ) 00232 9998 FORMAT( ' M=', I4, ' P=', I4, ', Q=', I4, ', type ', I2, 00233 $ ', test ', I2, ', ratio=', G13.6 ) 00234 RETURN 00235 * 00236 * End of DCKCSD 00237 * 00238 END 00239 * 00240 * 00241 * 00242 SUBROUTINE DLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK ) 00243 IMPLICIT NONE 00244 * 00245 INTEGER LDX, M, P, Q 00246 INTEGER ISEED( 4 ) 00247 DOUBLE PRECISION THETA( * ) 00248 DOUBLE PRECISION WORK( * ), X( LDX, * ) 00249 * 00250 DOUBLE PRECISION ONE, ZERO 00251 PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 ) 00252 * 00253 INTEGER I, INFO, R 00254 * 00255 R = MIN( P, M-P, Q, M-Q ) 00256 * 00257 CALL DLASET( 'Full', M, M, ZERO, ZERO, X, LDX ) 00258 * 00259 DO I = 1, MIN(P,Q)-R 00260 X(I,I) = ONE 00261 END DO 00262 DO I = 1, R 00263 X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) = COS(THETA(I)) 00264 END DO 00265 DO I = 1, MIN(P,M-Q)-R 00266 X(P-I+1,M-I+1) = -ONE 00267 END DO 00268 DO I = 1, R 00269 X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) = 00270 $ -SIN(THETA(R-I+1)) 00271 END DO 00272 DO I = 1, MIN(M-P,Q)-R 00273 X(M-I+1,Q-I+1) = ONE 00274 END DO 00275 DO I = 1, R 00276 X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) = 00277 $ SIN(THETA(R-I+1)) 00278 END DO 00279 DO I = 1, MIN(M-P,M-Q)-R 00280 X(P+I,Q+I) = ONE 00281 END DO 00282 DO I = 1, R 00283 X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) = 00284 $ COS(THETA(I)) 00285 END DO 00286 CALL DLAROR( 'Left', 'No init', P, M, X, LDX, ISEED, WORK, INFO ) 00287 CALL DLAROR( 'Left', 'No init', M-P, M, X(P+1,1), LDX, 00288 $ ISEED, WORK, INFO ) 00289 CALL DLAROR( 'Right', 'No init', M, Q, X, LDX, ISEED, 00290 $ WORK, INFO ) 00291 CALL DLAROR( 'Right', 'No init', M, M-Q, 00292 $ X(1,Q+1), LDX, ISEED, WORK, INFO ) 00293 * 00294 END 00295
1.7.3 
