1 SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
2 $ NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
3 $ COPYB, C, S, COPYS, WORK, IWORK, 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 LOGICAL TSTERR
11 INTEGER NM, NN, NNB, NNS, NOUT
12 DOUBLE PRECISION THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL DOTYPE( * )
16 INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
17 $ NVAL( * ), NXVAL( * )
18 DOUBLE PRECISION A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
19 $ COPYS( * ), S( * ), WORK( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DDRVLS tests the least squares driver routines DGELS, DGELSS, DGELSX,
26 * DGELSY and DGELSD.
27 *
28 * Arguments
29 * =========
30 *
31 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
32 * The matrix types to be used for testing. Matrices of type j
33 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
34 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
35 * The matrix of type j is generated as follows:
36 * j=1: A = U*D*V where U and V are random orthogonal matrices
37 * and D has random entries (> 0.1) taken from a uniform
38 * distribution (0,1). A is full rank.
39 * j=2: The same of 1, but A is scaled up.
40 * j=3: The same of 1, but A is scaled down.
41 * j=4: A = U*D*V where U and V are random orthogonal matrices
42 * and D has 3*min(M,N)/4 random entries (> 0.1) taken
43 * from a uniform distribution (0,1) and the remaining
44 * entries set to 0. A is rank-deficient.
45 * j=5: The same of 4, but A is scaled up.
46 * j=6: The same of 5, but A is scaled down.
47 *
48 * NM (input) INTEGER
49 * The number of values of M contained in the vector MVAL.
50 *
51 * MVAL (input) INTEGER array, dimension (NM)
52 * The values of the matrix row dimension M.
53 *
54 * NN (input) INTEGER
55 * The number of values of N contained in the vector NVAL.
56 *
57 * NVAL (input) INTEGER array, dimension (NN)
58 * The values of the matrix column dimension N.
59 *
60 * NNS (input) INTEGER
61 * The number of values of NRHS contained in the vector NSVAL.
62 *
63 * NSVAL (input) INTEGER array, dimension (NNS)
64 * The values of the number of right hand sides NRHS.
65 *
66 * NNB (input) INTEGER
67 * The number of values of NB and NX contained in the
68 * vectors NBVAL and NXVAL. The blocking parameters are used
69 * in pairs (NB,NX).
70 *
71 * NBVAL (input) INTEGER array, dimension (NNB)
72 * The values of the blocksize NB.
73 *
74 * NXVAL (input) INTEGER array, dimension (NNB)
75 * The values of the crossover point NX.
76 *
77 * THRESH (input) DOUBLE PRECISION
78 * The threshold value for the test ratios. A result is
79 * included in the output file if RESULT >= THRESH. To have
80 * every test ratio printed, use THRESH = 0.
81 *
82 * TSTERR (input) LOGICAL
83 * Flag that indicates whether error exits are to be tested.
84 *
85 * A (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
86 * where MMAX is the maximum value of M in MVAL and NMAX is the
87 * maximum value of N in NVAL.
88 *
89 * COPYA (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
90 *
91 * B (workspace) DOUBLE PRECISION array, dimension (MMAX*NSMAX)
92 * where MMAX is the maximum value of M in MVAL and NSMAX is the
93 * maximum value of NRHS in NSVAL.
94 *
95 * COPYB (workspace) DOUBLE PRECISION array, dimension (MMAX*NSMAX)
96 *
97 * C (workspace) DOUBLE PRECISION array, dimension (MMAX*NSMAX)
98 *
99 * S (workspace) DOUBLE PRECISION array, dimension
100 * (min(MMAX,NMAX))
101 *
102 * COPYS (workspace) DOUBLE PRECISION array, dimension
103 * (min(MMAX,NMAX))
104 *
105 * WORK (workspace) DOUBLE PRECISION array,
106 * dimension (MMAX*NMAX + 4*NMAX + MMAX).
107 *
108 * IWORK (workspace) INTEGER array, dimension (15*NMAX)
109 *
110 * NOUT (input) INTEGER
111 * The unit number for output.
112 *
113 * =====================================================================
114 *
115 * .. Parameters ..
116 INTEGER NTESTS
117 PARAMETER ( NTESTS = 18 )
118 INTEGER SMLSIZ
119 PARAMETER ( SMLSIZ = 25 )
120 DOUBLE PRECISION ONE, TWO, ZERO
121 PARAMETER ( ONE = 1.0D0, TWO = 2.0D0, ZERO = 0.0D0 )
122 * ..
123 * .. Local Scalars ..
124 CHARACTER TRANS
125 CHARACTER*3 PATH
126 INTEGER CRANK, I, IM, IN, INB, INFO, INS, IRANK,
127 $ ISCALE, ITRAN, ITYPE, J, K, LDA, LDB, LDWORK,
128 $ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS,
129 $ NFAIL, NLVL, NRHS, NROWS, NRUN, RANK
130 DOUBLE PRECISION EPS, NORMA, NORMB, RCOND
131 * ..
132 * .. Local Arrays ..
133 INTEGER ISEED( 4 ), ISEEDY( 4 )
134 DOUBLE PRECISION RESULT( NTESTS )
135 * ..
136 * .. External Functions ..
137 DOUBLE PRECISION DASUM, DLAMCH, DQRT12, DQRT14, DQRT17
138 EXTERNAL DASUM, DLAMCH, DQRT12, DQRT14, DQRT17
139 * ..
140 * .. External Subroutines ..
141 EXTERNAL ALAERH, ALAHD, ALASVM, DAXPY, DERRLS, DGELS,
142 $ DGELSD, DGELSS, DGELSX, DGELSY, DGEMM, DLACPY,
143 $ DLARNV, DLASRT, DQRT13, DQRT15, DQRT16, DSCAL,
144 $ XLAENV
145 * ..
146 * .. Intrinsic Functions ..
147 INTRINSIC DBLE, INT, LOG, MAX, MIN, SQRT
148 * ..
149 * .. Scalars in Common ..
150 LOGICAL LERR, OK
151 CHARACTER*32 SRNAMT
152 INTEGER INFOT, IOUNIT
153 * ..
154 * .. Common blocks ..
155 COMMON / INFOC / INFOT, IOUNIT, OK, LERR
156 COMMON / SRNAMC / SRNAMT
157 * ..
158 * .. Data statements ..
159 DATA ISEEDY / 1988, 1989, 1990, 1991 /
160 * ..
161 * .. Executable Statements ..
162 *
163 * Initialize constants and the random number seed.
164 *
165 PATH( 1: 1 ) = 'Double precision'
166 PATH( 2: 3 ) = 'LS'
167 NRUN = 0
168 NFAIL = 0
169 NERRS = 0
170 DO 10 I = 1, 4
171 ISEED( I ) = ISEEDY( I )
172 10 CONTINUE
173 EPS = DLAMCH( 'Epsilon' )
174 *
175 * Threshold for rank estimation
176 *
177 RCOND = SQRT( EPS ) - ( SQRT( EPS )-EPS ) / 2
178 *
179 * Test the error exits
180 *
181 CALL XLAENV( 2, 2 )
182 CALL XLAENV( 9, SMLSIZ )
183 IF( TSTERR )
184 $ CALL DERRLS( PATH, NOUT )
185 *
186 * Print the header if NM = 0 or NN = 0 and THRESH = 0.
187 *
188 IF( ( NM.EQ.0 .OR. NN.EQ.0 ) .AND. THRESH.EQ.ZERO )
189 $ CALL ALAHD( NOUT, PATH )
190 INFOT = 0
191 CALL XLAENV( 2, 2 )
192 CALL XLAENV( 9, SMLSIZ )
193 *
194 DO 150 IM = 1, NM
195 M = MVAL( IM )
196 LDA = MAX( 1, M )
197 *
198 DO 140 IN = 1, NN
199 N = NVAL( IN )
200 MNMIN = MIN( M, N )
201 LDB = MAX( 1, M, N )
202 *
203 DO 130 INS = 1, NNS
204 NRHS = NSVAL( INS )
205 NLVL = MAX( INT( LOG( MAX( ONE, DBLE( MNMIN ) ) /
206 $ DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1, 0 )
207 LWORK = MAX( 1, ( M+NRHS )*( N+2 ), ( N+NRHS )*( M+2 ),
208 $ M*N+4*MNMIN+MAX( M, N ), 12*MNMIN+2*MNMIN*SMLSIZ+
209 $ 8*MNMIN*NLVL+MNMIN*NRHS+(SMLSIZ+1)**2 )
210 *
211 DO 120 IRANK = 1, 2
212 DO 110 ISCALE = 1, 3
213 ITYPE = ( IRANK-1 )*3 + ISCALE
214 IF( .NOT.DOTYPE( ITYPE ) )
215 $ GO TO 110
216 *
217 IF( IRANK.EQ.1 ) THEN
218 *
219 * Test DGELS
220 *
221 * Generate a matrix of scaling type ISCALE
222 *
223 CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
224 $ ISEED )
225 DO 40 INB = 1, NNB
226 NB = NBVAL( INB )
227 CALL XLAENV( 1, NB )
228 CALL XLAENV( 3, NXVAL( INB ) )
229 *
230 DO 30 ITRAN = 1, 2
231 IF( ITRAN.EQ.1 ) THEN
232 TRANS = 'N'
233 NROWS = M
234 NCOLS = N
235 ELSE
236 TRANS = 'T'
237 NROWS = N
238 NCOLS = M
239 END IF
240 LDWORK = MAX( 1, NCOLS )
241 *
242 * Set up a consistent rhs
243 *
244 IF( NCOLS.GT.0 ) THEN
245 CALL DLARNV( 2, ISEED, NCOLS*NRHS,
246 $ WORK )
247 CALL DSCAL( NCOLS*NRHS,
248 $ ONE / DBLE( NCOLS ), WORK,
249 $ 1 )
250 END IF
251 CALL DGEMM( TRANS, 'No transpose', NROWS,
252 $ NRHS, NCOLS, ONE, COPYA, LDA,
253 $ WORK, LDWORK, ZERO, B, LDB )
254 CALL DLACPY( 'Full', NROWS, NRHS, B, LDB,
255 $ COPYB, LDB )
256 *
257 * Solve LS or overdetermined system
258 *
259 IF( M.GT.0 .AND. N.GT.0 ) THEN
260 CALL DLACPY( 'Full', M, N, COPYA, LDA,
261 $ A, LDA )
262 CALL DLACPY( 'Full', NROWS, NRHS,
263 $ COPYB, LDB, B, LDB )
264 END IF
265 SRNAMT = 'DGELS '
266 CALL DGELS( TRANS, M, N, NRHS, A, LDA, B,
267 $ LDB, WORK, LWORK, INFO )
268 IF( INFO.NE.0 )
269 $ CALL ALAERH( PATH, 'DGELS ', INFO, 0,
270 $ TRANS, M, N, NRHS, -1, NB,
271 $ ITYPE, NFAIL, NERRS,
272 $ NOUT )
273 *
274 * Check correctness of results
275 *
276 LDWORK = MAX( 1, NROWS )
277 IF( NROWS.GT.0 .AND. NRHS.GT.0 )
278 $ CALL DLACPY( 'Full', NROWS, NRHS,
279 $ COPYB, LDB, C, LDB )
280 CALL DQRT16( TRANS, M, N, NRHS, COPYA,
281 $ LDA, B, LDB, C, LDB, WORK,
282 $ RESULT( 1 ) )
283 *
284 IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
285 $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
286 *
287 * Solving LS system
288 *
289 RESULT( 2 ) = DQRT17( TRANS, 1, M, N,
290 $ NRHS, COPYA, LDA, B, LDB,
291 $ COPYB, LDB, C, WORK,
292 $ LWORK )
293 ELSE
294 *
295 * Solving overdetermined system
296 *
297 RESULT( 2 ) = DQRT14( TRANS, M, N,
298 $ NRHS, COPYA, LDA, B, LDB,
299 $ WORK, LWORK )
300 END IF
301 *
302 * Print information about the tests that
303 * did not pass the threshold.
304 *
305 DO 20 K = 1, 2
306 IF( RESULT( K ).GE.THRESH ) THEN
307 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
308 $ CALL ALAHD( NOUT, PATH )
309 WRITE( NOUT, FMT = 9999 )TRANS, M,
310 $ N, NRHS, NB, ITYPE, K,
311 $ RESULT( K )
312 NFAIL = NFAIL + 1
313 END IF
314 20 CONTINUE
315 NRUN = NRUN + 2
316 30 CONTINUE
317 40 CONTINUE
318 END IF
319 *
320 * Generate a matrix of scaling type ISCALE and rank
321 * type IRANK.
322 *
323 CALL DQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA,
324 $ COPYB, LDB, COPYS, RANK, NORMA, NORMB,
325 $ ISEED, WORK, LWORK )
326 *
327 * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
328 *
329 * Initialize vector IWORK.
330 *
331 DO 50 J = 1, N
332 IWORK( J ) = 0
333 50 CONTINUE
334 LDWORK = MAX( 1, M )
335 *
336 * Test DGELSX
337 *
338 * DGELSX: Compute the minimum-norm solution X
339 * to min( norm( A * X - B ) ) using a complete
340 * orthogonal factorization.
341 *
342 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
343 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B, LDB )
344 *
345 SRNAMT = 'DGELSX'
346 CALL DGELSX( M, N, NRHS, A, LDA, B, LDB, IWORK,
347 $ RCOND, CRANK, WORK, INFO )
348 IF( INFO.NE.0 )
349 $ CALL ALAERH( PATH, 'DGELSX', INFO, 0, ' ', M, N,
350 $ NRHS, -1, NB, ITYPE, NFAIL, NERRS,
351 $ NOUT )
352 *
353 * workspace used: MAX( MNMIN+3*N, 2*MNMIN+NRHS )
354 *
355 * Test 3: Compute relative error in svd
356 * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
357 *
358 RESULT( 3 ) = DQRT12( CRANK, CRANK, A, LDA, COPYS,
359 $ WORK, LWORK )
360 *
361 * Test 4: Compute error in solution
362 * workspace: M*NRHS + M
363 *
364 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
365 $ LDWORK )
366 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
367 $ LDA, B, LDB, WORK, LDWORK,
368 $ WORK( M*NRHS+1 ), RESULT( 4 ) )
369 *
370 * Test 5: Check norm of r'*A
371 * workspace: NRHS*(M+N)
372 *
373 RESULT( 5 ) = ZERO
374 IF( M.GT.CRANK )
375 $ RESULT( 5 ) = DQRT17( 'No transpose', 1, M, N,
376 $ NRHS, COPYA, LDA, B, LDB, COPYB,
377 $ LDB, C, WORK, LWORK )
378 *
379 * Test 6: Check if x is in the rowspace of A
380 * workspace: (M+NRHS)*(N+2)
381 *
382 RESULT( 6 ) = ZERO
383 *
384 IF( N.GT.CRANK )
385 $ RESULT( 6 ) = DQRT14( 'No transpose', M, N,
386 $ NRHS, COPYA, LDA, B, LDB, WORK,
387 $ LWORK )
388 *
389 * Print information about the tests that did not
390 * pass the threshold.
391 *
392 DO 60 K = 3, 6
393 IF( RESULT( K ).GE.THRESH ) THEN
394 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
395 $ CALL ALAHD( NOUT, PATH )
396 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
397 $ ITYPE, K, RESULT( K )
398 NFAIL = NFAIL + 1
399 END IF
400 60 CONTINUE
401 NRUN = NRUN + 4
402 *
403 * Loop for testing different block sizes.
404 *
405 DO 100 INB = 1, NNB
406 NB = NBVAL( INB )
407 CALL XLAENV( 1, NB )
408 CALL XLAENV( 3, NXVAL( INB ) )
409 *
410 * Test DGELSY
411 *
412 * DGELSY: Compute the minimum-norm solution X
413 * to min( norm( A * X - B ) )
414 * using the rank-revealing orthogonal
415 * factorization.
416 *
417 * Initialize vector IWORK.
418 *
419 DO 70 J = 1, N
420 IWORK( J ) = 0
421 70 CONTINUE
422 *
423 * Set LWLSY to the adequate value.
424 *
425 LWLSY = MAX( 1, MNMIN+2*N+NB*( N+1 ),
426 $ 2*MNMIN+NB*NRHS )
427 *
428 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
429 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B,
430 $ LDB )
431 *
432 SRNAMT = 'DGELSY'
433 CALL DGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
434 $ RCOND, CRANK, WORK, LWLSY, INFO )
435 IF( INFO.NE.0 )
436 $ CALL ALAERH( PATH, 'DGELSY', INFO, 0, ' ', M,
437 $ N, NRHS, -1, NB, ITYPE, NFAIL,
438 $ NERRS, NOUT )
439 *
440 * Test 7: Compute relative error in svd
441 * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
442 *
443 RESULT( 7 ) = DQRT12( CRANK, CRANK, A, LDA,
444 $ COPYS, WORK, LWORK )
445 *
446 * Test 8: Compute error in solution
447 * workspace: M*NRHS + M
448 *
449 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
450 $ LDWORK )
451 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
452 $ LDA, B, LDB, WORK, LDWORK,
453 $ WORK( M*NRHS+1 ), RESULT( 8 ) )
454 *
455 * Test 9: Check norm of r'*A
456 * workspace: NRHS*(M+N)
457 *
458 RESULT( 9 ) = ZERO
459 IF( M.GT.CRANK )
460 $ RESULT( 9 ) = DQRT17( 'No transpose', 1, M,
461 $ N, NRHS, COPYA, LDA, B, LDB,
462 $ COPYB, LDB, C, WORK, LWORK )
463 *
464 * Test 10: Check if x is in the rowspace of A
465 * workspace: (M+NRHS)*(N+2)
466 *
467 RESULT( 10 ) = ZERO
468 *
469 IF( N.GT.CRANK )
470 $ RESULT( 10 ) = DQRT14( 'No transpose', M, N,
471 $ NRHS, COPYA, LDA, B, LDB,
472 $ WORK, LWORK )
473 *
474 * Test DGELSS
475 *
476 * DGELSS: Compute the minimum-norm solution X
477 * to min( norm( A * X - B ) )
478 * using the SVD.
479 *
480 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
481 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B,
482 $ LDB )
483 SRNAMT = 'DGELSS'
484 CALL DGELSS( M, N, NRHS, A, LDA, B, LDB, S,
485 $ RCOND, CRANK, WORK, LWORK, INFO )
486 IF( INFO.NE.0 )
487 $ CALL ALAERH( PATH, 'DGELSS', INFO, 0, ' ', M,
488 $ N, NRHS, -1, NB, ITYPE, NFAIL,
489 $ NERRS, NOUT )
490 *
491 * workspace used: 3*min(m,n) +
492 * max(2*min(m,n),nrhs,max(m,n))
493 *
494 * Test 11: Compute relative error in svd
495 *
496 IF( RANK.GT.0 ) THEN
497 CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
498 RESULT( 11 ) = DASUM( MNMIN, S, 1 ) /
499 $ DASUM( MNMIN, COPYS, 1 ) /
500 $ ( EPS*DBLE( MNMIN ) )
501 ELSE
502 RESULT( 11 ) = ZERO
503 END IF
504 *
505 * Test 12: Compute error in solution
506 *
507 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
508 $ LDWORK )
509 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
510 $ LDA, B, LDB, WORK, LDWORK,
511 $ WORK( M*NRHS+1 ), RESULT( 12 ) )
512 *
513 * Test 13: Check norm of r'*A
514 *
515 RESULT( 13 ) = ZERO
516 IF( M.GT.CRANK )
517 $ RESULT( 13 ) = DQRT17( 'No transpose', 1, M,
518 $ N, NRHS, COPYA, LDA, B, LDB,
519 $ COPYB, LDB, C, WORK, LWORK )
520 *
521 * Test 14: Check if x is in the rowspace of A
522 *
523 RESULT( 14 ) = ZERO
524 IF( N.GT.CRANK )
525 $ RESULT( 14 ) = DQRT14( 'No transpose', M, N,
526 $ NRHS, COPYA, LDA, B, LDB,
527 $ WORK, LWORK )
528 *
529 * Test DGELSD
530 *
531 * DGELSD: Compute the minimum-norm solution X
532 * to min( norm( A * X - B ) ) using a
533 * divide and conquer SVD.
534 *
535 * Initialize vector IWORK.
536 *
537 DO 80 J = 1, N
538 IWORK( J ) = 0
539 80 CONTINUE
540 *
541 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
542 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B,
543 $ LDB )
544 *
545 SRNAMT = 'DGELSD'
546 CALL DGELSD( M, N, NRHS, A, LDA, B, LDB, S,
547 $ RCOND, CRANK, WORK, LWORK, IWORK,
548 $ INFO )
549 IF( INFO.NE.0 )
550 $ CALL ALAERH( PATH, 'DGELSD', INFO, 0, ' ', M,
551 $ N, NRHS, -1, NB, ITYPE, NFAIL,
552 $ NERRS, NOUT )
553 *
554 * Test 15: Compute relative error in svd
555 *
556 IF( RANK.GT.0 ) THEN
557 CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
558 RESULT( 15 ) = DASUM( MNMIN, S, 1 ) /
559 $ DASUM( MNMIN, COPYS, 1 ) /
560 $ ( EPS*DBLE( MNMIN ) )
561 ELSE
562 RESULT( 15 ) = ZERO
563 END IF
564 *
565 * Test 16: Compute error in solution
566 *
567 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
568 $ LDWORK )
569 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
570 $ LDA, B, LDB, WORK, LDWORK,
571 $ WORK( M*NRHS+1 ), RESULT( 16 ) )
572 *
573 * Test 17: Check norm of r'*A
574 *
575 RESULT( 17 ) = ZERO
576 IF( M.GT.CRANK )
577 $ RESULT( 17 ) = DQRT17( 'No transpose', 1, M,
578 $ N, NRHS, COPYA, LDA, B, LDB,
579 $ COPYB, LDB, C, WORK, LWORK )
580 *
581 * Test 18: Check if x is in the rowspace of A
582 *
583 RESULT( 18 ) = ZERO
584 IF( N.GT.CRANK )
585 $ RESULT( 18 ) = DQRT14( 'No transpose', M, N,
586 $ NRHS, COPYA, LDA, B, LDB,
587 $ WORK, LWORK )
588 *
589 * Print information about the tests that did not
590 * pass the threshold.
591 *
592 DO 90 K = 7, NTESTS
593 IF( RESULT( K ).GE.THRESH ) THEN
594 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
595 $ CALL ALAHD( NOUT, PATH )
596 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
597 $ ITYPE, K, RESULT( K )
598 NFAIL = NFAIL + 1
599 END IF
600 90 CONTINUE
601 NRUN = NRUN + 12
602 *
603 100 CONTINUE
604 110 CONTINUE
605 120 CONTINUE
606 130 CONTINUE
607 140 CONTINUE
608 150 CONTINUE
609 *
610 * Print a summary of the results.
611 *
612 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
613 *
614 9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4,
615 $ ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 )
616 9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4,
617 $ ', type', I2, ', test(', I2, ')=', G12.5 )
618 RETURN
619 *
620 * End of DDRVLS
621 *
622 END
2 $ NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
3 $ COPYB, C, S, COPYS, WORK, IWORK, 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 LOGICAL TSTERR
11 INTEGER NM, NN, NNB, NNS, NOUT
12 DOUBLE PRECISION THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL DOTYPE( * )
16 INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
17 $ NVAL( * ), NXVAL( * )
18 DOUBLE PRECISION A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
19 $ COPYS( * ), S( * ), WORK( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DDRVLS tests the least squares driver routines DGELS, DGELSS, DGELSX,
26 * DGELSY and DGELSD.
27 *
28 * Arguments
29 * =========
30 *
31 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
32 * The matrix types to be used for testing. Matrices of type j
33 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
34 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
35 * The matrix of type j is generated as follows:
36 * j=1: A = U*D*V where U and V are random orthogonal matrices
37 * and D has random entries (> 0.1) taken from a uniform
38 * distribution (0,1). A is full rank.
39 * j=2: The same of 1, but A is scaled up.
40 * j=3: The same of 1, but A is scaled down.
41 * j=4: A = U*D*V where U and V are random orthogonal matrices
42 * and D has 3*min(M,N)/4 random entries (> 0.1) taken
43 * from a uniform distribution (0,1) and the remaining
44 * entries set to 0. A is rank-deficient.
45 * j=5: The same of 4, but A is scaled up.
46 * j=6: The same of 5, but A is scaled down.
47 *
48 * NM (input) INTEGER
49 * The number of values of M contained in the vector MVAL.
50 *
51 * MVAL (input) INTEGER array, dimension (NM)
52 * The values of the matrix row dimension M.
53 *
54 * NN (input) INTEGER
55 * The number of values of N contained in the vector NVAL.
56 *
57 * NVAL (input) INTEGER array, dimension (NN)
58 * The values of the matrix column dimension N.
59 *
60 * NNS (input) INTEGER
61 * The number of values of NRHS contained in the vector NSVAL.
62 *
63 * NSVAL (input) INTEGER array, dimension (NNS)
64 * The values of the number of right hand sides NRHS.
65 *
66 * NNB (input) INTEGER
67 * The number of values of NB and NX contained in the
68 * vectors NBVAL and NXVAL. The blocking parameters are used
69 * in pairs (NB,NX).
70 *
71 * NBVAL (input) INTEGER array, dimension (NNB)
72 * The values of the blocksize NB.
73 *
74 * NXVAL (input) INTEGER array, dimension (NNB)
75 * The values of the crossover point NX.
76 *
77 * THRESH (input) DOUBLE PRECISION
78 * The threshold value for the test ratios. A result is
79 * included in the output file if RESULT >= THRESH. To have
80 * every test ratio printed, use THRESH = 0.
81 *
82 * TSTERR (input) LOGICAL
83 * Flag that indicates whether error exits are to be tested.
84 *
85 * A (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
86 * where MMAX is the maximum value of M in MVAL and NMAX is the
87 * maximum value of N in NVAL.
88 *
89 * COPYA (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
90 *
91 * B (workspace) DOUBLE PRECISION array, dimension (MMAX*NSMAX)
92 * where MMAX is the maximum value of M in MVAL and NSMAX is the
93 * maximum value of NRHS in NSVAL.
94 *
95 * COPYB (workspace) DOUBLE PRECISION array, dimension (MMAX*NSMAX)
96 *
97 * C (workspace) DOUBLE PRECISION array, dimension (MMAX*NSMAX)
98 *
99 * S (workspace) DOUBLE PRECISION array, dimension
100 * (min(MMAX,NMAX))
101 *
102 * COPYS (workspace) DOUBLE PRECISION array, dimension
103 * (min(MMAX,NMAX))
104 *
105 * WORK (workspace) DOUBLE PRECISION array,
106 * dimension (MMAX*NMAX + 4*NMAX + MMAX).
107 *
108 * IWORK (workspace) INTEGER array, dimension (15*NMAX)
109 *
110 * NOUT (input) INTEGER
111 * The unit number for output.
112 *
113 * =====================================================================
114 *
115 * .. Parameters ..
116 INTEGER NTESTS
117 PARAMETER ( NTESTS = 18 )
118 INTEGER SMLSIZ
119 PARAMETER ( SMLSIZ = 25 )
120 DOUBLE PRECISION ONE, TWO, ZERO
121 PARAMETER ( ONE = 1.0D0, TWO = 2.0D0, ZERO = 0.0D0 )
122 * ..
123 * .. Local Scalars ..
124 CHARACTER TRANS
125 CHARACTER*3 PATH
126 INTEGER CRANK, I, IM, IN, INB, INFO, INS, IRANK,
127 $ ISCALE, ITRAN, ITYPE, J, K, LDA, LDB, LDWORK,
128 $ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS,
129 $ NFAIL, NLVL, NRHS, NROWS, NRUN, RANK
130 DOUBLE PRECISION EPS, NORMA, NORMB, RCOND
131 * ..
132 * .. Local Arrays ..
133 INTEGER ISEED( 4 ), ISEEDY( 4 )
134 DOUBLE PRECISION RESULT( NTESTS )
135 * ..
136 * .. External Functions ..
137 DOUBLE PRECISION DASUM, DLAMCH, DQRT12, DQRT14, DQRT17
138 EXTERNAL DASUM, DLAMCH, DQRT12, DQRT14, DQRT17
139 * ..
140 * .. External Subroutines ..
141 EXTERNAL ALAERH, ALAHD, ALASVM, DAXPY, DERRLS, DGELS,
142 $ DGELSD, DGELSS, DGELSX, DGELSY, DGEMM, DLACPY,
143 $ DLARNV, DLASRT, DQRT13, DQRT15, DQRT16, DSCAL,
144 $ XLAENV
145 * ..
146 * .. Intrinsic Functions ..
147 INTRINSIC DBLE, INT, LOG, MAX, MIN, SQRT
148 * ..
149 * .. Scalars in Common ..
150 LOGICAL LERR, OK
151 CHARACTER*32 SRNAMT
152 INTEGER INFOT, IOUNIT
153 * ..
154 * .. Common blocks ..
155 COMMON / INFOC / INFOT, IOUNIT, OK, LERR
156 COMMON / SRNAMC / SRNAMT
157 * ..
158 * .. Data statements ..
159 DATA ISEEDY / 1988, 1989, 1990, 1991 /
160 * ..
161 * .. Executable Statements ..
162 *
163 * Initialize constants and the random number seed.
164 *
165 PATH( 1: 1 ) = 'Double precision'
166 PATH( 2: 3 ) = 'LS'
167 NRUN = 0
168 NFAIL = 0
169 NERRS = 0
170 DO 10 I = 1, 4
171 ISEED( I ) = ISEEDY( I )
172 10 CONTINUE
173 EPS = DLAMCH( 'Epsilon' )
174 *
175 * Threshold for rank estimation
176 *
177 RCOND = SQRT( EPS ) - ( SQRT( EPS )-EPS ) / 2
178 *
179 * Test the error exits
180 *
181 CALL XLAENV( 2, 2 )
182 CALL XLAENV( 9, SMLSIZ )
183 IF( TSTERR )
184 $ CALL DERRLS( PATH, NOUT )
185 *
186 * Print the header if NM = 0 or NN = 0 and THRESH = 0.
187 *
188 IF( ( NM.EQ.0 .OR. NN.EQ.0 ) .AND. THRESH.EQ.ZERO )
189 $ CALL ALAHD( NOUT, PATH )
190 INFOT = 0
191 CALL XLAENV( 2, 2 )
192 CALL XLAENV( 9, SMLSIZ )
193 *
194 DO 150 IM = 1, NM
195 M = MVAL( IM )
196 LDA = MAX( 1, M )
197 *
198 DO 140 IN = 1, NN
199 N = NVAL( IN )
200 MNMIN = MIN( M, N )
201 LDB = MAX( 1, M, N )
202 *
203 DO 130 INS = 1, NNS
204 NRHS = NSVAL( INS )
205 NLVL = MAX( INT( LOG( MAX( ONE, DBLE( MNMIN ) ) /
206 $ DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1, 0 )
207 LWORK = MAX( 1, ( M+NRHS )*( N+2 ), ( N+NRHS )*( M+2 ),
208 $ M*N+4*MNMIN+MAX( M, N ), 12*MNMIN+2*MNMIN*SMLSIZ+
209 $ 8*MNMIN*NLVL+MNMIN*NRHS+(SMLSIZ+1)**2 )
210 *
211 DO 120 IRANK = 1, 2
212 DO 110 ISCALE = 1, 3
213 ITYPE = ( IRANK-1 )*3 + ISCALE
214 IF( .NOT.DOTYPE( ITYPE ) )
215 $ GO TO 110
216 *
217 IF( IRANK.EQ.1 ) THEN
218 *
219 * Test DGELS
220 *
221 * Generate a matrix of scaling type ISCALE
222 *
223 CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
224 $ ISEED )
225 DO 40 INB = 1, NNB
226 NB = NBVAL( INB )
227 CALL XLAENV( 1, NB )
228 CALL XLAENV( 3, NXVAL( INB ) )
229 *
230 DO 30 ITRAN = 1, 2
231 IF( ITRAN.EQ.1 ) THEN
232 TRANS = 'N'
233 NROWS = M
234 NCOLS = N
235 ELSE
236 TRANS = 'T'
237 NROWS = N
238 NCOLS = M
239 END IF
240 LDWORK = MAX( 1, NCOLS )
241 *
242 * Set up a consistent rhs
243 *
244 IF( NCOLS.GT.0 ) THEN
245 CALL DLARNV( 2, ISEED, NCOLS*NRHS,
246 $ WORK )
247 CALL DSCAL( NCOLS*NRHS,
248 $ ONE / DBLE( NCOLS ), WORK,
249 $ 1 )
250 END IF
251 CALL DGEMM( TRANS, 'No transpose', NROWS,
252 $ NRHS, NCOLS, ONE, COPYA, LDA,
253 $ WORK, LDWORK, ZERO, B, LDB )
254 CALL DLACPY( 'Full', NROWS, NRHS, B, LDB,
255 $ COPYB, LDB )
256 *
257 * Solve LS or overdetermined system
258 *
259 IF( M.GT.0 .AND. N.GT.0 ) THEN
260 CALL DLACPY( 'Full', M, N, COPYA, LDA,
261 $ A, LDA )
262 CALL DLACPY( 'Full', NROWS, NRHS,
263 $ COPYB, LDB, B, LDB )
264 END IF
265 SRNAMT = 'DGELS '
266 CALL DGELS( TRANS, M, N, NRHS, A, LDA, B,
267 $ LDB, WORK, LWORK, INFO )
268 IF( INFO.NE.0 )
269 $ CALL ALAERH( PATH, 'DGELS ', INFO, 0,
270 $ TRANS, M, N, NRHS, -1, NB,
271 $ ITYPE, NFAIL, NERRS,
272 $ NOUT )
273 *
274 * Check correctness of results
275 *
276 LDWORK = MAX( 1, NROWS )
277 IF( NROWS.GT.0 .AND. NRHS.GT.0 )
278 $ CALL DLACPY( 'Full', NROWS, NRHS,
279 $ COPYB, LDB, C, LDB )
280 CALL DQRT16( TRANS, M, N, NRHS, COPYA,
281 $ LDA, B, LDB, C, LDB, WORK,
282 $ RESULT( 1 ) )
283 *
284 IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
285 $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
286 *
287 * Solving LS system
288 *
289 RESULT( 2 ) = DQRT17( TRANS, 1, M, N,
290 $ NRHS, COPYA, LDA, B, LDB,
291 $ COPYB, LDB, C, WORK,
292 $ LWORK )
293 ELSE
294 *
295 * Solving overdetermined system
296 *
297 RESULT( 2 ) = DQRT14( TRANS, M, N,
298 $ NRHS, COPYA, LDA, B, LDB,
299 $ WORK, LWORK )
300 END IF
301 *
302 * Print information about the tests that
303 * did not pass the threshold.
304 *
305 DO 20 K = 1, 2
306 IF( RESULT( K ).GE.THRESH ) THEN
307 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
308 $ CALL ALAHD( NOUT, PATH )
309 WRITE( NOUT, FMT = 9999 )TRANS, M,
310 $ N, NRHS, NB, ITYPE, K,
311 $ RESULT( K )
312 NFAIL = NFAIL + 1
313 END IF
314 20 CONTINUE
315 NRUN = NRUN + 2
316 30 CONTINUE
317 40 CONTINUE
318 END IF
319 *
320 * Generate a matrix of scaling type ISCALE and rank
321 * type IRANK.
322 *
323 CALL DQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA,
324 $ COPYB, LDB, COPYS, RANK, NORMA, NORMB,
325 $ ISEED, WORK, LWORK )
326 *
327 * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
328 *
329 * Initialize vector IWORK.
330 *
331 DO 50 J = 1, N
332 IWORK( J ) = 0
333 50 CONTINUE
334 LDWORK = MAX( 1, M )
335 *
336 * Test DGELSX
337 *
338 * DGELSX: Compute the minimum-norm solution X
339 * to min( norm( A * X - B ) ) using a complete
340 * orthogonal factorization.
341 *
342 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
343 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B, LDB )
344 *
345 SRNAMT = 'DGELSX'
346 CALL DGELSX( M, N, NRHS, A, LDA, B, LDB, IWORK,
347 $ RCOND, CRANK, WORK, INFO )
348 IF( INFO.NE.0 )
349 $ CALL ALAERH( PATH, 'DGELSX', INFO, 0, ' ', M, N,
350 $ NRHS, -1, NB, ITYPE, NFAIL, NERRS,
351 $ NOUT )
352 *
353 * workspace used: MAX( MNMIN+3*N, 2*MNMIN+NRHS )
354 *
355 * Test 3: Compute relative error in svd
356 * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
357 *
358 RESULT( 3 ) = DQRT12( CRANK, CRANK, A, LDA, COPYS,
359 $ WORK, LWORK )
360 *
361 * Test 4: Compute error in solution
362 * workspace: M*NRHS + M
363 *
364 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
365 $ LDWORK )
366 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
367 $ LDA, B, LDB, WORK, LDWORK,
368 $ WORK( M*NRHS+1 ), RESULT( 4 ) )
369 *
370 * Test 5: Check norm of r'*A
371 * workspace: NRHS*(M+N)
372 *
373 RESULT( 5 ) = ZERO
374 IF( M.GT.CRANK )
375 $ RESULT( 5 ) = DQRT17( 'No transpose', 1, M, N,
376 $ NRHS, COPYA, LDA, B, LDB, COPYB,
377 $ LDB, C, WORK, LWORK )
378 *
379 * Test 6: Check if x is in the rowspace of A
380 * workspace: (M+NRHS)*(N+2)
381 *
382 RESULT( 6 ) = ZERO
383 *
384 IF( N.GT.CRANK )
385 $ RESULT( 6 ) = DQRT14( 'No transpose', M, N,
386 $ NRHS, COPYA, LDA, B, LDB, WORK,
387 $ LWORK )
388 *
389 * Print information about the tests that did not
390 * pass the threshold.
391 *
392 DO 60 K = 3, 6
393 IF( RESULT( K ).GE.THRESH ) THEN
394 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
395 $ CALL ALAHD( NOUT, PATH )
396 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
397 $ ITYPE, K, RESULT( K )
398 NFAIL = NFAIL + 1
399 END IF
400 60 CONTINUE
401 NRUN = NRUN + 4
402 *
403 * Loop for testing different block sizes.
404 *
405 DO 100 INB = 1, NNB
406 NB = NBVAL( INB )
407 CALL XLAENV( 1, NB )
408 CALL XLAENV( 3, NXVAL( INB ) )
409 *
410 * Test DGELSY
411 *
412 * DGELSY: Compute the minimum-norm solution X
413 * to min( norm( A * X - B ) )
414 * using the rank-revealing orthogonal
415 * factorization.
416 *
417 * Initialize vector IWORK.
418 *
419 DO 70 J = 1, N
420 IWORK( J ) = 0
421 70 CONTINUE
422 *
423 * Set LWLSY to the adequate value.
424 *
425 LWLSY = MAX( 1, MNMIN+2*N+NB*( N+1 ),
426 $ 2*MNMIN+NB*NRHS )
427 *
428 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
429 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B,
430 $ LDB )
431 *
432 SRNAMT = 'DGELSY'
433 CALL DGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
434 $ RCOND, CRANK, WORK, LWLSY, INFO )
435 IF( INFO.NE.0 )
436 $ CALL ALAERH( PATH, 'DGELSY', INFO, 0, ' ', M,
437 $ N, NRHS, -1, NB, ITYPE, NFAIL,
438 $ NERRS, NOUT )
439 *
440 * Test 7: Compute relative error in svd
441 * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
442 *
443 RESULT( 7 ) = DQRT12( CRANK, CRANK, A, LDA,
444 $ COPYS, WORK, LWORK )
445 *
446 * Test 8: Compute error in solution
447 * workspace: M*NRHS + M
448 *
449 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
450 $ LDWORK )
451 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
452 $ LDA, B, LDB, WORK, LDWORK,
453 $ WORK( M*NRHS+1 ), RESULT( 8 ) )
454 *
455 * Test 9: Check norm of r'*A
456 * workspace: NRHS*(M+N)
457 *
458 RESULT( 9 ) = ZERO
459 IF( M.GT.CRANK )
460 $ RESULT( 9 ) = DQRT17( 'No transpose', 1, M,
461 $ N, NRHS, COPYA, LDA, B, LDB,
462 $ COPYB, LDB, C, WORK, LWORK )
463 *
464 * Test 10: Check if x is in the rowspace of A
465 * workspace: (M+NRHS)*(N+2)
466 *
467 RESULT( 10 ) = ZERO
468 *
469 IF( N.GT.CRANK )
470 $ RESULT( 10 ) = DQRT14( 'No transpose', M, N,
471 $ NRHS, COPYA, LDA, B, LDB,
472 $ WORK, LWORK )
473 *
474 * Test DGELSS
475 *
476 * DGELSS: Compute the minimum-norm solution X
477 * to min( norm( A * X - B ) )
478 * using the SVD.
479 *
480 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
481 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B,
482 $ LDB )
483 SRNAMT = 'DGELSS'
484 CALL DGELSS( M, N, NRHS, A, LDA, B, LDB, S,
485 $ RCOND, CRANK, WORK, LWORK, INFO )
486 IF( INFO.NE.0 )
487 $ CALL ALAERH( PATH, 'DGELSS', INFO, 0, ' ', M,
488 $ N, NRHS, -1, NB, ITYPE, NFAIL,
489 $ NERRS, NOUT )
490 *
491 * workspace used: 3*min(m,n) +
492 * max(2*min(m,n),nrhs,max(m,n))
493 *
494 * Test 11: Compute relative error in svd
495 *
496 IF( RANK.GT.0 ) THEN
497 CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
498 RESULT( 11 ) = DASUM( MNMIN, S, 1 ) /
499 $ DASUM( MNMIN, COPYS, 1 ) /
500 $ ( EPS*DBLE( MNMIN ) )
501 ELSE
502 RESULT( 11 ) = ZERO
503 END IF
504 *
505 * Test 12: Compute error in solution
506 *
507 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
508 $ LDWORK )
509 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
510 $ LDA, B, LDB, WORK, LDWORK,
511 $ WORK( M*NRHS+1 ), RESULT( 12 ) )
512 *
513 * Test 13: Check norm of r'*A
514 *
515 RESULT( 13 ) = ZERO
516 IF( M.GT.CRANK )
517 $ RESULT( 13 ) = DQRT17( 'No transpose', 1, M,
518 $ N, NRHS, COPYA, LDA, B, LDB,
519 $ COPYB, LDB, C, WORK, LWORK )
520 *
521 * Test 14: Check if x is in the rowspace of A
522 *
523 RESULT( 14 ) = ZERO
524 IF( N.GT.CRANK )
525 $ RESULT( 14 ) = DQRT14( 'No transpose', M, N,
526 $ NRHS, COPYA, LDA, B, LDB,
527 $ WORK, LWORK )
528 *
529 * Test DGELSD
530 *
531 * DGELSD: Compute the minimum-norm solution X
532 * to min( norm( A * X - B ) ) using a
533 * divide and conquer SVD.
534 *
535 * Initialize vector IWORK.
536 *
537 DO 80 J = 1, N
538 IWORK( J ) = 0
539 80 CONTINUE
540 *
541 CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
542 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B,
543 $ LDB )
544 *
545 SRNAMT = 'DGELSD'
546 CALL DGELSD( M, N, NRHS, A, LDA, B, LDB, S,
547 $ RCOND, CRANK, WORK, LWORK, IWORK,
548 $ INFO )
549 IF( INFO.NE.0 )
550 $ CALL ALAERH( PATH, 'DGELSD', INFO, 0, ' ', M,
551 $ N, NRHS, -1, NB, ITYPE, NFAIL,
552 $ NERRS, NOUT )
553 *
554 * Test 15: Compute relative error in svd
555 *
556 IF( RANK.GT.0 ) THEN
557 CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
558 RESULT( 15 ) = DASUM( MNMIN, S, 1 ) /
559 $ DASUM( MNMIN, COPYS, 1 ) /
560 $ ( EPS*DBLE( MNMIN ) )
561 ELSE
562 RESULT( 15 ) = ZERO
563 END IF
564 *
565 * Test 16: Compute error in solution
566 *
567 CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
568 $ LDWORK )
569 CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
570 $ LDA, B, LDB, WORK, LDWORK,
571 $ WORK( M*NRHS+1 ), RESULT( 16 ) )
572 *
573 * Test 17: Check norm of r'*A
574 *
575 RESULT( 17 ) = ZERO
576 IF( M.GT.CRANK )
577 $ RESULT( 17 ) = DQRT17( 'No transpose', 1, M,
578 $ N, NRHS, COPYA, LDA, B, LDB,
579 $ COPYB, LDB, C, WORK, LWORK )
580 *
581 * Test 18: Check if x is in the rowspace of A
582 *
583 RESULT( 18 ) = ZERO
584 IF( N.GT.CRANK )
585 $ RESULT( 18 ) = DQRT14( 'No transpose', M, N,
586 $ NRHS, COPYA, LDA, B, LDB,
587 $ WORK, LWORK )
588 *
589 * Print information about the tests that did not
590 * pass the threshold.
591 *
592 DO 90 K = 7, NTESTS
593 IF( RESULT( K ).GE.THRESH ) THEN
594 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
595 $ CALL ALAHD( NOUT, PATH )
596 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
597 $ ITYPE, K, RESULT( K )
598 NFAIL = NFAIL + 1
599 END IF
600 90 CONTINUE
601 NRUN = NRUN + 12
602 *
603 100 CONTINUE
604 110 CONTINUE
605 120 CONTINUE
606 130 CONTINUE
607 140 CONTINUE
608 150 CONTINUE
609 *
610 * Print a summary of the results.
611 *
612 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
613 *
614 9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4,
615 $ ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 )
616 9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4,
617 $ ', type', I2, ', test(', I2, ')=', G12.5 )
618 RETURN
619 *
620 * End of DDRVLS
621 *
622 END