LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cdrvls.f
Go to the documentation of this file.
1*> \brief \b CDRVLS
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 CDRVLS( 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* REAL THRESH
19* ..
20* .. Array Arguments ..
21* LOGICAL DOTYPE( * )
22* INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
23* $ NVAL( * ), NXVAL( * )
24* REAL COPYS( * ), S( * )
25* COMPLEX A( * ), B( * ), C( * ), COPYA( * ), COPYB( * )
26* ..
27*
28*
29*> \par Purpose:
30* =============
31*>
32*> \verbatim
33*>
34*> CDRVLS tests the least squares driver routines CGELS, CGELST,
35*> CGETSLS, CGELSS, CGELSY
36*> and CGELSD.
37*> \endverbatim
38*
39* Arguments:
40* ==========
41*
42*> \param[in] DOTYPE
43*> \verbatim
44*> DOTYPE is LOGICAL array, dimension (NTYPES)
45*> The matrix types to be used for testing. Matrices of type j
46*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
47*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
48*> The matrix of type j is generated as follows:
49*> j=1: A = U*D*V where U and V are random unitary matrices
50*> and D has random entries (> 0.1) taken from a uniform
51*> distribution (0,1). A is full rank.
52*> j=2: The same of 1, but A is scaled up.
53*> j=3: The same of 1, but A is scaled down.
54*> j=4: A = U*D*V where U and V are random unitary matrices
55*> and D has 3*min(M,N)/4 random entries (> 0.1) taken
56*> from a uniform distribution (0,1) and the remaining
57*> entries set to 0. A is rank-deficient.
58*> j=5: The same of 4, but A is scaled up.
59*> j=6: The same of 5, but A is scaled down.
60*> \endverbatim
61*>
62*> \param[in] NM
63*> \verbatim
64*> NM is INTEGER
65*> The number of values of M contained in the vector MVAL.
66*> \endverbatim
67*>
68*> \param[in] MVAL
69*> \verbatim
70*> MVAL is INTEGER array, dimension (NM)
71*> The values of the matrix row dimension M.
72*> \endverbatim
73*>
74*> \param[in] NN
75*> \verbatim
76*> NN is INTEGER
77*> The number of values of N contained in the vector NVAL.
78*> \endverbatim
79*>
80*> \param[in] NVAL
81*> \verbatim
82*> NVAL is INTEGER array, dimension (NN)
83*> The values of the matrix column dimension N.
84*> \endverbatim
85*>
86*> \param[in] NNB
87*> \verbatim
88*> NNB is INTEGER
89*> The number of values of NB and NX contained in the
90*> vectors NBVAL and NXVAL. The blocking parameters are used
91*> in pairs (NB,NX).
92*> \endverbatim
93*>
94*> \param[in] NBVAL
95*> \verbatim
96*> NBVAL is INTEGER array, dimension (NNB)
97*> The values of the blocksize NB.
98*> \endverbatim
99*>
100*> \param[in] NXVAL
101*> \verbatim
102*> NXVAL is INTEGER array, dimension (NNB)
103*> The values of the crossover point NX.
104*> \endverbatim
105*>
106*> \param[in] NNS
107*> \verbatim
108*> NNS is INTEGER
109*> The number of values of NRHS contained in the vector NSVAL.
110*> \endverbatim
111*>
112*> \param[in] NSVAL
113*> \verbatim
114*> NSVAL is INTEGER array, dimension (NNS)
115*> The values of the number of right hand sides NRHS.
116*> \endverbatim
117*>
118*> \param[in] THRESH
119*> \verbatim
120*> THRESH is REAL
121*> The threshold value for the test ratios. A result is
122*> included in the output file if RESULT >= THRESH. To have
123*> every test ratio printed, use THRESH = 0.
124*> \endverbatim
125*>
126*> \param[in] TSTERR
127*> \verbatim
128*> TSTERR is LOGICAL
129*> Flag that indicates whether error exits are to be tested.
130*> \endverbatim
131*>
132*> \param[out] A
133*> \verbatim
134*> A is COMPLEX array, dimension (MMAX*NMAX)
135*> where MMAX is the maximum value of M in MVAL and NMAX is the
136*> maximum value of N in NVAL.
137*> \endverbatim
138*>
139*> \param[out] COPYA
140*> \verbatim
141*> COPYA is COMPLEX array, dimension (MMAX*NMAX)
142*> \endverbatim
143*>
144*> \param[out] B
145*> \verbatim
146*> B is COMPLEX array, dimension (MMAX*NSMAX)
147*> where MMAX is the maximum value of M in MVAL and NSMAX is the
148*> maximum value of NRHS in NSVAL.
149*> \endverbatim
150*>
151*> \param[out] COPYB
152*> \verbatim
153*> COPYB is COMPLEX array, dimension (MMAX*NSMAX)
154*> \endverbatim
155*>
156*> \param[out] C
157*> \verbatim
158*> C is COMPLEX array, dimension (MMAX*NSMAX)
159*> \endverbatim
160*>
161*> \param[out] S
162*> \verbatim
163*> S is REAL array, dimension
164*> (min(MMAX,NMAX))
165*> \endverbatim
166*>
167*> \param[out] COPYS
168*> \verbatim
169*> COPYS is REAL array, dimension
170*> (min(MMAX,NMAX))
171*> \endverbatim
172*>
173*> \param[in] NOUT
174*> \verbatim
175*> NOUT is INTEGER
176*> The unit number for output.
177*> \endverbatim
178*
179* Authors:
180* ========
181*
182*> \author Univ. of Tennessee
183*> \author Univ. of California Berkeley
184*> \author Univ. of Colorado Denver
185*> \author NAG Ltd.
186*
187*> \ingroup complex_lin
188*
189* =====================================================================
190 SUBROUTINE cdrvls( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
191 $ NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
192 $ COPYB, C, S, COPYS, NOUT )
193*
194* -- LAPACK test routine --
195* -- LAPACK is a software package provided by Univ. of Tennessee, --
196* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
197*
198* .. Scalar Arguments ..
199 LOGICAL TSTERR
200 INTEGER NM, NN, NNB, NNS, NOUT
201 REAL THRESH
202* ..
203* .. Array Arguments ..
204 LOGICAL DOTYPE( * )
205 INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
206 $ nval( * ), nxval( * )
207 REAL COPYS( * ), S( * )
208 COMPLEX A( * ), B( * ), C( * ), COPYA( * ), COPYB( * )
209* ..
210*
211* =====================================================================
212*
213* .. Parameters ..
214 INTEGER NTESTS
215 PARAMETER ( NTESTS = 18 )
216 INTEGER SMLSIZ
217 parameter( smlsiz = 25 )
218 REAL ONE, ZERO
219 parameter( one = 1.0e+0, zero = 0.0e+0 )
220 COMPLEX CONE, CZERO
221 parameter( cone = ( 1.0e+0, 0.0e+0 ),
222 $ czero = ( 0.0e+0, 0.0e+0 ) )
223* ..
224* .. Local Scalars ..
225 CHARACTER TRANS
226 CHARACTER*3 PATH
227 INTEGER CRANK, I, IM, IMB, IN, INB, INFO, INS, IRANK,
228 $ iscale, itran, itype, j, k, lda, ldb, ldwork,
229 $ lwlsy, lwork, m, mnmin, n, nb, ncols, nerrs,
230 $ nfail, nrhs, nrows, nrun, rank, mb,
231 $ mmax, nmax, nsmax, liwork, lrwork,
232 $ lwork_cgels, lwork_cgelst, lwork_cgetsls,
233 $ lwork_cgelss, lwork_cgelsy, lwork_cgelsd,
234 $ lrwork_cgelsy, lrwork_cgelss, lrwork_cgelsd
235 REAL EPS, NORMA, NORMB, RCOND
236* ..
237* .. Local Arrays ..
238 INTEGER ISEED( 4 ), ISEEDY( 4 ), IWQ( 1 )
239 REAL RESULT( NTESTS ), RWQ( 1 )
240 COMPLEX WQ( 1 )
241* ..
242* .. Allocatable Arrays ..
243 COMPLEX, ALLOCATABLE :: WORK (:)
244 REAL, ALLOCATABLE :: RWORK (:), WORK2 (:)
245 INTEGER, ALLOCATABLE :: IWORK (:)
246* ..
247* .. External Functions ..
248 REAL CQRT12, CQRT14, CQRT17, SASUM, SLAMCH
249 EXTERNAL CQRT12, CQRT14, CQRT17, SASUM, SLAMCH
250* ..
251* .. External Subroutines ..
252 EXTERNAL alaerh, alahd, alasvm, cerrls, cgels, cgelsd,
253 $ cgelss, cgelst, cgelsy, cgemm, cgetsls, clacpy,
254 $ clarnv, cqrt13, cqrt15, cqrt16, csscal,
255 $ saxpy, xlaenv
256* ..
257* .. Intrinsic Functions ..
258 INTRINSIC max, min, int, real, sqrt
259* ..
260* .. Scalars in Common ..
261 LOGICAL LERR, OK
262 CHARACTER*32 SRNAMT
263 INTEGER INFOT, IOUNIT
264* ..
265* .. Common blocks ..
266 COMMON / infoc / infot, iounit, ok, lerr
267 COMMON / srnamc / srnamt
268* ..
269* .. Data statements ..
270 DATA iseedy / 1988, 1989, 1990, 1991 /
271* ..
272* .. Executable Statements ..
273*
274* Initialize constants and the random number seed.
275*
276 path( 1: 1 ) = 'Complex precision'
277 path( 2: 3 ) = 'LS'
278 nrun = 0
279 nfail = 0
280 nerrs = 0
281 DO 10 i = 1, 4
282 iseed( i ) = iseedy( i )
283 10 CONTINUE
284 eps = slamch( 'Epsilon' )
285*
286* Threshold for rank estimation
287*
288 rcond = sqrt( eps ) - ( sqrt( eps )-eps ) / 2
289*
290* Test the error exits
291*
292 CALL xlaenv( 9, smlsiz )
293 IF( tsterr )
294 $ CALL cerrls( path, nout )
295*
296* Print the header if NM = 0 or NN = 0 and THRESH = 0.
297*
298 IF( ( nm.EQ.0 .OR. nn.EQ.0 ) .AND. thresh.EQ.zero )
299 $ CALL alahd( nout, path )
300 infot = 0
301*
302* Compute maximal workspace needed for all routines
303*
304 nmax = 0
305 mmax = 0
306 nsmax = 0
307 DO i = 1, nm
308 IF ( mval( i ).GT.mmax ) THEN
309 mmax = mval( i )
310 END IF
311 ENDDO
312 DO i = 1, nn
313 IF ( nval( i ).GT.nmax ) THEN
314 nmax = nval( i )
315 END IF
316 ENDDO
317 DO i = 1, nns
318 IF ( nsval( i ).GT.nsmax ) THEN
319 nsmax = nsval( i )
320 END IF
321 ENDDO
322 m = mmax
323 n = nmax
324 nrhs = nsmax
325 mnmin = max( min( m, n ), 1 )
326*
327* Compute workspace needed for routines
328* CQRT14, CQRT17 (two side cases), CQRT15 and CQRT12
329*
330 lwork = max( 1, ( m+n )*nrhs,
331 $ ( n+nrhs )*( m+2 ), ( m+nrhs )*( n+2 ),
332 $ max( m+mnmin, nrhs*mnmin,2*n+m ),
333 $ max( m*n+4*mnmin+max(m,n), m*n+2*mnmin+4*n ) )
334 lrwork = 1
335 liwork = 1
336*
337* Iterate through all test cases and compute necessary workspace
338* sizes for ?GELS, ?GELST, ?GETSLS, ?GELSY, ?GELSS and ?GELSD
339* routines.
340*
341 DO im = 1, nm
342 m = mval( im )
343 lda = max( 1, m )
344 DO in = 1, nn
345 n = nval( in )
346 mnmin = max(min( m, n ),1)
347 ldb = max( 1, m, n )
348 DO ins = 1, nns
349 nrhs = nsval( ins )
350 DO irank = 1, 2
351 DO iscale = 1, 3
352 itype = ( irank-1 )*3 + iscale
353 IF( dotype( itype ) ) THEN
354 IF( irank.EQ.1 ) THEN
355 DO itran = 1, 2
356 IF( itran.EQ.1 ) THEN
357 trans = 'N'
358 ELSE
359 trans = 'C'
360 END IF
361*
362* Compute workspace needed for CGELS
363 CALL cgels( trans, m, n, nrhs, a, lda,
364 $ b, ldb, wq, -1, info )
365 lwork_cgels = int( wq( 1 ) )
366* Compute workspace needed for CGELST
367 CALL cgelst( trans, m, n, nrhs, a, lda,
368 $ b, ldb, wq, -1, info )
369 lwork_cgelst = int( wq( 1 ) )
370* Compute workspace needed for CGETSLS
371 CALL cgetsls( trans, m, n, nrhs, a, lda,
372 $ b, ldb, wq, -1, info )
373 lwork_cgetsls = int( wq( 1 ) )
374 ENDDO
375 END IF
376* Compute workspace needed for CGELSY
377 CALL cgelsy( m, n, nrhs, a, lda, b, ldb,
378 $ iwq, rcond, crank, wq, -1, rwq,
379 $ info )
380 lwork_cgelsy = int( wq( 1 ) )
381 lrwork_cgelsy = 2*n
382* Compute workspace needed for CGELSS
383 CALL cgelss( m, n, nrhs, a, lda, b, ldb, s,
384 $ rcond, crank, wq, -1, rwq, info )
385 lwork_cgelss = int( wq( 1 ) )
386 lrwork_cgelss = 5*mnmin
387* Compute workspace needed for CGELSD
388 CALL cgelsd( m, n, nrhs, a, lda, b, ldb, s,
389 $ rcond, crank, wq, -1, rwq, iwq,
390 $ info )
391 lwork_cgelsd = int( wq( 1 ) )
392 lrwork_cgelsd = int( rwq( 1 ) )
393* Compute LIWORK workspace needed for CGELSY and CGELSD
394 liwork = max( liwork, n, iwq( 1 ) )
395* Compute LRWORK workspace needed for CGELSY, CGELSS and CGELSD
396 lrwork = max( lrwork, lrwork_cgelsy,
397 $ lrwork_cgelss, lrwork_cgelsd )
398* Compute LWORK workspace needed for all functions
399 lwork = max( lwork, lwork_cgels, lwork_cgetsls,
400 $ lwork_cgelsy, lwork_cgelss,
401 $ lwork_cgelsd )
402 END IF
403 ENDDO
404 ENDDO
405 ENDDO
406 ENDDO
407 ENDDO
408*
409 lwlsy = lwork
410*
411 ALLOCATE( work( lwork ) )
412 ALLOCATE( iwork( liwork ) )
413 ALLOCATE( rwork( lrwork ) )
414 ALLOCATE( work2( 2 * lwork ) )
415*
416 DO 140 im = 1, nm
417 m = mval( im )
418 lda = max( 1, m )
419*
420 DO 130 in = 1, nn
421 n = nval( in )
422 mnmin = max(min( m, n ),1)
423 ldb = max( 1, m, n )
424 mb = (mnmin+1)
425*
426 DO 120 ins = 1, nns
427 nrhs = nsval( ins )
428*
429 DO 110 irank = 1, 2
430 DO 100 iscale = 1, 3
431 itype = ( irank-1 )*3 + iscale
432 IF( .NOT.dotype( itype ) )
433 $ GO TO 100
434* =====================================================
435* Begin test CGELS
436* =====================================================
437 IF( irank.EQ.1 ) THEN
438*
439* Generate a matrix of scaling type ISCALE
440*
441 CALL cqrt13( iscale, m, n, copya, lda, norma,
442 $ iseed )
443*
444* Loop for testing different block sizes.
445*
446 DO inb = 1, nnb
447 nb = nbval( inb )
448 CALL xlaenv( 1, nb )
449 CALL xlaenv( 3, nxval( inb ) )
450*
451* Loop for testing non-transposed and transposed.
452*
453 DO itran = 1, 2
454 IF( itran.EQ.1 ) THEN
455 trans = 'N'
456 nrows = m
457 ncols = n
458 ELSE
459 trans = 'C'
460 nrows = n
461 ncols = m
462 END IF
463 ldwork = max( 1, ncols )
464*
465* Set up a consistent rhs
466*
467 IF( ncols.GT.0 ) THEN
468 CALL clarnv( 2, iseed, ncols*nrhs,
469 $ work )
470 CALL csscal( ncols*nrhs,
471 $ one / real( ncols ), work,
472 $ 1 )
473 END IF
474 CALL cgemm( trans, 'No transpose', nrows,
475 $ nrhs, ncols, cone, copya, lda,
476 $ work, ldwork, czero, b, ldb )
477 CALL clacpy( 'Full', nrows, nrhs, b, ldb,
478 $ copyb, ldb )
479*
480* Solve LS or overdetermined system
481*
482 IF( m.GT.0 .AND. n.GT.0 ) THEN
483 CALL clacpy( 'Full', m, n, copya, lda,
484 $ a, lda )
485 CALL clacpy( 'Full', nrows, nrhs,
486 $ copyb, ldb, b, ldb )
487 END IF
488 srnamt = 'CGELS '
489 CALL cgels( trans, m, n, nrhs, a, lda, b,
490 $ ldb, work, lwork, info )
491*
492 IF( info.NE.0 )
493 $ CALL alaerh( path, 'CGELS ', info, 0,
494 $ trans, m, n, nrhs, -1, nb,
495 $ itype, nfail, nerrs,
496 $ nout )
497*
498* Test 1: Check correctness of results
499* for CGELS, compute the residual:
500* RESID = norm(B - A*X) /
501* / ( max(m,n) * norm(A) * norm(X) * EPS )
502*
503 IF( nrows.GT.0 .AND. nrhs.GT.0 )
504 $ CALL clacpy( 'Full', nrows, nrhs,
505 $ copyb, ldb, c, ldb )
506 CALL cqrt16( trans, m, n, nrhs, copya,
507 $ lda, b, ldb, c, ldb, rwork,
508 $ result( 1 ) )
509*
510* Test 2: Check correctness of results
511* for CGELS.
512*
513 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
514 $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
515*
516* Solving LS system
517*
518 result( 2 ) = cqrt17( trans, 1, m, n,
519 $ nrhs, copya, lda, b, ldb,
520 $ copyb, ldb, c, work,
521 $ lwork )
522 ELSE
523*
524* Solving overdetermined system
525*
526 result( 2 ) = cqrt14( trans, m, n,
527 $ nrhs, copya, lda, b, ldb,
528 $ work, lwork )
529 END IF
530*
531* Print information about the tests that
532* did not pass the threshold.
533*
534 DO k = 1, 2
535 IF( result( k ).GE.thresh ) THEN
536 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
537 $ CALL alahd( nout, path )
538 WRITE( nout, fmt = 9999 )trans, m,
539 $ n, nrhs, nb, itype, k,
540 $ result( k )
541 nfail = nfail + 1
542 END IF
543 END DO
544 nrun = nrun + 2
545 END DO
546 END DO
547 END IF
548* =====================================================
549* End test CGELS
550* =====================================================
551* =====================================================
552* Begin test CGELST
553* =====================================================
554 IF( irank.EQ.1 ) THEN
555*
556* Generate a matrix of scaling type ISCALE
557*
558 CALL cqrt13( iscale, m, n, copya, lda, norma,
559 $ iseed )
560*
561* Loop for testing different block sizes.
562*
563 DO inb = 1, nnb
564 nb = nbval( inb )
565 CALL xlaenv( 1, nb )
566 CALL xlaenv( 3, nxval( inb ) )
567*
568* Loop for testing non-transposed and transposed.
569*
570 DO itran = 1, 2
571 IF( itran.EQ.1 ) THEN
572 trans = 'N'
573 nrows = m
574 ncols = n
575 ELSE
576 trans = 'C'
577 nrows = n
578 ncols = m
579 END IF
580 ldwork = max( 1, ncols )
581*
582* Set up a consistent rhs
583*
584 IF( ncols.GT.0 ) THEN
585 CALL clarnv( 2, iseed, ncols*nrhs,
586 $ work )
587 CALL csscal( ncols*nrhs,
588 $ one / real( ncols ), work,
589 $ 1 )
590 END IF
591 CALL cgemm( trans, 'No transpose', nrows,
592 $ nrhs, ncols, cone, copya, lda,
593 $ work, ldwork, czero, b, ldb )
594 CALL clacpy( 'Full', nrows, nrhs, b, ldb,
595 $ copyb, ldb )
596*
597* Solve LS or overdetermined system
598*
599 IF( m.GT.0 .AND. n.GT.0 ) THEN
600 CALL clacpy( 'Full', m, n, copya, lda,
601 $ a, lda )
602 CALL clacpy( 'Full', nrows, nrhs,
603 $ copyb, ldb, b, ldb )
604 END IF
605 srnamt = 'CGELST'
606 CALL cgelst( trans, m, n, nrhs, a, lda, b,
607 $ ldb, work, lwork, info )
608*
609 IF( info.NE.0 )
610 $ CALL alaerh( path, 'CGELST', info, 0,
611 $ trans, m, n, nrhs, -1, nb,
612 $ itype, nfail, nerrs,
613 $ nout )
614*
615* Test 3: Check correctness of results
616* for CGELST, compute the residual:
617* RESID = norm(B - A*X) /
618* / ( max(m,n) * norm(A) * norm(X) * EPS )
619*
620 IF( nrows.GT.0 .AND. nrhs.GT.0 )
621 $ CALL clacpy( 'Full', nrows, nrhs,
622 $ copyb, ldb, c, ldb )
623 CALL cqrt16( trans, m, n, nrhs, copya,
624 $ lda, b, ldb, c, ldb, rwork,
625 $ result( 3 ) )
626*
627* Test 4: Check correctness of results
628* for CGELST.
629*
630 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
631 $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
632*
633* Solving LS system
634*
635 result( 4 ) = cqrt17( trans, 1, m, n,
636 $ nrhs, copya, lda, b, ldb,
637 $ copyb, ldb, c, work,
638 $ lwork )
639 ELSE
640*
641* Solving overdetermined system
642*
643 result( 4 ) = cqrt14( trans, m, n,
644 $ nrhs, copya, lda, b, ldb,
645 $ work, lwork )
646 END IF
647*
648* Print information about the tests that
649* did not pass the threshold.
650*
651 DO k = 3, 4
652 IF( result( k ).GE.thresh ) THEN
653 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
654 $ CALL alahd( nout, path )
655 WRITE( nout, fmt = 9999 )trans, m,
656 $ n, nrhs, nb, itype, k,
657 $ result( k )
658 nfail = nfail + 1
659 END IF
660 END DO
661 nrun = nrun + 2
662 END DO
663 END DO
664 END IF
665* =====================================================
666* End test CGELST
667* =====================================================
668* =====================================================
669* Begin test CGELSTSLS
670* =====================================================
671 IF( irank.EQ.1 ) THEN
672*
673* Generate a matrix of scaling type ISCALE
674*
675 CALL cqrt13( iscale, m, n, copya, lda, norma,
676 $ iseed )
677*
678* Loop for testing different block sizes MB.
679*
680 DO inb = 1, nnb
681 mb = nbval( inb )
682 CALL xlaenv( 1, mb )
683*
684* Loop for testing different block sizes NB.
685*
686 DO imb = 1, nnb
687 nb = nbval( imb )
688 CALL xlaenv( 2, nb )
689*
690* Loop for testing non-transposed
691* and transposed.
692*
693 DO itran = 1, 2
694 IF( itran.EQ.1 ) THEN
695 trans = 'N'
696 nrows = m
697 ncols = n
698 ELSE
699 trans = 'C'
700 nrows = n
701 ncols = m
702 END IF
703 ldwork = max( 1, ncols )
704*
705* Set up a consistent rhs
706*
707 IF( ncols.GT.0 ) THEN
708 CALL clarnv( 2, iseed, ncols*nrhs,
709 $ work )
710 CALL cscal( ncols*nrhs,
711 $ cone / real( ncols ),
712 $ work, 1 )
713 END IF
714 CALL cgemm( trans, 'No transpose',
715 $ nrows, nrhs, ncols, cone,
716 $ copya, lda, work, ldwork,
717 $ czero, b, ldb )
718 CALL clacpy( 'Full', nrows, nrhs,
719 $ b, ldb, copyb, ldb )
720*
721* Solve LS or overdetermined system
722*
723 IF( m.GT.0 .AND. n.GT.0 ) THEN
724 CALL clacpy( 'Full', m, n,
725 $ copya, lda, a, lda )
726 CALL clacpy( 'Full', nrows, nrhs,
727 $ copyb, ldb, b, ldb )
728 END IF
729 srnamt = 'CGETSLS '
730 CALL cgetsls( trans, m, n, nrhs, a,
731 $ lda, b, ldb, work, lwork,
732 $ info )
733 IF( info.NE.0 )
734 $ CALL alaerh( path, 'CGETSLS ', info,
735 $ 0, trans, m, n, nrhs,
736 $ -1, nb, itype, nfail,
737 $ nerrs, nout )
738*
739* Test 5: Check correctness of results
740* for CGETSLS, compute the residual:
741* RESID = norm(B - A*X) /
742* / ( max(m,n) * norm(A) * norm(X) * EPS )
743*
744 IF( nrows.GT.0 .AND. nrhs.GT.0 )
745 $ CALL clacpy( 'Full', nrows, nrhs,
746 $ copyb, ldb, c, ldb )
747 CALL cqrt16( trans, m, n, nrhs,
748 $ copya, lda, b, ldb,
749 $ c, ldb, work2,
750 $ result( 5 ) )
751*
752* Test 6: Check correctness of results
753* for CGETSLS.
754*
755 IF( ( itran.EQ.1 .AND. m.GE.n ) .OR.
756 $ ( itran.EQ.2 .AND. m.LT.n ) ) THEN
757*
758* Solving LS system, compute:
759* r = norm((B- A*X)**T * A) /
760* / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
761*
762 result( 6 ) = cqrt17( trans, 1, m,
763 $ n, nrhs, copya, lda,
764 $ b, ldb, copyb, ldb,
765 $ c, work, lwork )
766 ELSE
767*
768* Solving overdetermined system
769*
770 result( 6 ) = cqrt14( trans, m, n,
771 $ nrhs, copya, lda, b,
772 $ ldb, work, lwork )
773 END IF
774*
775* Print information about the tests that
776* did not pass the threshold.
777*
778 DO k = 5, 6
779 IF( result( k ).GE.thresh ) THEN
780 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
781 $ CALL alahd( nout, path )
782 WRITE( nout, fmt = 9997 )trans,
783 $ m, n, nrhs, mb, nb, itype, k,
784 $ result( k )
785 nfail = nfail + 1
786 END IF
787 END DO
788 nrun = nrun + 2
789 END DO
790 END DO
791 END DO
792 END IF
793* =====================================================
794* End test CGELSTSLS
795* ====================================================
796*
797* Generate a matrix of scaling type ISCALE and rank
798* type IRANK.
799*
800 CALL cqrt15( iscale, irank, m, n, nrhs, copya, lda,
801 $ copyb, ldb, copys, rank, norma, normb,
802 $ iseed, work, lwork )
803*
804* workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
805*
806 ldwork = max( 1, m )
807*
808* Loop for testing different block sizes.
809*
810 DO 90 inb = 1, nnb
811 nb = nbval( inb )
812 CALL xlaenv( 1, nb )
813 CALL xlaenv( 3, nxval( inb ) )
814*
815* Test CGELSY
816*
817* CGELSY: Compute the minimum-norm solution
818* X to min( norm( A * X - B ) )
819* using the rank-revealing orthogonal
820* factorization.
821*
822 CALL clacpy( 'Full', m, n, copya, lda, a, lda )
823 CALL clacpy( 'Full', m, nrhs, copyb, ldb, b,
824 $ ldb )
825*
826* Initialize vector IWORK.
827*
828 DO 70 j = 1, n
829 iwork( j ) = 0
830 70 CONTINUE
831*
832 srnamt = 'CGELSY'
833 CALL cgelsy( m, n, nrhs, a, lda, b, ldb, iwork,
834 $ rcond, crank, work, lwlsy, rwork,
835 $ info )
836 IF( info.NE.0 )
837 $ CALL alaerh( path, 'CGELSY', info, 0, ' ', m,
838 $ n, nrhs, -1, nb, itype, nfail,
839 $ nerrs, nout )
840*
841* workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS)
842*
843* Test 7: Compute relative error in svd
844* workspace: M*N + 4*MIN(M,N) + MAX(M,N)
845*
846 result( 7 ) = cqrt12( crank, crank, a, lda,
847 $ copys, work, lwork, rwork )
848*
849* Test 8: Compute error in solution
850* workspace: M*NRHS + M
851*
852 CALL clacpy( 'Full', m, nrhs, copyb, ldb, work,
853 $ ldwork )
854 CALL cqrt16( 'No transpose', m, n, nrhs, copya,
855 $ lda, b, ldb, work, ldwork, rwork,
856 $ result( 8 ) )
857*
858* Test 9: Check norm of r'*A
859* workspace: NRHS*(M+N)
860*
861 result( 9 ) = zero
862 IF( m.GT.crank )
863 $ result( 9 ) = cqrt17( 'No transpose', 1, m,
864 $ n, nrhs, copya, lda, b, ldb,
865 $ copyb, ldb, c, work, lwork )
866*
867* Test 10: Check if x is in the rowspace of A
868* workspace: (M+NRHS)*(N+2)
869*
870 result( 10 ) = zero
871*
872 IF( n.GT.crank )
873 $ result( 10 ) = cqrt14( 'No transpose', m, n,
874 $ nrhs, copya, lda, b, ldb,
875 $ work, lwork )
876*
877* Test CGELSS
878*
879* CGELSS: Compute the minimum-norm solution
880* X to min( norm( A * X - B ) )
881* using the SVD.
882*
883 CALL clacpy( 'Full', m, n, copya, lda, a, lda )
884 CALL clacpy( 'Full', m, nrhs, copyb, ldb, b,
885 $ ldb )
886 srnamt = 'CGELSS'
887 CALL cgelss( m, n, nrhs, a, lda, b, ldb, s,
888 $ rcond, crank, work, lwork, rwork,
889 $ info )
890*
891 IF( info.NE.0 )
892 $ CALL alaerh( path, 'CGELSS', info, 0, ' ', m,
893 $ n, nrhs, -1, nb, itype, nfail,
894 $ nerrs, nout )
895*
896* workspace used: 3*min(m,n) +
897* max(2*min(m,n),nrhs,max(m,n))
898*
899* Test 11: Compute relative error in svd
900*
901 IF( rank.GT.0 ) THEN
902 CALL saxpy( mnmin, -one, copys, 1, s, 1 )
903 result( 11 ) = sasum( mnmin, s, 1 ) /
904 $ sasum( mnmin, copys, 1 ) /
905 $ ( eps*real( mnmin ) )
906 ELSE
907 result( 11 ) = zero
908 END IF
909*
910* Test 12: Compute error in solution
911*
912 CALL clacpy( 'Full', m, nrhs, copyb, ldb, work,
913 $ ldwork )
914 CALL cqrt16( 'No transpose', m, n, nrhs, copya,
915 $ lda, b, ldb, work, ldwork, rwork,
916 $ result( 12 ) )
917*
918* Test 13: Check norm of r'*A
919*
920 result( 13 ) = zero
921 IF( m.GT.crank )
922 $ result( 13 ) = cqrt17( 'No transpose', 1, m,
923 $ n, nrhs, copya, lda, b, ldb,
924 $ copyb, ldb, c, work, lwork )
925*
926* Test 14: Check if x is in the rowspace of A
927*
928 result( 14 ) = zero
929 IF( n.GT.crank )
930 $ result( 14 ) = cqrt14( 'No transpose', m, n,
931 $ nrhs, copya, lda, b, ldb,
932 $ work, lwork )
933*
934* Test CGELSD
935*
936* CGELSD: Compute the minimum-norm solution X
937* to min( norm( A * X - B ) ) using a
938* divide and conquer SVD.
939*
940 CALL xlaenv( 9, 25 )
941*
942 CALL clacpy( 'Full', m, n, copya, lda, a, lda )
943 CALL clacpy( 'Full', m, nrhs, copyb, ldb, b,
944 $ ldb )
945*
946 srnamt = 'CGELSD'
947 CALL cgelsd( m, n, nrhs, a, lda, b, ldb, s,
948 $ rcond, crank, work, lwork, rwork,
949 $ iwork, info )
950 IF( info.NE.0 )
951 $ CALL alaerh( path, 'CGELSD', info, 0, ' ', m,
952 $ n, nrhs, -1, nb, itype, nfail,
953 $ nerrs, nout )
954*
955* Test 15: Compute relative error in svd
956*
957 IF( rank.GT.0 ) THEN
958 CALL saxpy( mnmin, -one, copys, 1, s, 1 )
959 result( 15 ) = sasum( mnmin, s, 1 ) /
960 $ sasum( mnmin, copys, 1 ) /
961 $ ( eps*real( mnmin ) )
962 ELSE
963 result( 15 ) = zero
964 END IF
965*
966* Test 16: Compute error in solution
967*
968 CALL clacpy( 'Full', m, nrhs, copyb, ldb, work,
969 $ ldwork )
970 CALL cqrt16( 'No transpose', m, n, nrhs, copya,
971 $ lda, b, ldb, work, ldwork, rwork,
972 $ result( 16 ) )
973*
974* Test 17: Check norm of r'*A
975*
976 result( 17 ) = zero
977 IF( m.GT.crank )
978 $ result( 17 ) = cqrt17( 'No transpose', 1, m,
979 $ n, nrhs, copya, lda, b, ldb,
980 $ copyb, ldb, c, work, lwork )
981*
982* Test 18: Check if x is in the rowspace of A
983*
984 result( 18 ) = zero
985 IF( n.GT.crank )
986 $ result( 18 ) = cqrt14( 'No transpose', m, n,
987 $ nrhs, copya, lda, b, ldb,
988 $ work, lwork )
989*
990* Print information about the tests that did not
991* pass the threshold.
992*
993 DO 80 k = 7, 18
994 IF( result( k ).GE.thresh ) THEN
995 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
996 $ CALL alahd( nout, path )
997 WRITE( nout, fmt = 9998 )m, n, nrhs, nb,
998 $ itype, k, result( k )
999 nfail = nfail + 1
1000 END IF
1001 80 CONTINUE
1002 nrun = nrun + 12
1003*
1004 90 CONTINUE
1005 100 CONTINUE
1006 110 CONTINUE
1007 120 CONTINUE
1008 130 CONTINUE
1009 140 CONTINUE
1010*
1011* Print a summary of the results.
1012*
1013 CALL alasvm( path, nout, nfail, nrun, nerrs )
1014*
1015 9999 FORMAT( ' TRANS=''', a1, ''', M=', i5, ', N=', i5, ', NRHS=', i4,
1016 $ ', NB=', i4, ', type', i2, ', test(', i2, ')=', g12.5 )
1017 9998 FORMAT( ' M=', i5, ', N=', i5, ', NRHS=', i4, ', NB=', i4,
1018 $ ', type', i2, ', test(', i2, ')=', g12.5 )
1019 9997 FORMAT( ' TRANS=''', a1,' M=', i5, ', N=', i5, ', NRHS=', i4,
1020 $ ', MB=', i4,', NB=', i4,', type', i2,
1021 $ ', test(', i2, ')=', g12.5 )
1022*
1023 DEALLOCATE( work )
1024 DEALLOCATE( rwork )
1025 DEALLOCATE( iwork )
1026 RETURN
1027*
1028* End of CDRVLS
1029*
1030 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
cdrvls
subroutine cdrvls(dotype, nm, mval, nn, nval, nns, nsval, nnb, nbval, nxval, thresh, tsterr, a, copya, b, copyb, c, s, copys, nout)
CDRVLS
Definition cdrvls.f:193
cerrls
subroutine cerrls(path, nunit)
CERRLS
Definition cerrls.f:55
cqrt13
subroutine cqrt13(scale, m, n, a, lda, norma, iseed)
CQRT13
Definition cqrt13.f:91
cqrt15
subroutine cqrt15(scale, rksel, m, n, nrhs, a, lda, b, ldb, s, rank, norma, normb, iseed, work, lwork)
CQRT15
Definition cqrt15.f:149
cqrt16
subroutine cqrt16(trans, m, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
CQRT16
Definition cqrt16.f:133
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
cgels
subroutine cgels(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
CGELS solves overdetermined or underdetermined systems for GE matrices
Definition cgels.f:191
cgelsd
subroutine cgelsd(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, rwork, iwork, info)
CGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices
Definition cgelsd.f:217
cgelss
subroutine cgelss(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, rwork, info)
CGELSS solves overdetermined or underdetermined systems for GE matrices
Definition cgelss.f:176
cgelst
subroutine cgelst(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
CGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization ...
Definition cgelst.f:202
cgelsy
subroutine cgelsy(m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, lwork, rwork, info)
CGELSY solves overdetermined or underdetermined systems for GE matrices
Definition cgelsy.f:211
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
cgetsls
subroutine cgetsls(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
CGETSLS
Definition cgetsls.f:171
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:101
clarnv
subroutine clarnv(idist, iseed, n, x)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition clarnv.f:97
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78