LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
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,
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
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine xlaenv(ispec, nvalue)
XLAENV
Definition xlaenv.f:81
subroutine alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
Definition alaerh.f:147
subroutine alahd(iounit, path)
ALAHD
Definition alahd.f:107
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
subroutine cerrls(path, nunit)
CERRLS
Definition cerrls.f:55
subroutine cqrt13(scale, m, n, a, lda, norma, iseed)
CQRT13
Definition cqrt13.f:91
subroutine cqrt15(scale, rksel, m, n, nrhs, a, lda, b, ldb, s, rank, norma, normb, iseed, work, lwork)
CQRT15
Definition cqrt15.f:149
subroutine cqrt16(trans, m, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
CQRT16
Definition cqrt16.f:133
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
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:182
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:219
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:178
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:194
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:212
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
subroutine cgetsls(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
CGETSLS
Definition cgetsls.f:162
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
subroutine clarnv(idist, iseed, n, x)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition clarnv.f:99
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78