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.EQ..OR..EQ..AND..EQ. IF( ( NM0 NN0 ) THRESHZERO )
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.GT. IF ( MVAL( I )MMAX ) THEN
308 MMAX = MVAL( I )
309 END IF
310 ENDDO
311 DO I = 1, NN
312.GT. IF ( NVAL( I )NMAX ) THEN
313 NMAX = NVAL( I )
314 END IF
315 ENDDO
316 DO I = 1, NNS
317.GT. IF ( NSVAL( I )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.EQ. IF( IRANK1 ) THEN
353 DO ITRAN = 1, 2
354.EQ. IF( ITRAN1 ) 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.NOT. IF( DOTYPE( ITYPE ) )
427 $ GO TO 100
428*
429.EQ. IF( IRANK1 ) 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.EQ. IF( ITRAN1 ) 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.GT. IF( NCOLS0 ) 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.GT..AND..GT. IF( M0 N0 ) 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