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 34 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873
  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. *> \date June 2017
  187. *
  188. *> \ingroup complex_lin
  189. *
  190. * =====================================================================
  191. SUBROUTINE CDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
  192. $ NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
  193. $ COPYB, C, S, COPYS, NOUT )
  194. *
  195. * -- LAPACK test routine (version 3.7.1) --
  196. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  197. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  198. * June 2017
  199. *
  200. * .. Scalar Arguments ..
  201. LOGICAL TSTERR
  202. INTEGER NM, NN, NNB, NNS, NOUT
  203. REAL THRESH
  204. * ..
  205. * .. Array Arguments ..
  206. LOGICAL DOTYPE( * )
  207. INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
  208. $ NVAL( * ), NXVAL( * )
  209. REAL COPYS( * ), S( * )
  210. COMPLEX A( * ), B( * ), C( * ), COPYA( * ), COPYB( * )
  211. * ..
  212. *
  213. * =====================================================================
  214. *
  215. * .. Parameters ..
  216. INTEGER NTESTS
  217. PARAMETER ( NTESTS = 16 )
  218. INTEGER SMLSIZ
  219. PARAMETER ( SMLSIZ = 25 )
  220. REAL ONE, ZERO
  221. PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
  222. COMPLEX CONE, CZERO
  223. PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ),
  224. $ CZERO = ( 0.0E+0, 0.0E+0 ) )
  225. * ..
  226. * .. Local Scalars ..
  227. CHARACTER TRANS
  228. CHARACTER*3 PATH
  229. INTEGER CRANK, I, IM, IMB, IN, INB, INFO, INS, IRANK,
  230. $ ISCALE, ITRAN, ITYPE, J, K, LDA, LDB, LDWORK,
  231. $ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS,
  232. $ NFAIL, NRHS, NROWS, NRUN, RANK, MB,
  233. $ MMAX, NMAX, NSMAX, LIWORK, LRWORK,
  234. $ LWORK_CGELS, LWORK_CGETSLS, LWORK_CGELSS,
  235. $ LWORK_CGELSY, LWORK_CGELSD,
  236. $ LRWORK_CGELSY, LRWORK_CGELSS, LRWORK_CGELSD
  237. REAL EPS, NORMA, NORMB, RCOND
  238. * ..
  239. * .. Local Arrays ..
  240. INTEGER ISEED( 4 ), ISEEDY( 4 ), IWQ( 1 )
  241. REAL RESULT( NTESTS ), RWQ( 1 )
  242. COMPLEX WQ( 1 )
  243. * ..
  244. * .. Allocatable Arrays ..
  245. COMPLEX, ALLOCATABLE :: WORK (:)
  246. REAL, ALLOCATABLE :: RWORK (:), WORK2 (:)
  247. INTEGER, ALLOCATABLE :: IWORK (:)
  248. * ..
  249. * .. External Functions ..
  250. REAL CQRT12, CQRT14, CQRT17, SASUM, SLAMCH
  251. EXTERNAL CQRT12, CQRT14, CQRT17, SASUM, SLAMCH
  252. * ..
  253. * .. External Subroutines ..
  254. EXTERNAL ALAERH, ALAHD, ALASVM, CERRLS, CGELS, CGELSD,
  255. $ CGELSS, CGELSY, CGEMM, CGETSLS, CLACPY,
  256. $ CLARNV, CQRT13, CQRT15, CQRT16, CSSCAL,
  257. $ SAXPY, XLAENV
  258. * ..
  259. * .. Intrinsic Functions ..
  260. INTRINSIC MAX, MIN, INT, REAL, SQRT
  261. * ..
  262. * .. Scalars in Common ..
  263. LOGICAL LERR, OK
  264. CHARACTER*32 SRNAMT
  265. INTEGER INFOT, IOUNIT
  266. * ..
  267. * .. Common blocks ..
  268. COMMON / INFOC / INFOT, IOUNIT, OK, LERR
  269. COMMON / SRNAMC / SRNAMT
  270. * ..
  271. * .. Data statements ..
  272. DATA ISEEDY / 1988, 1989, 1990, 1991 /
  273. * ..
  274. * .. Executable Statements ..
  275. *
  276. * Initialize constants and the random number seed.
  277. *
  278. PATH( 1: 1 ) = 'Complex precision'
  279. PATH( 2: 3 ) = 'LS'
  280. NRUN = 0
  281. NFAIL = 0
  282. NERRS = 0
  283. DO 10 I = 1, 4
  284. ISEED( I ) = ISEEDY( I )
  285. 10 CONTINUE
  286. EPS = SLAMCH( 'Epsilon' )
  287. *
  288. * Threshold for rank estimation
  289. *
  290. RCOND = SQRT( EPS ) - ( SQRT( EPS )-EPS ) / 2
  291. *
  292. * Test the error exits
  293. *
  294. CALL XLAENV( 9, SMLSIZ )
  295. IF( TSTERR )
  296. $ CALL CERRLS( PATH, NOUT )
  297. *
  298. * Print the header if NM = 0 or NN = 0 and THRESH = 0.
  299. *
  300. IF( ( NM.EQ.0 .OR. NN.EQ.0 ) .AND. THRESH.EQ.ZERO )
  301. $ CALL ALAHD( NOUT, PATH )
  302. INFOT = 0
  303. *
  304. * Compute maximal workspace needed for all routines
  305. *
  306. NMAX = 0
  307. MMAX = 0
  308. NSMAX = 0
  309. DO I = 1, NM
  310. IF ( MVAL( I ).GT.MMAX ) THEN
  311. MMAX = MVAL( I )
  312. END IF
  313. ENDDO
  314. DO I = 1, NN
  315. IF ( NVAL( I ).GT.NMAX ) THEN
  316. NMAX = NVAL( I )
  317. END IF
  318. ENDDO
  319. DO I = 1, NNS
  320. IF ( NSVAL( I ).GT.NSMAX ) THEN
  321. NSMAX = NSVAL( I )
  322. END IF
  323. ENDDO
  324. M = MMAX
  325. N = NMAX
  326. NRHS = NSMAX
  327. MNMIN = MAX( MIN( M, N ), 1 )
  328. *
  329. * Compute workspace needed for routines
  330. * CQRT14, CQRT17 (two side cases), CQRT15 and CQRT12
  331. *
  332. LWORK = MAX( 1, ( M+N )*NRHS,
  333. $ ( N+NRHS )*( M+2 ), ( M+NRHS )*( N+2 ),
  334. $ MAX( M+MNMIN, NRHS*MNMIN,2*N+M ),
  335. $ MAX( M*N+4*MNMIN+MAX(M,N), M*N+2*MNMIN+4*N ) )
  336. LRWORK = 1
  337. LIWORK = 1
  338. *
  339. * Iterate through all test cases and compute necessary workspace
  340. * sizes for ?GELS, ?GETSLS, ?GELSY, ?GELSS and ?GELSD routines.
  341. *
  342. DO IM = 1, NM
  343. M = MVAL( IM )
  344. LDA = MAX( 1, M )
  345. DO IN = 1, NN
  346. N = NVAL( IN )
  347. MNMIN = MAX(MIN( M, N ),1)
  348. LDB = MAX( 1, M, N )
  349. DO INS = 1, NNS
  350. NRHS = NSVAL( INS )
  351. DO IRANK = 1, 2
  352. DO ISCALE = 1, 3
  353. ITYPE = ( IRANK-1 )*3 + ISCALE
  354. IF( DOTYPE( ITYPE ) ) THEN
  355. IF( IRANK.EQ.1 ) THEN
  356. DO ITRAN = 1, 2
  357. IF( ITRAN.EQ.1 ) THEN
  358. TRANS = 'N'
  359. ELSE
  360. TRANS = 'C'
  361. END IF
  362. *
  363. * Compute workspace needed for CGELS
  364. CALL CGELS( TRANS, M, N, NRHS, A, LDA,
  365. $ B, LDB, WQ, -1, INFO )
  366. LWORK_CGELS = INT( WQ( 1 ) )
  367. * Compute workspace needed for CGETSLS
  368. CALL CGETSLS( TRANS, M, N, NRHS, A, LDA,
  369. $ B, LDB, WQ, -1, INFO )
  370. LWORK_CGETSLS = INT( WQ( 1 ) )
  371. ENDDO
  372. END IF
  373. * Compute workspace needed for CGELSY
  374. CALL CGELSY( M, N, NRHS, A, LDA, B, LDB,
  375. $ IWQ, RCOND, CRANK, WQ, -1, RWORK,
  376. $ INFO )
  377. LWORK_CGELSY = INT( WQ( 1 ) )
  378. LRWORK_CGELSY = 2*N
  379. * Compute workspace needed for CGELSS
  380. CALL CGELSS( M, N, NRHS, A, LDA, B, LDB, S,
  381. $ RCOND, CRANK, WQ, -1, RWORK, INFO )
  382. LWORK_CGELSS = INT( WQ( 1 ) )
  383. LRWORK_CGELSS = 5*MNMIN
  384. * Compute workspace needed for CGELSD
  385. CALL CGELSD( M, N, NRHS, A, LDA, B, LDB, S,
  386. $ RCOND, CRANK, WQ, -1, RWQ, IWQ,
  387. $ INFO )
  388. LWORK_CGELSD = INT( WQ( 1 ) )
  389. LRWORK_CGELSD = INT( RWQ ( 1 ) )
  390. * Compute LIWORK workspace needed for CGELSY and CGELSD
  391. LIWORK = MAX( LIWORK, N, IWQ ( 1 ) )
  392. * Compute LRWORK workspace needed for CGELSY, CGELSS and CGELSD
  393. LRWORK = MAX( LRWORK, LRWORK_CGELSY,
  394. $ LRWORK_CGELSS, LRWORK_CGELSD )
  395. * Compute LWORK workspace needed for all functions
  396. LWORK = MAX( LWORK, LWORK_CGELS, LWORK_CGETSLS,
  397. $ LWORK_CGELSY, LWORK_CGELSS,
  398. $ LWORK_CGELSD )
  399. END IF
  400. ENDDO
  401. ENDDO
  402. ENDDO
  403. ENDDO
  404. ENDDO
  405. *
  406. LWLSY = LWORK
  407. *
  408. ALLOCATE( WORK( LWORK ) )
  409. ALLOCATE( IWORK( LIWORK ) )
  410. ALLOCATE( RWORK( LRWORK ) )
  411. ALLOCATE( WORK2( 2 * LWORK ) )
  412. *
  413. DO 140 IM = 1, NM
  414. M = MVAL( IM )
  415. LDA = MAX( 1, M )
  416. *
  417. DO 130 IN = 1, NN
  418. N = NVAL( IN )
  419. MNMIN = MAX(MIN( M, N ),1)
  420. LDB = MAX( 1, M, N )
  421. MB = (MNMIN+1)
  422. *
  423. DO 120 INS = 1, NNS
  424. NRHS = NSVAL( INS )
  425. *
  426. DO 110 IRANK = 1, 2
  427. DO 100 ISCALE = 1, 3
  428. ITYPE = ( IRANK-1 )*3 + ISCALE
  429. IF( .NOT.DOTYPE( ITYPE ) )
  430. $ GO TO 100
  431. *
  432. IF( IRANK.EQ.1 ) THEN
  433. *
  434. * Test CGELS
  435. *
  436. * Generate a matrix of scaling type ISCALE
  437. *
  438. CALL CQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
  439. $ ISEED )
  440. DO 40 INB = 1, NNB
  441. NB = NBVAL( INB )
  442. CALL XLAENV( 1, NB )
  443. CALL XLAENV( 3, NXVAL( INB ) )
  444. *
  445. DO 30 ITRAN = 1, 2
  446. IF( ITRAN.EQ.1 ) THEN
  447. TRANS = 'N'
  448. NROWS = M
  449. NCOLS = N
  450. ELSE
  451. TRANS = 'C'
  452. NROWS = N
  453. NCOLS = M
  454. END IF
  455. LDWORK = MAX( 1, NCOLS )
  456. *
  457. * Set up a consistent rhs
  458. *
  459. IF( NCOLS.GT.0 ) THEN
  460. CALL CLARNV( 2, ISEED, NCOLS*NRHS,
  461. $ WORK )
  462. CALL CSSCAL( NCOLS*NRHS,
  463. $ ONE / REAL( NCOLS ), WORK,
  464. $ 1 )
  465. END IF
  466. CALL CGEMM( TRANS, 'No transpose', NROWS,
  467. $ NRHS, NCOLS, CONE, COPYA, LDA,
  468. $ WORK, LDWORK, CZERO, B, LDB )
  469. CALL CLACPY( 'Full', NROWS, NRHS, B, LDB,
  470. $ COPYB, LDB )
  471. *
  472. * Solve LS or overdetermined system
  473. *
  474. IF( M.GT.0 .AND. N.GT.0 ) THEN
  475. CALL CLACPY( 'Full', M, N, COPYA, LDA,
  476. $ A, LDA )
  477. CALL CLACPY( 'Full', NROWS, NRHS,
  478. $ COPYB, LDB, B, LDB )
  479. END IF
  480. SRNAMT = 'CGELS '
  481. CALL CGELS( TRANS, M, N, NRHS, A, LDA, B,
  482. $ LDB, WORK, LWORK, INFO )
  483. *
  484. IF( INFO.NE.0 )
  485. $ CALL ALAERH( PATH, 'CGELS ', INFO, 0,
  486. $ TRANS, M, N, NRHS, -1, NB,
  487. $ ITYPE, NFAIL, NERRS,
  488. $ NOUT )
  489. *
  490. * Check correctness of results
  491. *
  492. LDWORK = MAX( 1, NROWS )
  493. IF( NROWS.GT.0 .AND. NRHS.GT.0 )
  494. $ CALL CLACPY( 'Full', NROWS, NRHS,
  495. $ COPYB, LDB, C, LDB )
  496. CALL CQRT16( TRANS, M, N, NRHS, COPYA,
  497. $ LDA, B, LDB, C, LDB, RWORK,
  498. $ RESULT( 1 ) )
  499. *
  500. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
  501. $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
  502. *
  503. * Solving LS system
  504. *
  505. RESULT( 2 ) = CQRT17( TRANS, 1, M, N,
  506. $ NRHS, COPYA, LDA, B, LDB,
  507. $ COPYB, LDB, C, WORK,
  508. $ LWORK )
  509. ELSE
  510. *
  511. * Solving overdetermined system
  512. *
  513. RESULT( 2 ) = CQRT14( TRANS, M, N,
  514. $ NRHS, COPYA, LDA, B, LDB,
  515. $ WORK, LWORK )
  516. END IF
  517. *
  518. * Print information about the tests that
  519. * did not pass the threshold.
  520. *
  521. DO 20 K = 1, 2
  522. IF( RESULT( K ).GE.THRESH ) THEN
  523. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  524. $ CALL ALAHD( NOUT, PATH )
  525. WRITE( NOUT, FMT = 9999 )TRANS, M,
  526. $ N, NRHS, NB, ITYPE, K,
  527. $ RESULT( K )
  528. NFAIL = NFAIL + 1
  529. END IF
  530. 20 CONTINUE
  531. NRUN = NRUN + 2
  532. 30 CONTINUE
  533. 40 CONTINUE
  534. *
  535. *
  536. * Test CGETSLS
  537. *
  538. * Generate a matrix of scaling type ISCALE
  539. *
  540. CALL CQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
  541. $ ISEED )
  542. DO 65 INB = 1, NNB
  543. MB = NBVAL( INB )
  544. CALL XLAENV( 1, MB )
  545. DO 62 IMB = 1, NNB
  546. NB = NBVAL( IMB )
  547. CALL XLAENV( 2, NB )
  548. *
  549. DO 60 ITRAN = 1, 2
  550. IF( ITRAN.EQ.1 ) THEN
  551. TRANS = 'N'
  552. NROWS = M
  553. NCOLS = N
  554. ELSE
  555. TRANS = 'C'
  556. NROWS = N
  557. NCOLS = M
  558. END IF
  559. LDWORK = MAX( 1, NCOLS )
  560. *
  561. * Set up a consistent rhs
  562. *
  563. IF( NCOLS.GT.0 ) THEN
  564. CALL CLARNV( 2, ISEED, NCOLS*NRHS,
  565. $ WORK )
  566. CALL CSCAL( NCOLS*NRHS,
  567. $ ONE / REAL( NCOLS ), WORK,
  568. $ 1 )
  569. END IF
  570. CALL CGEMM( TRANS, 'No transpose', NROWS,
  571. $ NRHS, NCOLS, CONE, COPYA, LDA,
  572. $ WORK, LDWORK, CZERO, B, LDB )
  573. CALL CLACPY( 'Full', NROWS, NRHS, B, LDB,
  574. $ COPYB, LDB )
  575. *
  576. * Solve LS or overdetermined system
  577. *
  578. IF( M.GT.0 .AND. N.GT.0 ) THEN
  579. CALL CLACPY( 'Full', M, N, COPYA, LDA,
  580. $ A, LDA )
  581. CALL CLACPY( 'Full', NROWS, NRHS,
  582. $ COPYB, LDB, B, LDB )
  583. END IF
  584. SRNAMT = 'CGETSLS '
  585. CALL CGETSLS( TRANS, M, N, NRHS, A,
  586. $ LDA, B, LDB, WORK, LWORK, INFO )
  587. IF( INFO.NE.0 )
  588. $ CALL ALAERH( PATH, 'CGETSLS ', INFO, 0,
  589. $ TRANS, M, N, NRHS, -1, NB,
  590. $ ITYPE, NFAIL, NERRS,
  591. $ NOUT )
  592. *
  593. * Check correctness of results
  594. *
  595. LDWORK = MAX( 1, NROWS )
  596. IF( NROWS.GT.0 .AND. NRHS.GT.0 )
  597. $ CALL CLACPY( 'Full', NROWS, NRHS,
  598. $ COPYB, LDB, C, LDB )
  599. CALL CQRT16( TRANS, M, N, NRHS, COPYA,
  600. $ LDA, B, LDB, C, LDB, WORK2,
  601. $ RESULT( 15 ) )
  602. *
  603. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
  604. $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
  605. *
  606. * Solving LS system
  607. *
  608. RESULT( 16 ) = CQRT17( TRANS, 1, M, N,
  609. $ NRHS, COPYA, LDA, B, LDB,
  610. $ COPYB, LDB, C, WORK,
  611. $ LWORK )
  612. ELSE
  613. *
  614. * Solving overdetermined system
  615. *
  616. RESULT( 16 ) = CQRT14( TRANS, M, N,
  617. $ NRHS, COPYA, LDA, B, LDB,
  618. $ WORK, LWORK )
  619. END IF
  620. *
  621. * Print information about the tests that
  622. * did not pass the threshold.
  623. *
  624. DO 50 K = 15, 16
  625. IF( RESULT( K ).GE.THRESH ) THEN
  626. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  627. $ CALL ALAHD( NOUT, PATH )
  628. WRITE( NOUT, FMT = 9997 )TRANS, M,
  629. $ N, NRHS, MB, NB, ITYPE, K,
  630. $ RESULT( K )
  631. NFAIL = NFAIL + 1
  632. END IF
  633. 50 CONTINUE
  634. NRUN = NRUN + 2
  635. 60 CONTINUE
  636. 62 CONTINUE
  637. 65 CONTINUE
  638. END IF
  639. *
  640. * Generate a matrix of scaling type ISCALE and rank
  641. * type IRANK.
  642. *
  643. CALL CQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA,
  644. $ COPYB, LDB, COPYS, RANK, NORMA, NORMB,
  645. $ ISEED, WORK, LWORK )
  646. *
  647. * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
  648. *
  649. LDWORK = MAX( 1, M )
  650. *
  651. * Loop for testing different block sizes.
  652. *
  653. DO 90 INB = 1, NNB
  654. NB = NBVAL( INB )
  655. CALL XLAENV( 1, NB )
  656. CALL XLAENV( 3, NXVAL( INB ) )
  657. *
  658. * Test CGELSY
  659. *
  660. * CGELSY: Compute the minimum-norm solution
  661. * X to min( norm( A * X - B ) )
  662. * using the rank-revealing orthogonal
  663. * factorization.
  664. *
  665. CALL CLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  666. CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  667. $ LDB )
  668. *
  669. * Initialize vector IWORK.
  670. *
  671. DO 70 J = 1, N
  672. IWORK( J ) = 0
  673. 70 CONTINUE
  674. *
  675. SRNAMT = 'CGELSY'
  676. CALL CGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
  677. $ RCOND, CRANK, WORK, LWLSY, RWORK,
  678. $ INFO )
  679. IF( INFO.NE.0 )
  680. $ CALL ALAERH( PATH, 'CGELSY', INFO, 0, ' ', M,
  681. $ N, NRHS, -1, NB, ITYPE, NFAIL,
  682. $ NERRS, NOUT )
  683. *
  684. * workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS)
  685. *
  686. * Test 3: Compute relative error in svd
  687. * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
  688. *
  689. RESULT( 3 ) = CQRT12( CRANK, CRANK, A, LDA,
  690. $ COPYS, WORK, LWORK, RWORK )
  691. *
  692. * Test 4: Compute error in solution
  693. * workspace: M*NRHS + M
  694. *
  695. CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  696. $ LDWORK )
  697. CALL CQRT16( 'No transpose', M, N, NRHS, COPYA,
  698. $ LDA, B, LDB, WORK, LDWORK, RWORK,
  699. $ RESULT( 4 ) )
  700. *
  701. * Test 5: Check norm of r'*A
  702. * workspace: NRHS*(M+N)
  703. *
  704. RESULT( 5 ) = ZERO
  705. IF( M.GT.CRANK )
  706. $ RESULT( 5 ) = CQRT17( 'No transpose', 1, M,
  707. $ N, NRHS, COPYA, LDA, B, LDB,
  708. $ COPYB, LDB, C, WORK, LWORK )
  709. *
  710. * Test 6: Check if x is in the rowspace of A
  711. * workspace: (M+NRHS)*(N+2)
  712. *
  713. RESULT( 6 ) = ZERO
  714. *
  715. IF( N.GT.CRANK )
  716. $ RESULT( 6 ) = CQRT14( 'No transpose', M, N,
  717. $ NRHS, COPYA, LDA, B, LDB,
  718. $ WORK, LWORK )
  719. *
  720. * Test CGELSS
  721. *
  722. * CGELSS: Compute the minimum-norm solution
  723. * X to min( norm( A * X - B ) )
  724. * using the SVD.
  725. *
  726. CALL CLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  727. CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  728. $ LDB )
  729. SRNAMT = 'CGELSS'
  730. CALL CGELSS( M, N, NRHS, A, LDA, B, LDB, S,
  731. $ RCOND, CRANK, WORK, LWORK, RWORK,
  732. $ INFO )
  733. *
  734. IF( INFO.NE.0 )
  735. $ CALL ALAERH( PATH, 'CGELSS', INFO, 0, ' ', M,
  736. $ N, NRHS, -1, NB, ITYPE, NFAIL,
  737. $ NERRS, NOUT )
  738. *
  739. * workspace used: 3*min(m,n) +
  740. * max(2*min(m,n),nrhs,max(m,n))
  741. *
  742. * Test 7: Compute relative error in svd
  743. *
  744. IF( RANK.GT.0 ) THEN
  745. CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
  746. RESULT( 7 ) = SASUM( MNMIN, S, 1 ) /
  747. $ SASUM( MNMIN, COPYS, 1 ) /
  748. $ ( EPS*REAL( MNMIN ) )
  749. ELSE
  750. RESULT( 7 ) = ZERO
  751. END IF
  752. *
  753. * Test 8: Compute error in solution
  754. *
  755. CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  756. $ LDWORK )
  757. CALL CQRT16( 'No transpose', M, N, NRHS, COPYA,
  758. $ LDA, B, LDB, WORK, LDWORK, RWORK,
  759. $ RESULT( 8 ) )
  760. *
  761. * Test 9: Check norm of r'*A
  762. *
  763. RESULT( 9 ) = ZERO
  764. IF( M.GT.CRANK )
  765. $ RESULT( 9 ) = CQRT17( 'No transpose', 1, M,
  766. $ N, NRHS, COPYA, LDA, B, LDB,
  767. $ COPYB, LDB, C, WORK, LWORK )
  768. *
  769. * Test 10: Check if x is in the rowspace of A
  770. *
  771. RESULT( 10 ) = ZERO
  772. IF( N.GT.CRANK )
  773. $ RESULT( 10 ) = CQRT14( 'No transpose', M, N,
  774. $ NRHS, COPYA, LDA, B, LDB,
  775. $ WORK, LWORK )
  776. *
  777. * Test CGELSD
  778. *
  779. * CGELSD: Compute the minimum-norm solution X
  780. * to min( norm( A * X - B ) ) using a
  781. * divide and conquer SVD.
  782. *
  783. CALL XLAENV( 9, 25 )
  784. *
  785. CALL CLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  786. CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  787. $ LDB )
  788. *
  789. SRNAMT = 'CGELSD'
  790. CALL CGELSD( M, N, NRHS, A, LDA, B, LDB, S,
  791. $ RCOND, CRANK, WORK, LWORK, RWORK,
  792. $ IWORK, INFO )
  793. IF( INFO.NE.0 )
  794. $ CALL ALAERH( PATH, 'CGELSD', INFO, 0, ' ', M,
  795. $ N, NRHS, -1, NB, ITYPE, NFAIL,
  796. $ NERRS, NOUT )
  797. *
  798. * Test 11: Compute relative error in svd
  799. *
  800. IF( RANK.GT.0 ) THEN
  801. CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
  802. RESULT( 11 ) = SASUM( MNMIN, S, 1 ) /
  803. $ SASUM( MNMIN, COPYS, 1 ) /
  804. $ ( EPS*REAL( MNMIN ) )
  805. ELSE
  806. RESULT( 11 ) = ZERO
  807. END IF
  808. *
  809. * Test 12: Compute error in solution
  810. *
  811. CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  812. $ LDWORK )
  813. CALL CQRT16( 'No transpose', M, N, NRHS, COPYA,
  814. $ LDA, B, LDB, WORK, LDWORK, RWORK,
  815. $ RESULT( 12 ) )
  816. *
  817. * Test 13: Check norm of r'*A
  818. *
  819. RESULT( 13 ) = ZERO
  820. IF( M.GT.CRANK )
  821. $ RESULT( 13 ) = CQRT17( 'No transpose', 1, M,
  822. $ N, NRHS, COPYA, LDA, B, LDB,
  823. $ COPYB, LDB, C, WORK, LWORK )
  824. *
  825. * Test 14: Check if x is in the rowspace of A
  826. *
  827. RESULT( 14 ) = ZERO
  828. IF( N.GT.CRANK )
  829. $ RESULT( 14 ) = CQRT14( 'No transpose', M, N,
  830. $ NRHS, COPYA, LDA, B, LDB,
  831. $ WORK, LWORK )
  832. *
  833. * Print information about the tests that did not
  834. * pass the threshold.
  835. *
  836. DO 80 K = 3, 14
  837. IF( RESULT( K ).GE.THRESH ) THEN
  838. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  839. $ CALL ALAHD( NOUT, PATH )
  840. WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
  841. $ ITYPE, K, RESULT( K )
  842. NFAIL = NFAIL + 1
  843. END IF
  844. 80 CONTINUE
  845. NRUN = NRUN + 12
  846. *
  847. 90 CONTINUE
  848. 100 CONTINUE
  849. 110 CONTINUE
  850. 120 CONTINUE
  851. 130 CONTINUE
  852. 140 CONTINUE
  853. *
  854. * Print a summary of the results.
  855. *
  856. CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
  857. *
  858. 9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4,
  859. $ ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 )
  860. 9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4,
  861. $ ', type', I2, ', test(', I2, ')=', G12.5 )
  862. 9997 FORMAT( ' TRANS=''', A1,' M=', I5, ', N=', I5, ', NRHS=', I4,
  863. $ ', MB=', I4,', NB=', I4,', type', I2,
  864. $ ', test(', I2, ')=', G12.5 )
  865. *
  866. DEALLOCATE( WORK )
  867. DEALLOCATE( RWORK )
  868. DEALLOCATE( IWORK )
  869. RETURN
  870. *
  871. * End of CDRVLS
  872. *
  873. END