You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

cdrvls.f 42 kB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030
  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