1 SUBROUTINE DQRT15( SCALE, RKSEL, M, N, NRHS, A, LDA, B, LDB, S,
2 $ RANK, NORMA, NORMB, ISEED, WORK, LWORK )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER LDA, LDB, LWORK, M, N, NRHS, RANK, RKSEL, SCALE
10 DOUBLE PRECISION NORMA, NORMB
11 * ..
12 * .. Array Arguments ..
13 INTEGER ISEED( 4 )
14 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( LWORK )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DQRT15 generates a matrix with full or deficient rank and of various
21 * norms.
22 *
23 * Arguments
24 * =========
25 *
26 * SCALE (input) INTEGER
27 * SCALE = 1: normally scaled matrix
28 * SCALE = 2: matrix scaled up
29 * SCALE = 3: matrix scaled down
30 *
31 * RKSEL (input) INTEGER
32 * RKSEL = 1: full rank matrix
33 * RKSEL = 2: rank-deficient matrix
34 *
35 * M (input) INTEGER
36 * The number of rows of the matrix A.
37 *
38 * N (input) INTEGER
39 * The number of columns of A.
40 *
41 * NRHS (input) INTEGER
42 * The number of columns of B.
43 *
44 * A (output) DOUBLE PRECISION array, dimension (LDA,N)
45 * The M-by-N matrix A.
46 *
47 * LDA (input) INTEGER
48 * The leading dimension of the array A.
49 *
50 * B (output) DOUBLE PRECISION array, dimension (LDB, NRHS)
51 * A matrix that is in the range space of matrix A.
52 *
53 * LDB (input) INTEGER
54 * The leading dimension of the array B.
55 *
56 * S (output) DOUBLE PRECISION array, dimension MIN(M,N)
57 * Singular values of A.
58 *
59 * RANK (output) INTEGER
60 * number of nonzero singular values of A.
61 *
62 * NORMA (output) DOUBLE PRECISION
63 * one-norm of A.
64 *
65 * NORMB (output) DOUBLE PRECISION
66 * one-norm of B.
67 *
68 * ISEED (input/output) integer array, dimension (4)
69 * seed for random number generator.
70 *
71 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK)
72 *
73 * LWORK (input) INTEGER
74 * length of work space required.
75 * LWORK >= MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
76 *
77 * =====================================================================
78 *
79 * .. Parameters ..
80 DOUBLE PRECISION ZERO, ONE, TWO, SVMIN
81 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
82 $ SVMIN = 0.1D0 )
83 * ..
84 * .. Local Scalars ..
85 INTEGER INFO, J, MN
86 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TEMP
87 * ..
88 * .. Local Arrays ..
89 DOUBLE PRECISION DUMMY( 1 )
90 * ..
91 * .. External Functions ..
92 DOUBLE PRECISION DASUM, DLAMCH, DLANGE, DLARND, DNRM2
93 EXTERNAL DASUM, DLAMCH, DLANGE, DLARND, DNRM2
94 * ..
95 * .. External Subroutines ..
96 EXTERNAL DGEMM, DLAORD, DLARF, DLARNV, DLAROR, DLASCL,
97 $ DLASET, DSCAL, XERBLA
98 * ..
99 * .. Intrinsic Functions ..
100 INTRINSIC ABS, MAX, MIN
101 * ..
102 * .. Executable Statements ..
103 *
104 MN = MIN( M, N )
105 IF( LWORK.LT.MAX( M+MN, MN*NRHS, 2*N+M ) ) THEN
106 CALL XERBLA( 'DQRT15', 16 )
107 RETURN
108 END IF
109 *
110 SMLNUM = DLAMCH( 'Safe minimum' )
111 BIGNUM = ONE / SMLNUM
112 EPS = DLAMCH( 'Epsilon' )
113 SMLNUM = ( SMLNUM / EPS ) / EPS
114 BIGNUM = ONE / SMLNUM
115 *
116 * Determine rank and (unscaled) singular values
117 *
118 IF( RKSEL.EQ.1 ) THEN
119 RANK = MN
120 ELSE IF( RKSEL.EQ.2 ) THEN
121 RANK = ( 3*MN ) / 4
122 DO 10 J = RANK + 1, MN
123 S( J ) = ZERO
124 10 CONTINUE
125 ELSE
126 CALL XERBLA( 'DQRT15', 2 )
127 END IF
128 *
129 IF( RANK.GT.0 ) THEN
130 *
131 * Nontrivial case
132 *
133 S( 1 ) = ONE
134 DO 30 J = 2, RANK
135 20 CONTINUE
136 TEMP = DLARND( 1, ISEED )
137 IF( TEMP.GT.SVMIN ) THEN
138 S( J ) = ABS( TEMP )
139 ELSE
140 GO TO 20
141 END IF
142 30 CONTINUE
143 CALL DLAORD( 'Decreasing', RANK, S, 1 )
144 *
145 * Generate 'rank' columns of a random orthogonal matrix in A
146 *
147 CALL DLARNV( 2, ISEED, M, WORK )
148 CALL DSCAL( M, ONE / DNRM2( M, WORK, 1 ), WORK, 1 )
149 CALL DLASET( 'Full', M, RANK, ZERO, ONE, A, LDA )
150 CALL DLARF( 'Left', M, RANK, WORK, 1, TWO, A, LDA,
151 $ WORK( M+1 ) )
152 *
153 * workspace used: m+mn
154 *
155 * Generate consistent rhs in the range space of A
156 *
157 CALL DLARNV( 2, ISEED, RANK*NRHS, WORK )
158 CALL DGEMM( 'No transpose', 'No transpose', M, NRHS, RANK, ONE,
159 $ A, LDA, WORK, RANK, ZERO, B, LDB )
160 *
161 * work space used: <= mn *nrhs
162 *
163 * generate (unscaled) matrix A
164 *
165 DO 40 J = 1, RANK
166 CALL DSCAL( M, S( J ), A( 1, J ), 1 )
167 40 CONTINUE
168 IF( RANK.LT.N )
169 $ CALL DLASET( 'Full', M, N-RANK, ZERO, ZERO, A( 1, RANK+1 ),
170 $ LDA )
171 CALL DLAROR( 'Right', 'No initialization', M, N, A, LDA, ISEED,
172 $ WORK, INFO )
173 *
174 ELSE
175 *
176 * work space used 2*n+m
177 *
178 * Generate null matrix and rhs
179 *
180 DO 50 J = 1, MN
181 S( J ) = ZERO
182 50 CONTINUE
183 CALL DLASET( 'Full', M, N, ZERO, ZERO, A, LDA )
184 CALL DLASET( 'Full', M, NRHS, ZERO, ZERO, B, LDB )
185 *
186 END IF
187 *
188 * Scale the matrix
189 *
190 IF( SCALE.NE.1 ) THEN
191 NORMA = DLANGE( 'Max', M, N, A, LDA, DUMMY )
192 IF( NORMA.NE.ZERO ) THEN
193 IF( SCALE.EQ.2 ) THEN
194 *
195 * matrix scaled up
196 *
197 CALL DLASCL( 'General', 0, 0, NORMA, BIGNUM, M, N, A,
198 $ LDA, INFO )
199 CALL DLASCL( 'General', 0, 0, NORMA, BIGNUM, MN, 1, S,
200 $ MN, INFO )
201 CALL DLASCL( 'General', 0, 0, NORMA, BIGNUM, M, NRHS, B,
202 $ LDB, INFO )
203 ELSE IF( SCALE.EQ.3 ) THEN
204 *
205 * matrix scaled down
206 *
207 CALL DLASCL( 'General', 0, 0, NORMA, SMLNUM, M, N, A,
208 $ LDA, INFO )
209 CALL DLASCL( 'General', 0, 0, NORMA, SMLNUM, MN, 1, S,
210 $ MN, INFO )
211 CALL DLASCL( 'General', 0, 0, NORMA, SMLNUM, M, NRHS, B,
212 $ LDB, INFO )
213 ELSE
214 CALL XERBLA( 'DQRT15', 1 )
215 RETURN
216 END IF
217 END IF
218 END IF
219 *
220 NORMA = DASUM( MN, S, 1 )
221 NORMB = DLANGE( 'One-norm', M, NRHS, B, LDB, DUMMY )
222 *
223 RETURN
224 *
225 * End of DQRT15
226 *
227 END
2 $ RANK, NORMA, NORMB, ISEED, WORK, LWORK )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER LDA, LDB, LWORK, M, N, NRHS, RANK, RKSEL, SCALE
10 DOUBLE PRECISION NORMA, NORMB
11 * ..
12 * .. Array Arguments ..
13 INTEGER ISEED( 4 )
14 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( LWORK )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DQRT15 generates a matrix with full or deficient rank and of various
21 * norms.
22 *
23 * Arguments
24 * =========
25 *
26 * SCALE (input) INTEGER
27 * SCALE = 1: normally scaled matrix
28 * SCALE = 2: matrix scaled up
29 * SCALE = 3: matrix scaled down
30 *
31 * RKSEL (input) INTEGER
32 * RKSEL = 1: full rank matrix
33 * RKSEL = 2: rank-deficient matrix
34 *
35 * M (input) INTEGER
36 * The number of rows of the matrix A.
37 *
38 * N (input) INTEGER
39 * The number of columns of A.
40 *
41 * NRHS (input) INTEGER
42 * The number of columns of B.
43 *
44 * A (output) DOUBLE PRECISION array, dimension (LDA,N)
45 * The M-by-N matrix A.
46 *
47 * LDA (input) INTEGER
48 * The leading dimension of the array A.
49 *
50 * B (output) DOUBLE PRECISION array, dimension (LDB, NRHS)
51 * A matrix that is in the range space of matrix A.
52 *
53 * LDB (input) INTEGER
54 * The leading dimension of the array B.
55 *
56 * S (output) DOUBLE PRECISION array, dimension MIN(M,N)
57 * Singular values of A.
58 *
59 * RANK (output) INTEGER
60 * number of nonzero singular values of A.
61 *
62 * NORMA (output) DOUBLE PRECISION
63 * one-norm of A.
64 *
65 * NORMB (output) DOUBLE PRECISION
66 * one-norm of B.
67 *
68 * ISEED (input/output) integer array, dimension (4)
69 * seed for random number generator.
70 *
71 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK)
72 *
73 * LWORK (input) INTEGER
74 * length of work space required.
75 * LWORK >= MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
76 *
77 * =====================================================================
78 *
79 * .. Parameters ..
80 DOUBLE PRECISION ZERO, ONE, TWO, SVMIN
81 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
82 $ SVMIN = 0.1D0 )
83 * ..
84 * .. Local Scalars ..
85 INTEGER INFO, J, MN
86 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TEMP
87 * ..
88 * .. Local Arrays ..
89 DOUBLE PRECISION DUMMY( 1 )
90 * ..
91 * .. External Functions ..
92 DOUBLE PRECISION DASUM, DLAMCH, DLANGE, DLARND, DNRM2
93 EXTERNAL DASUM, DLAMCH, DLANGE, DLARND, DNRM2
94 * ..
95 * .. External Subroutines ..
96 EXTERNAL DGEMM, DLAORD, DLARF, DLARNV, DLAROR, DLASCL,
97 $ DLASET, DSCAL, XERBLA
98 * ..
99 * .. Intrinsic Functions ..
100 INTRINSIC ABS, MAX, MIN
101 * ..
102 * .. Executable Statements ..
103 *
104 MN = MIN( M, N )
105 IF( LWORK.LT.MAX( M+MN, MN*NRHS, 2*N+M ) ) THEN
106 CALL XERBLA( 'DQRT15', 16 )
107 RETURN
108 END IF
109 *
110 SMLNUM = DLAMCH( 'Safe minimum' )
111 BIGNUM = ONE / SMLNUM
112 EPS = DLAMCH( 'Epsilon' )
113 SMLNUM = ( SMLNUM / EPS ) / EPS
114 BIGNUM = ONE / SMLNUM
115 *
116 * Determine rank and (unscaled) singular values
117 *
118 IF( RKSEL.EQ.1 ) THEN
119 RANK = MN
120 ELSE IF( RKSEL.EQ.2 ) THEN
121 RANK = ( 3*MN ) / 4
122 DO 10 J = RANK + 1, MN
123 S( J ) = ZERO
124 10 CONTINUE
125 ELSE
126 CALL XERBLA( 'DQRT15', 2 )
127 END IF
128 *
129 IF( RANK.GT.0 ) THEN
130 *
131 * Nontrivial case
132 *
133 S( 1 ) = ONE
134 DO 30 J = 2, RANK
135 20 CONTINUE
136 TEMP = DLARND( 1, ISEED )
137 IF( TEMP.GT.SVMIN ) THEN
138 S( J ) = ABS( TEMP )
139 ELSE
140 GO TO 20
141 END IF
142 30 CONTINUE
143 CALL DLAORD( 'Decreasing', RANK, S, 1 )
144 *
145 * Generate 'rank' columns of a random orthogonal matrix in A
146 *
147 CALL DLARNV( 2, ISEED, M, WORK )
148 CALL DSCAL( M, ONE / DNRM2( M, WORK, 1 ), WORK, 1 )
149 CALL DLASET( 'Full', M, RANK, ZERO, ONE, A, LDA )
150 CALL DLARF( 'Left', M, RANK, WORK, 1, TWO, A, LDA,
151 $ WORK( M+1 ) )
152 *
153 * workspace used: m+mn
154 *
155 * Generate consistent rhs in the range space of A
156 *
157 CALL DLARNV( 2, ISEED, RANK*NRHS, WORK )
158 CALL DGEMM( 'No transpose', 'No transpose', M, NRHS, RANK, ONE,
159 $ A, LDA, WORK, RANK, ZERO, B, LDB )
160 *
161 * work space used: <= mn *nrhs
162 *
163 * generate (unscaled) matrix A
164 *
165 DO 40 J = 1, RANK
166 CALL DSCAL( M, S( J ), A( 1, J ), 1 )
167 40 CONTINUE
168 IF( RANK.LT.N )
169 $ CALL DLASET( 'Full', M, N-RANK, ZERO, ZERO, A( 1, RANK+1 ),
170 $ LDA )
171 CALL DLAROR( 'Right', 'No initialization', M, N, A, LDA, ISEED,
172 $ WORK, INFO )
173 *
174 ELSE
175 *
176 * work space used 2*n+m
177 *
178 * Generate null matrix and rhs
179 *
180 DO 50 J = 1, MN
181 S( J ) = ZERO
182 50 CONTINUE
183 CALL DLASET( 'Full', M, N, ZERO, ZERO, A, LDA )
184 CALL DLASET( 'Full', M, NRHS, ZERO, ZERO, B, LDB )
185 *
186 END IF
187 *
188 * Scale the matrix
189 *
190 IF( SCALE.NE.1 ) THEN
191 NORMA = DLANGE( 'Max', M, N, A, LDA, DUMMY )
192 IF( NORMA.NE.ZERO ) THEN
193 IF( SCALE.EQ.2 ) THEN
194 *
195 * matrix scaled up
196 *
197 CALL DLASCL( 'General', 0, 0, NORMA, BIGNUM, M, N, A,
198 $ LDA, INFO )
199 CALL DLASCL( 'General', 0, 0, NORMA, BIGNUM, MN, 1, S,
200 $ MN, INFO )
201 CALL DLASCL( 'General', 0, 0, NORMA, BIGNUM, M, NRHS, B,
202 $ LDB, INFO )
203 ELSE IF( SCALE.EQ.3 ) THEN
204 *
205 * matrix scaled down
206 *
207 CALL DLASCL( 'General', 0, 0, NORMA, SMLNUM, M, N, A,
208 $ LDA, INFO )
209 CALL DLASCL( 'General', 0, 0, NORMA, SMLNUM, MN, 1, S,
210 $ MN, INFO )
211 CALL DLASCL( 'General', 0, 0, NORMA, SMLNUM, M, NRHS, B,
212 $ LDB, INFO )
213 ELSE
214 CALL XERBLA( 'DQRT15', 1 )
215 RETURN
216 END IF
217 END IF
218 END IF
219 *
220 NORMA = DASUM( MN, S, 1 )
221 NORMB = DLANGE( 'One-norm', M, NRHS, B, LDB, DUMMY )
222 *
223 RETURN
224 *
225 * End of DQRT15
226 *
227 END