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.

zdrvls.f 34 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872
  1. *> \brief \b ZDRVLS
  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 ZDRVLS( 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. * DOUBLE PRECISION THRESH
  19. * ..
  20. * .. Array Arguments ..
  21. * LOGICAL DOTYPE( * )
  22. * INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
  23. * $ NVAL( * ), NXVAL( * )
  24. * DOUBLE PRECISION COPYS( * ), S( * )
  25. * COMPLEX*16 A( * ), B( * ), C( * ), COPYA( * ), COPYB( * )
  26. * ..
  27. *
  28. *
  29. *> \par Purpose:
  30. * =============
  31. *>
  32. *> \verbatim
  33. *>
  34. *> ZDRVLS tests the least squares driver routines ZGELS, ZGETSLS, ZGELSS, ZGELSY
  35. *> and ZGELSD.
  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 DOUBLE PRECISION
  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*16 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*16 array, dimension (MMAX*NMAX)
  141. *> \endverbatim
  142. *>
  143. *> \param[out] B
  144. *> \verbatim
  145. *> B is COMPLEX*16 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*16 array, dimension (MMAX*NSMAX)
  153. *> \endverbatim
  154. *>
  155. *> \param[out] C
  156. *> \verbatim
  157. *> C is COMPLEX*16 array, dimension (MMAX*NSMAX)
  158. *> \endverbatim
  159. *>
  160. *> \param[out] S
  161. *> \verbatim
  162. *> S is DOUBLE PRECISION array, dimension
  163. *> (min(MMAX,NMAX))
  164. *> \endverbatim
  165. *>
  166. *> \param[out] COPYS
  167. *> \verbatim
  168. *> COPYS is DOUBLE PRECISION 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 complex16_lin
  189. *
  190. * =====================================================================
  191. SUBROUTINE ZDRVLS( 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. DOUBLE PRECISION THRESH
  204. * ..
  205. * .. Array Arguments ..
  206. LOGICAL DOTYPE( * )
  207. INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
  208. $ NVAL( * ), NXVAL( * )
  209. DOUBLE PRECISION COPYS( * ), S( * )
  210. COMPLEX*16 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. DOUBLE PRECISION ONE, ZERO
  221. PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  222. COMPLEX*16 CONE, CZERO
  223. PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
  224. $ CZERO = ( 0.0D+0, 0.0D+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_ZGELS, LWORK_ZGETSLS, LWORK_ZGELSS,
  235. $ LWORK_ZGELSY, LWORK_ZGELSD,
  236. $ LRWORK_ZGELSY, LRWORK_ZGELSS, LRWORK_ZGELSD
  237. DOUBLE PRECISION EPS, NORMA, NORMB, RCOND
  238. * ..
  239. * .. Local Arrays ..
  240. INTEGER ISEED( 4 ), ISEEDY( 4 ), IWQ
  241. DOUBLE PRECISION RESULT( NTESTS ), RWQ
  242. COMPLEX*16 WQ
  243. * ..
  244. * .. Allocatable Arrays ..
  245. COMPLEX*16, ALLOCATABLE :: WORK (:)
  246. DOUBLE PRECISION, ALLOCATABLE :: RWORK (:)
  247. INTEGER, ALLOCATABLE :: IWORK (:)
  248. * ..
  249. * .. External Functions ..
  250. DOUBLE PRECISION DASUM, DLAMCH, ZQRT12, ZQRT14, ZQRT17
  251. EXTERNAL DASUM, DLAMCH, ZQRT12, ZQRT14, ZQRT17
  252. * ..
  253. * .. External Subroutines ..
  254. EXTERNAL ALAERH, ALAHD, ALASVM, DAXPY, DLASRT, XLAENV,
  255. $ ZDSCAL, ZERRLS, ZGELS, ZGELSD, ZGELSS,
  256. $ ZGELSY, ZGEMM, ZLACPY, ZLARNV, ZQRT13, ZQRT15,
  257. $ ZQRT16, ZGETSLS
  258. * ..
  259. * .. Intrinsic Functions ..
  260. INTRINSIC DBLE, MAX, MIN, INT, 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 ) = 'Zomplex 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 = DLAMCH( '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 ZERRLS( 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. * ZQRT14, ZQRT17 (two side cases), ZQRT15 and ZQRT12
  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 ZGELS
  364. CALL ZGELS( TRANS, M, N, NRHS, A, LDA,
  365. $ B, LDB, WQ, -1, INFO )
  366. LWORK_ZGELS = INT ( WQ )
  367. * Compute workspace needed for ZGETSLS
  368. CALL ZGETSLS( TRANS, M, N, NRHS, A, LDA,
  369. $ B, LDB, WQ, -1, INFO )
  370. LWORK_ZGETSLS = INT( WQ )
  371. ENDDO
  372. END IF
  373. * Compute workspace needed for ZGELSY
  374. CALL ZGELSY( M, N, NRHS, A, LDA, B, LDB, IWQ,
  375. $ RCOND, CRANK, WQ, -1, RWORK, INFO )
  376. LWORK_ZGELSY = INT( WQ )
  377. LRWORK_ZGELSY = 2*N
  378. * Compute workspace needed for ZGELSS
  379. CALL ZGELSS( M, N, NRHS, A, LDA, B, LDB, S,
  380. $ RCOND, CRANK, WQ, -1 , RWORK,
  381. $ INFO )
  382. LWORK_ZGELSS = INT( WQ )
  383. LRWORK_ZGELSS = 5*MNMIN
  384. * Compute workspace needed for ZGELSD
  385. CALL ZGELSD( M, N, NRHS, A, LDA, B, LDB, S,
  386. $ RCOND, CRANK, WQ, -1, RWQ, IWQ,
  387. $ INFO )
  388. LWORK_ZGELSD = INT( WQ )
  389. LRWORK_ZGELSD = INT( RWQ )
  390. * Compute LIWORK workspace needed for ZGELSY and ZGELSD
  391. LIWORK = MAX( LIWORK, N, IWQ )
  392. * Compute LRWORK workspace needed for ZGELSY, ZGELSS and ZGELSD
  393. LRWORK = MAX( LRWORK, LRWORK_ZGELSY,
  394. $ LRWORK_ZGELSS, LRWORK_ZGELSD )
  395. * Compute LWORK workspace needed for all functions
  396. LWORK = MAX( LWORK, LWORK_ZGELS, LWORK_ZGETSLS,
  397. $ LWORK_ZGELSY, LWORK_ZGELSS,
  398. $ LWORK_ZGELSD )
  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. *
  412. DO 140 IM = 1, NM
  413. M = MVAL( IM )
  414. LDA = MAX( 1, M )
  415. *
  416. DO 130 IN = 1, NN
  417. N = NVAL( IN )
  418. MNMIN = MAX(MIN( M, N ),1)
  419. LDB = MAX( 1, M, N )
  420. MB = (MNMIN+1)
  421. *
  422. DO 120 INS = 1, NNS
  423. NRHS = NSVAL( INS )
  424. *
  425. DO 110 IRANK = 1, 2
  426. DO 100 ISCALE = 1, 3
  427. ITYPE = ( IRANK-1 )*3 + ISCALE
  428. IF( .NOT.DOTYPE( ITYPE ) )
  429. $ GO TO 100
  430. *
  431. IF( IRANK.EQ.1 ) THEN
  432. *
  433. * Test ZGELS
  434. *
  435. * Generate a matrix of scaling type ISCALE
  436. *
  437. CALL ZQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
  438. $ ISEED )
  439. DO 40 INB = 1, NNB
  440. NB = NBVAL( INB )
  441. CALL XLAENV( 1, NB )
  442. CALL XLAENV( 3, NXVAL( INB ) )
  443. *
  444. DO 30 ITRAN = 1, 2
  445. IF( ITRAN.EQ.1 ) THEN
  446. TRANS = 'N'
  447. NROWS = M
  448. NCOLS = N
  449. ELSE
  450. TRANS = 'C'
  451. NROWS = N
  452. NCOLS = M
  453. END IF
  454. LDWORK = MAX( 1, NCOLS )
  455. *
  456. * Set up a consistent rhs
  457. *
  458. IF( NCOLS.GT.0 ) THEN
  459. CALL ZLARNV( 2, ISEED, NCOLS*NRHS,
  460. $ WORK )
  461. CALL ZDSCAL( NCOLS*NRHS,
  462. $ ONE / DBLE( NCOLS ), WORK,
  463. $ 1 )
  464. END IF
  465. CALL ZGEMM( TRANS, 'No transpose', NROWS,
  466. $ NRHS, NCOLS, CONE, COPYA, LDA,
  467. $ WORK, LDWORK, CZERO, B, LDB )
  468. CALL ZLACPY( 'Full', NROWS, NRHS, B, LDB,
  469. $ COPYB, LDB )
  470. *
  471. * Solve LS or overdetermined system
  472. *
  473. IF( M.GT.0 .AND. N.GT.0 ) THEN
  474. CALL ZLACPY( 'Full', M, N, COPYA, LDA,
  475. $ A, LDA )
  476. CALL ZLACPY( 'Full', NROWS, NRHS,
  477. $ COPYB, LDB, B, LDB )
  478. END IF
  479. SRNAMT = 'ZGELS '
  480. CALL ZGELS( TRANS, M, N, NRHS, A, LDA, B,
  481. $ LDB, WORK, LWORK, INFO )
  482. *
  483. IF( INFO.NE.0 )
  484. $ CALL ALAERH( PATH, 'ZGELS ', INFO, 0,
  485. $ TRANS, M, N, NRHS, -1, NB,
  486. $ ITYPE, NFAIL, NERRS,
  487. $ NOUT )
  488. *
  489. * Check correctness of results
  490. *
  491. LDWORK = MAX( 1, NROWS )
  492. IF( NROWS.GT.0 .AND. NRHS.GT.0 )
  493. $ CALL ZLACPY( 'Full', NROWS, NRHS,
  494. $ COPYB, LDB, C, LDB )
  495. CALL ZQRT16( TRANS, M, N, NRHS, COPYA,
  496. $ LDA, B, LDB, C, LDB, RWORK,
  497. $ RESULT( 1 ) )
  498. *
  499. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
  500. $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
  501. *
  502. * Solving LS system
  503. *
  504. RESULT( 2 ) = ZQRT17( TRANS, 1, M, N,
  505. $ NRHS, COPYA, LDA, B, LDB,
  506. $ COPYB, LDB, C, WORK,
  507. $ LWORK )
  508. ELSE
  509. *
  510. * Solving overdetermined system
  511. *
  512. RESULT( 2 ) = ZQRT14( TRANS, M, N,
  513. $ NRHS, COPYA, LDA, B, LDB,
  514. $ WORK, LWORK )
  515. END IF
  516. *
  517. * Print information about the tests that
  518. * did not pass the threshold.
  519. *
  520. DO 20 K = 1, 2
  521. IF( RESULT( K ).GE.THRESH ) THEN
  522. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  523. $ CALL ALAHD( NOUT, PATH )
  524. WRITE( NOUT, FMT = 9999 )TRANS, M,
  525. $ N, NRHS, NB, ITYPE, K,
  526. $ RESULT( K )
  527. NFAIL = NFAIL + 1
  528. END IF
  529. 20 CONTINUE
  530. NRUN = NRUN + 2
  531. 30 CONTINUE
  532. 40 CONTINUE
  533. *
  534. *
  535. * Test ZGETSLS
  536. *
  537. * Generate a matrix of scaling type ISCALE
  538. *
  539. CALL ZQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
  540. $ ISEED )
  541. DO 65 INB = 1, NNB
  542. MB = NBVAL( INB )
  543. CALL XLAENV( 1, MB )
  544. DO 62 IMB = 1, NNB
  545. NB = NBVAL( IMB )
  546. CALL XLAENV( 2, NB )
  547. *
  548. DO 60 ITRAN = 1, 2
  549. IF( ITRAN.EQ.1 ) THEN
  550. TRANS = 'N'
  551. NROWS = M
  552. NCOLS = N
  553. ELSE
  554. TRANS = 'C'
  555. NROWS = N
  556. NCOLS = M
  557. END IF
  558. LDWORK = MAX( 1, NCOLS )
  559. *
  560. * Set up a consistent rhs
  561. *
  562. IF( NCOLS.GT.0 ) THEN
  563. CALL ZLARNV( 2, ISEED, NCOLS*NRHS,
  564. $ WORK )
  565. CALL ZSCAL( NCOLS*NRHS,
  566. $ ONE / DBLE( NCOLS ), WORK,
  567. $ 1 )
  568. END IF
  569. CALL ZGEMM( TRANS, 'No transpose', NROWS,
  570. $ NRHS, NCOLS, CONE, COPYA, LDA,
  571. $ WORK, LDWORK, CZERO, B, LDB )
  572. CALL ZLACPY( 'Full', NROWS, NRHS, B, LDB,
  573. $ COPYB, LDB )
  574. *
  575. * Solve LS or overdetermined system
  576. *
  577. IF( M.GT.0 .AND. N.GT.0 ) THEN
  578. CALL ZLACPY( 'Full', M, N, COPYA, LDA,
  579. $ A, LDA )
  580. CALL ZLACPY( 'Full', NROWS, NRHS,
  581. $ COPYB, LDB, B, LDB )
  582. END IF
  583. SRNAMT = 'ZGETSLS '
  584. CALL ZGETSLS( TRANS, M, N, NRHS, A,
  585. $ LDA, B, LDB, WORK, LWORK, INFO )
  586. IF( INFO.NE.0 )
  587. $ CALL ALAERH( PATH, 'ZGETSLS ', INFO, 0,
  588. $ TRANS, M, N, NRHS, -1, NB,
  589. $ ITYPE, NFAIL, NERRS,
  590. $ NOUT )
  591. *
  592. * Check correctness of results
  593. *
  594. LDWORK = MAX( 1, NROWS )
  595. IF( NROWS.GT.0 .AND. NRHS.GT.0 )
  596. $ CALL ZLACPY( 'Full', NROWS, NRHS,
  597. $ COPYB, LDB, C, LDB )
  598. CALL ZQRT16( TRANS, M, N, NRHS, COPYA,
  599. $ LDA, B, LDB, C, LDB, WORK,
  600. $ RESULT( 15 ) )
  601. *
  602. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
  603. $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
  604. *
  605. * Solving LS system
  606. *
  607. RESULT( 16 ) = ZQRT17( TRANS, 1, M, N,
  608. $ NRHS, COPYA, LDA, B, LDB,
  609. $ COPYB, LDB, C, WORK,
  610. $ LWORK )
  611. ELSE
  612. *
  613. * Solving overdetermined system
  614. *
  615. RESULT( 16 ) = ZQRT14( TRANS, M, N,
  616. $ NRHS, COPYA, LDA, B, LDB,
  617. $ WORK, LWORK )
  618. END IF
  619. *
  620. * Print information about the tests that
  621. * did not pass the threshold.
  622. *
  623. DO 50 K = 15, 16
  624. IF( RESULT( K ).GE.THRESH ) THEN
  625. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  626. $ CALL ALAHD( NOUT, PATH )
  627. WRITE( NOUT, FMT = 9997 )TRANS, M,
  628. $ N, NRHS, MB, NB, ITYPE, K,
  629. $ RESULT( K )
  630. NFAIL = NFAIL + 1
  631. END IF
  632. 50 CONTINUE
  633. NRUN = NRUN + 2
  634. 60 CONTINUE
  635. 62 CONTINUE
  636. 65 CONTINUE
  637. END IF
  638. *
  639. * Generate a matrix of scaling type ISCALE and rank
  640. * type IRANK.
  641. *
  642. CALL ZQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA,
  643. $ COPYB, LDB, COPYS, RANK, NORMA, NORMB,
  644. $ ISEED, WORK, LWORK )
  645. *
  646. * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
  647. *
  648. LDWORK = MAX( 1, M )
  649. *
  650. * Loop for testing different block sizes.
  651. *
  652. DO 90 INB = 1, NNB
  653. NB = NBVAL( INB )
  654. CALL XLAENV( 1, NB )
  655. CALL XLAENV( 3, NXVAL( INB ) )
  656. *
  657. * Test ZGELSY
  658. *
  659. * ZGELSY: Compute the minimum-norm solution
  660. * X to min( norm( A * X - B ) )
  661. * using the rank-revealing orthogonal
  662. * factorization.
  663. *
  664. CALL ZLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  665. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  666. $ LDB )
  667. *
  668. * Initialize vector IWORK.
  669. *
  670. DO 70 J = 1, N
  671. IWORK( J ) = 0
  672. 70 CONTINUE
  673. *
  674. SRNAMT = 'ZGELSY'
  675. CALL ZGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
  676. $ RCOND, CRANK, WORK, LWLSY, RWORK,
  677. $ INFO )
  678. IF( INFO.NE.0 )
  679. $ CALL ALAERH( PATH, 'ZGELSY', INFO, 0, ' ', M,
  680. $ N, NRHS, -1, NB, ITYPE, NFAIL,
  681. $ NERRS, NOUT )
  682. *
  683. * workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS)
  684. *
  685. * Test 3: Compute relative error in svd
  686. * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
  687. *
  688. RESULT( 3 ) = ZQRT12( CRANK, CRANK, A, LDA,
  689. $ COPYS, WORK, LWORK, RWORK )
  690. *
  691. * Test 4: Compute error in solution
  692. * workspace: M*NRHS + M
  693. *
  694. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  695. $ LDWORK )
  696. CALL ZQRT16( 'No transpose', M, N, NRHS, COPYA,
  697. $ LDA, B, LDB, WORK, LDWORK, RWORK,
  698. $ RESULT( 4 ) )
  699. *
  700. * Test 5: Check norm of r'*A
  701. * workspace: NRHS*(M+N)
  702. *
  703. RESULT( 5 ) = ZERO
  704. IF( M.GT.CRANK )
  705. $ RESULT( 5 ) = ZQRT17( 'No transpose', 1, M,
  706. $ N, NRHS, COPYA, LDA, B, LDB,
  707. $ COPYB, LDB, C, WORK, LWORK )
  708. *
  709. * Test 6: Check if x is in the rowspace of A
  710. * workspace: (M+NRHS)*(N+2)
  711. *
  712. RESULT( 6 ) = ZERO
  713. *
  714. IF( N.GT.CRANK )
  715. $ RESULT( 6 ) = ZQRT14( 'No transpose', M, N,
  716. $ NRHS, COPYA, LDA, B, LDB,
  717. $ WORK, LWORK )
  718. *
  719. * Test ZGELSS
  720. *
  721. * ZGELSS: Compute the minimum-norm solution
  722. * X to min( norm( A * X - B ) )
  723. * using the SVD.
  724. *
  725. CALL ZLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  726. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  727. $ LDB )
  728. SRNAMT = 'ZGELSS'
  729. CALL ZGELSS( M, N, NRHS, A, LDA, B, LDB, S,
  730. $ RCOND, CRANK, WORK, LWORK, RWORK,
  731. $ INFO )
  732. *
  733. IF( INFO.NE.0 )
  734. $ CALL ALAERH( PATH, 'ZGELSS', INFO, 0, ' ', M,
  735. $ N, NRHS, -1, NB, ITYPE, NFAIL,
  736. $ NERRS, NOUT )
  737. *
  738. * workspace used: 3*min(m,n) +
  739. * max(2*min(m,n),nrhs,max(m,n))
  740. *
  741. * Test 7: Compute relative error in svd
  742. *
  743. IF( RANK.GT.0 ) THEN
  744. CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
  745. RESULT( 7 ) = DASUM( MNMIN, S, 1 ) /
  746. $ DASUM( MNMIN, COPYS, 1 ) /
  747. $ ( EPS*DBLE( MNMIN ) )
  748. ELSE
  749. RESULT( 7 ) = ZERO
  750. END IF
  751. *
  752. * Test 8: Compute error in solution
  753. *
  754. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  755. $ LDWORK )
  756. CALL ZQRT16( 'No transpose', M, N, NRHS, COPYA,
  757. $ LDA, B, LDB, WORK, LDWORK, RWORK,
  758. $ RESULT( 8 ) )
  759. *
  760. * Test 9: Check norm of r'*A
  761. *
  762. RESULT( 9 ) = ZERO
  763. IF( M.GT.CRANK )
  764. $ RESULT( 9 ) = ZQRT17( 'No transpose', 1, M,
  765. $ N, NRHS, COPYA, LDA, B, LDB,
  766. $ COPYB, LDB, C, WORK, LWORK )
  767. *
  768. * Test 10: Check if x is in the rowspace of A
  769. *
  770. RESULT( 10 ) = ZERO
  771. IF( N.GT.CRANK )
  772. $ RESULT( 10 ) = ZQRT14( 'No transpose', M, N,
  773. $ NRHS, COPYA, LDA, B, LDB,
  774. $ WORK, LWORK )
  775. *
  776. * Test ZGELSD
  777. *
  778. * ZGELSD: Compute the minimum-norm solution X
  779. * to min( norm( A * X - B ) ) using a
  780. * divide and conquer SVD.
  781. *
  782. CALL XLAENV( 9, 25 )
  783. *
  784. CALL ZLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  785. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  786. $ LDB )
  787. *
  788. SRNAMT = 'ZGELSD'
  789. CALL ZGELSD( M, N, NRHS, A, LDA, B, LDB, S,
  790. $ RCOND, CRANK, WORK, LWORK, RWORK,
  791. $ IWORK, INFO )
  792. IF( INFO.NE.0 )
  793. $ CALL ALAERH( PATH, 'ZGELSD', INFO, 0, ' ', M,
  794. $ N, NRHS, -1, NB, ITYPE, NFAIL,
  795. $ NERRS, NOUT )
  796. *
  797. * Test 11: Compute relative error in svd
  798. *
  799. IF( RANK.GT.0 ) THEN
  800. CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
  801. RESULT( 11 ) = DASUM( MNMIN, S, 1 ) /
  802. $ DASUM( MNMIN, COPYS, 1 ) /
  803. $ ( EPS*DBLE( MNMIN ) )
  804. ELSE
  805. RESULT( 11 ) = ZERO
  806. END IF
  807. *
  808. * Test 12: Compute error in solution
  809. *
  810. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  811. $ LDWORK )
  812. CALL ZQRT16( 'No transpose', M, N, NRHS, COPYA,
  813. $ LDA, B, LDB, WORK, LDWORK, RWORK,
  814. $ RESULT( 12 ) )
  815. *
  816. * Test 13: Check norm of r'*A
  817. *
  818. RESULT( 13 ) = ZERO
  819. IF( M.GT.CRANK )
  820. $ RESULT( 13 ) = ZQRT17( 'No transpose', 1, M,
  821. $ N, NRHS, COPYA, LDA, B, LDB,
  822. $ COPYB, LDB, C, WORK, LWORK )
  823. *
  824. * Test 14: Check if x is in the rowspace of A
  825. *
  826. RESULT( 14 ) = ZERO
  827. IF( N.GT.CRANK )
  828. $ RESULT( 14 ) = ZQRT14( 'No transpose', M, N,
  829. $ NRHS, COPYA, LDA, B, LDB,
  830. $ WORK, LWORK )
  831. *
  832. * Print information about the tests that did not
  833. * pass the threshold.
  834. *
  835. DO 80 K = 3, 14
  836. IF( RESULT( K ).GE.THRESH ) THEN
  837. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  838. $ CALL ALAHD( NOUT, PATH )
  839. WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
  840. $ ITYPE, K, RESULT( K )
  841. NFAIL = NFAIL + 1
  842. END IF
  843. 80 CONTINUE
  844. NRUN = NRUN + 12
  845. *
  846. 90 CONTINUE
  847. 100 CONTINUE
  848. 110 CONTINUE
  849. 120 CONTINUE
  850. 130 CONTINUE
  851. 140 CONTINUE
  852. *
  853. * Print a summary of the results.
  854. *
  855. CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
  856. *
  857. 9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4,
  858. $ ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 )
  859. 9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4,
  860. $ ', type', I2, ', test(', I2, ')=', G12.5 )
  861. 9997 FORMAT( ' TRANS=''', A1,' M=', I5, ', N=', I5, ', NRHS=', I4,
  862. $ ', MB=', I4,', NB=', I4,', type', I2,
  863. $ ', test(', I2, ')=', G12.5 )
  864. *
  865. DEALLOCATE( WORK )
  866. DEALLOCATE( IWORK )
  867. DEALLOCATE( RWORK )
  868. RETURN
  869. *
  870. * End of ZDRVLS
  871. *
  872. END