LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
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, 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 MVAL( * ), NBVAL( * ), NSVAL( * ),
23 * $ NVAL( * ), NXVAL( * )
24 * DOUBLE PRECISION COPYS( * ), S( * )
25 * COMPLEX*16 A( * ), B( * ), C( * ), COPYA( * ), COPYB( * )
26 * ..
27 *
28 *
29 *> \par Purpose:
30 * =============
31 *>
32 *> \verbatim
33 *>
34 *> ZDRVLS tests the least squares driver routines ZGELS, ZGETSLS, ZGELSS, ZGELSY
35 *> and ZGELSD.
36 *> \endverbatim
37 *
38 * Arguments:
39 * ==========
40 *
41 *> \param[in] DOTYPE
42 *> \verbatim
43 *> DOTYPE is LOGICAL array, dimension (NTYPES)
44 *> The matrix types to be used for testing. Matrices of type j
45 *> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
46 *> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
47 *> The matrix of type j is generated as follows:
48 *> j=1: A = U*D*V where U and V are random unitary matrices
49 *> and D has random entries (> 0.1) taken from a uniform
50 *> distribution (0,1). A is full rank.
51 *> j=2: The same of 1, but A is scaled up.
52 *> j=3: The same of 1, but A is scaled down.
53 *> j=4: A = U*D*V where U and V are random unitary matrices
54 *> and D has 3*min(M,N)/4 random entries (> 0.1) taken
55 *> from a uniform distribution (0,1) and the remaining
56 *> entries set to 0. A is rank-deficient.
57 *> j=5: The same of 4, but A is scaled up.
58 *> j=6: The same of 5, but A is scaled down.
59 *> \endverbatim
60 *>
61 *> \param[in] NM
62 *> \verbatim
63 *> NM is INTEGER
64 *> The number of values of M contained in the vector MVAL.
65 *> \endverbatim
66 *>
67 *> \param[in] MVAL
68 *> \verbatim
69 *> MVAL is INTEGER array, dimension (NM)
70 *> The values of the matrix row dimension M.
71 *> \endverbatim
72 *>
73 *> \param[in] NN
74 *> \verbatim
75 *> NN is INTEGER
76 *> The number of values of N contained in the vector NVAL.
77 *> \endverbatim
78 *>
79 *> \param[in] NVAL
80 *> \verbatim
81 *> NVAL is INTEGER array, dimension (NN)
82 *> The values of the matrix column dimension N.
83 *> \endverbatim
84 *>
85 *> \param[in] NNB
86 *> \verbatim
87 *> NNB is INTEGER
88 *> The number of values of NB and NX contained in the
89 *> vectors NBVAL and NXVAL. The blocking parameters are used
90 *> in pairs (NB,NX).
91 *> \endverbatim
92 *>
93 *> \param[in] NBVAL
94 *> \verbatim
95 *> NBVAL is INTEGER array, dimension (NNB)
96 *> The values of the blocksize NB.
97 *> \endverbatim
98 *>
99 *> \param[in] NXVAL
100 *> \verbatim
101 *> NXVAL is INTEGER array, dimension (NNB)
102 *> The values of the crossover point NX.
103 *> \endverbatim
104 *>
105 *> \param[in] NNS
106 *> \verbatim
107 *> NNS is INTEGER
108 *> The number of values of NRHS contained in the vector NSVAL.
109 *> \endverbatim
110 *>
111 *> \param[in] NSVAL
112 *> \verbatim
113 *> NSVAL is INTEGER array, dimension (NNS)
114 *> The values of the number of right hand sides NRHS.
115 *> \endverbatim
116 *>
117 *> \param[in] THRESH
118 *> \verbatim
119 *> THRESH is DOUBLE PRECISION
120 *> The threshold value for the test ratios. A result is
121 *> included in the output file if RESULT >= THRESH. To have
122 *> every test ratio printed, use THRESH = 0.
123 *> \endverbatim
124 *>
125 *> \param[in] TSTERR
126 *> \verbatim
127 *> TSTERR is LOGICAL
128 *> Flag that indicates whether error exits are to be tested.
129 *> \endverbatim
130 *>
131 *> \param[out] A
132 *> \verbatim
133 *> A is COMPLEX*16 array, dimension (MMAX*NMAX)
134 *> where MMAX is the maximum value of M in MVAL and NMAX is the
135 *> maximum value of N in NVAL.
136 *> \endverbatim
137 *>
138 *> \param[out] COPYA
139 *> \verbatim
140 *> COPYA is COMPLEX*16 array, dimension (MMAX*NMAX)
141 *> \endverbatim
142 *>
143 *> \param[out] B
144 *> \verbatim
145 *> B is COMPLEX*16 array, dimension (MMAX*NSMAX)
146 *> where MMAX is the maximum value of M in MVAL and NSMAX is the
147 *> maximum value of NRHS in NSVAL.
148 *> \endverbatim
149 *>
150 *> \param[out] COPYB
151 *> \verbatim
152 *> COPYB is COMPLEX*16 array, dimension (MMAX*NSMAX)
153 *> \endverbatim
154 *>
155 *> \param[out] C
156 *> \verbatim
157 *> C is COMPLEX*16 array, dimension (MMAX*NSMAX)
158 *> \endverbatim
159 *>
160 *> \param[out] S
161 *> \verbatim
162 *> S is DOUBLE PRECISION array, dimension
163 *> (min(MMAX,NMAX))
164 *> \endverbatim
165 *>
166 *> \param[out] COPYS
167 *> \verbatim
168 *> COPYS is DOUBLE PRECISION array, dimension
169 *> (min(MMAX,NMAX))
170 *> \endverbatim
171 *>
172 *> \param[in] NOUT
173 *> \verbatim
174 *> NOUT is INTEGER
175 *> The unit number for output.
176 *> \endverbatim
177 *
178 * Authors:
179 * ========
180 *
181 *> \author Univ. of Tennessee
182 *> \author Univ. of California Berkeley
183 *> \author Univ. of Colorado Denver
184 *> \author NAG Ltd.
185 *
186 *> \date June 2017
187 *
188 *> \ingroup complex16_lin
189 *
190 * =====================================================================
191  SUBROUTINE zdrvls( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
192  $ NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
193  $ COPYB, C, S, COPYS, NOUT )
194 *
195 * -- LAPACK test routine (version 3.7.1) --
196 * -- LAPACK is a software package provided by Univ. of Tennessee, --
197 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
198 * June 2017
199 *
200 * .. Scalar Arguments ..
201  LOGICAL tsterr
202  INTEGER nm, nn, nnb, nns, nout
203  DOUBLE PRECISION thresh
204 * ..
205 * .. Array Arguments ..
206  LOGICAL dotype( * )
207  INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
208  $ nval( * ), nxval( * )
209  DOUBLE PRECISION COPYS( * ), S( * )
210  COMPLEX*16 A( * ), B( * ), C( * ), COPYA( * ), COPYB( * )
211 * ..
212 *
213 * =====================================================================
214 *
215 * .. Parameters ..
216  INTEGER NTESTS
217  PARAMETER ( NTESTS = 16 )
218  INTEGER SMLSIZ
219  parameter( smlsiz = 25 )
220  DOUBLE PRECISION ONE, ZERO
221  parameter( one = 1.0d+0, zero = 0.0d+0 )
222  COMPLEX*16 CONE, CZERO
223  parameter( cone = ( 1.0d+0, 0.0d+0 ),
224  $ czero = ( 0.0d+0, 0.0d+0 ) )
225 * ..
226 * .. Local Scalars ..
227  CHARACTER TRANS
228  CHARACTER*3 PATH
229  INTEGER CRANK, I, IM, IMB, IN, INB, INFO, INS, IRANK,
230  $ iscale, itran, itype, j, k, lda, ldb, ldwork,
231  $ lwlsy, lwork, m, mnmin, n, nb, ncols, nerrs,
232  $ nfail, nrhs, nrows, nrun, rank, mb,
233  $ mmax, nmax, nsmax, liwork, lrwork,
234  $ lwork_zgels, lwork_zgetsls, lwork_zgelss,
235  $ lwork_zgelsy, lwork_zgelsd,
236  $ lrwork_zgelsy, lrwork_zgelss, lrwork_zgelsd
237  DOUBLE PRECISION EPS, NORMA, NORMB, RCOND
238 * ..
239 * .. Local Arrays ..
240  INTEGER ISEED( 4 ), ISEEDY( 4 ), IWQ( 1 )
241  DOUBLE PRECISION RESULT( NTESTS ), RWQ( 1 )
242  COMPLEX*16 WQ( 1 )
243 * ..
244 * .. Allocatable Arrays ..
245  COMPLEX*16, ALLOCATABLE :: WORK (:)
246  DOUBLE PRECISION, ALLOCATABLE :: RWORK (:), WORK2 (:)
247  INTEGER, ALLOCATABLE :: IWORK (:)
248 * ..
249 * .. External Functions ..
250  DOUBLE PRECISION DASUM, DLAMCH, ZQRT12, ZQRT14, ZQRT17
251  EXTERNAL DASUM, DLAMCH, ZQRT12, ZQRT14, ZQRT17
252 * ..
253 * .. External Subroutines ..
254  EXTERNAL alaerh, alahd, alasvm, daxpy, dlasrt, xlaenv,
255  $ zdscal, zerrls, zgels, zgelsd, zgelss,
256  $ zgelsy, zgemm, zlacpy, zlarnv, zqrt13, zqrt15,
257  $ zqrt16, zgetsls
258 * ..
259 * .. Intrinsic Functions ..
260  INTRINSIC dble, max, min, int, sqrt
261 * ..
262 * .. Scalars in Common ..
263  LOGICAL LERR, OK
264  CHARACTER*32 SRNAMT
265  INTEGER INFOT, IOUNIT
266 * ..
267 * .. Common blocks ..
268  COMMON / infoc / infot, iounit, ok, lerr
269  COMMON / srnamc / srnamt
270 * ..
271 * .. Data statements ..
272  DATA iseedy / 1988, 1989, 1990, 1991 /
273 * ..
274 * .. Executable Statements ..
275 *
276 * Initialize constants and the random number seed.
277 *
278  path( 1: 1 ) = 'Zomplex precision'
279  path( 2: 3 ) = 'LS'
280  nrun = 0
281  nfail = 0
282  nerrs = 0
283  DO 10 i = 1, 4
284  iseed( i ) = iseedy( i )
285  10 CONTINUE
286  eps = dlamch( 'Epsilon' )
287 *
288 * Threshold for rank estimation
289 *
290  rcond = sqrt( eps ) - ( sqrt( eps )-eps ) / 2
291 *
292 * Test the error exits
293 *
294  CALL xlaenv( 9, smlsiz )
295  IF( tsterr )
296  $ CALL zerrls( path, nout )
297 *
298 * Print the header if NM = 0 or NN = 0 and THRESH = 0.
299 *
300  IF( ( nm.EQ.0 .OR. nn.EQ.0 ) .AND. thresh.EQ.zero )
301  $ CALL alahd( nout, path )
302  infot = 0
303 *
304 * Compute maximal workspace needed for all routines
305 *
306  nmax = 0
307  mmax = 0
308  nsmax = 0
309  DO i = 1, nm
310  IF ( mval( i ).GT.mmax ) THEN
311  mmax = mval( i )
312  END IF
313  ENDDO
314  DO i = 1, nn
315  IF ( nval( i ).GT.nmax ) THEN
316  nmax = nval( i )
317  END IF
318  ENDDO
319  DO i = 1, nns
320  IF ( nsval( i ).GT.nsmax ) THEN
321  nsmax = nsval( i )
322  END IF
323  ENDDO
324  m = mmax
325  n = nmax
326  nrhs = nsmax
327  mnmin = max( min( m, n ), 1 )
328 *
329 * Compute workspace needed for routines
330 * ZQRT14, ZQRT17 (two side cases), ZQRT15 and ZQRT12
331 *
332  lwork = max( 1, ( m+n )*nrhs,
333  $ ( n+nrhs )*( m+2 ), ( m+nrhs )*( n+2 ),
334  $ max( m+mnmin, nrhs*mnmin,2*n+m ),
335  $ max( m*n+4*mnmin+max(m,n), m*n+2*mnmin+4*n ) )
336  lrwork = 1
337  liwork = 1
338 *
339 * Iterate through all test cases and compute necessary workspace
340 * sizes for ?GELS, ?GETSLS, ?GELSY, ?GELSS and ?GELSD routines.
341 *
342  DO im = 1, nm
343  m = mval( im )
344  lda = max( 1, m )
345  DO in = 1, nn
346  n = nval( in )
347  mnmin = max(min( m, n ),1)
348  ldb = max( 1, m, n )
349  DO ins = 1, nns
350  nrhs = nsval( ins )
351  DO irank = 1, 2
352  DO iscale = 1, 3
353  itype = ( irank-1 )*3 + iscale
354  IF( dotype( itype ) ) THEN
355  IF( irank.EQ.1 ) THEN
356  DO itran = 1, 2
357  IF( itran.EQ.1 ) THEN
358  trans = 'N'
359  ELSE
360  trans = 'C'
361  END IF
362 *
363 * Compute workspace needed for ZGELS
364  CALL zgels( trans, m, n, nrhs, a, lda,
365  $ b, ldb, wq, -1, info )
366  lwork_zgels = int( wq( 1 ) )
367 * Compute workspace needed for ZGETSLS
368  CALL zgetsls( trans, m, n, nrhs, a, lda,
369  $ b, ldb, wq, -1, info )
370  lwork_zgetsls = int( wq( 1 ) )
371  ENDDO
372  END IF
373 * Compute workspace needed for ZGELSY
374  CALL zgelsy( m, n, nrhs, a, lda, b, ldb, iwq,
375  $ rcond, crank, wq, -1, rwork, info )
376  lwork_zgelsy = int( wq( 1 ) )
377  lrwork_zgelsy = 2*n
378 * Compute workspace needed for ZGELSS
379  CALL zgelss( m, n, nrhs, a, lda, b, ldb, s,
380  $ rcond, crank, wq, -1 , rwork,
381  $ info )
382  lwork_zgelss = int( wq( 1 ) )
383  lrwork_zgelss = 5*mnmin
384 * Compute workspace needed for ZGELSD
385  CALL zgelsd( m, n, nrhs, a, lda, b, ldb, s,
386  $ rcond, crank, wq, -1, rwq, iwq,
387  $ info )
388  lwork_zgelsd = int( wq( 1 ) )
389  lrwork_zgelsd = int( rwq( 1 ) )
390 * Compute LIWORK workspace needed for ZGELSY and ZGELSD
391  liwork = max( liwork, n, iwq( 1 ) )
392 * Compute LRWORK workspace needed for ZGELSY, ZGELSS and ZGELSD
393  lrwork = max( lrwork, lrwork_zgelsy,
394  $ lrwork_zgelss, lrwork_zgelsd )
395 * Compute LWORK workspace needed for all functions
396  lwork = max( lwork, lwork_zgels, lwork_zgetsls,
397  $ lwork_zgelsy, lwork_zgelss,
398  $ lwork_zgelsd )
399  END IF
400  ENDDO
401  ENDDO
402  ENDDO
403  ENDDO
404  ENDDO
405 *
406  lwlsy = lwork
407 *
408  ALLOCATE( work( lwork ) )
409  ALLOCATE( work2( 2 * lwork ) )
410  ALLOCATE( iwork( liwork ) )
411  ALLOCATE( rwork( lrwork ) )
412 *
413  DO 140 im = 1, nm
414  m = mval( im )
415  lda = max( 1, m )
416 *
417  DO 130 in = 1, nn
418  n = nval( in )
419  mnmin = max(min( m, n ),1)
420  ldb = max( 1, m, n )
421  mb = (mnmin+1)
422 *
423  DO 120 ins = 1, nns
424  nrhs = nsval( ins )
425 *
426  DO 110 irank = 1, 2
427  DO 100 iscale = 1, 3
428  itype = ( irank-1 )*3 + iscale
429  IF( .NOT.dotype( itype ) )
430  $ GO TO 100
431 *
432  IF( irank.EQ.1 ) THEN
433 *
434 * Test ZGELS
435 *
436 * Generate a matrix of scaling type ISCALE
437 *
438  CALL zqrt13( iscale, m, n, copya, lda, norma,
439  $ iseed )
440  DO 40 inb = 1, nnb
441  nb = nbval( inb )
442  CALL xlaenv( 1, nb )
443  CALL xlaenv( 3, nxval( inb ) )
444 *
445  DO 30 itran = 1, 2
446  IF( itran.EQ.1 ) THEN
447  trans = 'N'
448  nrows = m
449  ncols = n
450  ELSE
451  trans = 'C'
452  nrows = n
453  ncols = m
454  END IF
455  ldwork = max( 1, ncols )
456 *
457 * Set up a consistent rhs
458 *
459  IF( ncols.GT.0 ) THEN
460  CALL zlarnv( 2, iseed, ncols*nrhs,
461  $ work )
462  CALL zdscal( ncols*nrhs,
463  $ one / dble( ncols ), work,
464  $ 1 )
465  END IF
466  CALL zgemm( trans, 'No transpose', nrows,
467  $ nrhs, ncols, cone, copya, lda,
468  $ work, ldwork, czero, b, ldb )
469  CALL zlacpy( 'Full', nrows, nrhs, b, ldb,
470  $ copyb, ldb )
471 *
472 * Solve LS or overdetermined system
473 *
474  IF( m.GT.0 .AND. n.GT.0 ) THEN
475  CALL zlacpy( 'Full', m, n, copya, lda,
476  $ a, lda )
477  CALL zlacpy( 'Full', nrows, nrhs,
478  $ copyb, ldb, b, ldb )
479  END IF
480  srnamt = 'ZGELS '
481  CALL zgels( trans, m, n, nrhs, a, lda, b,
482  $ ldb, work, lwork, info )
483 *
484  IF( info.NE.0 )
485  $ CALL alaerh( path, 'ZGELS ', info, 0,
486  $ trans, m, n, nrhs, -1, nb,
487  $ itype, nfail, nerrs,
488  $ nout )
489 *
490 * Check correctness of results
491 *
492  ldwork = max( 1, nrows )
493  IF( nrows.GT.0 .AND. nrhs.GT.0 )
494  $ CALL zlacpy( 'Full', nrows, nrhs,
495  $ copyb, ldb, c, ldb )
496  CALL zqrt16( trans, m, n, nrhs, copya,
497  $ lda, b, ldb, c, ldb, rwork,
498  $ result( 1 ) )
499 *
500  IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
501  $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
502 *
503 * Solving LS system
504 *
505  result( 2 ) = zqrt17( trans, 1, m, n,
506  $ nrhs, copya, lda, b, ldb,
507  $ copyb, ldb, c, work,
508  $ lwork )
509  ELSE
510 *
511 * Solving overdetermined system
512 *
513  result( 2 ) = zqrt14( trans, m, n,
514  $ nrhs, copya, lda, b, ldb,
515  $ work, lwork )
516  END IF
517 *
518 * Print information about the tests that
519 * did not pass the threshold.
520 *
521  DO 20 k = 1, 2
522  IF( result( k ).GE.thresh ) THEN
523  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
524  $ CALL alahd( nout, path )
525  WRITE( nout, fmt = 9999 )trans, m,
526  $ n, nrhs, nb, itype, k,
527  $ result( k )
528  nfail = nfail + 1
529  END IF
530  20 CONTINUE
531  nrun = nrun + 2
532  30 CONTINUE
533  40 CONTINUE
534 *
535 *
536 * Test ZGETSLS
537 *
538 * Generate a matrix of scaling type ISCALE
539 *
540  CALL zqrt13( iscale, m, n, copya, lda, norma,
541  $ iseed )
542  DO 65 inb = 1, nnb
543  mb = nbval( inb )
544  CALL xlaenv( 1, mb )
545  DO 62 imb = 1, nnb
546  nb = nbval( imb )
547  CALL xlaenv( 2, nb )
548 *
549  DO 60 itran = 1, 2
550  IF( itran.EQ.1 ) THEN
551  trans = 'N'
552  nrows = m
553  ncols = n
554  ELSE
555  trans = 'C'
556  nrows = n
557  ncols = m
558  END IF
559  ldwork = max( 1, ncols )
560 *
561 * Set up a consistent rhs
562 *
563  IF( ncols.GT.0 ) THEN
564  CALL zlarnv( 2, iseed, ncols*nrhs,
565  $ work )
566  CALL zscal( ncols*nrhs,
567  $ one / dble( ncols ), work,
568  $ 1 )
569  END IF
570  CALL zgemm( trans, 'No transpose', nrows,
571  $ nrhs, ncols, cone, copya, lda,
572  $ work, ldwork, czero, b, ldb )
573  CALL zlacpy( 'Full', nrows, nrhs, b, ldb,
574  $ copyb, ldb )
575 *
576 * Solve LS or overdetermined system
577 *
578  IF( m.GT.0 .AND. n.GT.0 ) THEN
579  CALL zlacpy( 'Full', m, n, copya, lda,
580  $ a, lda )
581  CALL zlacpy( 'Full', nrows, nrhs,
582  $ copyb, ldb, b, ldb )
583  END IF
584  srnamt = 'ZGETSLS '
585  CALL zgetsls( trans, m, n, nrhs, a,
586  $ lda, b, ldb, work, lwork, info )
587  IF( info.NE.0 )
588  $ CALL alaerh( path, 'ZGETSLS ', info, 0,
589  $ trans, m, n, nrhs, -1, nb,
590  $ itype, nfail, nerrs,
591  $ nout )
592 *
593 * Check correctness of results
594 *
595  ldwork = max( 1, nrows )
596  IF( nrows.GT.0 .AND. nrhs.GT.0 )
597  $ CALL zlacpy( 'Full', nrows, nrhs,
598  $ copyb, ldb, c, ldb )
599  CALL zqrt16( trans, m, n, nrhs, copya,
600  $ lda, b, ldb, c, ldb, work2,
601  $ result( 15 ) )
602 *
603  IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
604  $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
605 *
606 * Solving LS system
607 *
608  result( 16 ) = zqrt17( trans, 1, m, n,
609  $ nrhs, copya, lda, b, ldb,
610  $ copyb, ldb, c, work,
611  $ lwork )
612  ELSE
613 *
614 * Solving overdetermined system
615 *
616  result( 16 ) = zqrt14( trans, m, n,
617  $ nrhs, copya, lda, b, ldb,
618  $ work, lwork )
619  END IF
620 *
621 * Print information about the tests that
622 * did not pass the threshold.
623 *
624  DO 50 k = 15, 16
625  IF( result( k ).GE.thresh ) THEN
626  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
627  $ CALL alahd( nout, path )
628  WRITE( nout, fmt = 9997 )trans, m,
629  $ n, nrhs, mb, nb, itype, k,
630  $ result( k )
631  nfail = nfail + 1
632  END IF
633  50 CONTINUE
634  nrun = nrun + 2
635  60 CONTINUE
636  62 CONTINUE
637  65 CONTINUE
638  END IF
639 *
640 * Generate a matrix of scaling type ISCALE and rank
641 * type IRANK.
642 *
643  CALL zqrt15( iscale, irank, m, n, nrhs, copya, lda,
644  $ copyb, ldb, copys, rank, norma, normb,
645  $ iseed, work, lwork )
646 *
647 * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
648 *
649  ldwork = max( 1, m )
650 *
651 * Loop for testing different block sizes.
652 *
653  DO 90 inb = 1, nnb
654  nb = nbval( inb )
655  CALL xlaenv( 1, nb )
656  CALL xlaenv( 3, nxval( inb ) )
657 *
658 * Test ZGELSY
659 *
660 * ZGELSY: Compute the minimum-norm solution
661 * X to min( norm( A * X - B ) )
662 * using the rank-revealing orthogonal
663 * factorization.
664 *
665  CALL zlacpy( 'Full', m, n, copya, lda, a, lda )
666  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, b,
667  $ ldb )
668 *
669 * Initialize vector IWORK.
670 *
671  DO 70 j = 1, n
672  iwork( j ) = 0
673  70 CONTINUE
674 *
675  srnamt = 'ZGELSY'
676  CALL zgelsy( m, n, nrhs, a, lda, b, ldb, iwork,
677  $ rcond, crank, work, lwlsy, rwork,
678  $ info )
679  IF( info.NE.0 )
680  $ CALL alaerh( path, 'ZGELSY', info, 0, ' ', m,
681  $ n, nrhs, -1, nb, itype, nfail,
682  $ nerrs, nout )
683 *
684 * workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS)
685 *
686 * Test 3: Compute relative error in svd
687 * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
688 *
689  result( 3 ) = zqrt12( crank, crank, a, lda,
690  $ copys, work, lwork, rwork )
691 *
692 * Test 4: Compute error in solution
693 * workspace: M*NRHS + M
694 *
695  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, work,
696  $ ldwork )
697  CALL zqrt16( 'No transpose', m, n, nrhs, copya,
698  $ lda, b, ldb, work, ldwork, rwork,
699  $ result( 4 ) )
700 *
701 * Test 5: Check norm of r'*A
702 * workspace: NRHS*(M+N)
703 *
704  result( 5 ) = zero
705  IF( m.GT.crank )
706  $ result( 5 ) = zqrt17( 'No transpose', 1, m,
707  $ n, nrhs, copya, lda, b, ldb,
708  $ copyb, ldb, c, work, lwork )
709 *
710 * Test 6: Check if x is in the rowspace of A
711 * workspace: (M+NRHS)*(N+2)
712 *
713  result( 6 ) = zero
714 *
715  IF( n.GT.crank )
716  $ result( 6 ) = zqrt14( 'No transpose', m, n,
717  $ nrhs, copya, lda, b, ldb,
718  $ work, lwork )
719 *
720 * Test ZGELSS
721 *
722 * ZGELSS: Compute the minimum-norm solution
723 * X to min( norm( A * X - B ) )
724 * using the SVD.
725 *
726  CALL zlacpy( 'Full', m, n, copya, lda, a, lda )
727  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, b,
728  $ ldb )
729  srnamt = 'ZGELSS'
730  CALL zgelss( m, n, nrhs, a, lda, b, ldb, s,
731  $ rcond, crank, work, lwork, rwork,
732  $ info )
733 *
734  IF( info.NE.0 )
735  $ CALL alaerh( path, 'ZGELSS', info, 0, ' ', m,
736  $ n, nrhs, -1, nb, itype, nfail,
737  $ nerrs, nout )
738 *
739 * workspace used: 3*min(m,n) +
740 * max(2*min(m,n),nrhs,max(m,n))
741 *
742 * Test 7: Compute relative error in svd
743 *
744  IF( rank.GT.0 ) THEN
745  CALL daxpy( mnmin, -one, copys, 1, s, 1 )
746  result( 7 ) = dasum( mnmin, s, 1 ) /
747  $ dasum( mnmin, copys, 1 ) /
748  $ ( eps*dble( mnmin ) )
749  ELSE
750  result( 7 ) = zero
751  END IF
752 *
753 * Test 8: Compute error in solution
754 *
755  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, work,
756  $ ldwork )
757  CALL zqrt16( 'No transpose', m, n, nrhs, copya,
758  $ lda, b, ldb, work, ldwork, rwork,
759  $ result( 8 ) )
760 *
761 * Test 9: Check norm of r'*A
762 *
763  result( 9 ) = zero
764  IF( m.GT.crank )
765  $ result( 9 ) = zqrt17( 'No transpose', 1, m,
766  $ n, nrhs, copya, lda, b, ldb,
767  $ copyb, ldb, c, work, lwork )
768 *
769 * Test 10: Check if x is in the rowspace of A
770 *
771  result( 10 ) = zero
772  IF( n.GT.crank )
773  $ result( 10 ) = zqrt14( 'No transpose', m, n,
774  $ nrhs, copya, lda, b, ldb,
775  $ work, lwork )
776 *
777 * Test ZGELSD
778 *
779 * ZGELSD: Compute the minimum-norm solution X
780 * to min( norm( A * X - B ) ) using a
781 * divide and conquer SVD.
782 *
783  CALL xlaenv( 9, 25 )
784 *
785  CALL zlacpy( 'Full', m, n, copya, lda, a, lda )
786  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, b,
787  $ ldb )
788 *
789  srnamt = 'ZGELSD'
790  CALL zgelsd( m, n, nrhs, a, lda, b, ldb, s,
791  $ rcond, crank, work, lwork, rwork,
792  $ iwork, info )
793  IF( info.NE.0 )
794  $ CALL alaerh( path, 'ZGELSD', info, 0, ' ', m,
795  $ n, nrhs, -1, nb, itype, nfail,
796  $ nerrs, nout )
797 *
798 * Test 11: Compute relative error in svd
799 *
800  IF( rank.GT.0 ) THEN
801  CALL daxpy( mnmin, -one, copys, 1, s, 1 )
802  result( 11 ) = dasum( mnmin, s, 1 ) /
803  $ dasum( mnmin, copys, 1 ) /
804  $ ( eps*dble( mnmin ) )
805  ELSE
806  result( 11 ) = zero
807  END IF
808 *
809 * Test 12: Compute error in solution
810 *
811  CALL zlacpy( 'Full', m, nrhs, copyb, ldb, work,
812  $ ldwork )
813  CALL zqrt16( 'No transpose', m, n, nrhs, copya,
814  $ lda, b, ldb, work, ldwork, rwork,
815  $ result( 12 ) )
816 *
817 * Test 13: Check norm of r'*A
818 *
819  result( 13 ) = zero
820  IF( m.GT.crank )
821  $ result( 13 ) = zqrt17( 'No transpose', 1, m,
822  $ n, nrhs, copya, lda, b, ldb,
823  $ copyb, ldb, c, work, lwork )
824 *
825 * Test 14: Check if x is in the rowspace of A
826 *
827  result( 14 ) = zero
828  IF( n.GT.crank )
829  $ result( 14 ) = zqrt14( 'No transpose', m, n,
830  $ nrhs, copya, lda, b, ldb,
831  $ work, lwork )
832 *
833 * Print information about the tests that did not
834 * pass the threshold.
835 *
836  DO 80 k = 3, 14
837  IF( result( k ).GE.thresh ) THEN
838  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
839  $ CALL alahd( nout, path )
840  WRITE( nout, fmt = 9998 )m, n, nrhs, nb,
841  $ itype, k, result( k )
842  nfail = nfail + 1
843  END IF
844  80 CONTINUE
845  nrun = nrun + 12
846 *
847  90 CONTINUE
848  100 CONTINUE
849  110 CONTINUE
850  120 CONTINUE
851  130 CONTINUE
852  140 CONTINUE
853 *
854 * Print a summary of the results.
855 *
856  CALL alasvm( path, nout, nfail, nrun, nerrs )
857 *
858  9999 FORMAT( ' TRANS=''', a1, ''', M=', i5, ', N=', i5, ', NRHS=', i4,
859  $ ', NB=', i4, ', type', i2, ', test(', i2, ')=', g12.5 )
860  9998 FORMAT( ' M=', i5, ', N=', i5, ', NRHS=', i4, ', NB=', i4,
861  $ ', type', i2, ', test(', i2, ')=', g12.5 )
862  9997 FORMAT( ' TRANS=''', a1,' M=', i5, ', N=', i5, ', NRHS=', i4,
863  $ ', MB=', i4,', NB=', i4,', type', i2,
864  $ ', test(', i2, ')=', g12.5 )
865 *
866  DEALLOCATE( work )
867  DEALLOCATE( iwork )
868  DEALLOCATE( rwork )
869  RETURN
870 *
871 * End of ZDRVLS
872 *
873  END
zerrls
subroutine zerrls(PATH, NUNIT)
ZERRLS
Definition: zerrls.f:57
zgels
subroutine zgels(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
ZGELS solves overdetermined or underdetermined systems for GE matrices
Definition: zgels.f:184
zqrt13
subroutine zqrt13(SCALE, M, N, A, LDA, NORMA, ISEED)
ZQRT13
Definition: zqrt13.f:93
alasvm
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
zdrvls
subroutine zdrvls(DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB, NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B, COPYB, C, S, COPYS, NOUT)
ZDRVLS
Definition: zdrvls.f:194
zgelsy
subroutine zgelsy(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, LWORK, RWORK, INFO)
ZGELSY solves overdetermined or underdetermined systems for GE matrices
Definition: zgelsy.f:212
alahd
subroutine alahd(IOUNIT, PATH)
ALAHD
Definition: alahd.f:109
zdscal
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:80
zgelsd
subroutine zgelsd(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, IWORK, INFO)
ZGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices
Definition: zgelsd.f:227
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
zgetsls
subroutine zgetsls(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
ZGETSLS
Definition: zgetsls.f:164
zgelss
subroutine zgelss(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO)
ZGELSS solves overdetermined or underdetermined systems for GE matrices
Definition: zgelss.f:180
zlacpy
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
alaerh
subroutine alaerh(PATH, SUBNAM, INFO, INFOE, OPTS, M, N, KL, KU, N5, IMAT, NFAIL, NERRS, NOUT)
ALAERH
Definition: alaerh.f:149
zqrt15
subroutine zqrt15(SCALE, RKSEL, M, N, NRHS, A, LDA, B, LDB, S, RANK, NORMA, NORMB, ISEED, WORK, LWORK)
ZQRT15
Definition: zqrt15.f:151
zqrt16
subroutine zqrt16(TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, RWORK, RESID)
ZQRT16
Definition: zqrt16.f:135
zlarnv
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: zlarnv.f:101
zscal
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:80
dlasrt
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:90
xlaenv
subroutine xlaenv(ISPEC, NVALUE)
XLAENV
Definition: xlaenv.f:83
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:91