 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DQRT15( SCALE, RKSEL, M, N, NRHS, A, LDA, B, LDB, S, 00002 $ RANK, NORMA, NORMB, ISEED, WORK, LWORK ) 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 INTEGER LDA, LDB, LWORK, M, N, NRHS, RANK, RKSEL, SCALE 00010 DOUBLE PRECISION NORMA, NORMB 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER ISEED( 4 ) 00014 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( LWORK ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DQRT15 generates a matrix with full or deficient rank and of various 00021 * norms. 00022 * 00023 * Arguments 00024 * ========= 00025 * 00026 * SCALE (input) INTEGER 00027 * SCALE = 1: normally scaled matrix 00028 * SCALE = 2: matrix scaled up 00029 * SCALE = 3: matrix scaled down 00030 * 00031 * RKSEL (input) INTEGER 00032 * RKSEL = 1: full rank matrix 00033 * RKSEL = 2: rank-deficient matrix 00034 * 00035 * M (input) INTEGER 00036 * The number of rows of the matrix A. 00037 * 00038 * N (input) INTEGER 00039 * The number of columns of A. 00040 * 00041 * NRHS (input) INTEGER 00042 * The number of columns of B. 00043 * 00044 * A (output) DOUBLE PRECISION array, dimension (LDA,N) 00045 * The M-by-N matrix A. 00046 * 00047 * LDA (input) INTEGER 00048 * The leading dimension of the array A. 00049 * 00050 * B (output) DOUBLE PRECISION array, dimension (LDB, NRHS) 00051 * A matrix that is in the range space of matrix A. 00052 * 00053 * LDB (input) INTEGER 00054 * The leading dimension of the array B. 00055 * 00056 * S (output) DOUBLE PRECISION array, dimension MIN(M,N) 00057 * Singular values of A. 00058 * 00059 * RANK (output) INTEGER 00060 * number of nonzero singular values of A. 00061 * 00062 * NORMA (output) DOUBLE PRECISION 00063 * one-norm of A. 00064 * 00065 * NORMB (output) DOUBLE PRECISION 00066 * one-norm of B. 00067 * 00068 * ISEED (input/output) integer array, dimension (4) 00069 * seed for random number generator. 00070 * 00071 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 00072 * 00073 * LWORK (input) INTEGER 00074 * length of work space required. 00075 * LWORK >= MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M) 00076 * 00077 * ===================================================================== 00078 * 00079 * .. Parameters .. 00080 DOUBLE PRECISION ZERO, ONE, TWO, SVMIN 00081 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 00082 $ SVMIN = 0.1D0 ) 00083 * .. 00084 * .. Local Scalars .. 00085 INTEGER INFO, J, MN 00086 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TEMP 00087 * .. 00088 * .. Local Arrays .. 00089 DOUBLE PRECISION DUMMY( 1 ) 00090 * .. 00091 * .. External Functions .. 00092 DOUBLE PRECISION DASUM, DLAMCH, DLANGE, DLARND, DNRM2 00093 EXTERNAL DASUM, DLAMCH, DLANGE, DLARND, DNRM2 00094 * .. 00095 * .. External Subroutines .. 00096 EXTERNAL DGEMM, DLAORD, DLARF, DLARNV, DLAROR, DLASCL, 00097 $ DLASET, DSCAL, XERBLA 00098 * .. 00099 * .. Intrinsic Functions .. 00100 INTRINSIC ABS, MAX, MIN 00101 * .. 00102 * .. Executable Statements .. 00103 * 00104 MN = MIN( M, N ) 00105 IF( LWORK.LT.MAX( M+MN, MN*NRHS, 2*N+M ) ) THEN 00106 CALL XERBLA( 'DQRT15', 16 ) 00107 RETURN 00108 END IF 00109 * 00110 SMLNUM = DLAMCH( 'Safe minimum' ) 00111 BIGNUM = ONE / SMLNUM 00112 EPS = DLAMCH( 'Epsilon' ) 00113 SMLNUM = ( SMLNUM / EPS ) / EPS 00114 BIGNUM = ONE / SMLNUM 00115 * 00116 * Determine rank and (unscaled) singular values 00117 * 00118 IF( RKSEL.EQ.1 ) THEN 00119 RANK = MN 00120 ELSE IF( RKSEL.EQ.2 ) THEN 00121 RANK = ( 3*MN ) / 4 00122 DO 10 J = RANK + 1, MN 00123 S( J ) = ZERO 00124 10 CONTINUE 00125 ELSE 00126 CALL XERBLA( 'DQRT15', 2 ) 00127 END IF 00128 * 00129 IF( RANK.GT.0 ) THEN 00130 * 00131 * Nontrivial case 00132 * 00133 S( 1 ) = ONE 00134 DO 30 J = 2, RANK 00135 20 CONTINUE 00136 TEMP = DLARND( 1, ISEED ) 00137 IF( TEMP.GT.SVMIN ) THEN 00138 S( J ) = ABS( TEMP ) 00139 ELSE 00140 GO TO 20 00141 END IF 00142 30 CONTINUE 00143 CALL DLAORD( 'Decreasing', RANK, S, 1 ) 00144 * 00145 * Generate 'rank' columns of a random orthogonal matrix in A 00146 * 00147 CALL DLARNV( 2, ISEED, M, WORK ) 00148 CALL DSCAL( M, ONE / DNRM2( M, WORK, 1 ), WORK, 1 ) 00149 CALL DLASET( 'Full', M, RANK, ZERO, ONE, A, LDA ) 00150 CALL DLARF( 'Left', M, RANK, WORK, 1, TWO, A, LDA, 00151 $ WORK( M+1 ) ) 00152 * 00153 * workspace used: m+mn 00154 * 00155 * Generate consistent rhs in the range space of A 00156 * 00157 CALL DLARNV( 2, ISEED, RANK*NRHS, WORK ) 00158 CALL DGEMM( 'No transpose', 'No transpose', M, NRHS, RANK, ONE, 00159 $ A, LDA, WORK, RANK, ZERO, B, LDB ) 00160 * 00161 * work space used: <= mn *nrhs 00162 * 00163 * generate (unscaled) matrix A 00164 * 00165 DO 40 J = 1, RANK 00166 CALL DSCAL( M, S( J ), A( 1, J ), 1 ) 00167 40 CONTINUE 00168 IF( RANK.LT.N ) 00169 $ CALL DLASET( 'Full', M, N-RANK, ZERO, ZERO, A( 1, RANK+1 ), 00170 $ LDA ) 00171 CALL DLAROR( 'Right', 'No initialization', M, N, A, LDA, ISEED, 00172 $ WORK, INFO ) 00173 * 00174 ELSE 00175 * 00176 * work space used 2*n+m 00177 * 00178 * Generate null matrix and rhs 00179 * 00180 DO 50 J = 1, MN 00181 S( J ) = ZERO 00182 50 CONTINUE 00183 CALL DLASET( 'Full', M, N, ZERO, ZERO, A, LDA ) 00184 CALL DLASET( 'Full', M, NRHS, ZERO, ZERO, B, LDB ) 00185 * 00186 END IF 00187 * 00188 * Scale the matrix 00189 * 00190 IF( SCALE.NE.1 ) THEN 00191 NORMA = DLANGE( 'Max', M, N, A, LDA, DUMMY ) 00192 IF( NORMA.NE.ZERO ) THEN 00193 IF( SCALE.EQ.2 ) THEN 00194 * 00195 * matrix scaled up 00196 * 00197 CALL DLASCL( 'General', 0, 0, NORMA, BIGNUM, M, N, A, 00198 $ LDA, INFO ) 00199 CALL DLASCL( 'General', 0, 0, NORMA, BIGNUM, MN, 1, S, 00200 $ MN, INFO ) 00201 CALL DLASCL( 'General', 0, 0, NORMA, BIGNUM, M, NRHS, B, 00202 $ LDB, INFO ) 00203 ELSE IF( SCALE.EQ.3 ) THEN 00204 * 00205 * matrix scaled down 00206 * 00207 CALL DLASCL( 'General', 0, 0, NORMA, SMLNUM, M, N, A, 00208 $ LDA, INFO ) 00209 CALL DLASCL( 'General', 0, 0, NORMA, SMLNUM, MN, 1, S, 00210 $ MN, INFO ) 00211 CALL DLASCL( 'General', 0, 0, NORMA, SMLNUM, M, NRHS, B, 00212 $ LDB, INFO ) 00213 ELSE 00214 CALL XERBLA( 'DQRT15', 1 ) 00215 RETURN 00216 END IF 00217 END IF 00218 END IF 00219 * 00220 NORMA = DASUM( MN, S, 1 ) 00221 NORMB = DLANGE( 'One-norm', M, NRHS, B, LDB, DUMMY ) 00222 * 00223 RETURN 00224 * 00225 * End of DQRT15 00226 * 00227 END
1.7.3 
