LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
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,
252 $ zgelsd, zgelss, zgelst, zgelsy, zgemm,
253 $ zgetsls, zlacpy, zlarnv, zqrt13, zqrt15,
254 $ zqrt16, zdscal, xlaenv
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
alasvm
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
xlaenv
subroutine xlaenv(ispec, nvalue)
XLAENV
Definition xlaenv.f:81
alaerh
subroutine alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
Definition alaerh.f:147
alahd
subroutine alahd(iounit, path)
ALAHD
Definition alahd.f:107
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
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:182
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:219
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:178
zgelst
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
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
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
zgetsls
subroutine zgetsls(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
ZGETSLS
Definition zgetsls.f:162
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:103
zlarnv
subroutine zlarnv(idist, iseed, n, x)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition zlarnv.f:99
zdscal
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
zscal
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
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:192
zerrls
subroutine zerrls(path, nunit)
ZERRLS
Definition zerrls.f:55
zqrt13
subroutine zqrt13(scale, m, n, a, lda, norma, iseed)
ZQRT13
Definition zqrt13.f:91
zqrt15
subroutine zqrt15(scale, rksel, m, n, nrhs, a, lda, b, ldb, s, rank, norma, normb, iseed, work, lwork)
ZQRT15
Definition zqrt15.f:149
zqrt16
subroutine zqrt16(trans, m, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
ZQRT16
Definition zqrt16.f:133