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

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029
  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, ZGELST,
  35. *> ZGETSLS, ZGELSS, ZGELSY 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 = 18 )
  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_ZGELST, LWORK_ZGETSLS,
  232. $ LWORK_ZGELSS, 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, ZERRLS, ZGELS,
  252. $ ZGELSD, ZGELSS, ZGELST, ZGELSY, ZGEMM,
  253. $ ZGETSLS, ZLACPY, ZLARNV, ZQRT13, ZQRT15,
  254. $ ZQRT16, ZDSCAL, XLAENV
  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, ?GELST, ?GETSLS, ?GELSY, ?GELSS and ?GELSD
  338. * routines.
  339. *
  340. DO IM = 1, NM
  341. M = MVAL( IM )
  342. LDA = MAX( 1, M )
  343. DO IN = 1, NN
  344. N = NVAL( IN )
  345. MNMIN = MAX(MIN( M, N ),1)
  346. LDB = MAX( 1, M, N )
  347. DO INS = 1, NNS
  348. NRHS = NSVAL( INS )
  349. DO IRANK = 1, 2
  350. DO ISCALE = 1, 3
  351. ITYPE = ( IRANK-1 )*3 + ISCALE
  352. IF( DOTYPE( ITYPE ) ) THEN
  353. IF( IRANK.EQ.1 ) THEN
  354. DO ITRAN = 1, 2
  355. IF( ITRAN.EQ.1 ) THEN
  356. TRANS = 'N'
  357. ELSE
  358. TRANS = 'C'
  359. END IF
  360. *
  361. * Compute workspace needed for ZGELS
  362. CALL ZGELS( TRANS, M, N, NRHS, A, LDA,
  363. $ B, LDB, WQ, -1, INFO )
  364. LWORK_ZGELS = INT ( WQ( 1 ) )
  365. * Compute workspace needed for ZGELST
  366. CALL ZGELST( TRANS, M, N, NRHS, A, LDA,
  367. $ B, LDB, WQ, -1, INFO )
  368. LWORK_ZGELST = INT ( WQ ( 1 ) )
  369. * Compute workspace needed for ZGETSLS
  370. CALL ZGETSLS( TRANS, M, N, NRHS, A, LDA,
  371. $ B, LDB, WQ, -1, INFO )
  372. LWORK_ZGETSLS = INT( WQ( 1 ) )
  373. ENDDO
  374. END IF
  375. * Compute workspace needed for ZGELSY
  376. CALL ZGELSY( M, N, NRHS, A, LDA, B, LDB, IWQ,
  377. $ RCOND, CRANK, WQ, -1, RWQ, INFO )
  378. LWORK_ZGELSY = INT( WQ( 1 ) )
  379. LRWORK_ZGELSY = 2*N
  380. * Compute workspace needed for ZGELSS
  381. CALL ZGELSS( M, N, NRHS, A, LDA, B, LDB, S,
  382. $ RCOND, CRANK, WQ, -1 , RWQ,
  383. $ INFO )
  384. LWORK_ZGELSS = INT( WQ( 1 ) )
  385. LRWORK_ZGELSS = 5*MNMIN
  386. * Compute workspace needed for ZGELSD
  387. CALL ZGELSD( M, N, NRHS, A, LDA, B, LDB, S,
  388. $ RCOND, CRANK, WQ, -1, RWQ, IWQ,
  389. $ INFO )
  390. LWORK_ZGELSD = INT( WQ( 1 ) )
  391. LRWORK_ZGELSD = INT( RWQ ( 1 ) )
  392. * Compute LIWORK workspace needed for ZGELSY and ZGELSD
  393. LIWORK = MAX( LIWORK, N, IWQ( 1 ) )
  394. * Compute LRWORK workspace needed for ZGELSY, ZGELSS and ZGELSD
  395. LRWORK = MAX( LRWORK, LRWORK_ZGELSY,
  396. $ LRWORK_ZGELSS, LRWORK_ZGELSD )
  397. * Compute LWORK workspace needed for all functions
  398. LWORK = MAX( LWORK, LWORK_ZGELS, LWORK_ZGELST,
  399. $ LWORK_ZGETSLS, LWORK_ZGELSY,
  400. $ LWORK_ZGELSS, LWORK_ZGELSD )
  401. END IF
  402. ENDDO
  403. ENDDO
  404. ENDDO
  405. ENDDO
  406. ENDDO
  407. *
  408. LWLSY = LWORK
  409. *
  410. ALLOCATE( WORK( LWORK ) )
  411. ALLOCATE( WORK2( 2 * LWORK ) )
  412. ALLOCATE( IWORK( LIWORK ) )
  413. ALLOCATE( RWORK( LRWORK ) )
  414. *
  415. DO 140 IM = 1, NM
  416. M = MVAL( IM )
  417. LDA = MAX( 1, M )
  418. *
  419. DO 130 IN = 1, NN
  420. N = NVAL( IN )
  421. MNMIN = MAX(MIN( M, N ),1)
  422. LDB = MAX( 1, M, N )
  423. MB = (MNMIN+1)
  424. *
  425. DO 120 INS = 1, NNS
  426. NRHS = NSVAL( INS )
  427. *
  428. DO 110 IRANK = 1, 2
  429. DO 100 ISCALE = 1, 3
  430. ITYPE = ( IRANK-1 )*3 + ISCALE
  431. IF( .NOT.DOTYPE( ITYPE ) )
  432. $ GO TO 100
  433. * =====================================================
  434. * Begin test ZGELS
  435. * =====================================================
  436. IF( IRANK.EQ.1 ) THEN
  437. *
  438. * Generate a matrix of scaling type ISCALE
  439. *
  440. CALL ZQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
  441. $ ISEED )
  442. *
  443. * Loop for testing different block sizes.
  444. *
  445. DO INB = 1, NNB
  446. NB = NBVAL( INB )
  447. CALL XLAENV( 1, NB )
  448. CALL XLAENV( 3, NXVAL( INB ) )
  449. *
  450. * Loop for testing non-transposed and transposed.
  451. *
  452. DO ITRAN = 1, 2
  453. IF( ITRAN.EQ.1 ) THEN
  454. TRANS = 'N'
  455. NROWS = M
  456. NCOLS = N
  457. ELSE
  458. TRANS = 'C'
  459. NROWS = N
  460. NCOLS = M
  461. END IF
  462. LDWORK = MAX( 1, NCOLS )
  463. *
  464. * Set up a consistent rhs
  465. *
  466. IF( NCOLS.GT.0 ) THEN
  467. CALL ZLARNV( 2, ISEED, NCOLS*NRHS,
  468. $ WORK )
  469. CALL ZDSCAL( NCOLS*NRHS,
  470. $ ONE / DBLE( NCOLS ), WORK,
  471. $ 1 )
  472. END IF
  473. CALL ZGEMM( TRANS, 'No transpose', NROWS,
  474. $ NRHS, NCOLS, CONE, COPYA, LDA,
  475. $ WORK, LDWORK, CZERO, B, LDB )
  476. CALL ZLACPY( 'Full', NROWS, NRHS, B, LDB,
  477. $ COPYB, LDB )
  478. *
  479. * Solve LS or overdetermined system
  480. *
  481. IF( M.GT.0 .AND. N.GT.0 ) THEN
  482. CALL ZLACPY( 'Full', M, N, COPYA, LDA,
  483. $ A, LDA )
  484. CALL ZLACPY( 'Full', NROWS, NRHS,
  485. $ COPYB, LDB, B, LDB )
  486. END IF
  487. SRNAMT = 'ZGELS '
  488. CALL ZGELS( TRANS, M, N, NRHS, A, LDA, B,
  489. $ LDB, WORK, LWORK, INFO )
  490. *
  491. IF( INFO.NE.0 )
  492. $ CALL ALAERH( PATH, 'ZGELS ', INFO, 0,
  493. $ TRANS, M, N, NRHS, -1, NB,
  494. $ ITYPE, NFAIL, NERRS,
  495. $ NOUT )
  496. *
  497. * Test 1: Check correctness of results
  498. * for ZGELS, compute the residual:
  499. * RESID = norm(B - A*X) /
  500. * / ( max(m,n) * norm(A) * norm(X) * EPS )
  501. *
  502. IF( NROWS.GT.0 .AND. NRHS.GT.0 )
  503. $ CALL ZLACPY( 'Full', NROWS, NRHS,
  504. $ COPYB, LDB, C, LDB )
  505. CALL ZQRT16( TRANS, M, N, NRHS, COPYA,
  506. $ LDA, B, LDB, C, LDB, RWORK,
  507. $ RESULT( 1 ) )
  508. *
  509. * Test 2: Check correctness of results
  510. * for ZGELS.
  511. *
  512. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
  513. $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
  514. *
  515. * Solving LS system
  516. *
  517. RESULT( 2 ) = ZQRT17( TRANS, 1, M, N,
  518. $ NRHS, COPYA, LDA, B, LDB,
  519. $ COPYB, LDB, C, WORK,
  520. $ LWORK )
  521. ELSE
  522. *
  523. * Solving overdetermined system
  524. *
  525. RESULT( 2 ) = ZQRT14( TRANS, M, N,
  526. $ NRHS, COPYA, LDA, B, LDB,
  527. $ WORK, LWORK )
  528. END IF
  529. *
  530. * Print information about the tests that
  531. * did not pass the threshold.
  532. *
  533. DO K = 1, 2
  534. IF( RESULT( K ).GE.THRESH ) THEN
  535. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  536. $ CALL ALAHD( NOUT, PATH )
  537. WRITE( NOUT, FMT = 9999 )TRANS, M,
  538. $ N, NRHS, NB, ITYPE, K,
  539. $ RESULT( K )
  540. NFAIL = NFAIL + 1
  541. END IF
  542. END DO
  543. NRUN = NRUN + 2
  544. END DO
  545. END DO
  546. END IF
  547. * =====================================================
  548. * End test ZGELS
  549. * =====================================================
  550. * =====================================================
  551. * Begin test ZGELST
  552. * =====================================================
  553. IF( IRANK.EQ.1 ) THEN
  554. *
  555. * Generate a matrix of scaling type ISCALE
  556. *
  557. CALL ZQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
  558. $ ISEED )
  559. *
  560. * Loop for testing different block sizes.
  561. *
  562. DO INB = 1, NNB
  563. NB = NBVAL( INB )
  564. CALL XLAENV( 1, NB )
  565. CALL XLAENV( 3, NXVAL( INB ) )
  566. *
  567. * Loop for testing non-transposed and transposed.
  568. *
  569. DO ITRAN = 1, 2
  570. IF( ITRAN.EQ.1 ) THEN
  571. TRANS = 'N'
  572. NROWS = M
  573. NCOLS = N
  574. ELSE
  575. TRANS = 'C'
  576. NROWS = N
  577. NCOLS = M
  578. END IF
  579. LDWORK = MAX( 1, NCOLS )
  580. *
  581. * Set up a consistent rhs
  582. *
  583. IF( NCOLS.GT.0 ) THEN
  584. CALL ZLARNV( 2, ISEED, NCOLS*NRHS,
  585. $ WORK )
  586. CALL ZDSCAL( NCOLS*NRHS,
  587. $ ONE / DBLE( NCOLS ), WORK,
  588. $ 1 )
  589. END IF
  590. CALL ZGEMM( TRANS, 'No transpose', NROWS,
  591. $ NRHS, NCOLS, CONE, COPYA, LDA,
  592. $ WORK, LDWORK, CZERO, B, LDB )
  593. CALL ZLACPY( 'Full', NROWS, NRHS, B, LDB,
  594. $ COPYB, LDB )
  595. *
  596. * Solve LS or overdetermined system
  597. *
  598. IF( M.GT.0 .AND. N.GT.0 ) THEN
  599. CALL ZLACPY( 'Full', M, N, COPYA, LDA,
  600. $ A, LDA )
  601. CALL ZLACPY( 'Full', NROWS, NRHS,
  602. $ COPYB, LDB, B, LDB )
  603. END IF
  604. SRNAMT = 'ZGELST'
  605. CALL ZGELST( TRANS, M, N, NRHS, A, LDA, B,
  606. $ LDB, WORK, LWORK, INFO )
  607. *
  608. IF( INFO.NE.0 )
  609. $ CALL ALAERH( PATH, 'ZGELST', INFO, 0,
  610. $ TRANS, M, N, NRHS, -1, NB,
  611. $ ITYPE, NFAIL, NERRS,
  612. $ NOUT )
  613. *
  614. * Test 3: Check correctness of results
  615. * for ZGELST, compute the residual:
  616. * RESID = norm(B - A*X) /
  617. * / ( max(m,n) * norm(A) * norm(X) * EPS )
  618. *
  619. IF( NROWS.GT.0 .AND. NRHS.GT.0 )
  620. $ CALL ZLACPY( 'Full', NROWS, NRHS,
  621. $ COPYB, LDB, C, LDB )
  622. CALL ZQRT16( TRANS, M, N, NRHS, COPYA,
  623. $ LDA, B, LDB, C, LDB, RWORK,
  624. $ RESULT( 3 ) )
  625. *
  626. * Test 4: Check correctness of results
  627. * for ZGELST.
  628. *
  629. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
  630. $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
  631. *
  632. * Solving LS system
  633. *
  634. RESULT( 4 ) = ZQRT17( TRANS, 1, M, N,
  635. $ NRHS, COPYA, LDA, B, LDB,
  636. $ COPYB, LDB, C, WORK,
  637. $ LWORK )
  638. ELSE
  639. *
  640. * Solving overdetermined system
  641. *
  642. RESULT( 4 ) = ZQRT14( TRANS, M, N,
  643. $ NRHS, COPYA, LDA, B, LDB,
  644. $ WORK, LWORK )
  645. END IF
  646. *
  647. * Print information about the tests that
  648. * did not pass the threshold.
  649. *
  650. DO K = 3, 4
  651. IF( RESULT( K ).GE.THRESH ) THEN
  652. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  653. $ CALL ALAHD( NOUT, PATH )
  654. WRITE( NOUT, FMT = 9999 )TRANS, M,
  655. $ N, NRHS, NB, ITYPE, K,
  656. $ RESULT( K )
  657. NFAIL = NFAIL + 1
  658. END IF
  659. END DO
  660. NRUN = NRUN + 2
  661. END DO
  662. END DO
  663. END IF
  664. * =====================================================
  665. * End test ZGELST
  666. * =====================================================
  667. * =====================================================
  668. * Begin test ZGELSTSLS
  669. * =====================================================
  670. IF( IRANK.EQ.1 ) THEN
  671. *
  672. * Generate a matrix of scaling type ISCALE
  673. *
  674. CALL ZQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
  675. $ ISEED )
  676. *
  677. * Loop for testing different block sizes MB.
  678. *
  679. DO INB = 1, NNB
  680. MB = NBVAL( INB )
  681. CALL XLAENV( 1, MB )
  682. *
  683. * Loop for testing different block sizes NB.
  684. *
  685. DO IMB = 1, NNB
  686. NB = NBVAL( IMB )
  687. CALL XLAENV( 2, NB )
  688. *
  689. * Loop for testing non-transposed
  690. * and transposed.
  691. *
  692. DO ITRAN = 1, 2
  693. IF( ITRAN.EQ.1 ) THEN
  694. TRANS = 'N'
  695. NROWS = M
  696. NCOLS = N
  697. ELSE
  698. TRANS = 'C'
  699. NROWS = N
  700. NCOLS = M
  701. END IF
  702. LDWORK = MAX( 1, NCOLS )
  703. *
  704. * Set up a consistent rhs
  705. *
  706. IF( NCOLS.GT.0 ) THEN
  707. CALL ZLARNV( 2, ISEED, NCOLS*NRHS,
  708. $ WORK )
  709. CALL ZSCAL( NCOLS*NRHS,
  710. $ CONE / DBLE( NCOLS ),
  711. $ WORK, 1 )
  712. END IF
  713. CALL ZGEMM( TRANS, 'No transpose',
  714. $ NROWS, NRHS, NCOLS, CONE,
  715. $ COPYA, LDA, WORK, LDWORK,
  716. $ CZERO, B, LDB )
  717. CALL ZLACPY( 'Full', NROWS, NRHS,
  718. $ B, LDB, COPYB, LDB )
  719. *
  720. * Solve LS or overdetermined system
  721. *
  722. IF( M.GT.0 .AND. N.GT.0 ) THEN
  723. CALL ZLACPY( 'Full', M, N,
  724. $ COPYA, LDA, A, LDA )
  725. CALL ZLACPY( 'Full', NROWS, NRHS,
  726. $ COPYB, LDB, B, LDB )
  727. END IF
  728. SRNAMT = 'ZGETSLS '
  729. CALL ZGETSLS( TRANS, M, N, NRHS, A,
  730. $ LDA, B, LDB, WORK, LWORK,
  731. $ INFO )
  732. IF( INFO.NE.0 )
  733. $ CALL ALAERH( PATH, 'ZGETSLS ', INFO,
  734. $ 0, TRANS, M, N, NRHS,
  735. $ -1, NB, ITYPE, NFAIL,
  736. $ NERRS, NOUT )
  737. *
  738. * Test 5: Check correctness of results
  739. * for ZGETSLS, compute the residual:
  740. * RESID = norm(B - A*X) /
  741. * / ( max(m,n) * norm(A) * norm(X) * EPS )
  742. *
  743. IF( NROWS.GT.0 .AND. NRHS.GT.0 )
  744. $ CALL ZLACPY( 'Full', NROWS, NRHS,
  745. $ COPYB, LDB, C, LDB )
  746. CALL ZQRT16( TRANS, M, N, NRHS,
  747. $ COPYA, LDA, B, LDB,
  748. $ C, LDB, WORK2,
  749. $ RESULT( 5 ) )
  750. *
  751. * Test 6: Check correctness of results
  752. * for ZGETSLS.
  753. *
  754. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
  755. $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
  756. *
  757. * Solving LS system, compute:
  758. * r = norm((B- A*X)**T * A) /
  759. * / (norm(A)*norm(B)*max(M,N,NRHS)*EPS)
  760. *
  761. RESULT( 6 ) = ZQRT17( TRANS, 1, M,
  762. $ N, NRHS, COPYA, LDA,
  763. $ B, LDB, COPYB, LDB,
  764. $ C, WORK, LWORK )
  765. ELSE
  766. *
  767. * Solving overdetermined system
  768. *
  769. RESULT( 6 ) = ZQRT14( TRANS, M, N,
  770. $ NRHS, COPYA, LDA, B,
  771. $ LDB, WORK, LWORK )
  772. END IF
  773. *
  774. * Print information about the tests that
  775. * did not pass the threshold.
  776. *
  777. DO K = 5, 6
  778. IF( RESULT( K ).GE.THRESH ) THEN
  779. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  780. $ CALL ALAHD( NOUT, PATH )
  781. WRITE( NOUT, FMT = 9997 )TRANS,
  782. $ M, N, NRHS, MB, NB, ITYPE, K,
  783. $ RESULT( K )
  784. NFAIL = NFAIL + 1
  785. END IF
  786. END DO
  787. NRUN = NRUN + 2
  788. END DO
  789. END DO
  790. END DO
  791. END IF
  792. * =====================================================
  793. * End test ZGELSTSLS
  794. * =====================================================
  795. *
  796. * Generate a matrix of scaling type ISCALE and rank
  797. * type IRANK.
  798. *
  799. CALL ZQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA,
  800. $ COPYB, LDB, COPYS, RANK, NORMA, NORMB,
  801. $ ISEED, WORK, LWORK )
  802. *
  803. * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
  804. *
  805. LDWORK = MAX( 1, M )
  806. *
  807. * Loop for testing different block sizes.
  808. *
  809. DO 90 INB = 1, NNB
  810. NB = NBVAL( INB )
  811. CALL XLAENV( 1, NB )
  812. CALL XLAENV( 3, NXVAL( INB ) )
  813. *
  814. * Test ZGELSY
  815. *
  816. * ZGELSY: Compute the minimum-norm solution
  817. * X to min( norm( A * X - B ) )
  818. * using the rank-revealing orthogonal
  819. * factorization.
  820. *
  821. CALL ZLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  822. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  823. $ LDB )
  824. *
  825. * Initialize vector IWORK.
  826. *
  827. DO 70 J = 1, N
  828. IWORK( J ) = 0
  829. 70 CONTINUE
  830. *
  831. SRNAMT = 'ZGELSY'
  832. CALL ZGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
  833. $ RCOND, CRANK, WORK, LWLSY, RWORK,
  834. $ INFO )
  835. IF( INFO.NE.0 )
  836. $ CALL ALAERH( PATH, 'ZGELSY', INFO, 0, ' ', M,
  837. $ N, NRHS, -1, NB, ITYPE, NFAIL,
  838. $ NERRS, NOUT )
  839. *
  840. * workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS)
  841. *
  842. * Test 7: Compute relative error in svd
  843. * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
  844. *
  845. RESULT( 7 ) = ZQRT12( CRANK, CRANK, A, LDA,
  846. $ COPYS, WORK, LWORK, RWORK )
  847. *
  848. * Test 8: Compute error in solution
  849. * workspace: M*NRHS + M
  850. *
  851. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  852. $ LDWORK )
  853. CALL ZQRT16( 'No transpose', M, N, NRHS, COPYA,
  854. $ LDA, B, LDB, WORK, LDWORK, RWORK,
  855. $ RESULT( 8 ) )
  856. *
  857. * Test 9: Check norm of r'*A
  858. * workspace: NRHS*(M+N)
  859. *
  860. RESULT( 9 ) = ZERO
  861. IF( M.GT.CRANK )
  862. $ RESULT( 9 ) = ZQRT17( 'No transpose', 1, M,
  863. $ N, NRHS, COPYA, LDA, B, LDB,
  864. $ COPYB, LDB, C, WORK, LWORK )
  865. *
  866. * Test 10: Check if x is in the rowspace of A
  867. * workspace: (M+NRHS)*(N+2)
  868. *
  869. RESULT( 10 ) = ZERO
  870. *
  871. IF( N.GT.CRANK )
  872. $ RESULT( 10 ) = ZQRT14( 'No transpose', M, N,
  873. $ NRHS, COPYA, LDA, B, LDB,
  874. $ WORK, LWORK )
  875. *
  876. * Test ZGELSS
  877. *
  878. * ZGELSS: Compute the minimum-norm solution
  879. * X to min( norm( A * X - B ) )
  880. * using the SVD.
  881. *
  882. CALL ZLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  883. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  884. $ LDB )
  885. SRNAMT = 'ZGELSS'
  886. CALL ZGELSS( M, N, NRHS, A, LDA, B, LDB, S,
  887. $ RCOND, CRANK, WORK, LWORK, RWORK,
  888. $ INFO )
  889. *
  890. IF( INFO.NE.0 )
  891. $ CALL ALAERH( PATH, 'ZGELSS', INFO, 0, ' ', M,
  892. $ N, NRHS, -1, NB, ITYPE, NFAIL,
  893. $ NERRS, NOUT )
  894. *
  895. * workspace used: 3*min(m,n) +
  896. * max(2*min(m,n),nrhs,max(m,n))
  897. *
  898. * Test 11: Compute relative error in svd
  899. *
  900. IF( RANK.GT.0 ) THEN
  901. CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
  902. RESULT( 11 ) = DASUM( MNMIN, S, 1 ) /
  903. $ DASUM( MNMIN, COPYS, 1 ) /
  904. $ ( EPS*DBLE( MNMIN ) )
  905. ELSE
  906. RESULT( 11 ) = ZERO
  907. END IF
  908. *
  909. * Test 12: Compute error in solution
  910. *
  911. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  912. $ LDWORK )
  913. CALL ZQRT16( 'No transpose', M, N, NRHS, COPYA,
  914. $ LDA, B, LDB, WORK, LDWORK, RWORK,
  915. $ RESULT( 12 ) )
  916. *
  917. * Test 13: Check norm of r'*A
  918. *
  919. RESULT( 13 ) = ZERO
  920. IF( M.GT.CRANK )
  921. $ RESULT( 13 ) = ZQRT17( 'No transpose', 1, M,
  922. $ N, NRHS, COPYA, LDA, B, LDB,
  923. $ COPYB, LDB, C, WORK, LWORK )
  924. *
  925. * Test 14: Check if x is in the rowspace of A
  926. *
  927. RESULT( 14 ) = ZERO
  928. IF( N.GT.CRANK )
  929. $ RESULT( 14 ) = ZQRT14( 'No transpose', M, N,
  930. $ NRHS, COPYA, LDA, B, LDB,
  931. $ WORK, LWORK )
  932. *
  933. * Test ZGELSD
  934. *
  935. * ZGELSD: Compute the minimum-norm solution X
  936. * to min( norm( A * X - B ) ) using a
  937. * divide and conquer SVD.
  938. *
  939. CALL XLAENV( 9, 25 )
  940. *
  941. CALL ZLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  942. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  943. $ LDB )
  944. *
  945. SRNAMT = 'ZGELSD'
  946. CALL ZGELSD( M, N, NRHS, A, LDA, B, LDB, S,
  947. $ RCOND, CRANK, WORK, LWORK, RWORK,
  948. $ IWORK, INFO )
  949. IF( INFO.NE.0 )
  950. $ CALL ALAERH( PATH, 'ZGELSD', INFO, 0, ' ', M,
  951. $ N, NRHS, -1, NB, ITYPE, NFAIL,
  952. $ NERRS, NOUT )
  953. *
  954. * Test 15: Compute relative error in svd
  955. *
  956. IF( RANK.GT.0 ) THEN
  957. CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
  958. RESULT( 15 ) = DASUM( MNMIN, S, 1 ) /
  959. $ DASUM( MNMIN, COPYS, 1 ) /
  960. $ ( EPS*DBLE( MNMIN ) )
  961. ELSE
  962. RESULT( 15 ) = ZERO
  963. END IF
  964. *
  965. * Test 16: Compute error in solution
  966. *
  967. CALL ZLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  968. $ LDWORK )
  969. CALL ZQRT16( 'No transpose', M, N, NRHS, COPYA,
  970. $ LDA, B, LDB, WORK, LDWORK, RWORK,
  971. $ RESULT( 16 ) )
  972. *
  973. * Test 17: Check norm of r'*A
  974. *
  975. RESULT( 17 ) = ZERO
  976. IF( M.GT.CRANK )
  977. $ RESULT( 17 ) = ZQRT17( 'No transpose', 1, M,
  978. $ N, NRHS, COPYA, LDA, B, LDB,
  979. $ COPYB, LDB, C, WORK, LWORK )
  980. *
  981. * Test 18: Check if x is in the rowspace of A
  982. *
  983. RESULT( 18 ) = ZERO
  984. IF( N.GT.CRANK )
  985. $ RESULT( 18 ) = ZQRT14( 'No transpose', M, N,
  986. $ NRHS, COPYA, LDA, B, LDB,
  987. $ WORK, LWORK )
  988. *
  989. * Print information about the tests that did not
  990. * pass the threshold.
  991. *
  992. DO 80 K = 7, 18
  993. IF( RESULT( K ).GE.THRESH ) THEN
  994. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  995. $ CALL ALAHD( NOUT, PATH )
  996. WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
  997. $ ITYPE, K, RESULT( K )
  998. NFAIL = NFAIL + 1
  999. END IF
  1000. 80 CONTINUE
  1001. NRUN = NRUN + 12
  1002. *
  1003. 90 CONTINUE
  1004. 100 CONTINUE
  1005. 110 CONTINUE
  1006. 120 CONTINUE
  1007. 130 CONTINUE
  1008. 140 CONTINUE
  1009. *
  1010. * Print a summary of the results.
  1011. *
  1012. CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
  1013. *
  1014. 9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4,
  1015. $ ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 )
  1016. 9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4,
  1017. $ ', type', I2, ', test(', I2, ')=', G12.5 )
  1018. 9997 FORMAT( ' TRANS=''', A1,' M=', I5, ', N=', I5, ', NRHS=', I4,
  1019. $ ', MB=', I4,', NB=', I4,', type', I2,
  1020. $ ', test(', I2, ')=', G12.5 )
  1021. *
  1022. DEALLOCATE( WORK )
  1023. DEALLOCATE( IWORK )
  1024. DEALLOCATE( RWORK )
  1025. RETURN
  1026. *
  1027. * End of ZDRVLS
  1028. *
  1029. END