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

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870
  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. *> \ingroup complex16_lin
  187. *
  188. * =====================================================================
  189. SUBROUTINE ZDRVLS( 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. DOUBLE PRECISION THRESH
  201. * ..
  202. * .. Array Arguments ..
  203. LOGICAL DOTYPE( * )
  204. INTEGER MVAL( * ), NBVAL( * ), NSVAL( * ),
  205. $ NVAL( * ), NXVAL( * )
  206. DOUBLE PRECISION COPYS( * ), S( * )
  207. COMPLEX*16 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. DOUBLE PRECISION ONE, ZERO
  218. PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  219. COMPLEX*16 CONE, CZERO
  220. PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
  221. $ CZERO = ( 0.0D+0, 0.0D+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_ZGELS, LWORK_ZGETSLS, LWORK_ZGELSS,
  232. $ LWORK_ZGELSY, LWORK_ZGELSD,
  233. $ LRWORK_ZGELSY, LRWORK_ZGELSS, LRWORK_ZGELSD
  234. DOUBLE PRECISION EPS, NORMA, NORMB, RCOND
  235. * ..
  236. * .. Local Arrays ..
  237. INTEGER ISEED( 4 ), ISEEDY( 4 ), IWQ( 1 )
  238. DOUBLE PRECISION RESULT( NTESTS ), RWQ( 1 )
  239. COMPLEX*16 WQ( 1 )
  240. * ..
  241. * .. Allocatable Arrays ..
  242. COMPLEX*16, ALLOCATABLE :: WORK (:)
  243. DOUBLE PRECISION, ALLOCATABLE :: RWORK (:), WORK2 (:)
  244. INTEGER, ALLOCATABLE :: IWORK (:)
  245. * ..
  246. * .. External Functions ..
  247. DOUBLE PRECISION DASUM, DLAMCH, ZQRT12, ZQRT14, ZQRT17
  248. EXTERNAL DASUM, DLAMCH, ZQRT12, ZQRT14, ZQRT17
  249. * ..
  250. * .. External Subroutines ..
  251. EXTERNAL ALAERH, ALAHD, ALASVM, DAXPY, DLASRT, XLAENV,
  252. $ ZDSCAL, ZERRLS, ZGELS, ZGELSD, ZGELSS,
  253. $ ZGELSY, ZGEMM, ZLACPY, ZLARNV, ZQRT13, ZQRT15,
  254. $ ZQRT16, ZGETSLS
  255. * ..
  256. * .. Intrinsic Functions ..
  257. INTRINSIC DBLE, MAX, MIN, INT, 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 ) = 'Zomplex 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 = DLAMCH( '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 ZERRLS( 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. * ZQRT14, ZQRT17 (two side cases), ZQRT15 and ZQRT12
  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 ZGELS
  361. CALL ZGELS( TRANS, M, N, NRHS, A, LDA,
  362. $ B, LDB, WQ, -1, INFO )
  363. LWORK_ZGELS = INT ( WQ( 1 ) )
  364. * Compute workspace needed for ZGETSLS
  365. CALL ZGETSLS( TRANS, M, N, NRHS, A, LDA,
  366. $ B, LDB, WQ, -1, INFO )
  367. LWORK_ZGETSLS = INT( WQ( 1 ) )
  368. ENDDO
  369. END IF
  370. * Compute workspace needed for ZGELSY
  371. CALL ZGELSY( M, N, NRHS, A, LDA, B, LDB, IWQ,
  372. $ RCOND, CRANK, WQ, -1, RWQ, INFO )
  373. LWORK_ZGELSY = INT( WQ( 1 ) )
  374. LRWORK_ZGELSY = 2*N
  375. * Compute workspace needed for ZGELSS
  376. CALL ZGELSS( M, N, NRHS, A, LDA, B, LDB, S,
  377. $ RCOND, CRANK, WQ, -1 , RWQ,
  378. $ INFO )
  379. LWORK_ZGELSS = INT( WQ( 1 ) )
  380. LRWORK_ZGELSS = 5*MNMIN
  381. * Compute workspace needed for ZGELSD
  382. CALL ZGELSD( M, N, NRHS, A, LDA, B, LDB, S,
  383. $ RCOND, CRANK, WQ, -1, RWQ, IWQ,
  384. $ INFO )
  385. LWORK_ZGELSD = INT( WQ( 1 ) )
  386. LRWORK_ZGELSD = INT( RWQ ( 1 ) )
  387. * Compute LIWORK workspace needed for ZGELSY and ZGELSD
  388. LIWORK = MAX( LIWORK, N, IWQ( 1 ) )
  389. * Compute LRWORK workspace needed for ZGELSY, ZGELSS and ZGELSD
  390. LRWORK = MAX( LRWORK, LRWORK_ZGELSY,
  391. $ LRWORK_ZGELSS, LRWORK_ZGELSD )
  392. * Compute LWORK workspace needed for all functions
  393. LWORK = MAX( LWORK, LWORK_ZGELS, LWORK_ZGETSLS,
  394. $ LWORK_ZGELSY, LWORK_ZGELSS,
  395. $ LWORK_ZGELSD )
  396. END IF
  397. ENDDO
  398. ENDDO
  399. ENDDO
  400. ENDDO
  401. ENDDO
  402. *
  403. LWLSY = LWORK
  404. *
  405. ALLOCATE( WORK( LWORK ) )
  406. ALLOCATE( WORK2( 2 * LWORK ) )
  407. ALLOCATE( IWORK( LIWORK ) )
  408. ALLOCATE( RWORK( LRWORK ) )
  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 ZGELS
  432. *
  433. * Generate a matrix of scaling type ISCALE
  434. *
  435. CALL ZQRT13( 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 ZLARNV( 2, ISEED, NCOLS*NRHS,
  458. $ WORK )
  459. CALL ZDSCAL( NCOLS*NRHS,
  460. $ ONE / DBLE( NCOLS ), WORK,
  461. $ 1 )
  462. END IF
  463. CALL ZGEMM( TRANS, 'No transpose', NROWS,
  464. $ NRHS, NCOLS, CONE, COPYA, LDA,
  465. $ WORK, LDWORK, CZERO, B, LDB )
  466. CALL ZLACPY( '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 ZLACPY( 'Full', M, N, COPYA, LDA,
  473. $ A, LDA )
  474. CALL ZLACPY( 'Full', NROWS, NRHS,
  475. $ COPYB, LDB, B, LDB )
  476. END IF
  477. SRNAMT = 'ZGELS '
  478. CALL ZGELS( TRANS, M, N, NRHS, A, LDA, B,
  479. $ LDB, WORK, LWORK, INFO )
  480. *
  481. IF( INFO.NE.0 )
  482. $ CALL ALAERH( PATH, 'ZGELS ', 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 ZLACPY( 'Full', NROWS, NRHS,
  492. $ COPYB, LDB, C, LDB )
  493. CALL ZQRT16( 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 ) = ZQRT17( 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 ) = ZQRT14( 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 ZGETSLS
  534. *
  535. * Generate a matrix of scaling type ISCALE
  536. *
  537. CALL ZQRT13( 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. IF( NCOLS.GT.0 ) THEN
  561. CALL ZLARNV( 2, ISEED, NCOLS*NRHS,
  562. $ WORK )
  563. CALL ZSCAL( NCOLS*NRHS,
  564. $ CONE / DBLE( NCOLS ), WORK,
  565. $ 1 )
  566. END IF
  567. CALL ZGEMM( TRANS, 'No transpose', NROWS,
  568. $ NRHS, NCOLS, CONE, COPYA, LDA,
  569. $ WORK, LDWORK, CZERO, B, LDB )
  570. CALL ZLACPY( 'Full', NROWS, NRHS, B, LDB,
  571. $ COPYB, LDB )
  572. *
  573. * Solve LS or overdetermined system
  574. *
  575. IF( M.GT.0 .AND. N.GT.0 ) THEN
  576. CALL ZLACPY( 'Full', M, N, COPYA, LDA,
  577. $ A, LDA )
  578. CALL ZLACPY( 'Full', NROWS, NRHS,
  579. $ COPYB, LDB, B, LDB )
  580. END IF
  581. SRNAMT = 'ZGETSLS '
  582. CALL ZGETSLS( TRANS, M, N, NRHS, A,
  583. $ LDA, B, LDB, WORK, LWORK, INFO )
  584. IF( INFO.NE.0 )
  585. $ CALL ALAERH( PATH, 'ZGETSLS ', 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. IF( NROWS.GT.0 .AND. NRHS.GT.0 )
  594. $ CALL ZLACPY( 'Full', NROWS, NRHS,
  595. $ COPYB, LDB, C, LDB )
  596. CALL ZQRT16( TRANS, M, N, NRHS, COPYA,
  597. $ LDA, B, LDB, C, LDB, WORK2,
  598. $ RESULT( 15 ) )
  599. *
  600. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
  601. $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
  602. *
  603. * Solving LS system
  604. *
  605. RESULT( 16 ) = ZQRT17( 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 ) = ZQRT14( 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. IF( RESULT( K ).GE.THRESH ) THEN
  623. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  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 ZQRT15( 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 ZGELSY
  656. *
  657. * ZGELSY: Compute the minimum-norm solution
  658. * X to min( norm( A * X - B ) )
  659. * using the rank-revealing orthogonal
  660. * factorization.
  661. *
  662. CALL ZLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  663. CALL ZLACPY( '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 = 'ZGELSY'
  673. CALL ZGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
  674. $ RCOND, CRANK, WORK, LWLSY, RWORK,
  675. $ INFO )
  676. IF( INFO.NE.0 )
  677. $ CALL ALAERH( PATH, 'ZGELSY', 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 ) = ZQRT12( 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 ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  693. $ LDWORK )
  694. CALL ZQRT16( '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. IF( M.GT.CRANK )
  703. $ RESULT( 5 ) = ZQRT17( '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. IF( N.GT.CRANK )
  713. $ RESULT( 6 ) = ZQRT14( 'No transpose', M, N,
  714. $ NRHS, COPYA, LDA, B, LDB,
  715. $ WORK, LWORK )
  716. *
  717. * Test ZGELSS
  718. *
  719. * ZGELSS: Compute the minimum-norm solution
  720. * X to min( norm( A * X - B ) )
  721. * using the SVD.
  722. *
  723. CALL ZLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  724. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  725. $ LDB )
  726. SRNAMT = 'ZGELSS'
  727. CALL ZGELSS( M, N, NRHS, A, LDA, B, LDB, S,
  728. $ RCOND, CRANK, WORK, LWORK, RWORK,
  729. $ INFO )
  730. *
  731. IF( INFO.NE.0 )
  732. $ CALL ALAERH( PATH, 'ZGELSS', 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. IF( RANK.GT.0 ) THEN
  742. CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
  743. RESULT( 7 ) = DASUM( MNMIN, S, 1 ) /
  744. $ DASUM( MNMIN, COPYS, 1 ) /
  745. $ ( EPS*DBLE( MNMIN ) )
  746. ELSE
  747. RESULT( 7 ) = ZERO
  748. END IF
  749. *
  750. * Test 8: Compute error in solution
  751. *
  752. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  753. $ LDWORK )
  754. CALL ZQRT16( '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. IF( M.GT.CRANK )
  762. $ RESULT( 9 ) = ZQRT17( '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. IF( N.GT.CRANK )
  770. $ RESULT( 10 ) = ZQRT14( 'No transpose', M, N,
  771. $ NRHS, COPYA, LDA, B, LDB,
  772. $ WORK, LWORK )
  773. *
  774. * Test ZGELSD
  775. *
  776. * ZGELSD: 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 ZLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  783. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  784. $ LDB )
  785. *
  786. SRNAMT = 'ZGELSD'
  787. CALL ZGELSD( M, N, NRHS, A, LDA, B, LDB, S,
  788. $ RCOND, CRANK, WORK, LWORK, RWORK,
  789. $ IWORK, INFO )
  790. IF( INFO.NE.0 )
  791. $ CALL ALAERH( PATH, 'ZGELSD', INFO, 0, ' ', M,
  792. $ N, NRHS, -1, NB, ITYPE, NFAIL,
  793. $ NERRS, NOUT )
  794. *
  795. * Test 11: Compute relative error in svd
  796. *
  797. IF( RANK.GT.0 ) THEN
  798. CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
  799. RESULT( 11 ) = DASUM( MNMIN, S, 1 ) /
  800. $ DASUM( MNMIN, COPYS, 1 ) /
  801. $ ( EPS*DBLE( MNMIN ) )
  802. ELSE
  803. RESULT( 11 ) = ZERO
  804. END IF
  805. *
  806. * Test 12: Compute error in solution
  807. *
  808. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  809. $ LDWORK )
  810. CALL ZQRT16( '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. IF( M.GT.CRANK )
  818. $ RESULT( 13 ) = ZQRT17( '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. IF( N.GT.CRANK )
  826. $ RESULT( 14 ) = ZQRT14( '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. IF( RESULT( K ).GE.THRESH ) THEN
  835. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  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( IWORK )
  865. DEALLOCATE( RWORK )
  866. RETURN
  867. *
  868. * End of ZDRVLS
  869. *
  870. END