LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
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, ZGELST,
35*> ZGETSLS, ZGELSS, ZGELSY 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*> \ingroup complex16_lin
187*
188* =====================================================================
189 SUBROUTINE zdrvls( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
190 \$ NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
191 \$ COPYB, C, S, COPYS, NOUT )
192*
193* -- LAPACK test routine --
194* -- LAPACK is a software package provided by Univ. of Tennessee, --
195* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
196*
197* .. Scalar Arguments ..
198 LOGICAL TSTERR
199 INTEGER NM, NN, NNB, NNS, NOUT
200 DOUBLE PRECISION THRESH
201* ..
202* .. Array Arguments ..
203 LOGICAL DOTYPE( * )
204 INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
205 \$ nval( * ), nxval( * )
206 DOUBLE PRECISION COPYS( * ), S( * )
207 COMPLEX*16 A( * ), B( * ), C( * ), COPYA( * ), COPYB( * )
208* ..
209*
210* =====================================================================
211*
212* .. Parameters ..
213 INTEGER NTESTS
214 PARAMETER ( NTESTS = 18 )
215 INTEGER SMLSIZ
216 parameter( smlsiz = 25 )
217 DOUBLE PRECISION ONE, ZERO
218 parameter( one = 1.0d+0, zero = 0.0d+0 )
219 COMPLEX*16 CONE, CZERO
220 parameter( cone = ( 1.0d+0, 0.0d+0 ),
221 \$ czero = ( 0.0d+0, 0.0d+0 ) )
222* ..
223* .. Local Scalars ..
224 CHARACTER TRANS
225 CHARACTER*3 PATH
226 INTEGER CRANK, I, IM, IMB, IN, INB, INFO, INS, IRANK,
227 \$ iscale, itran, itype, j, k, lda, ldb, ldwork,
228 \$ lwlsy, lwork, m, mnmin, n, nb, ncols, nerrs,
229 \$ nfail, nrhs, nrows, nrun, rank, mb,
230 \$ mmax, nmax, nsmax, liwork, lrwork,
231 \$ lwork_zgels, lwork_zgelst, lwork_zgetsls,
232 \$ lwork_zgelss, lwork_zgelsy, lwork_zgelsd,
233 \$ lrwork_zgelsy, lrwork_zgelss, lrwork_zgelsd
234 DOUBLE PRECISION EPS, NORMA, NORMB, RCOND
235* ..
236* .. Local Arrays ..
237 INTEGER ISEED( 4 ), ISEEDY( 4 ), IWQ( 1 )
238 DOUBLE PRECISION RESULT( NTESTS ), RWQ( 1 )
239 COMPLEX*16 WQ( 1 )
240* ..
241* .. Allocatable Arrays ..
242 COMPLEX*16, ALLOCATABLE :: WORK (:)
243 DOUBLE PRECISION, ALLOCATABLE :: RWORK (:), WORK2 (:)
244 INTEGER, ALLOCATABLE :: IWORK (:)
245* ..
246* .. External Functions ..
247 DOUBLE PRECISION DASUM, DLAMCH, ZQRT12, ZQRT14, ZQRT17
248 EXTERNAL DASUM, DLAMCH, ZQRT12, ZQRT14, ZQRT17
249* ..
250* .. External Subroutines ..
251 EXTERNAL alaerh, alahd, alasvm, daxpy, zerrls, zgels,
255* ..
256* .. Intrinsic Functions ..
257 INTRINSIC dble, max, min, int, sqrt
258* ..
259* .. Scalars in Common ..
260 LOGICAL LERR, OK
261 CHARACTER*32 SRNAMT
262 INTEGER INFOT, IOUNIT
263* ..
264* .. Common blocks ..
265 COMMON / infoc / infot, iounit, ok, lerr
266 COMMON / srnamc / srnamt
267* ..
268* .. Data statements ..
269 DATA iseedy / 1988, 1989, 1990, 1991 /
270* ..
271* .. Executable Statements ..
272*
273* Initialize constants and the random number seed.
274*
275 path( 1: 1 ) = 'Zomplex precision'
276 path( 2: 3 ) = 'LS'
277 nrun = 0
278 nfail = 0
279 nerrs = 0
280 DO 10 i = 1, 4
281 iseed( i ) = iseedy( i )
282 10 CONTINUE
283 eps = dlamch( 'Epsilon' )
284*
285* Threshold for rank estimation
286*
287 rcond = sqrt( eps ) - ( sqrt( eps )-eps ) / 2
288*
289* Test the error exits
290*
291 CALL xlaenv( 9, smlsiz )
292 IF( tsterr )
293 \$ CALL zerrls( path, nout )
294*
295* Print the header if NM = 0 or NN = 0 and THRESH = 0.
296*
297 IF( ( nm.EQ.0 .OR. nn.EQ.0 ) .AND. thresh.EQ.zero )
298 \$ CALL alahd( nout, path )
299 infot = 0
300*
301* Compute maximal workspace needed for all routines
302*
303 nmax = 0
304 mmax = 0
305 nsmax = 0
306 DO i = 1, nm
307 IF ( mval( i ).GT.mmax ) THEN
308 mmax = mval( i )
309 END IF
310 ENDDO
311 DO i = 1, nn
312 IF ( nval( i ).GT.nmax ) THEN
313 nmax = nval( i )
314 END IF
315 ENDDO
316 DO i = 1, nns
317 IF ( nsval( i ).GT.nsmax ) THEN
318 nsmax = nsval( i )
319 END IF
320 ENDDO
321 m = mmax
322 n = nmax
323 nrhs = nsmax
324 mnmin = max( min( m, n ), 1 )
325*
326* Compute workspace needed for routines
327* ZQRT14, ZQRT17 (two side cases), ZQRT15 and ZQRT12
328*
329 lwork = max( 1, ( m+n )*nrhs,
330 \$ ( n+nrhs )*( m+2 ), ( m+nrhs )*( n+2 ),
331 \$ max( m+mnmin, nrhs*mnmin,2*n+m ),
332 \$ max( m*n+4*mnmin+max(m,n), m*n+2*mnmin+4*n ) )
333 lrwork = 1
334 liwork = 1
335*
336* Iterate through all test cases and compute necessary workspace
337* sizes for ?GELS, ?GELST, ?GETSLS, ?GELSY, ?GELSS and ?GELSD
338* routines.
339*
340 DO im = 1, nm
341 m = mval( im )
342 lda = max( 1, m )
343 DO in = 1, nn
344 n = nval( in )
345 mnmin = max(min( m, n ),1)
346 ldb = max( 1, m, n )
347 DO ins = 1, nns
348 nrhs = nsval( ins )
349 DO irank = 1, 2
350 DO iscale = 1, 3
351 itype = ( irank-1 )*3 + iscale
352 IF( dotype( itype ) ) THEN
353 IF( irank.EQ.1 ) THEN
354 DO itran = 1, 2
355 IF( itran.EQ.1 ) THEN
356 trans = 'N'
357 ELSE
358 trans = 'C'
359 END IF
360*
361* Compute workspace needed for ZGELS
362 CALL zgels( trans, m, n, nrhs, a, lda,
363 \$ b, ldb, wq, -1, info )
364 lwork_zgels = int( wq( 1 ) )
365* Compute workspace needed for ZGELST
366 CALL zgelst( trans, m, n, nrhs, a, lda,
367 \$ b, ldb, wq, -1, info )
368 lwork_zgelst = int( wq( 1 ) )
369* Compute workspace needed for ZGETSLS
370 CALL zgetsls( trans, m, n, nrhs, a, lda,
371 \$ b, ldb, wq, -1, info )
372 lwork_zgetsls = int( wq( 1 ) )
373 ENDDO
374 END IF
375* Compute workspace needed for ZGELSY
376 CALL zgelsy( m, n, nrhs, a, lda, b, ldb, iwq,
377 \$ rcond, crank, wq, -1, rwq, info )
378 lwork_zgelsy = int( wq( 1 ) )
379 lrwork_zgelsy = 2*n
380* Compute workspace needed for ZGELSS
381 CALL zgelss( m, n, nrhs, a, lda, b, ldb, s,
382 \$ rcond, crank, wq, -1 , rwq,
383 \$ info )
384 lwork_zgelss = int( wq( 1 ) )
385 lrwork_zgelss = 5*mnmin
386* Compute workspace needed for ZGELSD
387 CALL zgelsd( m, n, nrhs, a, lda, b, ldb, s,
388 \$ rcond, crank, wq, -1, rwq, iwq,
389 \$ info )
390 lwork_zgelsd = int( wq( 1 ) )
391 lrwork_zgelsd = int( rwq( 1 ) )
392* Compute LIWORK workspace needed for ZGELSY and ZGELSD
393 liwork = max( liwork, n, iwq( 1 ) )
394* Compute LRWORK workspace needed for ZGELSY, ZGELSS and ZGELSD
395 lrwork = max( lrwork, lrwork_zgelsy,
396 \$ lrwork_zgelss, lrwork_zgelsd )
397* Compute LWORK workspace needed for all functions
398 lwork = max( lwork, lwork_zgels, lwork_zgelst,
399 \$ lwork_zgetsls, lwork_zgelsy,
400 \$ lwork_zgelss, lwork_zgelsd )
401 END IF
402 ENDDO
403 ENDDO
404 ENDDO
405 ENDDO
406 ENDDO
407*
408 lwlsy = lwork
409*
410 ALLOCATE( work( lwork ) )
411 ALLOCATE( work2( 2 * lwork ) )
412 ALLOCATE( iwork( liwork ) )
413 ALLOCATE( rwork( lrwork ) )
414*
415 DO 140 im = 1, nm
416 m = mval( im )
417 lda = max( 1, m )
418*
419 DO 130 in = 1, nn
420 n = nval( in )
421 mnmin = max(min( m, n ),1)
422 ldb = max( 1, m, n )
423 mb = (mnmin+1)
424*
425 DO 120 ins = 1, nns
426 nrhs = nsval( ins )
427*
428 DO 110 irank = 1, 2
429 DO 100 iscale = 1, 3
430 itype = ( irank-1 )*3 + iscale
431 IF( .NOT.dotype( itype ) )
432 \$ GO TO 100
433* =====================================================
434* Begin test ZGELS
435* =====================================================
436 IF( irank.EQ.1 ) THEN
437*
438* Generate a matrix of scaling type ISCALE
439*
440 CALL zqrt13( iscale, m, n, copya, lda, norma,
441 \$ iseed )
442*
443* Loop for testing different block sizes.
444*
445 DO inb = 1, nnb
446 nb = nbval( inb )
447 CALL xlaenv( 1, nb )
448 CALL xlaenv( 3, nxval( inb ) )
449*
450* Loop for testing non-transposed and transposed.
451*
452 DO itran = 1, 2
453 IF( itran.EQ.1 ) THEN
454 trans = 'N'
455 nrows = m
456 ncols = n
457 ELSE
458 trans = 'C'
459 nrows = n
460 ncols = m
461 END IF
462 ldwork = max( 1, ncols )
463*
464* Set up a consistent rhs
465*
466 IF( ncols.GT.0 ) THEN
467 CALL zlarnv( 2, iseed, ncols*nrhs,
468 \$ work )
469 CALL zdscal( ncols*nrhs,
470 \$ one / dble( ncols ), work,
471 \$ 1 )
472 END IF
473 CALL zgemm( trans, 'No transpose', nrows,
474 \$ nrhs, ncols, cone, copya, lda,
475 \$ work, ldwork, czero, b, ldb )
476 CALL zlacpy( 'Full', nrows, nrhs, b, ldb,
477 \$ copyb, ldb )
478*
479* Solve LS or overdetermined system
480*
481 IF( m.GT.0 .AND. n.GT.0 ) THEN
482 CALL zlacpy( 'Full', m, n, copya, lda,
483 \$ a, lda )
484 CALL zlacpy( 'Full', nrows, nrhs,
485 \$ copyb, ldb, b, ldb )
486 END IF
487 srnamt = 'ZGELS '
488 CALL zgels( trans, m, n, nrhs, a, lda, b,
489 \$ ldb, work, lwork, info )
490*
491 IF( info.NE.0 )
492 \$ CALL alaerh( path, 'ZGELS ', info, 0,
493 \$ trans, m, n, nrhs, -1, nb,
494 \$ itype, nfail, nerrs,
495 \$ nout )
496*
497* Test 1: Check correctness of results
498* for ZGELS, compute the residual:
499* RESID = norm(B - A*X) /
500* / ( max(m,n) * norm(A) * norm(X) * EPS )
501*
502 IF( nrows.GT.0 .AND. nrhs.GT.0 )
503 \$ CALL zlacpy( 'Full', nrows, nrhs,
504 \$ copyb, ldb, c, ldb )
505 CALL zqrt16( trans, m, n, nrhs, copya,
506 \$ lda, b, ldb, c, ldb, rwork,
507 \$ result( 1 ) )
508*
509* Test 2: Check correctness of results
510* for ZGELS.
511*
512 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
513 \$ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
514*
515* Solving LS system
516*
517 result( 2 ) = zqrt17( trans, 1, m, n,
518 \$ nrhs, copya, lda, b, ldb,
519 \$ copyb, ldb, c, work,
520 \$ lwork )
521 ELSE
522*
523* Solving overdetermined system
524*
525 result( 2 ) = zqrt14( trans, m, n,
526 \$ nrhs, copya, lda, b, ldb,
527 \$ work, lwork )
528 END IF
529*
530* Print information about the tests that
531* did not pass the threshold.
532*
533 DO k = 1, 2
534 IF( result( k ).GE.thresh ) THEN
535 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
536 \$ CALL alahd( nout, path )
537 WRITE( nout, fmt = 9999 )trans, m,
538 \$ n, nrhs, nb, itype, k,
539 \$ result( k )
540 nfail = nfail + 1
541 END IF
542 END DO
543 nrun = nrun + 2
544 END DO
545 END DO
546 END IF
547* =====================================================
548* End test ZGELS
549* =====================================================
550* =====================================================
551* Begin test ZGELST
552* =====================================================
553 IF( irank.EQ.1 ) THEN
554*
555* Generate a matrix of scaling type ISCALE
556*
557 CALL zqrt13( iscale, m, n, copya, lda, norma,
558 \$ iseed )
559*
560* Loop for testing different block sizes.
561*
562 DO inb = 1, nnb
563 nb = nbval( inb )
564 CALL xlaenv( 1, nb )
565 CALL xlaenv( 3, nxval( inb ) )
566*
567* Loop for testing non-transposed and transposed.
568*
569 DO itran = 1, 2
570 IF( itran.EQ.1 ) THEN
571 trans = 'N'
572 nrows = m
573 ncols = n
574 ELSE
575 trans = 'C'
576 nrows = n
577 ncols = m
578 END IF
579 ldwork = max( 1, ncols )
580*
581* Set up a consistent rhs
582*
583 IF( ncols.GT.0 ) THEN
584 CALL zlarnv( 2, iseed, ncols*nrhs,
585 \$ work )
586 CALL zdscal( ncols*nrhs,
587 \$ one / dble( ncols ), work,
588 \$ 1 )
589 END IF
590 CALL zgemm( trans, 'No transpose', nrows,
591 \$ nrhs, ncols, cone, copya, lda,
592 \$ work, ldwork, czero, b, ldb )
593 CALL zlacpy( 'Full', nrows, nrhs, b, ldb,
594 \$ copyb, ldb )
595*
596* Solve LS or overdetermined system
597*
598 IF( m.GT.0 .AND. n.GT.0 ) THEN
599 CALL zlacpy( 'Full', m, n, copya, lda,
600 \$ a, lda )
601 CALL zlacpy( 'Full', nrows, nrhs,
602 \$ copyb, ldb, b, ldb )
603 END IF
604 srnamt = 'ZGELST'
605 CALL zgelst( trans, m, n, nrhs, a, lda, b,
606 \$ ldb, work, lwork, info )
607*
608 IF( info.NE.0 )
609 \$ CALL alaerh( path, 'ZGELST', info, 0,
610 \$ trans, m, n, nrhs, -1, nb,
611 \$ itype, nfail, nerrs,
612 \$ nout )
613*
614* Test 3: Check correctness of results
615* for ZGELST, compute the residual:
616* RESID = norm(B - A*X) /
617* / ( max(m,n) * norm(A) * norm(X) * EPS )
618*
619 IF( nrows.GT.0 .AND. nrhs.GT.0 )
620 \$ CALL zlacpy( 'Full', nrows, nrhs,
621 \$ copyb, ldb, c, ldb )
622 CALL zqrt16( trans, m, n, nrhs, copya,
623 \$ lda, b, ldb, c, ldb, rwork,
624 \$ result( 3 ) )
625*
626* Test 4: Check correctness of results
627* for ZGELST.
628*
629 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
630 \$ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
631*
632* Solving LS system
633*
634 result( 4 ) = zqrt17( trans, 1, m, n,
635 \$ nrhs, copya, lda, b, ldb,
636 \$ copyb, ldb, c, work,
637 \$ lwork )
638 ELSE
639*
640* Solving overdetermined system
641*
642 result( 4 ) = zqrt14( trans, m, n,
643 \$ nrhs, copya, lda, b, ldb,
644 \$ work, lwork )
645 END IF
646*
647* Print information about the tests that
648* did not pass the threshold.
649*
650 DO k = 3, 4
651 IF( result( k ).GE.thresh ) THEN
652 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
653 \$ CALL alahd( nout, path )
654 WRITE( nout, fmt = 9999 )trans, m,
655 \$ n, nrhs, nb, itype, k,
656 \$ result( k )
657 nfail = nfail + 1
658 END IF
659 END DO
660 nrun = nrun + 2
661 END DO
662 END DO
663 END IF
664* =====================================================
665* End test ZGELST
666* =====================================================
667* =====================================================
668* Begin test ZGELSTSLS
669* =====================================================
670 IF( irank.EQ.1 ) THEN
671*
672* Generate a matrix of scaling type ISCALE
673*
674 CALL zqrt13( iscale, m, n, copya, lda, norma,
675 \$ iseed )
676*
677* Loop for testing different block sizes MB.
678*
679 DO inb = 1, nnb
680 mb = nbval( inb )
681 CALL xlaenv( 1, mb )
682*
683* Loop for testing different block sizes NB.
684*
685 DO imb = 1, nnb
686 nb = nbval( imb )
687 CALL xlaenv( 2, nb )
688*
689* Loop for testing non-transposed
690* and transposed.
691*
692 DO itran = 1, 2
693 IF( itran.EQ.1 ) THEN
694 trans = 'N'
695 nrows = m
696 ncols = n
697 ELSE
698 trans = 'C'
699 nrows = n
700 ncols = m
701 END IF
702 ldwork = max( 1, ncols )
703*
704* Set up a consistent rhs
705*
706 IF( ncols.GT.0 ) THEN
707 CALL zlarnv( 2, iseed, ncols*nrhs,
708 \$ work )
709 CALL zscal( ncols*nrhs,
710 \$ cone / dble( ncols ),
711 \$ work, 1 )
712 END IF
713 CALL zgemm( trans, 'No transpose',
714 \$ nrows, nrhs, ncols, cone,
715 \$ copya, lda, work, ldwork,
716 \$ czero, b, ldb )
717 CALL zlacpy( 'Full', nrows, nrhs,
718 \$ b, ldb, copyb, ldb )
719*
720* Solve LS or overdetermined system
721*
722 IF( m.GT.0 .AND. n.GT.0 ) THEN
723 CALL zlacpy( 'Full', m, n,
724 \$ copya, lda, a, lda )
725 CALL zlacpy( 'Full', nrows, nrhs,
726 \$ copyb, ldb, b, ldb )
727 END IF
728 srnamt = 'ZGETSLS '
729 CALL zgetsls( trans, m, n, nrhs, a,
730 \$ lda, b, ldb, work, lwork,
731 \$ info )
732 IF( info.NE.0 )
733 \$ CALL alaerh( path, 'ZGETSLS ', info,
734 \$ 0, trans, m, n, nrhs,
735 \$ -1, nb, itype, nfail,
736 \$ nerrs, nout )
737*
738* Test 5: Check correctness of results
739* for ZGETSLS, compute the residual:
740* RESID = norm(B - A*X) /
741* / ( max(m,n) * norm(A) * norm(X) * EPS )
742*
743 IF( nrows.GT.0 .AND. nrhs.GT.0 )
744 \$ CALL zlacpy( 'Full', nrows, nrhs,
745 \$ copyb, ldb, c, ldb )
746 CALL zqrt16( trans, m, n, nrhs,
747 \$ copya, lda, b, ldb,
748 \$ c, ldb, work2,
749 \$ result( 5 ) )
750*
751* Test 6: Check correctness of results
752* for ZGETSLS.
753*
754 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
755 \$ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
756*
757* Solving LS system, compute:
758* r = norm((B- A*X)**T * A) /
759* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
760*
761 result( 6 ) = zqrt17( trans, 1, m,
762 \$ n, nrhs, copya, lda,
763 \$ b, ldb, copyb, ldb,
764 \$ c, work, lwork )
765 ELSE
766*
767* Solving overdetermined system
768*
769 result( 6 ) = zqrt14( trans, m, n,
770 \$ nrhs, copya, lda, b,
771 \$ ldb, work, lwork )
772 END IF
773*
774* Print information about the tests that
775* did not pass the threshold.
776*
777 DO k = 5, 6
778 IF( result( k ).GE.thresh ) THEN
779 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
780 \$ CALL alahd( nout, path )
781 WRITE( nout, fmt = 9997 )trans,
782 \$ m, n, nrhs, mb, nb, itype, k,
783 \$ result( k )
784 nfail = nfail + 1
785 END IF
786 END DO
787 nrun = nrun + 2
788 END DO
789 END DO
790 END DO
791 END IF
792* =====================================================
793* End test ZGELSTSLS
794* =====================================================
795*
796* Generate a matrix of scaling type ISCALE and rank
797* type IRANK.
798*
799 CALL zqrt15( iscale, irank, m, n, nrhs, copya, lda,
800 \$ copyb, ldb, copys, rank, norma, normb,
801 \$ iseed, work, lwork )
802*
803* workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
804*
805 ldwork = max( 1, m )
806*
807* Loop for testing different block sizes.
808*
809 DO 90 inb = 1, nnb
810 nb = nbval( inb )
811 CALL xlaenv( 1, nb )
812 CALL xlaenv( 3, nxval( inb ) )
813*
814* Test ZGELSY
815*
816* ZGELSY: Compute the minimum-norm solution
817* X to min( norm( A * X - B ) )
818* using the rank-revealing orthogonal
819* factorization.
820*
821 CALL zlacpy( 'Full', m, n, copya, lda, a, lda )
822 CALL zlacpy( 'Full', m, nrhs, copyb, ldb, b,
823 \$ ldb )
824*
825* Initialize vector IWORK.
826*
827 DO 70 j = 1, n
828 iwork( j ) = 0
829 70 CONTINUE
830*
831 srnamt = 'ZGELSY'
832 CALL zgelsy( m, n, nrhs, a, lda, b, ldb, iwork,
833 \$ rcond, crank, work, lwlsy, rwork,
834 \$ info )
835 IF( info.NE.0 )
836 \$ CALL alaerh( path, 'ZGELSY', info, 0, ' ', m,
837 \$ n, nrhs, -1, nb, itype, nfail,
838 \$ nerrs, nout )
839*
840* workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS)
841*
842* Test 7: Compute relative error in svd
843* workspace: M*N + 4*MIN(M,N) + MAX(M,N)
844*
845 result( 7 ) = zqrt12( crank, crank, a, lda,
846 \$ copys, work, lwork, rwork )
847*
848* Test 8: Compute error in solution
849* workspace: M*NRHS + M
850*
851 CALL zlacpy( 'Full', m, nrhs, copyb, ldb, work,
852 \$ ldwork )
853 CALL zqrt16( 'No transpose', m, n, nrhs, copya,
854 \$ lda, b, ldb, work, ldwork, rwork,
855 \$ result( 8 ) )
856*
857* Test 9: Check norm of r'*A
858* workspace: NRHS*(M+N)
859*
860 result( 9 ) = zero
861 IF( m.GT.crank )
862 \$ result( 9 ) = zqrt17( 'No transpose', 1, m,
863 \$ n, nrhs, copya, lda, b, ldb,
864 \$ copyb, ldb, c, work, lwork )
865*
866* Test 10: Check if x is in the rowspace of A
867* workspace: (M+NRHS)*(N+2)
868*
869 result( 10 ) = zero
870*
871 IF( n.GT.crank )
872 \$ result( 10 ) = zqrt14( 'No transpose', m, n,
873 \$ nrhs, copya, lda, b, ldb,
874 \$ work, lwork )
875*
876* Test ZGELSS
877*
878* ZGELSS: Compute the minimum-norm solution
879* X to min( norm( A * X - B ) )
880* using the SVD.
881*
882 CALL zlacpy( 'Full', m, n, copya, lda, a, lda )
883 CALL zlacpy( 'Full', m, nrhs, copyb, ldb, b,
884 \$ ldb )
885 srnamt = 'ZGELSS'
886 CALL zgelss( m, n, nrhs, a, lda, b, ldb, s,
887 \$ rcond, crank, work, lwork, rwork,
888 \$ info )
889*
890 IF( info.NE.0 )
891 \$ CALL alaerh( path, 'ZGELSS', info, 0, ' ', m,
892 \$ n, nrhs, -1, nb, itype, nfail,
893 \$ nerrs, nout )
894*
895* workspace used: 3*min(m,n) +
896* max(2*min(m,n),nrhs,max(m,n))
897*
898* Test 11: Compute relative error in svd
899*
900 IF( rank.GT.0 ) THEN
901 CALL daxpy( mnmin, -one, copys, 1, s, 1 )
902 result( 11 ) = dasum( mnmin, s, 1 ) /
903 \$ dasum( mnmin, copys, 1 ) /
904 \$ ( eps*dble( mnmin ) )
905 ELSE
906 result( 11 ) = zero
907 END IF
908*
909* Test 12: Compute error in solution
910*
911 CALL zlacpy( 'Full', m, nrhs, copyb, ldb, work,
912 \$ ldwork )
913 CALL zqrt16( 'No transpose', m, n, nrhs, copya,
914 \$ lda, b, ldb, work, ldwork, rwork,
915 \$ result( 12 ) )
916*
917* Test 13: Check norm of r'*A
918*
919 result( 13 ) = zero
920 IF( m.GT.crank )
921 \$ result( 13 ) = zqrt17( 'No transpose', 1, m,
922 \$ n, nrhs, copya, lda, b, ldb,
923 \$ copyb, ldb, c, work, lwork )
924*
925* Test 14: Check if x is in the rowspace of A
926*
927 result( 14 ) = zero
928 IF( n.GT.crank )
929 \$ result( 14 ) = zqrt14( 'No transpose', m, n,
930 \$ nrhs, copya, lda, b, ldb,
931 \$ work, lwork )
932*
933* Test ZGELSD
934*
935* ZGELSD: Compute the minimum-norm solution X
936* to min( norm( A * X - B ) ) using a
937* divide and conquer SVD.
938*
939 CALL xlaenv( 9, 25 )
940*
941 CALL zlacpy( 'Full', m, n, copya, lda, a, lda )
942 CALL zlacpy( 'Full', m, nrhs, copyb, ldb, b,
943 \$ ldb )
944*
945 srnamt = 'ZGELSD'
946 CALL zgelsd( m, n, nrhs, a, lda, b, ldb, s,
947 \$ rcond, crank, work, lwork, rwork,
948 \$ iwork, info )
949 IF( info.NE.0 )
950 \$ CALL alaerh( path, 'ZGELSD', info, 0, ' ', m,
951 \$ n, nrhs, -1, nb, itype, nfail,
952 \$ nerrs, nout )
953*
954* Test 15: Compute relative error in svd
955*
956 IF( rank.GT.0 ) THEN
957 CALL daxpy( mnmin, -one, copys, 1, s, 1 )
958 result( 15 ) = dasum( mnmin, s, 1 ) /
959 \$ dasum( mnmin, copys, 1 ) /
960 \$ ( eps*dble( mnmin ) )
961 ELSE
962 result( 15 ) = zero
963 END IF
964*
965* Test 16: Compute error in solution
966*
967 CALL zlacpy( 'Full', m, nrhs, copyb, ldb, work,
968 \$ ldwork )
969 CALL zqrt16( 'No transpose', m, n, nrhs, copya,
970 \$ lda, b, ldb, work, ldwork, rwork,
971 \$ result( 16 ) )
972*
973* Test 17: Check norm of r'*A
974*
975 result( 17 ) = zero
976 IF( m.GT.crank )
977 \$ result( 17 ) = zqrt17( 'No transpose', 1, m,
978 \$ n, nrhs, copya, lda, b, ldb,
979 \$ copyb, ldb, c, work, lwork )
980*
981* Test 18: Check if x is in the rowspace of A
982*
983 result( 18 ) = zero
984 IF( n.GT.crank )
985 \$ result( 18 ) = zqrt14( 'No transpose', m, n,
986 \$ nrhs, copya, lda, b, ldb,
987 \$ work, lwork )
988*
989* Print information about the tests that did not
990* pass the threshold.
991*
992 DO 80 k = 7, 18
993 IF( result( k ).GE.thresh ) THEN
994 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
995 \$ CALL alahd( nout, path )
996 WRITE( nout, fmt = 9998 )m, n, nrhs, nb,
997 \$ itype, k, result( k )
998 nfail = nfail + 1
999 END IF
1000 80 CONTINUE
1001 nrun = nrun + 12
1002*
1003 90 CONTINUE
1004 100 CONTINUE
1005 110 CONTINUE
1006 120 CONTINUE
1007 130 CONTINUE
1008 140 CONTINUE
1009*
1010* Print a summary of the results.
1011*
1012 CALL alasvm( path, nout, nfail, nrun, nerrs )
1013*
1014 9999 FORMAT( ' TRANS=''', a1, ''', M=', i5, ', N=', i5, ', NRHS=', i4,
1015 \$ ', NB=', i4, ', type', i2, ', test(', i2, ')=', g12.5 )
1016 9998 FORMAT( ' M=', i5, ', N=', i5, ', NRHS=', i4, ', NB=', i4,
1017 \$ ', type', i2, ', test(', i2, ')=', g12.5 )
1018 9997 FORMAT( ' TRANS=''', a1,' M=', i5, ', N=', i5, ', NRHS=', i4,
1019 \$ ', MB=', i4,', NB=', i4,', type', i2,
1020 \$ ', test(', i2, ')=', g12.5 )
1021*
1022 DEALLOCATE( work )
1023 DEALLOCATE( iwork )
1024 DEALLOCATE( rwork )
1025 RETURN
1026*
1027* End of ZDRVLS
1028*
1029 END
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:73
subroutine xlaenv(ISPEC, NVALUE)
XLAENV
Definition: xlaenv.f:81
subroutine alahd(IOUNIT, PATH)
ALAHD
Definition: alahd.f:107
subroutine alaerh(PATH, SUBNAM, INFO, INFOE, OPTS, M, N, KL, KU, N5, IMAT, NFAIL, NERRS, NOUT)
ALAERH
Definition: alaerh.f:147
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zqrt16(TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, RWORK, RESID)
ZQRT16
Definition: zqrt16.f:133
subroutine zqrt15(SCALE, RKSEL, M, N, NRHS, A, LDA, B, LDB, S, RANK, NORMA, NORMB, ISEED, WORK, LWORK)
ZQRT15
Definition: zqrt15.f:149
subroutine zerrls(PATH, NUNIT)
ZERRLS
Definition: zerrls.f:55
subroutine zqrt13(SCALE, M, N, A, LDA, NORMA, ISEED)
ZQRT13
Definition: zqrt13.f:91
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:192
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:182
subroutine zgetsls(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
ZGETSLS
Definition: zgetsls.f:162
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:178
subroutine zgelst(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
ZGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization ...
Definition: zgelst.f:194
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:210
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:225
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: zlarnv.f:99
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89