1 SUBROUTINE SDRVLS( 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 REAL THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL DOTYPE( * )
16 INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
17 $ NVAL( * ), NXVAL( * )
18 REAL A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
19 $ COPYS( * ), S( * ), WORK( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * SDRVLS tests the least squares driver routines SGELS, SGELSS, SGELSX,
26 * SGELSY and SGELSD.
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) REAL
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) REAL 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) REAL array, dimension (MMAX*NMAX)
90 *
91 * B (workspace) REAL 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) REAL array, dimension (MMAX*NSMAX)
96 *
97 * C (workspace) REAL array, dimension (MMAX*NSMAX)
98 *
99 * S (workspace) REAL array, dimension
100 * (min(MMAX,NMAX))
101 *
102 * COPYS (workspace) REAL array, dimension
103 * (min(MMAX,NMAX))
104 *
105 * WORK (workspace) REAL 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 REAL ONE, TWO, ZERO
121 PARAMETER ( ONE = 1.0E0, TWO = 2.0E0, ZERO = 0.0E0 )
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 REAL EPS, NORMA, NORMB, RCOND
131 * ..
132 * .. Local Arrays ..
133 INTEGER ISEED( 4 ), ISEEDY( 4 )
134 REAL RESULT( NTESTS )
135 * ..
136 * .. External Functions ..
137 REAL SASUM, SLAMCH, SQRT12, SQRT14, SQRT17
138 EXTERNAL SASUM, SLAMCH, SQRT12, SQRT14, SQRT17
139 * ..
140 * .. External Subroutines ..
141 EXTERNAL ALAERH, ALAHD, ALASVM, SAXPY, SERRLS, SGELS,
142 $ SGELSD, SGELSS, SGELSX, SGELSY, SGEMM, SLACPY,
143 $ SLARNV, SQRT13, SQRT15, SQRT16, SSCAL,
144 $ XLAENV
145 * ..
146 * .. Intrinsic Functions ..
147 INTRINSIC INT, LOG, MAX, MIN, REAL, 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 ) = 'Single 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 = SLAMCH( '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 SERRLS( 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 *
192 DO 150 IM = 1, NM
193 M = MVAL( IM )
194 LDA = MAX( 1, M )
195 *
196 DO 140 IN = 1, NN
197 N = NVAL( IN )
198 MNMIN = MIN( M, N )
199 LDB = MAX( 1, M, N )
200 *
201 DO 130 INS = 1, NNS
202 NRHS = NSVAL( INS )
203 NLVL = MAX( INT( LOG( MAX( ONE, REAL( MNMIN ) ) /
204 $ REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1, 0 )
205 LWORK = MAX( 1, ( M+NRHS )*( N+2 ), ( N+NRHS )*( M+2 ),
206 $ M*N+4*MNMIN+MAX( M, N ), 12*MNMIN+2*MNMIN*SMLSIZ+
207 $ 8*MNMIN*NLVL+MNMIN*NRHS+(SMLSIZ+1)**2 )
208 *
209 DO 120 IRANK = 1, 2
210 DO 110 ISCALE = 1, 3
211 ITYPE = ( IRANK-1 )*3 + ISCALE
212 IF( .NOT.DOTYPE( ITYPE ) )
213 $ GO TO 110
214 *
215 IF( IRANK.EQ.1 ) THEN
216 *
217 * Test SGELS
218 *
219 * Generate a matrix of scaling type ISCALE
220 *
221 CALL SQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
222 $ ISEED )
223 DO 40 INB = 1, NNB
224 NB = NBVAL( INB )
225 CALL XLAENV( 1, NB )
226 CALL XLAENV( 3, NXVAL( INB ) )
227 *
228 DO 30 ITRAN = 1, 2
229 IF( ITRAN.EQ.1 ) THEN
230 TRANS = 'N'
231 NROWS = M
232 NCOLS = N
233 ELSE
234 TRANS = 'T'
235 NROWS = N
236 NCOLS = M
237 END IF
238 LDWORK = MAX( 1, NCOLS )
239 *
240 * Set up a consistent rhs
241 *
242 IF( NCOLS.GT.0 ) THEN
243 CALL SLARNV( 2, ISEED, NCOLS*NRHS,
244 $ WORK )
245 CALL SSCAL( NCOLS*NRHS,
246 $ ONE / REAL( NCOLS ), WORK,
247 $ 1 )
248 END IF
249 CALL SGEMM( TRANS, 'No transpose', NROWS,
250 $ NRHS, NCOLS, ONE, COPYA, LDA,
251 $ WORK, LDWORK, ZERO, B, LDB )
252 CALL SLACPY( 'Full', NROWS, NRHS, B, LDB,
253 $ COPYB, LDB )
254 *
255 * Solve LS or overdetermined system
256 *
257 IF( M.GT.0 .AND. N.GT.0 ) THEN
258 CALL SLACPY( 'Full', M, N, COPYA, LDA,
259 $ A, LDA )
260 CALL SLACPY( 'Full', NROWS, NRHS,
261 $ COPYB, LDB, B, LDB )
262 END IF
263 SRNAMT = 'SGELS '
264 CALL SGELS( TRANS, M, N, NRHS, A, LDA, B,
265 $ LDB, WORK, LWORK, INFO )
266 IF( INFO.NE.0 )
267 $ CALL ALAERH( PATH, 'SGELS ', INFO, 0,
268 $ TRANS, M, N, NRHS, -1, NB,
269 $ ITYPE, NFAIL, NERRS,
270 $ NOUT )
271 *
272 * Check correctness of results
273 *
274 LDWORK = MAX( 1, NROWS )
275 IF( NROWS.GT.0 .AND. NRHS.GT.0 )
276 $ CALL SLACPY( 'Full', NROWS, NRHS,
277 $ COPYB, LDB, C, LDB )
278 CALL SQRT16( TRANS, M, N, NRHS, COPYA,
279 $ LDA, B, LDB, C, LDB, WORK,
280 $ RESULT( 1 ) )
281 *
282 IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
283 $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
284 *
285 * Solving LS system
286 *
287 RESULT( 2 ) = SQRT17( TRANS, 1, M, N,
288 $ NRHS, COPYA, LDA, B, LDB,
289 $ COPYB, LDB, C, WORK,
290 $ LWORK )
291 ELSE
292 *
293 * Solving overdetermined system
294 *
295 RESULT( 2 ) = SQRT14( TRANS, M, N,
296 $ NRHS, COPYA, LDA, B, LDB,
297 $ WORK, LWORK )
298 END IF
299 *
300 * Print information about the tests that
301 * did not pass the threshold.
302 *
303 DO 20 K = 1, 2
304 IF( RESULT( K ).GE.THRESH ) THEN
305 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
306 $ CALL ALAHD( NOUT, PATH )
307 WRITE( NOUT, FMT = 9999 )TRANS, M,
308 $ N, NRHS, NB, ITYPE, K,
309 $ RESULT( K )
310 NFAIL = NFAIL + 1
311 END IF
312 20 CONTINUE
313 NRUN = NRUN + 2
314 30 CONTINUE
315 40 CONTINUE
316 END IF
317 *
318 * Generate a matrix of scaling type ISCALE and rank
319 * type IRANK.
320 *
321 CALL SQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA,
322 $ COPYB, LDB, COPYS, RANK, NORMA, NORMB,
323 $ ISEED, WORK, LWORK )
324 *
325 * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
326 *
327 * Initialize vector IWORK.
328 *
329 DO 50 J = 1, N
330 IWORK( J ) = 0
331 50 CONTINUE
332 LDWORK = MAX( 1, M )
333 *
334 * Test SGELSX
335 *
336 * SGELSX: Compute the minimum-norm solution X
337 * to min( norm( A * X - B ) ) using a complete
338 * orthogonal factorization.
339 *
340 CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
341 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B, LDB )
342 *
343 SRNAMT = 'SGELSX'
344 CALL SGELSX( M, N, NRHS, A, LDA, B, LDB, IWORK,
345 $ RCOND, CRANK, WORK, INFO )
346 IF( INFO.NE.0 )
347 $ CALL ALAERH( PATH, 'SGELSX', INFO, 0, ' ', M, N,
348 $ NRHS, -1, NB, ITYPE, NFAIL, NERRS,
349 $ NOUT )
350 *
351 * workspace used: MAX( MNMIN+3*N, 2*MNMIN+NRHS )
352 *
353 * Test 3: Compute relative error in svd
354 * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
355 *
356 RESULT( 3 ) = SQRT12( CRANK, CRANK, A, LDA, COPYS,
357 $ WORK, LWORK )
358 *
359 * Test 4: Compute error in solution
360 * workspace: M*NRHS + M
361 *
362 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
363 $ LDWORK )
364 CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
365 $ LDA, B, LDB, WORK, LDWORK,
366 $ WORK( M*NRHS+1 ), RESULT( 4 ) )
367 *
368 * Test 5: Check norm of r'*A
369 * workspace: NRHS*(M+N)
370 *
371 RESULT( 5 ) = ZERO
372 IF( M.GT.CRANK )
373 $ RESULT( 5 ) = SQRT17( 'No transpose', 1, M, N,
374 $ NRHS, COPYA, LDA, B, LDB, COPYB,
375 $ LDB, C, WORK, LWORK )
376 *
377 * Test 6: Check if x is in the rowspace of A
378 * workspace: (M+NRHS)*(N+2)
379 *
380 RESULT( 6 ) = ZERO
381 *
382 IF( N.GT.CRANK )
383 $ RESULT( 6 ) = SQRT14( 'No transpose', M, N,
384 $ NRHS, COPYA, LDA, B, LDB, WORK,
385 $ LWORK )
386 *
387 * Print information about the tests that did not
388 * pass the threshold.
389 *
390 DO 60 K = 3, 6
391 IF( RESULT( K ).GE.THRESH ) THEN
392 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
393 $ CALL ALAHD( NOUT, PATH )
394 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
395 $ ITYPE, K, RESULT( K )
396 NFAIL = NFAIL + 1
397 END IF
398 60 CONTINUE
399 NRUN = NRUN + 4
400 *
401 * Loop for testing different block sizes.
402 *
403 DO 100 INB = 1, NNB
404 NB = NBVAL( INB )
405 CALL XLAENV( 1, NB )
406 CALL XLAENV( 3, NXVAL( INB ) )
407 *
408 * Test SGELSY
409 *
410 * SGELSY: Compute the minimum-norm solution X
411 * to min( norm( A * X - B ) )
412 * using the rank-revealing orthogonal
413 * factorization.
414 *
415 * Initialize vector IWORK.
416 *
417 DO 70 J = 1, N
418 IWORK( J ) = 0
419 70 CONTINUE
420 *
421 * Set LWLSY to the adequate value.
422 *
423 LWLSY = MAX( 1, MNMIN+2*N+NB*( N+1 ),
424 $ 2*MNMIN+NB*NRHS )
425 *
426 CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
427 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B,
428 $ LDB )
429 *
430 SRNAMT = 'SGELSY'
431 CALL SGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
432 $ RCOND, CRANK, WORK, LWLSY, INFO )
433 IF( INFO.NE.0 )
434 $ CALL ALAERH( PATH, 'SGELSY', INFO, 0, ' ', M,
435 $ N, NRHS, -1, NB, ITYPE, NFAIL,
436 $ NERRS, NOUT )
437 *
438 * Test 7: Compute relative error in svd
439 * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
440 *
441 RESULT( 7 ) = SQRT12( CRANK, CRANK, A, LDA,
442 $ COPYS, WORK, LWORK )
443 *
444 * Test 8: Compute error in solution
445 * workspace: M*NRHS + M
446 *
447 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
448 $ LDWORK )
449 CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
450 $ LDA, B, LDB, WORK, LDWORK,
451 $ WORK( M*NRHS+1 ), RESULT( 8 ) )
452 *
453 * Test 9: Check norm of r'*A
454 * workspace: NRHS*(M+N)
455 *
456 RESULT( 9 ) = ZERO
457 IF( M.GT.CRANK )
458 $ RESULT( 9 ) = SQRT17( 'No transpose', 1, M,
459 $ N, NRHS, COPYA, LDA, B, LDB,
460 $ COPYB, LDB, C, WORK, LWORK )
461 *
462 * Test 10: Check if x is in the rowspace of A
463 * workspace: (M+NRHS)*(N+2)
464 *
465 RESULT( 10 ) = ZERO
466 *
467 IF( N.GT.CRANK )
468 $ RESULT( 10 ) = SQRT14( 'No transpose', M, N,
469 $ NRHS, COPYA, LDA, B, LDB,
470 $ WORK, LWORK )
471 *
472 * Test SGELSS
473 *
474 * SGELSS: Compute the minimum-norm solution X
475 * to min( norm( A * X - B ) )
476 * using the SVD.
477 *
478 CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
479 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B,
480 $ LDB )
481 SRNAMT = 'SGELSS'
482 CALL SGELSS( M, N, NRHS, A, LDA, B, LDB, S,
483 $ RCOND, CRANK, WORK, LWORK, INFO )
484 IF( INFO.NE.0 )
485 $ CALL ALAERH( PATH, 'SGELSS', INFO, 0, ' ', M,
486 $ N, NRHS, -1, NB, ITYPE, NFAIL,
487 $ NERRS, NOUT )
488 *
489 * workspace used: 3*min(m,n) +
490 * max(2*min(m,n),nrhs,max(m,n))
491 *
492 * Test 11: Compute relative error in svd
493 *
494 IF( RANK.GT.0 ) THEN
495 CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
496 RESULT( 11 ) = SASUM( MNMIN, S, 1 ) /
497 $ SASUM( MNMIN, COPYS, 1 ) /
498 $ ( EPS*REAL( MNMIN ) )
499 ELSE
500 RESULT( 11 ) = ZERO
501 END IF
502 *
503 * Test 12: Compute error in solution
504 *
505 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
506 $ LDWORK )
507 CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
508 $ LDA, B, LDB, WORK, LDWORK,
509 $ WORK( M*NRHS+1 ), RESULT( 12 ) )
510 *
511 * Test 13: Check norm of r'*A
512 *
513 RESULT( 13 ) = ZERO
514 IF( M.GT.CRANK )
515 $ RESULT( 13 ) = SQRT17( 'No transpose', 1, M,
516 $ N, NRHS, COPYA, LDA, B, LDB,
517 $ COPYB, LDB, C, WORK, LWORK )
518 *
519 * Test 14: Check if x is in the rowspace of A
520 *
521 RESULT( 14 ) = ZERO
522 IF( N.GT.CRANK )
523 $ RESULT( 14 ) = SQRT14( 'No transpose', M, N,
524 $ NRHS, COPYA, LDA, B, LDB,
525 $ WORK, LWORK )
526 *
527 * Test SGELSD
528 *
529 * SGELSD: Compute the minimum-norm solution X
530 * to min( norm( A * X - B ) ) using a
531 * divide and conquer SVD.
532 *
533 * Initialize vector IWORK.
534 *
535 DO 80 J = 1, N
536 IWORK( J ) = 0
537 80 CONTINUE
538 *
539 CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
540 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B,
541 $ LDB )
542 *
543 SRNAMT = 'SGELSD'
544 CALL SGELSD( M, N, NRHS, A, LDA, B, LDB, S,
545 $ RCOND, CRANK, WORK, LWORK, IWORK,
546 $ INFO )
547 IF( INFO.NE.0 )
548 $ CALL ALAERH( PATH, 'SGELSD', INFO, 0, ' ', M,
549 $ N, NRHS, -1, NB, ITYPE, NFAIL,
550 $ NERRS, NOUT )
551 *
552 * Test 15: Compute relative error in svd
553 *
554 IF( RANK.GT.0 ) THEN
555 CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
556 RESULT( 15 ) = SASUM( MNMIN, S, 1 ) /
557 $ SASUM( MNMIN, COPYS, 1 ) /
558 $ ( EPS*REAL( MNMIN ) )
559 ELSE
560 RESULT( 15 ) = ZERO
561 END IF
562 *
563 * Test 16: Compute error in solution
564 *
565 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
566 $ LDWORK )
567 CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
568 $ LDA, B, LDB, WORK, LDWORK,
569 $ WORK( M*NRHS+1 ), RESULT( 16 ) )
570 *
571 * Test 17: Check norm of r'*A
572 *
573 RESULT( 17 ) = ZERO
574 IF( M.GT.CRANK )
575 $ RESULT( 17 ) = SQRT17( 'No transpose', 1, M,
576 $ N, NRHS, COPYA, LDA, B, LDB,
577 $ COPYB, LDB, C, WORK, LWORK )
578 *
579 * Test 18: Check if x is in the rowspace of A
580 *
581 RESULT( 18 ) = ZERO
582 IF( N.GT.CRANK )
583 $ RESULT( 18 ) = SQRT14( 'No transpose', M, N,
584 $ NRHS, COPYA, LDA, B, LDB,
585 $ WORK, LWORK )
586 *
587 * Print information about the tests that did not
588 * pass the threshold.
589 *
590 DO 90 K = 7, NTESTS
591 IF( RESULT( K ).GE.THRESH ) THEN
592 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
593 $ CALL ALAHD( NOUT, PATH )
594 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
595 $ ITYPE, K, RESULT( K )
596 NFAIL = NFAIL + 1
597 END IF
598 90 CONTINUE
599 NRUN = NRUN + 12
600 *
601 100 CONTINUE
602 110 CONTINUE
603 120 CONTINUE
604 130 CONTINUE
605 140 CONTINUE
606 150 CONTINUE
607 *
608 * Print a summary of the results.
609 *
610 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
611 *
612 9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4,
613 $ ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 )
614 9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4,
615 $ ', type', I2, ', test(', I2, ')=', G12.5 )
616 RETURN
617 *
618 * End of SDRVLS
619 *
620 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 REAL THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL DOTYPE( * )
16 INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
17 $ NVAL( * ), NXVAL( * )
18 REAL A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
19 $ COPYS( * ), S( * ), WORK( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * SDRVLS tests the least squares driver routines SGELS, SGELSS, SGELSX,
26 * SGELSY and SGELSD.
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) REAL
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) REAL 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) REAL array, dimension (MMAX*NMAX)
90 *
91 * B (workspace) REAL 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) REAL array, dimension (MMAX*NSMAX)
96 *
97 * C (workspace) REAL array, dimension (MMAX*NSMAX)
98 *
99 * S (workspace) REAL array, dimension
100 * (min(MMAX,NMAX))
101 *
102 * COPYS (workspace) REAL array, dimension
103 * (min(MMAX,NMAX))
104 *
105 * WORK (workspace) REAL 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 REAL ONE, TWO, ZERO
121 PARAMETER ( ONE = 1.0E0, TWO = 2.0E0, ZERO = 0.0E0 )
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 REAL EPS, NORMA, NORMB, RCOND
131 * ..
132 * .. Local Arrays ..
133 INTEGER ISEED( 4 ), ISEEDY( 4 )
134 REAL RESULT( NTESTS )
135 * ..
136 * .. External Functions ..
137 REAL SASUM, SLAMCH, SQRT12, SQRT14, SQRT17
138 EXTERNAL SASUM, SLAMCH, SQRT12, SQRT14, SQRT17
139 * ..
140 * .. External Subroutines ..
141 EXTERNAL ALAERH, ALAHD, ALASVM, SAXPY, SERRLS, SGELS,
142 $ SGELSD, SGELSS, SGELSX, SGELSY, SGEMM, SLACPY,
143 $ SLARNV, SQRT13, SQRT15, SQRT16, SSCAL,
144 $ XLAENV
145 * ..
146 * .. Intrinsic Functions ..
147 INTRINSIC INT, LOG, MAX, MIN, REAL, 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 ) = 'Single 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 = SLAMCH( '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 SERRLS( 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 *
192 DO 150 IM = 1, NM
193 M = MVAL( IM )
194 LDA = MAX( 1, M )
195 *
196 DO 140 IN = 1, NN
197 N = NVAL( IN )
198 MNMIN = MIN( M, N )
199 LDB = MAX( 1, M, N )
200 *
201 DO 130 INS = 1, NNS
202 NRHS = NSVAL( INS )
203 NLVL = MAX( INT( LOG( MAX( ONE, REAL( MNMIN ) ) /
204 $ REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1, 0 )
205 LWORK = MAX( 1, ( M+NRHS )*( N+2 ), ( N+NRHS )*( M+2 ),
206 $ M*N+4*MNMIN+MAX( M, N ), 12*MNMIN+2*MNMIN*SMLSIZ+
207 $ 8*MNMIN*NLVL+MNMIN*NRHS+(SMLSIZ+1)**2 )
208 *
209 DO 120 IRANK = 1, 2
210 DO 110 ISCALE = 1, 3
211 ITYPE = ( IRANK-1 )*3 + ISCALE
212 IF( .NOT.DOTYPE( ITYPE ) )
213 $ GO TO 110
214 *
215 IF( IRANK.EQ.1 ) THEN
216 *
217 * Test SGELS
218 *
219 * Generate a matrix of scaling type ISCALE
220 *
221 CALL SQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
222 $ ISEED )
223 DO 40 INB = 1, NNB
224 NB = NBVAL( INB )
225 CALL XLAENV( 1, NB )
226 CALL XLAENV( 3, NXVAL( INB ) )
227 *
228 DO 30 ITRAN = 1, 2
229 IF( ITRAN.EQ.1 ) THEN
230 TRANS = 'N'
231 NROWS = M
232 NCOLS = N
233 ELSE
234 TRANS = 'T'
235 NROWS = N
236 NCOLS = M
237 END IF
238 LDWORK = MAX( 1, NCOLS )
239 *
240 * Set up a consistent rhs
241 *
242 IF( NCOLS.GT.0 ) THEN
243 CALL SLARNV( 2, ISEED, NCOLS*NRHS,
244 $ WORK )
245 CALL SSCAL( NCOLS*NRHS,
246 $ ONE / REAL( NCOLS ), WORK,
247 $ 1 )
248 END IF
249 CALL SGEMM( TRANS, 'No transpose', NROWS,
250 $ NRHS, NCOLS, ONE, COPYA, LDA,
251 $ WORK, LDWORK, ZERO, B, LDB )
252 CALL SLACPY( 'Full', NROWS, NRHS, B, LDB,
253 $ COPYB, LDB )
254 *
255 * Solve LS or overdetermined system
256 *
257 IF( M.GT.0 .AND. N.GT.0 ) THEN
258 CALL SLACPY( 'Full', M, N, COPYA, LDA,
259 $ A, LDA )
260 CALL SLACPY( 'Full', NROWS, NRHS,
261 $ COPYB, LDB, B, LDB )
262 END IF
263 SRNAMT = 'SGELS '
264 CALL SGELS( TRANS, M, N, NRHS, A, LDA, B,
265 $ LDB, WORK, LWORK, INFO )
266 IF( INFO.NE.0 )
267 $ CALL ALAERH( PATH, 'SGELS ', INFO, 0,
268 $ TRANS, M, N, NRHS, -1, NB,
269 $ ITYPE, NFAIL, NERRS,
270 $ NOUT )
271 *
272 * Check correctness of results
273 *
274 LDWORK = MAX( 1, NROWS )
275 IF( NROWS.GT.0 .AND. NRHS.GT.0 )
276 $ CALL SLACPY( 'Full', NROWS, NRHS,
277 $ COPYB, LDB, C, LDB )
278 CALL SQRT16( TRANS, M, N, NRHS, COPYA,
279 $ LDA, B, LDB, C, LDB, WORK,
280 $ RESULT( 1 ) )
281 *
282 IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
283 $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
284 *
285 * Solving LS system
286 *
287 RESULT( 2 ) = SQRT17( TRANS, 1, M, N,
288 $ NRHS, COPYA, LDA, B, LDB,
289 $ COPYB, LDB, C, WORK,
290 $ LWORK )
291 ELSE
292 *
293 * Solving overdetermined system
294 *
295 RESULT( 2 ) = SQRT14( TRANS, M, N,
296 $ NRHS, COPYA, LDA, B, LDB,
297 $ WORK, LWORK )
298 END IF
299 *
300 * Print information about the tests that
301 * did not pass the threshold.
302 *
303 DO 20 K = 1, 2
304 IF( RESULT( K ).GE.THRESH ) THEN
305 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
306 $ CALL ALAHD( NOUT, PATH )
307 WRITE( NOUT, FMT = 9999 )TRANS, M,
308 $ N, NRHS, NB, ITYPE, K,
309 $ RESULT( K )
310 NFAIL = NFAIL + 1
311 END IF
312 20 CONTINUE
313 NRUN = NRUN + 2
314 30 CONTINUE
315 40 CONTINUE
316 END IF
317 *
318 * Generate a matrix of scaling type ISCALE and rank
319 * type IRANK.
320 *
321 CALL SQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA,
322 $ COPYB, LDB, COPYS, RANK, NORMA, NORMB,
323 $ ISEED, WORK, LWORK )
324 *
325 * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
326 *
327 * Initialize vector IWORK.
328 *
329 DO 50 J = 1, N
330 IWORK( J ) = 0
331 50 CONTINUE
332 LDWORK = MAX( 1, M )
333 *
334 * Test SGELSX
335 *
336 * SGELSX: Compute the minimum-norm solution X
337 * to min( norm( A * X - B ) ) using a complete
338 * orthogonal factorization.
339 *
340 CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
341 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B, LDB )
342 *
343 SRNAMT = 'SGELSX'
344 CALL SGELSX( M, N, NRHS, A, LDA, B, LDB, IWORK,
345 $ RCOND, CRANK, WORK, INFO )
346 IF( INFO.NE.0 )
347 $ CALL ALAERH( PATH, 'SGELSX', INFO, 0, ' ', M, N,
348 $ NRHS, -1, NB, ITYPE, NFAIL, NERRS,
349 $ NOUT )
350 *
351 * workspace used: MAX( MNMIN+3*N, 2*MNMIN+NRHS )
352 *
353 * Test 3: Compute relative error in svd
354 * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
355 *
356 RESULT( 3 ) = SQRT12( CRANK, CRANK, A, LDA, COPYS,
357 $ WORK, LWORK )
358 *
359 * Test 4: Compute error in solution
360 * workspace: M*NRHS + M
361 *
362 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
363 $ LDWORK )
364 CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
365 $ LDA, B, LDB, WORK, LDWORK,
366 $ WORK( M*NRHS+1 ), RESULT( 4 ) )
367 *
368 * Test 5: Check norm of r'*A
369 * workspace: NRHS*(M+N)
370 *
371 RESULT( 5 ) = ZERO
372 IF( M.GT.CRANK )
373 $ RESULT( 5 ) = SQRT17( 'No transpose', 1, M, N,
374 $ NRHS, COPYA, LDA, B, LDB, COPYB,
375 $ LDB, C, WORK, LWORK )
376 *
377 * Test 6: Check if x is in the rowspace of A
378 * workspace: (M+NRHS)*(N+2)
379 *
380 RESULT( 6 ) = ZERO
381 *
382 IF( N.GT.CRANK )
383 $ RESULT( 6 ) = SQRT14( 'No transpose', M, N,
384 $ NRHS, COPYA, LDA, B, LDB, WORK,
385 $ LWORK )
386 *
387 * Print information about the tests that did not
388 * pass the threshold.
389 *
390 DO 60 K = 3, 6
391 IF( RESULT( K ).GE.THRESH ) THEN
392 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
393 $ CALL ALAHD( NOUT, PATH )
394 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
395 $ ITYPE, K, RESULT( K )
396 NFAIL = NFAIL + 1
397 END IF
398 60 CONTINUE
399 NRUN = NRUN + 4
400 *
401 * Loop for testing different block sizes.
402 *
403 DO 100 INB = 1, NNB
404 NB = NBVAL( INB )
405 CALL XLAENV( 1, NB )
406 CALL XLAENV( 3, NXVAL( INB ) )
407 *
408 * Test SGELSY
409 *
410 * SGELSY: Compute the minimum-norm solution X
411 * to min( norm( A * X - B ) )
412 * using the rank-revealing orthogonal
413 * factorization.
414 *
415 * Initialize vector IWORK.
416 *
417 DO 70 J = 1, N
418 IWORK( J ) = 0
419 70 CONTINUE
420 *
421 * Set LWLSY to the adequate value.
422 *
423 LWLSY = MAX( 1, MNMIN+2*N+NB*( N+1 ),
424 $ 2*MNMIN+NB*NRHS )
425 *
426 CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
427 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B,
428 $ LDB )
429 *
430 SRNAMT = 'SGELSY'
431 CALL SGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
432 $ RCOND, CRANK, WORK, LWLSY, INFO )
433 IF( INFO.NE.0 )
434 $ CALL ALAERH( PATH, 'SGELSY', INFO, 0, ' ', M,
435 $ N, NRHS, -1, NB, ITYPE, NFAIL,
436 $ NERRS, NOUT )
437 *
438 * Test 7: Compute relative error in svd
439 * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
440 *
441 RESULT( 7 ) = SQRT12( CRANK, CRANK, A, LDA,
442 $ COPYS, WORK, LWORK )
443 *
444 * Test 8: Compute error in solution
445 * workspace: M*NRHS + M
446 *
447 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
448 $ LDWORK )
449 CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
450 $ LDA, B, LDB, WORK, LDWORK,
451 $ WORK( M*NRHS+1 ), RESULT( 8 ) )
452 *
453 * Test 9: Check norm of r'*A
454 * workspace: NRHS*(M+N)
455 *
456 RESULT( 9 ) = ZERO
457 IF( M.GT.CRANK )
458 $ RESULT( 9 ) = SQRT17( 'No transpose', 1, M,
459 $ N, NRHS, COPYA, LDA, B, LDB,
460 $ COPYB, LDB, C, WORK, LWORK )
461 *
462 * Test 10: Check if x is in the rowspace of A
463 * workspace: (M+NRHS)*(N+2)
464 *
465 RESULT( 10 ) = ZERO
466 *
467 IF( N.GT.CRANK )
468 $ RESULT( 10 ) = SQRT14( 'No transpose', M, N,
469 $ NRHS, COPYA, LDA, B, LDB,
470 $ WORK, LWORK )
471 *
472 * Test SGELSS
473 *
474 * SGELSS: Compute the minimum-norm solution X
475 * to min( norm( A * X - B ) )
476 * using the SVD.
477 *
478 CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
479 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B,
480 $ LDB )
481 SRNAMT = 'SGELSS'
482 CALL SGELSS( M, N, NRHS, A, LDA, B, LDB, S,
483 $ RCOND, CRANK, WORK, LWORK, INFO )
484 IF( INFO.NE.0 )
485 $ CALL ALAERH( PATH, 'SGELSS', INFO, 0, ' ', M,
486 $ N, NRHS, -1, NB, ITYPE, NFAIL,
487 $ NERRS, NOUT )
488 *
489 * workspace used: 3*min(m,n) +
490 * max(2*min(m,n),nrhs,max(m,n))
491 *
492 * Test 11: Compute relative error in svd
493 *
494 IF( RANK.GT.0 ) THEN
495 CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
496 RESULT( 11 ) = SASUM( MNMIN, S, 1 ) /
497 $ SASUM( MNMIN, COPYS, 1 ) /
498 $ ( EPS*REAL( MNMIN ) )
499 ELSE
500 RESULT( 11 ) = ZERO
501 END IF
502 *
503 * Test 12: Compute error in solution
504 *
505 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
506 $ LDWORK )
507 CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
508 $ LDA, B, LDB, WORK, LDWORK,
509 $ WORK( M*NRHS+1 ), RESULT( 12 ) )
510 *
511 * Test 13: Check norm of r'*A
512 *
513 RESULT( 13 ) = ZERO
514 IF( M.GT.CRANK )
515 $ RESULT( 13 ) = SQRT17( 'No transpose', 1, M,
516 $ N, NRHS, COPYA, LDA, B, LDB,
517 $ COPYB, LDB, C, WORK, LWORK )
518 *
519 * Test 14: Check if x is in the rowspace of A
520 *
521 RESULT( 14 ) = ZERO
522 IF( N.GT.CRANK )
523 $ RESULT( 14 ) = SQRT14( 'No transpose', M, N,
524 $ NRHS, COPYA, LDA, B, LDB,
525 $ WORK, LWORK )
526 *
527 * Test SGELSD
528 *
529 * SGELSD: Compute the minimum-norm solution X
530 * to min( norm( A * X - B ) ) using a
531 * divide and conquer SVD.
532 *
533 * Initialize vector IWORK.
534 *
535 DO 80 J = 1, N
536 IWORK( J ) = 0
537 80 CONTINUE
538 *
539 CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
540 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B,
541 $ LDB )
542 *
543 SRNAMT = 'SGELSD'
544 CALL SGELSD( M, N, NRHS, A, LDA, B, LDB, S,
545 $ RCOND, CRANK, WORK, LWORK, IWORK,
546 $ INFO )
547 IF( INFO.NE.0 )
548 $ CALL ALAERH( PATH, 'SGELSD', INFO, 0, ' ', M,
549 $ N, NRHS, -1, NB, ITYPE, NFAIL,
550 $ NERRS, NOUT )
551 *
552 * Test 15: Compute relative error in svd
553 *
554 IF( RANK.GT.0 ) THEN
555 CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
556 RESULT( 15 ) = SASUM( MNMIN, S, 1 ) /
557 $ SASUM( MNMIN, COPYS, 1 ) /
558 $ ( EPS*REAL( MNMIN ) )
559 ELSE
560 RESULT( 15 ) = ZERO
561 END IF
562 *
563 * Test 16: Compute error in solution
564 *
565 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
566 $ LDWORK )
567 CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
568 $ LDA, B, LDB, WORK, LDWORK,
569 $ WORK( M*NRHS+1 ), RESULT( 16 ) )
570 *
571 * Test 17: Check norm of r'*A
572 *
573 RESULT( 17 ) = ZERO
574 IF( M.GT.CRANK )
575 $ RESULT( 17 ) = SQRT17( 'No transpose', 1, M,
576 $ N, NRHS, COPYA, LDA, B, LDB,
577 $ COPYB, LDB, C, WORK, LWORK )
578 *
579 * Test 18: Check if x is in the rowspace of A
580 *
581 RESULT( 18 ) = ZERO
582 IF( N.GT.CRANK )
583 $ RESULT( 18 ) = SQRT14( 'No transpose', M, N,
584 $ NRHS, COPYA, LDA, B, LDB,
585 $ WORK, LWORK )
586 *
587 * Print information about the tests that did not
588 * pass the threshold.
589 *
590 DO 90 K = 7, NTESTS
591 IF( RESULT( K ).GE.THRESH ) THEN
592 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
593 $ CALL ALAHD( NOUT, PATH )
594 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
595 $ ITYPE, K, RESULT( K )
596 NFAIL = NFAIL + 1
597 END IF
598 90 CONTINUE
599 NRUN = NRUN + 12
600 *
601 100 CONTINUE
602 110 CONTINUE
603 120 CONTINUE
604 130 CONTINUE
605 140 CONTINUE
606 150 CONTINUE
607 *
608 * Print a summary of the results.
609 *
610 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
611 *
612 9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4,
613 $ ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 )
614 9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4,
615 $ ', type', I2, ', test(', I2, ')=', G12.5 )
616 RETURN
617 *
618 * End of SDRVLS
619 *
620 END