1 SUBROUTINE DCHKQ3( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL,
2 $ THRESH, A, COPYA, S, COPYS, TAU, WORK, IWORK,
3 $ NOUT )
4 *
5 * -- LAPACK test routine (version 3.1.1) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * January 2007
8 *
9 * .. Scalar Arguments ..
10 INTEGER NM, NN, NNB, NOUT
11 DOUBLE PRECISION THRESH
12 * ..
13 * .. Array Arguments ..
14 LOGICAL DOTYPE( * )
15 INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NVAL( * ),
16 $ NXVAL( * )
17 DOUBLE PRECISION A( * ), COPYA( * ), COPYS( * ), S( * ),
18 $ TAU( * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DCHKQ3 tests DGEQP3.
25 *
26 * Arguments
27 * =========
28 *
29 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
30 * The matrix types to be used for testing. Matrices of type j
31 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
32 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
33 *
34 * NM (input) INTEGER
35 * The number of values of M contained in the vector MVAL.
36 *
37 * MVAL (input) INTEGER array, dimension (NM)
38 * The values of the matrix row dimension M.
39 *
40 * NN (input) INTEGER
41 * The number of values of N contained in the vector NVAL.
42 *
43 * NVAL (input) INTEGER array, dimension (NN)
44 * The values of the matrix column dimension N.
45 *
46 * NNB (input) INTEGER
47 * The number of values of NB and NX contained in the
48 * vectors NBVAL and NXVAL. The blocking parameters are used
49 * in pairs (NB,NX).
50 *
51 * NBVAL (input) INTEGER array, dimension (NNB)
52 * The values of the blocksize NB.
53 *
54 * NXVAL (input) INTEGER array, dimension (NNB)
55 * The values of the crossover point NX.
56 *
57 * THRESH (input) DOUBLE PRECISION
58 * The threshold value for the test ratios. A result is
59 * included in the output file if RESULT >= THRESH. To have
60 * every test ratio printed, use THRESH = 0.
61 *
62 * A (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
63 * where MMAX is the maximum value of M in MVAL and NMAX is the
64 * maximum value of N in NVAL.
65 *
66 * COPYA (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
67 *
68 * S (workspace) DOUBLE PRECISION array, dimension
69 * (min(MMAX,NMAX))
70 *
71 * COPYS (workspace) DOUBLE PRECISION array, dimension
72 * (min(MMAX,NMAX))
73 *
74 * TAU (workspace) DOUBLE PRECISION array, dimension (MMAX)
75 *
76 * WORK (workspace) DOUBLE PRECISION array, dimension
77 * (MMAX*NMAX + 4*NMAX + MMAX)
78 *
79 * IWORK (workspace) INTEGER array, dimension (2*NMAX)
80 *
81 * NOUT (input) INTEGER
82 * The unit number for output.
83 *
84 * =====================================================================
85 *
86 * .. Parameters ..
87 INTEGER NTYPES
88 PARAMETER ( NTYPES = 6 )
89 INTEGER NTESTS
90 PARAMETER ( NTESTS = 3 )
91 DOUBLE PRECISION ONE, ZERO
92 PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 )
93 * ..
94 * .. Local Scalars ..
95 CHARACTER*3 PATH
96 INTEGER I, IHIGH, ILOW, IM, IMODE, IN, INB, INFO,
97 $ ISTEP, K, LDA, LW, LWORK, M, MNMIN, MODE, N,
98 $ NB, NERRS, NFAIL, NRUN, NX
99 DOUBLE PRECISION EPS
100 * ..
101 * .. Local Arrays ..
102 INTEGER ISEED( 4 ), ISEEDY( 4 )
103 DOUBLE PRECISION RESULT( NTESTS )
104 * ..
105 * .. External Functions ..
106 DOUBLE PRECISION DLAMCH, DQPT01, DQRT11, DQRT12
107 EXTERNAL DLAMCH, DQPT01, DQRT11, DQRT12
108 * ..
109 * .. External Subroutines ..
110 EXTERNAL ALAHD, ALASUM, DGEQP3, DLACPY, DLAORD, DLASET,
111 $ DLATMS, ICOPY, XLAENV
112 * ..
113 * .. Intrinsic Functions ..
114 INTRINSIC MAX, MIN
115 * ..
116 * .. Scalars in Common ..
117 LOGICAL LERR, OK
118 CHARACTER*32 SRNAMT
119 INTEGER INFOT, IOUNIT
120 * ..
121 * .. Common blocks ..
122 COMMON / INFOC / INFOT, IOUNIT, OK, LERR
123 COMMON / SRNAMC / SRNAMT
124 * ..
125 * .. Data statements ..
126 DATA ISEEDY / 1988, 1989, 1990, 1991 /
127 * ..
128 * .. Executable Statements ..
129 *
130 * Initialize constants and the random number seed.
131 *
132 PATH( 1: 1 ) = 'Double precision'
133 PATH( 2: 3 ) = 'Q3'
134 NRUN = 0
135 NFAIL = 0
136 NERRS = 0
137 DO 10 I = 1, 4
138 ISEED( I ) = ISEEDY( I )
139 10 CONTINUE
140 EPS = DLAMCH( 'Epsilon' )
141 INFOT = 0
142 *
143 DO 90 IM = 1, NM
144 *
145 * Do for each value of M in MVAL.
146 *
147 M = MVAL( IM )
148 LDA = MAX( 1, M )
149 *
150 DO 80 IN = 1, NN
151 *
152 * Do for each value of N in NVAL.
153 *
154 N = NVAL( IN )
155 MNMIN = MIN( M, N )
156 LWORK = MAX( 1, M*MAX( M, N )+4*MNMIN+MAX( M, N ),
157 $ M*N + 2*MNMIN + 4*N )
158 *
159 DO 70 IMODE = 1, NTYPES
160 IF( .NOT.DOTYPE( IMODE ) )
161 $ GO TO 70
162 *
163 * Do for each type of matrix
164 * 1: zero matrix
165 * 2: one small singular value
166 * 3: geometric distribution of singular values
167 * 4: first n/2 columns fixed
168 * 5: last n/2 columns fixed
169 * 6: every second column fixed
170 *
171 MODE = IMODE
172 IF( IMODE.GT.3 )
173 $ MODE = 1
174 *
175 * Generate test matrix of size m by n using
176 * singular value distribution indicated by `mode'.
177 *
178 DO 20 I = 1, N
179 IWORK( I ) = 0
180 20 CONTINUE
181 IF( IMODE.EQ.1 ) THEN
182 CALL DLASET( 'Full', M, N, ZERO, ZERO, COPYA, LDA )
183 DO 30 I = 1, MNMIN
184 COPYS( I ) = ZERO
185 30 CONTINUE
186 ELSE
187 CALL DLATMS( M, N, 'Uniform', ISEED, 'Nonsymm', COPYS,
188 $ MODE, ONE / EPS, ONE, M, N, 'No packing',
189 $ COPYA, LDA, WORK, INFO )
190 IF( IMODE.GE.4 ) THEN
191 IF( IMODE.EQ.4 ) THEN
192 ILOW = 1
193 ISTEP = 1
194 IHIGH = MAX( 1, N / 2 )
195 ELSE IF( IMODE.EQ.5 ) THEN
196 ILOW = MAX( 1, N / 2 )
197 ISTEP = 1
198 IHIGH = N
199 ELSE IF( IMODE.EQ.6 ) THEN
200 ILOW = 1
201 ISTEP = 2
202 IHIGH = N
203 END IF
204 DO 40 I = ILOW, IHIGH, ISTEP
205 IWORK( I ) = 1
206 40 CONTINUE
207 END IF
208 CALL DLAORD( 'Decreasing', MNMIN, COPYS, 1 )
209 END IF
210 *
211 DO 60 INB = 1, NNB
212 *
213 * Do for each pair of values (NB,NX) in NBVAL and NXVAL.
214 *
215 NB = NBVAL( INB )
216 CALL XLAENV( 1, NB )
217 NX = NXVAL( INB )
218 CALL XLAENV( 3, NX )
219 *
220 * Get a working copy of COPYA into A and a copy of
221 * vector IWORK.
222 *
223 CALL DLACPY( 'All', M, N, COPYA, LDA, A, LDA )
224 CALL ICOPY( N, IWORK( 1 ), 1, IWORK( N+1 ), 1 )
225 *
226 * Compute the QR factorization with pivoting of A
227 *
228 LW = MAX( 1, 2*N+NB*( N+1 ) )
229 *
230 * Compute the QP3 factorization of A
231 *
232 SRNAMT = 'DGEQP3'
233 CALL DGEQP3( M, N, A, LDA, IWORK( N+1 ), TAU, WORK,
234 $ LW, INFO )
235 *
236 * Compute norm(svd(a) - svd(r))
237 *
238 RESULT( 1 ) = DQRT12( M, N, A, LDA, COPYS, WORK,
239 $ LWORK )
240 *
241 * Compute norm( A*P - Q*R )
242 *
243 RESULT( 2 ) = DQPT01( M, N, MNMIN, COPYA, A, LDA, TAU,
244 $ IWORK( N+1 ), WORK, LWORK )
245 *
246 * Compute Q'*Q
247 *
248 RESULT( 3 ) = DQRT11( M, MNMIN, A, LDA, TAU, WORK,
249 $ LWORK )
250 *
251 * Print information about the tests that did not pass
252 * the threshold.
253 *
254 DO 50 K = 1, NTESTS
255 IF( RESULT( K ).GE.THRESH ) THEN
256 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
257 $ CALL ALAHD( NOUT, PATH )
258 WRITE( NOUT, FMT = 9999 )'DGEQP3', M, N, NB,
259 $ IMODE, K, RESULT( K )
260 NFAIL = NFAIL + 1
261 END IF
262 50 CONTINUE
263 NRUN = NRUN + NTESTS
264 *
265 60 CONTINUE
266 70 CONTINUE
267 80 CONTINUE
268 90 CONTINUE
269 *
270 * Print a summary of the results.
271 *
272 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
273 *
274 9999 FORMAT( 1X, A, ' M =', I5, ', N =', I5, ', NB =', I4, ', type ',
275 $ I2, ', test ', I2, ', ratio =', G12.5 )
276 *
277 * End of DCHKQ3
278 *
279 END
2 $ THRESH, A, COPYA, S, COPYS, TAU, WORK, IWORK,
3 $ NOUT )
4 *
5 * -- LAPACK test routine (version 3.1.1) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * January 2007
8 *
9 * .. Scalar Arguments ..
10 INTEGER NM, NN, NNB, NOUT
11 DOUBLE PRECISION THRESH
12 * ..
13 * .. Array Arguments ..
14 LOGICAL DOTYPE( * )
15 INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NVAL( * ),
16 $ NXVAL( * )
17 DOUBLE PRECISION A( * ), COPYA( * ), COPYS( * ), S( * ),
18 $ TAU( * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DCHKQ3 tests DGEQP3.
25 *
26 * Arguments
27 * =========
28 *
29 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
30 * The matrix types to be used for testing. Matrices of type j
31 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
32 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
33 *
34 * NM (input) INTEGER
35 * The number of values of M contained in the vector MVAL.
36 *
37 * MVAL (input) INTEGER array, dimension (NM)
38 * The values of the matrix row dimension M.
39 *
40 * NN (input) INTEGER
41 * The number of values of N contained in the vector NVAL.
42 *
43 * NVAL (input) INTEGER array, dimension (NN)
44 * The values of the matrix column dimension N.
45 *
46 * NNB (input) INTEGER
47 * The number of values of NB and NX contained in the
48 * vectors NBVAL and NXVAL. The blocking parameters are used
49 * in pairs (NB,NX).
50 *
51 * NBVAL (input) INTEGER array, dimension (NNB)
52 * The values of the blocksize NB.
53 *
54 * NXVAL (input) INTEGER array, dimension (NNB)
55 * The values of the crossover point NX.
56 *
57 * THRESH (input) DOUBLE PRECISION
58 * The threshold value for the test ratios. A result is
59 * included in the output file if RESULT >= THRESH. To have
60 * every test ratio printed, use THRESH = 0.
61 *
62 * A (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
63 * where MMAX is the maximum value of M in MVAL and NMAX is the
64 * maximum value of N in NVAL.
65 *
66 * COPYA (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
67 *
68 * S (workspace) DOUBLE PRECISION array, dimension
69 * (min(MMAX,NMAX))
70 *
71 * COPYS (workspace) DOUBLE PRECISION array, dimension
72 * (min(MMAX,NMAX))
73 *
74 * TAU (workspace) DOUBLE PRECISION array, dimension (MMAX)
75 *
76 * WORK (workspace) DOUBLE PRECISION array, dimension
77 * (MMAX*NMAX + 4*NMAX + MMAX)
78 *
79 * IWORK (workspace) INTEGER array, dimension (2*NMAX)
80 *
81 * NOUT (input) INTEGER
82 * The unit number for output.
83 *
84 * =====================================================================
85 *
86 * .. Parameters ..
87 INTEGER NTYPES
88 PARAMETER ( NTYPES = 6 )
89 INTEGER NTESTS
90 PARAMETER ( NTESTS = 3 )
91 DOUBLE PRECISION ONE, ZERO
92 PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 )
93 * ..
94 * .. Local Scalars ..
95 CHARACTER*3 PATH
96 INTEGER I, IHIGH, ILOW, IM, IMODE, IN, INB, INFO,
97 $ ISTEP, K, LDA, LW, LWORK, M, MNMIN, MODE, N,
98 $ NB, NERRS, NFAIL, NRUN, NX
99 DOUBLE PRECISION EPS
100 * ..
101 * .. Local Arrays ..
102 INTEGER ISEED( 4 ), ISEEDY( 4 )
103 DOUBLE PRECISION RESULT( NTESTS )
104 * ..
105 * .. External Functions ..
106 DOUBLE PRECISION DLAMCH, DQPT01, DQRT11, DQRT12
107 EXTERNAL DLAMCH, DQPT01, DQRT11, DQRT12
108 * ..
109 * .. External Subroutines ..
110 EXTERNAL ALAHD, ALASUM, DGEQP3, DLACPY, DLAORD, DLASET,
111 $ DLATMS, ICOPY, XLAENV
112 * ..
113 * .. Intrinsic Functions ..
114 INTRINSIC MAX, MIN
115 * ..
116 * .. Scalars in Common ..
117 LOGICAL LERR, OK
118 CHARACTER*32 SRNAMT
119 INTEGER INFOT, IOUNIT
120 * ..
121 * .. Common blocks ..
122 COMMON / INFOC / INFOT, IOUNIT, OK, LERR
123 COMMON / SRNAMC / SRNAMT
124 * ..
125 * .. Data statements ..
126 DATA ISEEDY / 1988, 1989, 1990, 1991 /
127 * ..
128 * .. Executable Statements ..
129 *
130 * Initialize constants and the random number seed.
131 *
132 PATH( 1: 1 ) = 'Double precision'
133 PATH( 2: 3 ) = 'Q3'
134 NRUN = 0
135 NFAIL = 0
136 NERRS = 0
137 DO 10 I = 1, 4
138 ISEED( I ) = ISEEDY( I )
139 10 CONTINUE
140 EPS = DLAMCH( 'Epsilon' )
141 INFOT = 0
142 *
143 DO 90 IM = 1, NM
144 *
145 * Do for each value of M in MVAL.
146 *
147 M = MVAL( IM )
148 LDA = MAX( 1, M )
149 *
150 DO 80 IN = 1, NN
151 *
152 * Do for each value of N in NVAL.
153 *
154 N = NVAL( IN )
155 MNMIN = MIN( M, N )
156 LWORK = MAX( 1, M*MAX( M, N )+4*MNMIN+MAX( M, N ),
157 $ M*N + 2*MNMIN + 4*N )
158 *
159 DO 70 IMODE = 1, NTYPES
160 IF( .NOT.DOTYPE( IMODE ) )
161 $ GO TO 70
162 *
163 * Do for each type of matrix
164 * 1: zero matrix
165 * 2: one small singular value
166 * 3: geometric distribution of singular values
167 * 4: first n/2 columns fixed
168 * 5: last n/2 columns fixed
169 * 6: every second column fixed
170 *
171 MODE = IMODE
172 IF( IMODE.GT.3 )
173 $ MODE = 1
174 *
175 * Generate test matrix of size m by n using
176 * singular value distribution indicated by `mode'.
177 *
178 DO 20 I = 1, N
179 IWORK( I ) = 0
180 20 CONTINUE
181 IF( IMODE.EQ.1 ) THEN
182 CALL DLASET( 'Full', M, N, ZERO, ZERO, COPYA, LDA )
183 DO 30 I = 1, MNMIN
184 COPYS( I ) = ZERO
185 30 CONTINUE
186 ELSE
187 CALL DLATMS( M, N, 'Uniform', ISEED, 'Nonsymm', COPYS,
188 $ MODE, ONE / EPS, ONE, M, N, 'No packing',
189 $ COPYA, LDA, WORK, INFO )
190 IF( IMODE.GE.4 ) THEN
191 IF( IMODE.EQ.4 ) THEN
192 ILOW = 1
193 ISTEP = 1
194 IHIGH = MAX( 1, N / 2 )
195 ELSE IF( IMODE.EQ.5 ) THEN
196 ILOW = MAX( 1, N / 2 )
197 ISTEP = 1
198 IHIGH = N
199 ELSE IF( IMODE.EQ.6 ) THEN
200 ILOW = 1
201 ISTEP = 2
202 IHIGH = N
203 END IF
204 DO 40 I = ILOW, IHIGH, ISTEP
205 IWORK( I ) = 1
206 40 CONTINUE
207 END IF
208 CALL DLAORD( 'Decreasing', MNMIN, COPYS, 1 )
209 END IF
210 *
211 DO 60 INB = 1, NNB
212 *
213 * Do for each pair of values (NB,NX) in NBVAL and NXVAL.
214 *
215 NB = NBVAL( INB )
216 CALL XLAENV( 1, NB )
217 NX = NXVAL( INB )
218 CALL XLAENV( 3, NX )
219 *
220 * Get a working copy of COPYA into A and a copy of
221 * vector IWORK.
222 *
223 CALL DLACPY( 'All', M, N, COPYA, LDA, A, LDA )
224 CALL ICOPY( N, IWORK( 1 ), 1, IWORK( N+1 ), 1 )
225 *
226 * Compute the QR factorization with pivoting of A
227 *
228 LW = MAX( 1, 2*N+NB*( N+1 ) )
229 *
230 * Compute the QP3 factorization of A
231 *
232 SRNAMT = 'DGEQP3'
233 CALL DGEQP3( M, N, A, LDA, IWORK( N+1 ), TAU, WORK,
234 $ LW, INFO )
235 *
236 * Compute norm(svd(a) - svd(r))
237 *
238 RESULT( 1 ) = DQRT12( M, N, A, LDA, COPYS, WORK,
239 $ LWORK )
240 *
241 * Compute norm( A*P - Q*R )
242 *
243 RESULT( 2 ) = DQPT01( M, N, MNMIN, COPYA, A, LDA, TAU,
244 $ IWORK( N+1 ), WORK, LWORK )
245 *
246 * Compute Q'*Q
247 *
248 RESULT( 3 ) = DQRT11( M, MNMIN, A, LDA, TAU, WORK,
249 $ LWORK )
250 *
251 * Print information about the tests that did not pass
252 * the threshold.
253 *
254 DO 50 K = 1, NTESTS
255 IF( RESULT( K ).GE.THRESH ) THEN
256 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
257 $ CALL ALAHD( NOUT, PATH )
258 WRITE( NOUT, FMT = 9999 )'DGEQP3', M, N, NB,
259 $ IMODE, K, RESULT( K )
260 NFAIL = NFAIL + 1
261 END IF
262 50 CONTINUE
263 NRUN = NRUN + NTESTS
264 *
265 60 CONTINUE
266 70 CONTINUE
267 80 CONTINUE
268 90 CONTINUE
269 *
270 * Print a summary of the results.
271 *
272 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
273 *
274 9999 FORMAT( 1X, A, ' M =', I5, ', N =', I5, ', NB =', I4, ', type ',
275 $ I2, ', test ', I2, ', ratio =', G12.5 )
276 *
277 * End of DCHKQ3
278 *
279 END