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.

ddrvls.f 25 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661
  1. *> \brief \b DDRVLS
  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 DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
  12. * NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
  13. * COPYB, C, S, COPYS, WORK, IWORK, 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 IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
  23. * $ NVAL( * ), NXVAL( * )
  24. * DOUBLE PRECISION A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
  25. * $ COPYS( * ), S( * ), WORK( * )
  26. * ..
  27. *
  28. *
  29. *> \par Purpose:
  30. * =============
  31. *>
  32. *> \verbatim
  33. *>
  34. *> DDRVLS tests the least squares driver routines DGELS, DGELSS, DGELSY,
  35. *> and DGELSD.
  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 orthogonal 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 orthogonal 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] NNS
  86. *> \verbatim
  87. *> NNS is INTEGER
  88. *> The number of values of NRHS contained in the vector NSVAL.
  89. *> \endverbatim
  90. *>
  91. *> \param[in] NSVAL
  92. *> \verbatim
  93. *> NSVAL is INTEGER array, dimension (NNS)
  94. *> The values of the number of right hand sides NRHS.
  95. *> \endverbatim
  96. *>
  97. *> \param[in] NNB
  98. *> \verbatim
  99. *> NNB is INTEGER
  100. *> The number of values of NB and NX contained in the
  101. *> vectors NBVAL and NXVAL. The blocking parameters are used
  102. *> in pairs (NB,NX).
  103. *> \endverbatim
  104. *>
  105. *> \param[in] NBVAL
  106. *> \verbatim
  107. *> NBVAL is INTEGER array, dimension (NNB)
  108. *> The values of the blocksize NB.
  109. *> \endverbatim
  110. *>
  111. *> \param[in] NXVAL
  112. *> \verbatim
  113. *> NXVAL is INTEGER array, dimension (NNB)
  114. *> The values of the crossover point NX.
  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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (MMAX*NMAX)
  141. *> \endverbatim
  142. *>
  143. *> \param[out] B
  144. *> \verbatim
  145. *> B is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (MMAX*NSMAX)
  153. *> \endverbatim
  154. *>
  155. *> \param[out] C
  156. *> \verbatim
  157. *> C is DOUBLE PRECISION 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[out] WORK
  173. *> \verbatim
  174. *> WORK is DOUBLE PRECISION array,
  175. *> dimension (MMAX*NMAX + 4*NMAX + MMAX).
  176. *> \endverbatim
  177. *>
  178. *> \param[out] IWORK
  179. *> \verbatim
  180. *> IWORK is INTEGER array, dimension (15*NMAX)
  181. *> \endverbatim
  182. *>
  183. *> \param[in] NOUT
  184. *> \verbatim
  185. *> NOUT is INTEGER
  186. *> The unit number for output.
  187. *> \endverbatim
  188. *
  189. * Authors:
  190. * ========
  191. *
  192. *> \author Univ. of Tennessee
  193. *> \author Univ. of California Berkeley
  194. *> \author Univ. of Colorado Denver
  195. *> \author NAG Ltd.
  196. *
  197. *> \date November 2015
  198. *
  199. *> \ingroup double_lin
  200. *
  201. * =====================================================================
  202. SUBROUTINE DDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB,
  203. $ NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B,
  204. $ COPYB, C, S, COPYS, WORK, IWORK, NOUT )
  205. *
  206. * -- LAPACK test routine (version 3.6.0) --
  207. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  208. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  209. * November 2015
  210. *
  211. * .. Scalar Arguments ..
  212. LOGICAL TSTERR
  213. INTEGER NM, NN, NNB, NNS, NOUT
  214. DOUBLE PRECISION THRESH
  215. * ..
  216. * .. Array Arguments ..
  217. LOGICAL DOTYPE( * )
  218. INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
  219. $ NVAL( * ), NXVAL( * )
  220. DOUBLE PRECISION A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ),
  221. $ COPYS( * ), S( * ), WORK( * )
  222. * ..
  223. *
  224. * =====================================================================
  225. *
  226. * .. Parameters ..
  227. INTEGER NTESTS
  228. PARAMETER ( NTESTS = 14 )
  229. INTEGER SMLSIZ
  230. PARAMETER ( SMLSIZ = 25 )
  231. DOUBLE PRECISION ONE, TWO, ZERO
  232. PARAMETER ( ONE = 1.0D0, TWO = 2.0D0, ZERO = 0.0D0 )
  233. * ..
  234. * .. Local Scalars ..
  235. CHARACTER TRANS
  236. CHARACTER*3 PATH
  237. INTEGER CRANK, I, IM, IN, INB, INFO, INS, IRANK,
  238. $ ISCALE, ITRAN, ITYPE, J, K, LDA, LDB, LDWORK,
  239. $ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS,
  240. $ NFAIL, NLVL, NRHS, NROWS, NRUN, RANK
  241. DOUBLE PRECISION EPS, NORMA, NORMB, RCOND
  242. * ..
  243. * .. Local Arrays ..
  244. INTEGER ISEED( 4 ), ISEEDY( 4 )
  245. DOUBLE PRECISION RESULT( NTESTS )
  246. * ..
  247. * .. External Functions ..
  248. DOUBLE PRECISION DASUM, DLAMCH, DQRT12, DQRT14, DQRT17
  249. EXTERNAL DASUM, DLAMCH, DQRT12, DQRT14, DQRT17
  250. * ..
  251. * .. External Subroutines ..
  252. EXTERNAL ALAERH, ALAHD, ALASVM, DAXPY, DERRLS, DGELS,
  253. $ DGELSD, DGELSS, DGELSY, DGEMM, DLACPY,
  254. $ DLARNV, DLASRT, DQRT13, DQRT15, DQRT16, DSCAL,
  255. $ XLAENV
  256. * ..
  257. * .. Intrinsic Functions ..
  258. INTRINSIC DBLE, INT, LOG, MAX, MIN, SQRT
  259. * ..
  260. * .. Scalars in Common ..
  261. LOGICAL LERR, OK
  262. CHARACTER*32 SRNAMT
  263. INTEGER INFOT, IOUNIT
  264. * ..
  265. * .. Common blocks ..
  266. COMMON / INFOC / INFOT, IOUNIT, OK, LERR
  267. COMMON / SRNAMC / SRNAMT
  268. * ..
  269. * .. Data statements ..
  270. DATA ISEEDY / 1988, 1989, 1990, 1991 /
  271. * ..
  272. * .. Executable Statements ..
  273. *
  274. * Initialize constants and the random number seed.
  275. *
  276. PATH( 1: 1 ) = 'Double precision'
  277. PATH( 2: 3 ) = 'LS'
  278. NRUN = 0
  279. NFAIL = 0
  280. NERRS = 0
  281. DO 10 I = 1, 4
  282. ISEED( I ) = ISEEDY( I )
  283. 10 CONTINUE
  284. EPS = DLAMCH( 'Epsilon' )
  285. *
  286. * Threshold for rank estimation
  287. *
  288. RCOND = SQRT( EPS ) - ( SQRT( EPS )-EPS ) / 2
  289. *
  290. * Test the error exits
  291. *
  292. CALL XLAENV( 2, 2 )
  293. CALL XLAENV( 9, SMLSIZ )
  294. IF( TSTERR )
  295. $ CALL DERRLS( PATH, NOUT )
  296. *
  297. * Print the header if NM = 0 or NN = 0 and THRESH = 0.
  298. *
  299. IF( ( NM.EQ.0 .OR. NN.EQ.0 ) .AND. THRESH.EQ.ZERO )
  300. $ CALL ALAHD( NOUT, PATH )
  301. INFOT = 0
  302. CALL XLAENV( 2, 2 )
  303. CALL XLAENV( 9, SMLSIZ )
  304. *
  305. DO 150 IM = 1, NM
  306. M = MVAL( IM )
  307. LDA = MAX( 1, M )
  308. *
  309. DO 140 IN = 1, NN
  310. N = NVAL( IN )
  311. MNMIN = MIN( M, N )
  312. LDB = MAX( 1, M, N )
  313. *
  314. DO 130 INS = 1, NNS
  315. NRHS = NSVAL( INS )
  316. NLVL = MAX( INT( LOG( MAX( ONE, DBLE( MNMIN ) ) /
  317. $ DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1, 0 )
  318. LWORK = MAX( 1, ( M+NRHS )*( N+2 ), ( N+NRHS )*( M+2 ),
  319. $ M*N+4*MNMIN+MAX( M, N ), 12*MNMIN+2*MNMIN*SMLSIZ+
  320. $ 8*MNMIN*NLVL+MNMIN*NRHS+(SMLSIZ+1)**2 )
  321. *
  322. DO 120 IRANK = 1, 2
  323. DO 110 ISCALE = 1, 3
  324. ITYPE = ( IRANK-1 )*3 + ISCALE
  325. IF( .NOT.DOTYPE( ITYPE ) )
  326. $ GO TO 110
  327. *
  328. IF( IRANK.EQ.1 ) THEN
  329. *
  330. * Test DGELS
  331. *
  332. * Generate a matrix of scaling type ISCALE
  333. *
  334. CALL DQRT13( ISCALE, M, N, COPYA, LDA, NORMA,
  335. $ ISEED )
  336. DO 40 INB = 1, NNB
  337. NB = NBVAL( INB )
  338. CALL XLAENV( 1, NB )
  339. CALL XLAENV( 3, NXVAL( INB ) )
  340. *
  341. DO 30 ITRAN = 1, 2
  342. IF( ITRAN.EQ.1 ) THEN
  343. TRANS = 'N'
  344. NROWS = M
  345. NCOLS = N
  346. ELSE
  347. TRANS = 'T'
  348. NROWS = N
  349. NCOLS = M
  350. END IF
  351. LDWORK = MAX( 1, NCOLS )
  352. *
  353. * Set up a consistent rhs
  354. *
  355. IF( NCOLS.GT.0 ) THEN
  356. CALL DLARNV( 2, ISEED, NCOLS*NRHS,
  357. $ WORK )
  358. CALL DSCAL( NCOLS*NRHS,
  359. $ ONE / DBLE( NCOLS ), WORK,
  360. $ 1 )
  361. END IF
  362. CALL DGEMM( TRANS, 'No transpose', NROWS,
  363. $ NRHS, NCOLS, ONE, COPYA, LDA,
  364. $ WORK, LDWORK, ZERO, B, LDB )
  365. CALL DLACPY( 'Full', NROWS, NRHS, B, LDB,
  366. $ COPYB, LDB )
  367. *
  368. * Solve LS or overdetermined system
  369. *
  370. IF( M.GT.0 .AND. N.GT.0 ) THEN
  371. CALL DLACPY( 'Full', M, N, COPYA, LDA,
  372. $ A, LDA )
  373. CALL DLACPY( 'Full', NROWS, NRHS,
  374. $ COPYB, LDB, B, LDB )
  375. END IF
  376. SRNAMT = 'DGELS '
  377. CALL DGELS( TRANS, M, N, NRHS, A, LDA, B,
  378. $ LDB, WORK, LWORK, INFO )
  379. IF( INFO.NE.0 )
  380. $ CALL ALAERH( PATH, 'DGELS ', INFO, 0,
  381. $ TRANS, M, N, NRHS, -1, NB,
  382. $ ITYPE, NFAIL, NERRS,
  383. $ NOUT )
  384. *
  385. * Check correctness of results
  386. *
  387. LDWORK = MAX( 1, NROWS )
  388. IF( NROWS.GT.0 .AND. NRHS.GT.0 )
  389. $ CALL DLACPY( 'Full', NROWS, NRHS,
  390. $ COPYB, LDB, C, LDB )
  391. CALL DQRT16( TRANS, M, N, NRHS, COPYA,
  392. $ LDA, B, LDB, C, LDB, WORK,
  393. $ RESULT( 1 ) )
  394. *
  395. IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR.
  396. $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN
  397. *
  398. * Solving LS system
  399. *
  400. RESULT( 2 ) = DQRT17( TRANS, 1, M, N,
  401. $ NRHS, COPYA, LDA, B, LDB,
  402. $ COPYB, LDB, C, WORK,
  403. $ LWORK )
  404. ELSE
  405. *
  406. * Solving overdetermined system
  407. *
  408. RESULT( 2 ) = DQRT14( TRANS, M, N,
  409. $ NRHS, COPYA, LDA, B, LDB,
  410. $ WORK, LWORK )
  411. END IF
  412. *
  413. * Print information about the tests that
  414. * did not pass the threshold.
  415. *
  416. DO 20 K = 1, 2
  417. IF( RESULT( K ).GE.THRESH ) THEN
  418. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  419. $ CALL ALAHD( NOUT, PATH )
  420. WRITE( NOUT, FMT = 9999 )TRANS, M,
  421. $ N, NRHS, NB, ITYPE, K,
  422. $ RESULT( K )
  423. NFAIL = NFAIL + 1
  424. END IF
  425. 20 CONTINUE
  426. NRUN = NRUN + 2
  427. 30 CONTINUE
  428. 40 CONTINUE
  429. END IF
  430. *
  431. * Generate a matrix of scaling type ISCALE and rank
  432. * type IRANK.
  433. *
  434. CALL DQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA,
  435. $ COPYB, LDB, COPYS, RANK, NORMA, NORMB,
  436. $ ISEED, WORK, LWORK )
  437. *
  438. * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
  439. *
  440. LDWORK = MAX( 1, M )
  441. *
  442. * Loop for testing different block sizes.
  443. *
  444. DO 100 INB = 1, NNB
  445. NB = NBVAL( INB )
  446. CALL XLAENV( 1, NB )
  447. CALL XLAENV( 3, NXVAL( INB ) )
  448. *
  449. * Test DGELSY
  450. *
  451. * DGELSY: Compute the minimum-norm solution X
  452. * to min( norm( A * X - B ) )
  453. * using the rank-revealing orthogonal
  454. * factorization.
  455. *
  456. * Initialize vector IWORK.
  457. *
  458. DO 70 J = 1, N
  459. IWORK( J ) = 0
  460. 70 CONTINUE
  461. *
  462. * Set LWLSY to the adequate value.
  463. *
  464. LWLSY = MAX( 1, MNMIN+2*N+NB*( N+1 ),
  465. $ 2*MNMIN+NB*NRHS )
  466. *
  467. CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  468. CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  469. $ LDB )
  470. *
  471. SRNAMT = 'DGELSY'
  472. CALL DGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
  473. $ RCOND, CRANK, WORK, LWLSY, INFO )
  474. IF( INFO.NE.0 )
  475. $ CALL ALAERH( PATH, 'DGELSY', INFO, 0, ' ', M,
  476. $ N, NRHS, -1, NB, ITYPE, NFAIL,
  477. $ NERRS, NOUT )
  478. *
  479. * Test 3: Compute relative error in svd
  480. * workspace: M*N + 4*MIN(M,N) + MAX(M,N)
  481. *
  482. RESULT( 3 ) = DQRT12( CRANK, CRANK, A, LDA,
  483. $ COPYS, WORK, LWORK )
  484. *
  485. * Test 4: Compute error in solution
  486. * workspace: M*NRHS + M
  487. *
  488. CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  489. $ LDWORK )
  490. CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
  491. $ LDA, B, LDB, WORK, LDWORK,
  492. $ WORK( M*NRHS+1 ), RESULT( 4 ) )
  493. *
  494. * Test 5: Check norm of r'*A
  495. * workspace: NRHS*(M+N)
  496. *
  497. RESULT( 5 ) = ZERO
  498. IF( M.GT.CRANK )
  499. $ RESULT( 5 ) = DQRT17( 'No transpose', 1, M,
  500. $ N, NRHS, COPYA, LDA, B, LDB,
  501. $ COPYB, LDB, C, WORK, LWORK )
  502. *
  503. * Test 6: Check if x is in the rowspace of A
  504. * workspace: (M+NRHS)*(N+2)
  505. *
  506. RESULT( 6 ) = ZERO
  507. *
  508. IF( N.GT.CRANK )
  509. $ RESULT( 6 ) = DQRT14( 'No transpose', M, N,
  510. $ NRHS, COPYA, LDA, B, LDB,
  511. $ WORK, LWORK )
  512. *
  513. * Test DGELSS
  514. *
  515. * DGELSS: Compute the minimum-norm solution X
  516. * to min( norm( A * X - B ) )
  517. * using the SVD.
  518. *
  519. CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  520. CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  521. $ LDB )
  522. SRNAMT = 'DGELSS'
  523. CALL DGELSS( M, N, NRHS, A, LDA, B, LDB, S,
  524. $ RCOND, CRANK, WORK, LWORK, INFO )
  525. IF( INFO.NE.0 )
  526. $ CALL ALAERH( PATH, 'DGELSS', INFO, 0, ' ', M,
  527. $ N, NRHS, -1, NB, ITYPE, NFAIL,
  528. $ NERRS, NOUT )
  529. *
  530. * workspace used: 3*min(m,n) +
  531. * max(2*min(m,n),nrhs,max(m,n))
  532. *
  533. * Test 7: Compute relative error in svd
  534. *
  535. IF( RANK.GT.0 ) THEN
  536. CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
  537. RESULT( 7 ) = DASUM( MNMIN, S, 1 ) /
  538. $ DASUM( MNMIN, COPYS, 1 ) /
  539. $ ( EPS*DBLE( MNMIN ) )
  540. ELSE
  541. RESULT( 7 ) = ZERO
  542. END IF
  543. *
  544. * Test 8: Compute error in solution
  545. *
  546. CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  547. $ LDWORK )
  548. CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
  549. $ LDA, B, LDB, WORK, LDWORK,
  550. $ WORK( M*NRHS+1 ), RESULT( 8 ) )
  551. *
  552. * Test 9: Check norm of r'*A
  553. *
  554. RESULT( 9 ) = ZERO
  555. IF( M.GT.CRANK )
  556. $ RESULT( 9 ) = DQRT17( 'No transpose', 1, M,
  557. $ N, NRHS, COPYA, LDA, B, LDB,
  558. $ COPYB, LDB, C, WORK, LWORK )
  559. *
  560. * Test 10: Check if x is in the rowspace of A
  561. *
  562. RESULT( 10 ) = ZERO
  563. IF( N.GT.CRANK )
  564. $ RESULT( 10 ) = DQRT14( 'No transpose', M, N,
  565. $ NRHS, COPYA, LDA, B, LDB,
  566. $ WORK, LWORK )
  567. *
  568. * Test DGELSD
  569. *
  570. * DGELSD: Compute the minimum-norm solution X
  571. * to min( norm( A * X - B ) ) using a
  572. * divide and conquer SVD.
  573. *
  574. * Initialize vector IWORK.
  575. *
  576. DO 80 J = 1, N
  577. IWORK( J ) = 0
  578. 80 CONTINUE
  579. *
  580. CALL DLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
  581. CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, B,
  582. $ LDB )
  583. *
  584. SRNAMT = 'DGELSD'
  585. CALL DGELSD( M, N, NRHS, A, LDA, B, LDB, S,
  586. $ RCOND, CRANK, WORK, LWORK, IWORK,
  587. $ INFO )
  588. IF( INFO.NE.0 )
  589. $ CALL ALAERH( PATH, 'DGELSD', INFO, 0, ' ', M,
  590. $ N, NRHS, -1, NB, ITYPE, NFAIL,
  591. $ NERRS, NOUT )
  592. *
  593. * Test 11: Compute relative error in svd
  594. *
  595. IF( RANK.GT.0 ) THEN
  596. CALL DAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
  597. RESULT( 11 ) = DASUM( MNMIN, S, 1 ) /
  598. $ DASUM( MNMIN, COPYS, 1 ) /
  599. $ ( EPS*DBLE( MNMIN ) )
  600. ELSE
  601. RESULT( 11 ) = ZERO
  602. END IF
  603. *
  604. * Test 12: Compute error in solution
  605. *
  606. CALL DLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
  607. $ LDWORK )
  608. CALL DQRT16( 'No transpose', M, N, NRHS, COPYA,
  609. $ LDA, B, LDB, WORK, LDWORK,
  610. $ WORK( M*NRHS+1 ), RESULT( 12 ) )
  611. *
  612. * Test 13: Check norm of r'*A
  613. *
  614. RESULT( 13 ) = ZERO
  615. IF( M.GT.CRANK )
  616. $ RESULT( 13 ) = DQRT17( 'No transpose', 1, M,
  617. $ N, NRHS, COPYA, LDA, B, LDB,
  618. $ COPYB, LDB, C, WORK, LWORK )
  619. *
  620. * Test 14: Check if x is in the rowspace of A
  621. *
  622. RESULT( 14 ) = ZERO
  623. IF( N.GT.CRANK )
  624. $ RESULT( 14 ) = DQRT14( 'No transpose', M, N,
  625. $ NRHS, COPYA, LDA, B, LDB,
  626. $ WORK, LWORK )
  627. *
  628. * Print information about the tests that did not
  629. * pass the threshold.
  630. *
  631. DO 90 K = 3, NTESTS
  632. IF( RESULT( K ).GE.THRESH ) THEN
  633. IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
  634. $ CALL ALAHD( NOUT, PATH )
  635. WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
  636. $ ITYPE, K, RESULT( K )
  637. NFAIL = NFAIL + 1
  638. END IF
  639. 90 CONTINUE
  640. NRUN = NRUN + 12
  641. *
  642. 100 CONTINUE
  643. 110 CONTINUE
  644. 120 CONTINUE
  645. 130 CONTINUE
  646. 140 CONTINUE
  647. 150 CONTINUE
  648. *
  649. * Print a summary of the results.
  650. *
  651. CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
  652. *
  653. 9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4,
  654. $ ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 )
  655. 9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4,
  656. $ ', type', I2, ', test(', I2, ')=', G12.5 )
  657. RETURN
  658. *
  659. * End of DDRVLS
  660. *
  661. END