OpenRadioss 2025.1.11
OpenRadioss project
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, CGETSLS, CGELSS, CGELSY
35*> and CGELSD.
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 REAL
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 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 array, dimension (MMAX*NMAX)
141*> \endverbatim
142*>
143*> \param[out] B
144*> \verbatim
145*> B is COMPLEX 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 array, dimension (MMAX*NSMAX)
153*> \endverbatim
154*>
155*> \param[out] C
156*> \verbatim
157*> C is COMPLEX array, dimension (MMAX*NSMAX)
158*> \endverbatim
159*>
160*> \param[out] S
161*> \verbatim
162*> S is REAL array, dimension
163*> (min(MMAX,NMAX))
164*> \endverbatim
165*>
166*> \param[out] COPYS
167*> \verbatim
168*> COPYS is REAL 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 complex_lin
187*
188* =====================================================================
189 SUBROUTINE cdrvls( 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 REAL THRESH
201* ..
202* .. Array Arguments ..
203 LOGICAL DOTYPE( * )
204 INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
205 $ nval( * ), nxval( * )
206 REAL COPYS( * ), S( * )
207 COMPLEX 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 REAL ONE, ZERO
218 parameter( one = 1.0e+0, zero = 0.0e+0 )
219 COMPLEX CONE, CZERO
220 parameter( cone = ( 1.0e+0, 0.0e+0 ),
221 $ czero = ( 0.0e+0, 0.0e+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_cgels, lwork_cgetsls, lwork_cgelss,
232 $ lwork_cgelsy, lwork_cgelsd,
233 $ lrwork_cgelsy, lrwork_cgelss, lrwork_cgelsd
234 REAL EPS, NORMA, NORMB, RCOND
235* ..
236* .. Local Arrays ..
237 INTEGER ISEED( 4 ), ISEEDY( 4 ), IWQ( 1 )
238 REAL RESULT( NTESTS ), RWQ( 1 )
239 COMPLEX WQ( 1 )
240* ..
241* .. Allocatable Arrays ..
242 COMPLEX, ALLOCATABLE :: WORK (:)
243 REAL, ALLOCATABLE :: RWORK (:), WORK2 (:)
244 INTEGER, ALLOCATABLE :: IWORK (:)
245* ..
246* .. External Functions ..
247 REAL CQRT12, CQRT14, CQRT17, SASUM, SLAMCH
248 EXTERNAL CQRT12, CQRT14, CQRT17, SASUM, SLAMCH
249* ..
250* .. External Subroutines ..
251 EXTERNAL alaerh, alahd, alasvm, cerrls, cgels, cgelsd,
252 $ cgelss, cgelsy, cgemm, cgetsls, clacpy,
253 $ clarnv, cqrt13, cqrt15, cqrt16, csscal,
254 $ saxpy, xlaenv
255* ..
256* .. Intrinsic Functions ..
257 INTRINSIC max, min, int, real, 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 ) = 'Complex 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 = slamch( '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 cerrls( 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* CQRT14, CQRT17 (two side cases), CQRT15 and CQRT12
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 CGELS
361 CALL cgels( trans, m, n, nrhs, a, lda,
362 $ b, ldb, wq, -1, info )
363 lwork_cgels = int( wq( 1 ) )
364* Compute workspace needed for CGETSLS
365 CALL cgetsls( trans, m, n, nrhs, a, lda,
366 $ b, ldb, wq, -1, info )
367 lwork_cgetsls = int( wq( 1 ) )
368 ENDDO
369 END IF
370* Compute workspace needed for CGELSY
371 CALL cgelsy( m, n, nrhs, a, lda, b, ldb,
372 $ iwq, rcond, crank, wq, -1, rwq,
373 $ info )
374 lwork_cgelsy = int( wq( 1 ) )
375 lrwork_cgelsy = 2*n
376* Compute workspace needed for CGELSS
377 CALL cgelss( m, n, nrhs, a, lda, b, ldb, s,
378 $ rcond, crank, wq, -1, rwq, info )
379 lwork_cgelss = int( wq( 1 ) )
380 lrwork_cgelss = 5*mnmin
381* Compute workspace needed for CGELSD
382 CALL cgelsd( m, n, nrhs, a, lda, b, ldb, s,
383 $ rcond, crank, wq, -1, rwq, iwq,
384 $ info )
385 lwork_cgelsd = int( wq( 1 ) )
386 lrwork_cgelsd = int( rwq( 1 ) )
387* Compute LIWORK workspace needed for CGELSY and CGELSD
388 liwork = max( liwork, n, iwq( 1 ) )
389* Compute LRWORK workspace needed for CGELSY, CGELSS and CGELSD
390 lrwork = max( lrwork, lrwork_cgelsy,
391 $ lrwork_cgelss, lrwork_cgelsd )
392* Compute LWORK workspace needed for all functions
393 lwork = max( lwork, lwork_cgels, lwork_cgetsls,
394 $ lwork_cgelsy, lwork_cgelss,
395 $ lwork_cgelsd )
396 END IF
397 ENDDO
398 ENDDO
399 ENDDO
400 ENDDO
401 ENDDO
402*
403 lwlsy = lwork
404*
405 ALLOCATE( work( lwork ) )
406 ALLOCATE( iwork( liwork ) )
407 ALLOCATE( rwork( lrwork ) )
408 ALLOCATE( work2( 2 * lwork ) )
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 CGELS
432*
433* Generate a matrix of scaling type ISCALE
434*
435 CALL cqrt13( 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 clarnv( 2, iseed, ncols*nrhs,
458 $ work )
459 CALL csscal( ncols*nrhs,
460 $ one / real( ncols ), work,
461 $ 1 )
462 END IF
463 CALL cgemm( trans, 'No transpose', nrows,
464 $ nrhs, ncols, cone, copya, lda,
465 $ work, ldwork, czero, b, ldb )
466 CALL clacpy( '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 clacpy( 'Full', m, n, copya, lda,
473 $ a, lda )
474 CALL clacpy( 'Full', nrows, nrhs,
475 $ copyb, ldb, b, ldb )
476 END IF
477 srnamt = 'CGELS '
478 CALL cgels( trans, m, n, nrhs, a, lda, b,
479 $ ldb, work, lwork, info )
480*
481 IF( info.NE.0 )
482 $ CALL alaerh( path, 'CGELS ', 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 clacpy( 'Full', nrows, nrhs,
492 $ copyb, ldb, c, ldb )
493 CALL cqrt16( 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 ) = cqrt17( 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 ) = cqrt14( 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 CGETSLS
534*
535* Generate a matrix of scaling type ISCALE
536*
537 CALL cqrt13( 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.GT. IF( NCOLS0 ) THEN
561 CALL CLARNV( 2, ISEED, NCOLS*NRHS,
562 $ WORK )
563 CALL CSCAL( NCOLS*NRHS,
564 $ CONE / REAL( NCOLS ), WORK,
565 $ 1 )
566 END IF
567 CALL CGEMM( TRANS, 'no transpose', NROWS,
568 $ NRHS, NCOLS, CONE, COPYA, LDA,
569 $ WORK, LDWORK, CZERO, B, LDB )
570 CALL CLACPY( 'full', NROWS, NRHS, B, LDB,
571 $ COPYB, LDB )
572*
573* Solve LS or overdetermined system
574*
575.GT..AND..GT. IF( M0 N0 ) THEN
576 CALL CLACPY( 'full', M, N, COPYA, LDA,
577 $ A, LDA )
578 CALL CLACPY( 'full', NROWS, NRHS,
579 $ COPYB, LDB, B, LDB )
580 END IF
581 SRNAMT = 'cgetsls '
582 CALL CGETSLS( TRANS, M, N, NRHS, A,
583 $ LDA, B, LDB, WORK, LWORK, INFO )
584.NE. IF( INFO0 )
585 $ CALL ALAERH( PATH, 'cgetsls ', 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.GT..AND..GT. IF( NROWS0 NRHS0 )
594 $ CALL CLACPY( 'full', NROWS, NRHS,
595 $ COPYB, LDB, C, LDB )
596 CALL CQRT16( TRANS, M, N, NRHS, COPYA,
597 $ LDA, B, LDB, C, LDB, WORK2,
598 $ RESULT( 15 ) )
599*
600.EQ..AND..GE..OR. IF( ( ITRAN1 MN )
601.EQ..AND..LT. $ ( ITRAN2 MN ) ) THEN
602*
603* Solving LS system
604*
605 RESULT( 16 ) = CQRT17( 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 ) = CQRT14( 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.GE. IF( RESULT( K )THRESH ) THEN
623.EQ..AND..EQ. IF( NFAIL0 NERRS0 )
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 CQRT15( 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 CGELSY
656*
657* CGELSY: Compute the minimum-norm solution
658* X to min( norm( A * X - B ) )
659* using the rank-revealing orthogonal
660* factorization.
661*
662 CALL CLACPY( 'full', M, N, COPYA, LDA, A, LDA )
663 CALL CLACPY( '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 = 'cgelsy'
673 CALL CGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
674 $ RCOND, CRANK, WORK, LWLSY, RWORK,
675 $ INFO )
676.NE. IF( INFO0 )
677 $ CALL ALAERH( PATH, 'cgelsy', 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 ) = CQRT12( 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 CLACPY( 'full', M, NRHS, COPYB, LDB, WORK,
693 $ LDWORK )
694 CALL CQRT16( '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.GT. IF( MCRANK )
703 $ RESULT( 5 ) = CQRT17( '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.GT. IF( NCRANK )
713 $ RESULT( 6 ) = CQRT14( 'no transpose', M, N,
714 $ NRHS, COPYA, LDA, B, LDB,
715 $ WORK, LWORK )
716*
717* Test CGELSS
718*
719* CGELSS: Compute the minimum-norm solution
720* X to min( norm( A * X - B ) )
721* using the SVD.
722*
723 CALL CLACPY( 'full', M, N, COPYA, LDA, A, LDA )
724 CALL CLACPY( 'full', M, NRHS, COPYB, LDB, B,
725 $ LDB )
726 SRNAMT = 'cgelss'
727 CALL CGELSS( M, N, NRHS, A, LDA, B, LDB, S,
728 $ RCOND, CRANK, WORK, LWORK, RWORK,
729 $ INFO )
730*
731.NE. IF( INFO0 )
732 $ CALL ALAERH( PATH, 'cgelss', 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.GT. IF( RANK0 ) THEN
742 CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
743 RESULT( 7 ) = SASUM( MNMIN, S, 1 ) /
744 $ SASUM( MNMIN, COPYS, 1 ) /
745 $ ( EPS*REAL( MNMIN ) )
746 ELSE
747 RESULT( 7 ) = ZERO
748 END IF
749*
750* Test 8: Compute error in solution
751*
752 CALL CLACPY( 'full', M, NRHS, COPYB, LDB, WORK,
753 $ LDWORK )
754 CALL CQRT16( '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.GT. IF( MCRANK )
762 $ RESULT( 9 ) = CQRT17( '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.GT. IF( NCRANK )
770 $ RESULT( 10 ) = CQRT14( 'no transpose', M, N,
771 $ NRHS, COPYA, LDA, B, LDB,
772 $ WORK, LWORK )
773*
774* Test CGELSD
775*
776* CGELSD: 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 CLACPY( 'full', M, N, COPYA, LDA, A, LDA )
783 CALL CLACPY( 'full', M, NRHS, COPYB, LDB, B,
784 $ LDB )
785*
786 SRNAMT = 'cgelsd'
787 CALL CGELSD( M, N, NRHS, A, LDA, B, LDB, S,
788 $ RCOND, CRANK, WORK, LWORK, RWORK,
789 $ IWORK, INFO )
790.NE. IF( INFO0 )
791 $ CALL ALAERH( PATH, 'cgelsd', INFO, 0, ' ', M,
792 $ N, NRHS, -1, NB, ITYPE, NFAIL,
793 $ NERRS, NOUT )
794*
795* Test 11: Compute relative error in svd
796*
797.GT. IF( RANK0 ) THEN
798 CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
799 RESULT( 11 ) = SASUM( MNMIN, S, 1 ) /
800 $ SASUM( MNMIN, COPYS, 1 ) /
801 $ ( EPS*REAL( MNMIN ) )
802 ELSE
803 RESULT( 11 ) = ZERO
804 END IF
805*
806* Test 12: Compute error in solution
807*
808 CALL CLACPY( 'full', M, NRHS, COPYB, LDB, WORK,
809 $ LDWORK )
810 CALL CQRT16( '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.GT. IF( MCRANK )
818 $ RESULT( 13 ) = CQRT17( '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.GT. IF( NCRANK )
826 $ RESULT( 14 ) = CQRT14( '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.GE. IF( RESULT( K )THRESH ) THEN
835.EQ..AND..EQ. IF( NFAIL0 NERRS0 )
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( rwork )
865 DEALLOCATE( iwork )
866 RETURN
867*
868* End of CDRVLS
869*
870 END
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
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:182
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:210
cgetsls
subroutine cgetsls(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
CGETSLS
Definition cgetsls.f:162
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:178
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:225
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:103
clarnv
subroutine clarnv(idist, iseed, n, x)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition clarnv.f:99
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:187
cqrt13
subroutine cqrt13(scale, m, n, a, lda, norma, iseed)
CQRT13
Definition cqrt13.f:91
cqrt16
subroutine cqrt16(trans, m, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
CQRT16
Definition cqrt16.f:133
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:192
cerrls
subroutine cerrls(path, nunit)
CERRLS
Definition cerrls.f:55
cqrt15
subroutine cqrt15(scale, rksel, m, n, nrhs, a, lda, b, ldb, s, rank, norma, normb, iseed, work, lwork)
CQRT15
Definition cqrt15.f:149
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21