LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zdrvls.f
Go to the documentation of this file.
1 *> \brief \b ZDRVLS
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE ZDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
12 * NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
13 * COPYB, C, S, COPYS, WORK, RWORK, IWORK, NOUT )
14 *
15 * .. Scalar Arguments ..
16 * LOGICAL TSTERR
17 * INTEGER NM, NN, NNB, NNS, NOUT
18 * DOUBLE PRECISION THRESH
19 * ..
20 * .. Array Arguments ..
21 * LOGICAL DOTYPE( * )
22 * INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
23 * $ NVAL( * ), NXVAL( * )
24 * DOUBLE PRECISION COPYS( * ), RWORK( * ), S( * )
25 * COMPLEX*16 A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
26 * $ WORK( * )
27 * ..
28 *
29 *
30 *> \par Purpose:
31 * =============
32 *>
33 *> \verbatim
34 *>
35 *> ZDRVLS tests the least squares driver routines ZGELS, CGELSX, CGELSS,
36 *> ZGELSY and CGELSD.
37 *> \endverbatim
38 *
39 * Arguments:
40 * ==========
41 *
42 *> \param[in] DOTYPE
43 *> \verbatim
44 *> DOTYPE is LOGICAL array, dimension (NTYPES)
45 *> The matrix types to be used for testing. Matrices of type j
46 *> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
47 *> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
48 *> The matrix of type j is generated as follows:
49 *> j=1: A = U*D*V where U and V are random unitary matrices
50 *> and D has random entries (> 0.1) taken from a uniform
51 *> distribution (0,1). A is full rank.
52 *> j=2: The same of 1, but A is scaled up.
53 *> j=3: The same of 1, but A is scaled down.
54 *> j=4: A = U*D*V where U and V are random unitary matrices
55 *> and D has 3*min(M,N)/4 random entries (> 0.1) taken
56 *> from a uniform distribution (0,1) and the remaining
57 *> entries set to 0. A is rank-deficient.
58 *> j=5: The same of 4, but A is scaled up.
59 *> j=6: The same of 5, but A is scaled down.
60 *> \endverbatim
61 *>
62 *> \param[in] NM
63 *> \verbatim
64 *> NM is INTEGER
65 *> The number of values of M contained in the vector MVAL.
66 *> \endverbatim
67 *>
68 *> \param[in] MVAL
69 *> \verbatim
70 *> MVAL is INTEGER array, dimension (NM)
71 *> The values of the matrix row dimension M.
72 *> \endverbatim
73 *>
74 *> \param[in] NN
75 *> \verbatim
76 *> NN is INTEGER
77 *> The number of values of N contained in the vector NVAL.
78 *> \endverbatim
79 *>
80 *> \param[in] NVAL
81 *> \verbatim
82 *> NVAL is INTEGER array, dimension (NN)
83 *> The values of the matrix column dimension N.
84 *> \endverbatim
85 *>
86 *> \param[in] NNB
87 *> \verbatim
88 *> NNB is INTEGER
89 *> The number of values of NB and NX contained in the
90 *> vectors NBVAL and NXVAL. The blocking parameters are used
91 *> in pairs (NB,NX).
92 *> \endverbatim
93 *>
94 *> \param[in] NBVAL
95 *> \verbatim
96 *> NBVAL is INTEGER array, dimension (NNB)
97 *> The values of the blocksize NB.
98 *> \endverbatim
99 *>
100 *> \param[in] NXVAL
101 *> \verbatim
102 *> NXVAL is INTEGER array, dimension (NNB)
103 *> The values of the crossover point NX.
104 *> \endverbatim
105 *>
106 *> \param[in] NNS
107 *> \verbatim
108 *> NNS is INTEGER
109 *> The number of values of NRHS contained in the vector NSVAL.
110 *> \endverbatim
111 *>
112 *> \param[in] NSVAL
113 *> \verbatim
114 *> NSVAL is INTEGER array, dimension (NNS)
115 *> The values of the number of right hand sides NRHS.
116 *> \endverbatim
117 *>
118 *> \param[in] THRESH
119 *> \verbatim
120 *> THRESH is DOUBLE PRECISION
121 *> The threshold value for the test ratios. A result is
122 *> included in the output file if RESULT >= THRESH. To have
123 *> every test ratio printed, use THRESH = 0.
124 *> \endverbatim
125 *>
126 *> \param[in] TSTERR
127 *> \verbatim
128 *> TSTERR is LOGICAL
129 *> Flag that indicates whether error exits are to be tested.
130 *> \endverbatim
131 *>
132 *> \param[out] A
133 *> \verbatim
134 *> A is COMPLEX*16 array, dimension (MMAX*NMAX)
135 *> where MMAX is the maximum value of M in MVAL and NMAX is the
136 *> maximum value of N in NVAL.
137 *> \endverbatim
138 *>
139 *> \param[out] COPYA
140 *> \verbatim
141 *> COPYA is COMPLEX*16 array, dimension (MMAX*NMAX)
142 *> \endverbatim
143 *>
144 *> \param[out] B
145 *> \verbatim
146 *> B is COMPLEX*16 array, dimension (MMAX*NSMAX)
147 *> where MMAX is the maximum value of M in MVAL and NSMAX is the
148 *> maximum value of NRHS in NSVAL.
149 *> \endverbatim
150 *>
151 *> \param[out] COPYB
152 *> \verbatim
153 *> COPYB is COMPLEX*16 array, dimension (MMAX*NSMAX)
154 *> \endverbatim
155 *>
156 *> \param[out] C
157 *> \verbatim
158 *> C is COMPLEX*16 array, dimension (MMAX*NSMAX)
159 *> \endverbatim
160 *>
161 *> \param[out] S
162 *> \verbatim
163 *> S is DOUBLE PRECISION array, dimension
164 *> (min(MMAX,NMAX))
165 *> \endverbatim
166 *>
167 *> \param[out] COPYS
168 *> \verbatim
169 *> COPYS is DOUBLE PRECISION array, dimension
170 *> (min(MMAX,NMAX))
171 *> \endverbatim
172 *>
173 *> \param[out] WORK
174 *> \verbatim
175 *> WORK is COMPLEX*16 array, dimension
176 *> (MMAX*NMAX + 4*NMAX + MMAX).
177 *> \endverbatim
178 *>
179 *> \param[out] RWORK
180 *> \verbatim
181 *> RWORK is DOUBLE PRECISION array, dimension (5*NMAX-1)
182 *> \endverbatim
183 *>
184 *> \param[out] IWORK
185 *> \verbatim
186 *> IWORK is INTEGER array, dimension (15*NMAX)
187 *> \endverbatim
188 *>
189 *> \param[in] NOUT
190 *> \verbatim
191 *> NOUT is INTEGER
192 *> The unit number for output.
193 *> \endverbatim
194 *
195 * Authors:
196 * ========
197 *
198 *> \author Univ. of Tennessee
199 *> \author Univ. of California Berkeley
200 *> \author Univ. of Colorado Denver
201 *> \author NAG Ltd.
202 *
203 *> \date November 2011
204 *
205 *> \ingroup complex16_lin
206 *
207 * =====================================================================
208  SUBROUTINE zdrvls( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
209  $ nbval, nxval, thresh, tsterr, a, copya, b,
210  $ copyb, c, s, copys, work, rwork, iwork, nout )
211 *
212 * -- LAPACK test routine (version 3.4.0) --
213 * -- LAPACK is a software package provided by Univ. of Tennessee, --
214 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
215 * November 2011
216 *
217 * .. Scalar Arguments ..
218  LOGICAL tsterr
219  INTEGER nm, nn, nnb, nns, nout
220  DOUBLE PRECISION thresh
221 * ..
222 * .. Array Arguments ..
223  LOGICAL dotype( * )
224  INTEGER iwork( * ), mval( * ), nbval( * ), nsval( * ),
225  $ nval( * ), nxval( * )
226  DOUBLE PRECISION copys( * ), rwork( * ), s( * )
227  COMPLEX*16 a( * ), b( * ), c( * ), copya( * ), copyb( * ),
228  $ work( * )
229 * ..
230 *
231 * =====================================================================
232 *
233 * .. Parameters ..
234  INTEGER ntests
235  parameter( ntests = 18 )
236  INTEGER smlsiz
237  parameter( smlsiz = 25 )
238  DOUBLE PRECISION one, zero
239  parameter( one = 1.0d+0, zero = 0.0d+0 )
240  COMPLEX*16 cone, czero
241  parameter( cone = ( 1.0d+0, 0.0d+0 ),
242  $ czero = ( 0.0d+0, 0.0d+0 ) )
243 * ..
244 * .. Local Scalars ..
245  CHARACTER trans
246  CHARACTER*3 path
247  INTEGER crank, i, im, in, inb, info, ins, irank,
248  $ iscale, itran, itype, j, k, lda, ldb, ldwork,
249  $ lwlsy, lwork, m, mnmin, n, nb, ncols, nerrs,
250  $ nfail, nrhs, nrows, nrun, rank
251  DOUBLE PRECISION eps, norma, normb, rcond
252 * ..
253 * .. Local Arrays ..
254  INTEGER iseed( 4 ), iseedy( 4 )
255  DOUBLE PRECISION result( ntests )
256 * ..
257 * .. External Functions ..
258  DOUBLE PRECISION dasum, dlamch, zqrt12, zqrt14, zqrt17
259  EXTERNAL dasum, dlamch, zqrt12, zqrt14, zqrt17
260 * ..
261 * .. External Subroutines ..
262  EXTERNAL alaerh, alahd, alasvm, daxpy, dlasrt, xlaenv,
265  $ zqrt16
266 * ..
267 * .. Intrinsic Functions ..
268  INTRINSIC dble, max, min, sqrt
269 * ..
270 * .. Scalars in Common ..
271  LOGICAL lerr, ok
272  CHARACTER*32 srnamt
273  INTEGER infot, iounit
274 * ..
275 * .. Common blocks ..
276  common / infoc / infot, iounit, ok, lerr
277  common / srnamc / srnamt
278 * ..
279 * .. Data statements ..
280  DATA iseedy / 1988, 1989, 1990, 1991 /
281 * ..
282 * .. Executable Statements ..
283 *
284 * Initialize constants and the random number seed.
285 *
286  path( 1: 1 ) = 'Zomplex precision'
287  path( 2: 3 ) = 'LS'
288  nrun = 0
289  nfail = 0
290  nerrs = 0
291  DO 10 i = 1, 4
292  iseed( i ) = iseedy( i )
293  10 continue
294  eps = dlamch( 'Epsilon' )
295 *
296 * Threshold for rank estimation
297 *
298  rcond = sqrt( eps ) - ( sqrt( eps )-eps ) / 2
299 *
300 * Test the error exits
301 *
302  CALL xlaenv( 9, smlsiz )
303  IF( tsterr )
304  $ CALL zerrls( path, nout )
305 *
306 * Print the header if NM = 0 or NN = 0 and THRESH = 0.
307 *
308  IF( ( nm.EQ.0 .OR. nn.EQ.0 ) .AND. thresh.EQ.zero )
309  $ CALL alahd( nout, path )
310  infot = 0
311 *
312  DO 140 im = 1, nm
313  m = mval( im )
314  lda = max( 1, m )
315 *
316  DO 130 in = 1, nn
317  n = nval( in )
318  mnmin = min( m, n )
319  ldb = max( 1, m, n )
320 *
321  DO 120 ins = 1, nns
322  nrhs = nsval( ins )
323  lwork = max( 1, ( m+nrhs )*( n+2 ), ( n+nrhs )*( m+2 ),
324  $ m*n+4*mnmin+max( m, n ), 2*n+m )
325 *
326  DO 110 irank = 1, 2
327  DO 100 iscale = 1, 3
328  itype = ( irank-1 )*3 + iscale
329  IF( .NOT.dotype( itype ) )
330  $ go to 100
331 *
332  IF( irank.EQ.1 ) THEN
333 *
334 * Test ZGELS
335 *
336 * Generate a matrix of scaling type ISCALE
337 *
338  CALL zqrt13( iscale, m, n, copya, lda, norma,
339  $ iseed )
340  DO 40 inb = 1, nnb
341  nb = nbval( inb )
342  CALL xlaenv( 1, nb )
343  CALL xlaenv( 3, nxval( inb ) )
344 *
345  DO 30 itran = 1, 2
346  IF( itran.EQ.1 ) THEN
347  trans = 'N'
348  nrows = m
349  ncols = n
350  ELSE
351  trans = 'C'
352  nrows = n
353  ncols = m
354  END IF
355  ldwork = max( 1, ncols )
356 *
357 * Set up a consistent rhs
358 *
359  IF( ncols.GT.0 ) THEN
360  CALL zlarnv( 2, iseed, ncols*nrhs,
361  $ work )
362  CALL zdscal( ncols*nrhs,
363  $ one / dble( ncols ), work,
364  $ 1 )
365  END IF
366  CALL zgemm( trans, 'No transpose', nrows,
367  $ nrhs, ncols, cone, copya, lda,
368  $ work, ldwork, czero, b, ldb )
369  CALL zlacpy( 'Full', nrows, nrhs, b, ldb,
370  $ copyb, ldb )
371 *
372 * Solve LS or overdetermined system
373 *
374  IF( m.GT.0 .AND. n.GT.0 ) THEN
375  CALL zlacpy( 'Full', m, n, copya, lda,
376  $ a, lda )
377  CALL zlacpy( 'Full', nrows, nrhs,
378  $ copyb, ldb, b, ldb )
379  END IF
380  srnamt = 'ZGELS '
381  CALL zgels( trans, m, n, nrhs, a, lda, b,
382  $ ldb, work, lwork, info )
383 *
384  IF( info.NE.0 )
385  $ CALL alaerh( path, 'ZGELS ', info, 0,
386  $ trans, m, n, nrhs, -1, nb,
387  $ itype, nfail, nerrs,
388  $ nout )
389 *
390 * Check correctness of results
391 *
392  ldwork = max( 1, nrows )
393  IF( nrows.GT.0 .AND. nrhs.GT.0 )
394  $ CALL zlacpy( 'Full', nrows, nrhs,
395  $ copyb, ldb, c, ldb )
396  CALL zqrt16( trans, m, n, nrhs, copya,
397  $ lda, b, ldb, c, ldb, rwork,
398  $ result( 1 ) )
399 *
400  IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
401  $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
402 *
403 * Solving LS system
404 *
405  result( 2 ) = zqrt17( trans, 1, m, n,
406  $ nrhs, copya, lda, b, ldb,
407  $ copyb, ldb, c, work,
408  $ lwork )
409  ELSE
410 *
411 * Solving overdetermined system
412 *
413  result( 2 ) = zqrt14( trans, m, n,
414  $ nrhs, copya, lda, b, ldb,
415  $ work, lwork )
416  END IF
417 *
418 * Print information about the tests that
419 * did not pass the threshold.
420 *
421  DO 20 k = 1, 2
422  IF( result( k ).GE.thresh ) THEN
423  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
424  $ CALL alahd( nout, path )
425  WRITE( nout, fmt = 9999 )trans, m,
426  $ n, nrhs, nb, itype, k,
427  $ result( k )
428  nfail = nfail + 1
429  END IF
430  20 continue
431  nrun = nrun + 2
432  30 continue
433  40 continue
434  END IF
435 *
436 * Generate a matrix of scaling type ISCALE and rank
437 * type IRANK.
438 *
439  CALL zqrt15( iscale, irank, m, n, nrhs, copya, lda,
440  $ copyb, ldb, copys, rank, norma, normb,
441  $ iseed, work, lwork )
442 *
443 * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
444 *
445  DO 50 j = 1, n
446  iwork( j ) = 0
447  50 continue
448  ldwork = max( 1, m )
449 *
450 * Test ZGELSX
451 *
452 * ZGELSX: Compute the minimum-norm solution X
453 * to min( norm( A * X - B ) )
454 * using a complete orthogonal factorization.
455 *
456  CALL zlacpy( 'Full', m, n, copya, lda, a, lda )
457  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, b, ldb )
458 *
459  srnamt = 'ZGELSX'
460  CALL zgelsx( m, n, nrhs, a, lda, b, ldb, iwork,
461  $ rcond, crank, work, rwork, info )
462 *
463  IF( info.NE.0 )
464  $ CALL alaerh( path, 'ZGELSX', info, 0, ' ', m, n,
465  $ nrhs, -1, nb, itype, nfail, nerrs,
466  $ nout )
467 *
468 * workspace used: MAX( MNMIN+3*N, 2*MNMIN+NRHS )
469 *
470 * Test 3: Compute relative error in svd
471 * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
472 *
473  result( 3 ) = zqrt12( crank, crank, a, lda, copys,
474  $ work, lwork, rwork )
475 *
476 * Test 4: Compute error in solution
477 * workspace: M*NRHS + M
478 *
479  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, work,
480  $ ldwork )
481  CALL zqrt16( 'No transpose', m, n, nrhs, copya,
482  $ lda, b, ldb, work, ldwork, rwork,
483  $ result( 4 ) )
484 *
485 * Test 5: Check norm of r'*A
486 * workspace: NRHS*(M+N)
487 *
488  result( 5 ) = zero
489  IF( m.GT.crank )
490  $ result( 5 ) = zqrt17( 'No transpose', 1, m, n,
491  $ nrhs, copya, lda, b, ldb, copyb,
492  $ ldb, c, work, lwork )
493 *
494 * Test 6: Check if x is in the rowspace of A
495 * workspace: (M+NRHS)*(N+2)
496 *
497  result( 6 ) = zero
498 *
499  IF( n.GT.crank )
500  $ result( 6 ) = zqrt14( 'No transpose', m, n,
501  $ nrhs, copya, lda, b, ldb, work,
502  $ lwork )
503 *
504 * Print information about the tests that did not
505 * pass the threshold.
506 *
507  DO 60 k = 3, 6
508  IF( result( k ).GE.thresh ) THEN
509  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
510  $ CALL alahd( nout, path )
511  WRITE( nout, fmt = 9998 )m, n, nrhs, 0,
512  $ itype, k, result( k )
513  nfail = nfail + 1
514  END IF
515  60 continue
516  nrun = nrun + 4
517 *
518 * Loop for testing different block sizes.
519 *
520  DO 90 inb = 1, nnb
521  nb = nbval( inb )
522  CALL xlaenv( 1, nb )
523  CALL xlaenv( 3, nxval( inb ) )
524 *
525 * Test ZGELSY
526 *
527 * ZGELSY: Compute the minimum-norm solution
528 * X to min( norm( A * X - B ) )
529 * using the rank-revealing orthogonal
530 * factorization.
531 *
532  CALL zlacpy( 'Full', m, n, copya, lda, a, lda )
533  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, b,
534  $ ldb )
535 *
536 * Initialize vector IWORK.
537 *
538  DO 70 j = 1, n
539  iwork( j ) = 0
540  70 continue
541 *
542 * Set LWLSY to the adequate value.
543 *
544  lwlsy = mnmin + max( 2*mnmin, nb*( n+1 ),
545  $ mnmin+nb*nrhs )
546  lwlsy = max( 1, lwlsy )
547 *
548  srnamt = 'ZGELSY'
549  CALL zgelsy( m, n, nrhs, a, lda, b, ldb, iwork,
550  $ rcond, crank, work, lwlsy, rwork,
551  $ info )
552  IF( info.NE.0 )
553  $ CALL alaerh( path, 'ZGELSY', info, 0, ' ', m,
554  $ n, nrhs, -1, nb, itype, nfail,
555  $ nerrs, nout )
556 *
557 * workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS)
558 *
559 * Test 7: Compute relative error in svd
560 * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
561 *
562  result( 7 ) = zqrt12( crank, crank, a, lda,
563  $ copys, work, lwork, rwork )
564 *
565 * Test 8: Compute error in solution
566 * workspace: M*NRHS + M
567 *
568  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, work,
569  $ ldwork )
570  CALL zqrt16( 'No transpose', m, n, nrhs, copya,
571  $ lda, b, ldb, work, ldwork, rwork,
572  $ result( 8 ) )
573 *
574 * Test 9: Check norm of r'*A
575 * workspace: NRHS*(M+N)
576 *
577  result( 9 ) = zero
578  IF( m.GT.crank )
579  $ result( 9 ) = zqrt17( 'No transpose', 1, m,
580  $ n, nrhs, copya, lda, b, ldb,
581  $ copyb, ldb, c, work, lwork )
582 *
583 * Test 10: Check if x is in the rowspace of A
584 * workspace: (M+NRHS)*(N+2)
585 *
586  result( 10 ) = zero
587 *
588  IF( n.GT.crank )
589  $ result( 10 ) = zqrt14( 'No transpose', m, n,
590  $ nrhs, copya, lda, b, ldb,
591  $ work, lwork )
592 *
593 * Test ZGELSS
594 *
595 * ZGELSS: Compute the minimum-norm solution
596 * X to min( norm( A * X - B ) )
597 * using the SVD.
598 *
599  CALL zlacpy( 'Full', m, n, copya, lda, a, lda )
600  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, b,
601  $ ldb )
602  srnamt = 'ZGELSS'
603  CALL zgelss( m, n, nrhs, a, lda, b, ldb, s,
604  $ rcond, crank, work, lwork, rwork,
605  $ info )
606 *
607  IF( info.NE.0 )
608  $ CALL alaerh( path, 'ZGELSS', info, 0, ' ', m,
609  $ n, nrhs, -1, nb, itype, nfail,
610  $ nerrs, nout )
611 *
612 * workspace used: 3*min(m,n) +
613 * max(2*min(m,n),nrhs,max(m,n))
614 *
615 * Test 11: Compute relative error in svd
616 *
617  IF( rank.GT.0 ) THEN
618  CALL daxpy( mnmin, -one, copys, 1, s, 1 )
619  result( 11 ) = dasum( mnmin, s, 1 ) /
620  $ dasum( mnmin, copys, 1 ) /
621  $ ( eps*dble( mnmin ) )
622  ELSE
623  result( 11 ) = zero
624  END IF
625 *
626 * Test 12: Compute error in solution
627 *
628  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, work,
629  $ ldwork )
630  CALL zqrt16( 'No transpose', m, n, nrhs, copya,
631  $ lda, b, ldb, work, ldwork, rwork,
632  $ result( 12 ) )
633 *
634 * Test 13: Check norm of r'*A
635 *
636  result( 13 ) = zero
637  IF( m.GT.crank )
638  $ result( 13 ) = zqrt17( 'No transpose', 1, m,
639  $ n, nrhs, copya, lda, b, ldb,
640  $ copyb, ldb, c, work, lwork )
641 *
642 * Test 14: Check if x is in the rowspace of A
643 *
644  result( 14 ) = zero
645  IF( n.GT.crank )
646  $ result( 14 ) = zqrt14( 'No transpose', m, n,
647  $ nrhs, copya, lda, b, ldb,
648  $ work, lwork )
649 *
650 * Test ZGELSD
651 *
652 * ZGELSD: Compute the minimum-norm solution X
653 * to min( norm( A * X - B ) ) using a
654 * divide and conquer SVD.
655 *
656  CALL xlaenv( 9, 25 )
657 *
658  CALL zlacpy( 'Full', m, n, copya, lda, a, lda )
659  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, b,
660  $ ldb )
661 *
662  srnamt = 'ZGELSD'
663  CALL zgelsd( m, n, nrhs, a, lda, b, ldb, s,
664  $ rcond, crank, work, lwork, rwork,
665  $ iwork, info )
666  IF( info.NE.0 )
667  $ CALL alaerh( path, 'ZGELSD', info, 0, ' ', m,
668  $ n, nrhs, -1, nb, itype, nfail,
669  $ nerrs, nout )
670 *
671 * Test 15: Compute relative error in svd
672 *
673  IF( rank.GT.0 ) THEN
674  CALL daxpy( mnmin, -one, copys, 1, s, 1 )
675  result( 15 ) = dasum( mnmin, s, 1 ) /
676  $ dasum( mnmin, copys, 1 ) /
677  $ ( eps*dble( mnmin ) )
678  ELSE
679  result( 15 ) = zero
680  END IF
681 *
682 * Test 16: Compute error in solution
683 *
684  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, work,
685  $ ldwork )
686  CALL zqrt16( 'No transpose', m, n, nrhs, copya,
687  $ lda, b, ldb, work, ldwork, rwork,
688  $ result( 16 ) )
689 *
690 * Test 17: Check norm of r'*A
691 *
692  result( 17 ) = zero
693  IF( m.GT.crank )
694  $ result( 17 ) = zqrt17( 'No transpose', 1, m,
695  $ n, nrhs, copya, lda, b, ldb,
696  $ copyb, ldb, c, work, lwork )
697 *
698 * Test 18: Check if x is in the rowspace of A
699 *
700  result( 18 ) = zero
701  IF( n.GT.crank )
702  $ result( 18 ) = zqrt14( 'No transpose', m, n,
703  $ nrhs, copya, lda, b, ldb,
704  $ work, lwork )
705 *
706 * Print information about the tests that did not
707 * pass the threshold.
708 *
709  DO 80 k = 7, ntests
710  IF( result( k ).GE.thresh ) THEN
711  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
712  $ CALL alahd( nout, path )
713  WRITE( nout, fmt = 9998 )m, n, nrhs, nb,
714  $ itype, k, result( k )
715  nfail = nfail + 1
716  END IF
717  80 continue
718  nrun = nrun + 12
719 *
720  90 continue
721  100 continue
722  110 continue
723  120 continue
724  130 continue
725  140 continue
726 *
727 * Print a summary of the results.
728 *
729  CALL alasvm( path, nout, nfail, nrun, nerrs )
730 *
731  9999 format( ' TRANS=''', a1, ''', M=', i5, ', N=', i5, ', NRHS=', i4,
732  $ ', NB=', i4, ', type', i2, ', test(', i2, ')=', g12.5 )
733  9998 format( ' M=', i5, ', N=', i5, ', NRHS=', i4, ', NB=', i4,
734  $ ', type', i2, ', test(', i2, ')=', g12.5 )
735  return
736 *
737 * End of ZDRVLS
738 *
739  END