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