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