OpenRadioss 2025.1.11
OpenRadioss project
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, 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*> \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 = 16 )
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_zgetsls, lwork_zgelss,
232 $ 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, dlasrt, xlaenv,
252 $ zdscal, zerrls, zgels, zgelsd, zgelss,
253 $ zgelsy, zgemm, zlacpy, zlarnv, zqrt13, zqrt15,
254 $ zqrt16, zgetsls
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, ?GETSLS, ?GELSY, ?GELSS and ?GELSD routines.
338*
339 DO im = 1, nm
340 m = mval( im )
341 lda = max( 1, m )
342 DO in = 1, nn
343 n = nval( in )
344 mnmin = max(min( m, n ),1)
345 ldb = max( 1, m, n )
346 DO ins = 1, nns
347 nrhs = nsval( ins )
348 DO irank = 1, 2
349 DO iscale = 1, 3
350 itype = ( irank-1 )*3 + iscale
351 IF( dotype( itype ) ) THEN
352 IF( irank.EQ.1 ) THEN
353 DO itran = 1, 2
354 IF( itran.EQ.1 ) THEN
355 trans = 'N'
356 ELSE
357 trans = 'C'
358 END IF
359*
360* Compute workspace needed for ZGELS
361 CALL zgels( trans, m, n, nrhs, a, lda,
362 $ b, ldb, wq, -1, info )
363 lwork_zgels = int( wq( 1 ) )
364* Compute workspace needed for ZGETSLS
365 CALL zgetsls( trans, m, n, nrhs, a, lda,
366 $ b, ldb, wq, -1, info )
367 lwork_zgetsls = int( wq( 1 ) )
368 ENDDO
369 END IF
370* Compute workspace needed for ZGELSY
371 CALL zgelsy( m, n, nrhs, a, lda, b, ldb, iwq,
372 $ rcond, crank, wq, -1, rwq, info )
373 lwork_zgelsy = int( wq( 1 ) )
374 lrwork_zgelsy = 2*n
375* Compute workspace needed for ZGELSS
376 CALL zgelss( m, n, nrhs, a, lda, b, ldb, s,
377 $ rcond, crank, wq, -1 , rwq,
378 $ info )
379 lwork_zgelss = int( wq( 1 ) )
380 lrwork_zgelss = 5*mnmin
381* Compute workspace needed for ZGELSD
382 CALL zgelsd( m, n, nrhs, a, lda, b, ldb, s,
383 $ rcond, crank, wq, -1, rwq, iwq,
384 $ info )
385 lwork_zgelsd = int( wq( 1 ) )
386 lrwork_zgelsd = int( rwq( 1 ) )
387* Compute LIWORK workspace needed for ZGELSY and ZGELSD
388 liwork = max( liwork, n, iwq( 1 ) )
389* Compute LRWORK workspace needed for ZGELSY, ZGELSS and ZGELSD
390 lrwork = max( lrwork, lrwork_zgelsy,
391 $ lrwork_zgelss, lrwork_zgelsd )
392* Compute LWORK workspace needed for all functions
393 lwork = max( lwork, lwork_zgels, lwork_zgetsls,
394 $ lwork_zgelsy, lwork_zgelss,
395 $ lwork_zgelsd )
396 END IF
397 ENDDO
398 ENDDO
399 ENDDO
400 ENDDO
401 ENDDO
402*
403 lwlsy = lwork
404*
405 ALLOCATE( work( lwork ) )
406 ALLOCATE( work2( 2 * lwork ) )
407 ALLOCATE( iwork( liwork ) )
408 ALLOCATE( rwork( lrwork ) )
409*
410 DO 140 im = 1, nm
411 m = mval( im )
412 lda = max( 1, m )
413*
414 DO 130 in = 1, nn
415 n = nval( in )
416 mnmin = max(min( m, n ),1)
417 ldb = max( 1, m, n )
418 mb = (mnmin+1)
419*
420 DO 120 ins = 1, nns
421 nrhs = nsval( ins )
422*
423 DO 110 irank = 1, 2
424 DO 100 iscale = 1, 3
425 itype = ( irank-1 )*3 + iscale
426 IF( .NOT.dotype( itype ) )
427 $ GO TO 100
428*
429 IF( irank.EQ.1 ) THEN
430*
431* Test ZGELS
432*
433* Generate a matrix of scaling type ISCALE
434*
435 CALL zqrt13( iscale, m, n, copya, lda, norma,
436 $ iseed )
437 DO 40 inb = 1, nnb
438 nb = nbval( inb )
439 CALL xlaenv( 1, nb )
440 CALL xlaenv( 3, nxval( inb ) )
441*
442 DO 30 itran = 1, 2
443 IF( itran.EQ.1 ) THEN
444 trans = 'N'
445 nrows = m
446 ncols = n
447 ELSE
448 trans = 'C'
449 nrows = n
450 ncols = m
451 END IF
452 ldwork = max( 1, ncols )
453*
454* Set up a consistent rhs
455*
456 IF( ncols.GT.0 ) THEN
457 CALL zlarnv( 2, iseed, ncols*nrhs,
458 $ work )
459 CALL zdscal( ncols*nrhs,
460 $ one / dble( ncols ), work,
461 $ 1 )
462 END IF
463 CALL zgemm( trans, 'No transpose', nrows,
464 $ nrhs, ncols, cone, copya, lda,
465 $ work, ldwork, czero, b, ldb )
466 CALL zlacpy( 'Full', nrows, nrhs, b, ldb,
467 $ copyb, ldb )
468*
469* Solve LS or overdetermined system
470*
471 IF( m.GT.0 .AND. n.GT.0 ) THEN
472 CALL zlacpy( 'Full', m, n, copya, lda,
473 $ a, lda )
474 CALL zlacpy( 'Full', nrows, nrhs,
475 $ copyb, ldb, b, ldb )
476 END IF
477 srnamt = 'ZGELS '
478 CALL zgels( trans, m, n, nrhs, a, lda, b,
479 $ ldb, work, lwork, info )
480*
481 IF( info.NE.0 )
482 $ CALL alaerh( path, 'ZGELS ', info, 0,
483 $ trans, m, n, nrhs, -1, nb,
484 $ itype, nfail, nerrs,
485 $ nout )
486*
487* Check correctness of results
488*
489 ldwork = max( 1, nrows )
490 IF( nrows.GT.0 .AND. nrhs.GT.0 )
491 $ CALL zlacpy( 'Full', nrows, nrhs,
492 $ copyb, ldb, c, ldb )
493 CALL zqrt16( trans, m, n, nrhs, copya,
494 $ lda, b, ldb, c, ldb, rwork,
495 $ result( 1 ) )
496*
497 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
498 $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
499*
500* Solving LS system
501*
502 result( 2 ) = zqrt17( trans, 1, m, n,
503 $ nrhs, copya, lda, b, ldb,
504 $ copyb, ldb, c, work,
505 $ lwork )
506 ELSE
507*
508* Solving overdetermined system
509*
510 result( 2 ) = zqrt14( trans, m, n,
511 $ nrhs, copya, lda, b, ldb,
512 $ work, lwork )
513 END IF
514*
515* Print information about the tests that
516* did not pass the threshold.
517*
518 DO 20 k = 1, 2
519 IF( result( k ).GE.thresh ) THEN
520 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
521 $ CALL alahd( nout, path )
522 WRITE( nout, fmt = 9999 )trans, m,
523 $ n, nrhs, nb, itype, k,
524 $ result( k )
525 nfail = nfail + 1
526 END IF
527 20 CONTINUE
528 nrun = nrun + 2
529 30 CONTINUE
530 40 CONTINUE
531*
532*
533* Test ZGETSLS
534*
535* Generate a matrix of scaling type ISCALE
536*
537 CALL zqrt13( iscale, m, n, copya, lda, norma,
538 $ iseed )
539 DO 65 inb = 1, nnb
540 mb = nbval( inb )
541 CALL xlaenv( 1, mb )
542 DO 62 imb = 1, nnb
543 nb = nbval( imb )
544 CALL xlaenv( 2, nb )
545*
546 DO 60 itran = 1, 2
547 IF( itran.EQ.1 ) THEN
548 trans = 'N'
549 nrows = m
550 ncols = n
551 ELSE
552 trans = 'C'
553 nrows = n
554 ncols = m
555 END IF
556 ldwork = max( 1, ncols )
557*
558* Set up a consistent rhs
559*
560 IF( ncols.GT.0 ) THEN
561 CALL zlarnv( 2, iseed, ncols*nrhs,
562 $ work )
563 CALL zscal( ncols*nrhs,
564 $ cone / dble( ncols ), work,
565 $ 1 )
566 END IF
567 CALL zgemm( trans, 'No transpose', nrows,
568 $ nrhs, ncols, cone, copya, lda,
569 $ work, ldwork, czero, b, ldb )
570 CALL zlacpy( 'Full', nrows, nrhs, b, ldb,
571 $ copyb, ldb )
572*
573* Solve LS or overdetermined system
574*
575 IF( m.GT.0 .AND. n.GT.0 ) THEN
576 CALL zlacpy( 'Full', m, n, copya, lda,
577 $ a, lda )
578 CALL zlacpy( 'Full', nrows, nrhs,
579 $ copyb, ldb, b, ldb )
580 END IF
581 srnamt = 'ZGETSLS '
582 CALL zgetsls( trans, m, n, nrhs, a,
583 $ lda, b, ldb, work, lwork, info )
584 IF( info.NE.0 )
585 $ CALL alaerh( path, 'ZGETSLS ', info, 0,
586 $ trans, m, n, nrhs, -1, nb,
587 $ itype, nfail, nerrs,
588 $ nout )
589*
590* Check correctness of results
591*
592 ldwork = max( 1, nrows )
593 IF( nrows.GT.0 .AND. nrhs.GT.0 )
594 $ CALL zlacpy( 'Full', nrows, nrhs,
595 $ copyb, ldb, c, ldb )
596 CALL zqrt16( trans, m, n, nrhs, copya,
597 $ lda, b, ldb, c, ldb, work2,
598 $ result( 15 ) )
599*
600 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
601 $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
602*
603* Solving LS system
604*
605 result( 16 ) = zqrt17( trans, 1, m, n,
606 $ nrhs, copya, lda, b, ldb,
607 $ copyb, ldb, c, work,
608 $ lwork )
609 ELSE
610*
611* Solving overdetermined system
612*
613 result( 16 ) = zqrt14( trans, m, n,
614 $ nrhs, copya, lda, b, ldb,
615 $ work, lwork )
616 END IF
617*
618* Print information about the tests that
619* did not pass the threshold.
620*
621 DO 50 k = 15, 16
622 IF( result( k ).GE.thresh ) THEN
623 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
624 $ CALL alahd( nout, path )
625 WRITE( nout, fmt = 9997 )trans, m,
626 $ n, nrhs, mb, nb, itype, k,
627 $ result( k )
628 nfail = nfail + 1
629 END IF
630 50 CONTINUE
631 nrun = nrun + 2
632 60 CONTINUE
633 62 CONTINUE
634 65 CONTINUE
635 END IF
636*
637* Generate a matrix of scaling type ISCALE and rank
638* type IRANK.
639*
640 CALL zqrt15( iscale, irank, m, n, nrhs, copya, lda,
641 $ copyb, ldb, copys, rank, norma, normb,
642 $ iseed, work, lwork )
643*
644* workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
645*
646 ldwork = max( 1, m )
647*
648* Loop for testing different block sizes.
649*
650 DO 90 inb = 1, nnb
651 nb = nbval( inb )
652 CALL xlaenv( 1, nb )
653 CALL xlaenv( 3, nxval( inb ) )
654*
655* Test ZGELSY
656*
657* ZGELSY: Compute the minimum-norm solution
658* X to min( norm( A * X - B ) )
659* using the rank-revealing orthogonal
660* factorization.
661*
662 CALL zlacpy( 'Full', m, n, copya, lda, a, lda )
663 CALL zlacpy( 'Full', m, nrhs, copyb, ldb, b,
664 $ ldb )
665*
666* Initialize vector IWORK.
667*
668 DO 70 j = 1, n
669 iwork( j ) = 0
670 70 CONTINUE
671*
672 srnamt = 'ZGELSY'
673 CALL zgelsy( m, n, nrhs, a, lda, b, ldb, iwork,
674 $ rcond, crank, work, lwlsy, rwork,
675 $ info )
676 IF( info.NE.0 )
677 $ CALL alaerh( path, 'ZGELSY', info, 0, ' ', m,
678 $ n, nrhs, -1, nb, itype, nfail,
679 $ nerrs, nout )
680*
681* workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS)
682*
683* Test 3: Compute relative error in svd
684* workspace: M*N + 4*MIN(M,N) + MAX(M,N)
685*
686 result( 3 ) = zqrt12( crank, crank, a, lda,
687 $ copys, work, lwork, rwork )
688*
689* Test 4: Compute error in solution
690* workspace: M*NRHS + M
691*
692 CALL zlacpy( 'Full', m, nrhs, copyb, ldb, work,
693 $ ldwork )
694 CALL zqrt16( 'No transpose', m, n, nrhs, copya,
695 $ lda, b, ldb, work, ldwork, rwork,
696 $ result( 4 ) )
697*
698* Test 5: Check norm of r'*A
699* workspace: NRHS*(M+N)
700*
701 result( 5 ) = zero
702 IF( m.GT.crank )
703 $ result( 5 ) = zqrt17( 'No transpose', 1, m,
704 $ n, nrhs, copya, lda, b, ldb,
705 $ copyb, ldb, c, work, lwork )
706*
707* Test 6: Check if x is in the rowspace of A
708* workspace: (M+NRHS)*(N+2)
709*
710 result( 6 ) = zero
711*
712 IF( n.GT.crank )
713 $ result( 6 ) = zqrt14( 'No transpose', m, n,
714 $ nrhs, copya, lda, b, ldb,
715 $ work, lwork )
716*
717* Test ZGELSS
718*
719* ZGELSS: Compute the minimum-norm solution
720* X to min( norm( A * X - B ) )
721* using the SVD.
722*
723 CALL zlacpy( 'Full', m, n, copya, lda, a, lda )
724 CALL zlacpy( 'Full', m, nrhs, copyb, ldb, b,
725 $ ldb )
726 srnamt = 'ZGELSS'
727 CALL zgelss( m, n, nrhs, a, lda, b, ldb, s,
728 $ rcond, crank, work, lwork, rwork,
729 $ info )
730*
731 IF( info.NE.0 )
732 $ CALL alaerh( path, 'ZGELSS', info, 0, ' ', m,
733 $ n, nrhs, -1, nb, itype, nfail,
734 $ nerrs, nout )
735*
736* workspace used: 3*min(m,n) +
737* max(2*min(m,n),nrhs,max(m,n))
738*
739* Test 7: Compute relative error in svd
740*
741 IF( rank.GT.0 ) THEN
742 CALL daxpy( mnmin, -one, copys, 1, s, 1 )
743 result( 7 ) = dasum( mnmin, s, 1 ) /
744 $ dasum( mnmin, copys, 1 ) /
745 $ ( eps*dble( mnmin ) )
746 ELSE
747 result( 7 ) = zero
748 END IF
749*
750* Test 8: Compute error in solution
751*
752 CALL zlacpy( 'Full', m, nrhs, copyb, ldb, work,
753 $ ldwork )
754 CALL zqrt16( 'No transpose', m, n, nrhs, copya,
755 $ lda, b, ldb, work, ldwork, rwork,
756 $ result( 8 ) )
757*
758* Test 9: Check norm of r'*A
759*
760 result( 9 ) = zero
761 IF( m.GT.crank )
762 $ result( 9 ) = zqrt17( 'No transpose', 1, m,
763 $ n, nrhs, copya, lda, b, ldb,
764 $ copyb, ldb, c, work, lwork )
765*
766* Test 10: Check if x is in the rowspace of A
767*
768 result( 10 ) = zero
769 IF( n.GT.crank )
770 $ result( 10 ) = zqrt14( 'No transpose', m, n,
771 $ nrhs, copya, lda, b, ldb,
772 $ work, lwork )
773*
774* Test ZGELSD
775*
776* ZGELSD: Compute the minimum-norm solution X
777* to min( norm( A * X - B ) ) using a
778* divide and conquer SVD.
779*
780 CALL xlaenv( 9, 25 )
781*
782 CALL zlacpy( 'Full', m, n, copya, lda, a, lda )
783 CALL zlacpy( 'Full', m, nrhs, copyb, ldb, b,
784 $ ldb )
785*
786 srnamt = 'ZGELSD'
787 CALL zgelsd( m, n, nrhs, a, lda, b, ldb, s,
788 $ rcond, crank, work, lwork, rwork,
789 $ iwork, info )
790 IF( info.NE.0 )
791 $ CALL alaerh( path, 'ZGELSD', info, 0, ' ', m,
792 $ n, nrhs, -1, nb, itype, nfail,
793 $ nerrs, nout )
794*
795* Test 11: Compute relative error in svd
796*
797 IF( rank.GT.0 ) THEN
798 CALL daxpy( mnmin, -one, copys, 1, s, 1 )
799 result( 11 ) = dasum( mnmin, s, 1 ) /
800 $ dasum( mnmin, copys, 1 ) /
801 $ ( eps*dble( mnmin ) )
802 ELSE
803 result( 11 ) = zero
804 END IF
805*
806* Test 12: Compute error in solution
807*
808 CALL zlacpy( 'Full', m, nrhs, copyb, ldb, work,
809 $ ldwork )
810 CALL zqrt16( 'No transpose', m, n, nrhs, copya,
811 $ lda, b, ldb, work, ldwork, rwork,
812 $ result( 12 ) )
813*
814* Test 13: Check norm of r'*A
815*
816 result( 13 ) = zero
817 IF( m.GT.crank )
818 $ result( 13 ) = zqrt17( 'No transpose', 1, m,
819 $ n, nrhs, copya, lda, b, ldb,
820 $ copyb, ldb, c, work, lwork )
821*
822* Test 14: Check if x is in the rowspace of A
823*
824 result( 14 ) = zero
825 IF( n.GT.crank )
826 $ result( 14 ) = zqrt14( 'No transpose', m, n,
827 $ nrhs, copya, lda, b, ldb,
828 $ work, lwork )
829*
830* Print information about the tests that did not
831* pass the threshold.
832*
833 DO 80 k = 3, 14
834 IF( result( k ).GE.thresh ) THEN
835 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
836 $ CALL alahd( nout, path )
837 WRITE( nout, fmt = 9998 )m, n, nrhs, nb,
838 $ itype, k, result( k )
839 nfail = nfail + 1
840 END IF
841 80 CONTINUE
842 nrun = nrun + 12
843*
844 90 CONTINUE
845 100 CONTINUE
846 110 CONTINUE
847 120 CONTINUE
848 130 CONTINUE
849 140 CONTINUE
850*
851* Print a summary of the results.
852*
853 CALL alasvm( path, nout, nfail, nrun, nerrs )
854*
855 9999 FORMAT( ' TRANS=''', a1, ''', M=', i5, ', N=', i5, ', NRHS=', i4,
856 $ ', NB=', i4, ', type', i2, ', test(', i2, ')=', g12.5 )
857 9998 FORMAT( ' M=', i5, ', N=', i5, ', NRHS=', i4, ', NB=', i4,
858 $ ', type', i2, ', test(', i2, ')=', g12.5 )
859 9997 FORMAT( ' TRANS=''', a1,' M=', i5, ', N=', i5, ', NRHS=', i4,
860 $ ', MB=', i4,', NB=', i4,', type', i2,
861 $ ', test(', i2, ')=', g12.5 )
862*
863 DEALLOCATE( work )
864 DEALLOCATE( iwork )
865 DEALLOCATE( rwork )
866 RETURN
867*
868* End of ZDRVLS
869*
870 END
dlasrt
subroutine dlasrt(id, n, d, info)
DLASRT sorts numbers in increasing or decreasing order.
Definition dlasrt.f:88
xlaenv
subroutine xlaenv(ispec, nvalue)
XLAENV
Definition xlaenv.f:81
alasvm
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
alahd
subroutine alahd(iounit, path)
ALAHD
Definition alahd.f:107
alaerh
subroutine alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
Definition alaerh.f:147
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:210
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
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:225
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
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:187
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
zqrt13
subroutine zqrt13(scale, m, n, a, lda, norma, iseed)
ZQRT13
Definition zqrt13.f:91
zerrls
subroutine zerrls(path, nunit)
ZERRLS
Definition zerrls.f:55
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
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21