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.

ddrvbd.f 46 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254
  1. *> \brief \b DDRVBD
  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 DDRVBD( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH,
  12. * A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S,
  13. * SSAV, E, WORK, LWORK, IWORK, NOUT, INFO )
  14. *
  15. * .. Scalar Arguments ..
  16. * INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUT, NSIZES,
  17. * $ NTYPES
  18. * DOUBLE PRECISION THRESH
  19. * ..
  20. * .. Array Arguments ..
  21. * LOGICAL DOTYPE( * )
  22. * INTEGER ISEED( 4 ), IWORK( * ), MM( * ), NN( * )
  23. * DOUBLE PRECISION A( LDA, * ), ASAV( LDA, * ), E( * ), S( * ),
  24. * $ SSAV( * ), U( LDU, * ), USAV( LDU, * ),
  25. * $ VT( LDVT, * ), VTSAV( LDVT, * ), WORK( * )
  26. * ..
  27. *
  28. *
  29. *> \par Purpose:
  30. * =============
  31. *>
  32. *> \verbatim
  33. *>
  34. *> DDRVBD checks the singular value decomposition (SVD) drivers
  35. *> DGESVD, DGESDD, DGESVDQ, DGESVJ, DGEJSV, and DGESVDX.
  36. *>
  37. *> Both DGESVD and DGESDD factor A = U diag(S) VT, where U and VT are
  38. *> orthogonal and diag(S) is diagonal with the entries of the array S
  39. *> on its diagonal. The entries of S are the singular values,
  40. *> nonnegative and stored in decreasing order. U and VT can be
  41. *> optionally not computed, overwritten on A, or computed partially.
  42. *>
  43. *> A is M by N. Let MNMIN = min( M, N ). S has dimension MNMIN.
  44. *> U can be M by M or M by MNMIN. VT can be N by N or MNMIN by N.
  45. *>
  46. *> When DDRVBD is called, a number of matrix "sizes" (M's and N's)
  47. *> and a number of matrix "types" are specified. For each size (M,N)
  48. *> and each type of matrix, and for the minimal workspace as well as
  49. *> workspace adequate to permit blocking, an M x N matrix "A" will be
  50. *> generated and used to test the SVD routines. For each matrix, A will
  51. *> be factored as A = U diag(S) VT and the following 12 tests computed:
  52. *>
  53. *> Test for DGESVD:
  54. *>
  55. *> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
  56. *>
  57. *> (2) | I - U'U | / ( M ulp )
  58. *>
  59. *> (3) | I - VT VT' | / ( N ulp )
  60. *>
  61. *> (4) S contains MNMIN nonnegative values in decreasing order.
  62. *> (Return 0 if true, 1/ULP if false.)
  63. *>
  64. *> (5) | U - Upartial | / ( M ulp ) where Upartial is a partially
  65. *> computed U.
  66. *>
  67. *> (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
  68. *> computed VT.
  69. *>
  70. *> (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
  71. *> vector of singular values from the partial SVD
  72. *>
  73. *> Test for DGESDD:
  74. *>
  75. *> (8) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
  76. *>
  77. *> (9) | I - U'U | / ( M ulp )
  78. *>
  79. *> (10) | I - VT VT' | / ( N ulp )
  80. *>
  81. *> (11) S contains MNMIN nonnegative values in decreasing order.
  82. *> (Return 0 if true, 1/ULP if false.)
  83. *>
  84. *> (12) | U - Upartial | / ( M ulp ) where Upartial is a partially
  85. *> computed U.
  86. *>
  87. *> (13) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
  88. *> computed VT.
  89. *>
  90. *> (14) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
  91. *> vector of singular values from the partial SVD
  92. *>
  93. *> Test for DGESVDQ:
  94. *>
  95. *> (36) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
  96. *>
  97. *> (37) | I - U'U | / ( M ulp )
  98. *>
  99. *> (38) | I - VT VT' | / ( N ulp )
  100. *>
  101. *> (39) S contains MNMIN nonnegative values in decreasing order.
  102. *> (Return 0 if true, 1/ULP if false.)
  103. *>
  104. *> Test for DGESVJ:
  105. *>
  106. *> (15) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
  107. *>
  108. *> (16) | I - U'U | / ( M ulp )
  109. *>
  110. *> (17) | I - VT VT' | / ( N ulp )
  111. *>
  112. *> (18) S contains MNMIN nonnegative values in decreasing order.
  113. *> (Return 0 if true, 1/ULP if false.)
  114. *>
  115. *> Test for DGEJSV:
  116. *>
  117. *> (19) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
  118. *>
  119. *> (20) | I - U'U | / ( M ulp )
  120. *>
  121. *> (21) | I - VT VT' | / ( N ulp )
  122. *>
  123. *> (22) S contains MNMIN nonnegative values in decreasing order.
  124. *> (Return 0 if true, 1/ULP if false.)
  125. *>
  126. *> Test for DGESVDX( 'V', 'V', 'A' )/DGESVDX( 'N', 'N', 'A' )
  127. *>
  128. *> (23) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
  129. *>
  130. *> (24) | I - U'U | / ( M ulp )
  131. *>
  132. *> (25) | I - VT VT' | / ( N ulp )
  133. *>
  134. *> (26) S contains MNMIN nonnegative values in decreasing order.
  135. *> (Return 0 if true, 1/ULP if false.)
  136. *>
  137. *> (27) | U - Upartial | / ( M ulp ) where Upartial is a partially
  138. *> computed U.
  139. *>
  140. *> (28) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
  141. *> computed VT.
  142. *>
  143. *> (29) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
  144. *> vector of singular values from the partial SVD
  145. *>
  146. *> Test for DGESVDX( 'V', 'V', 'I' )
  147. *>
  148. *> (30) | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
  149. *>
  150. *> (31) | I - U'U | / ( M ulp )
  151. *>
  152. *> (32) | I - VT VT' | / ( N ulp )
  153. *>
  154. *> Test for DGESVDX( 'V', 'V', 'V' )
  155. *>
  156. *> (33) | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
  157. *>
  158. *> (34) | I - U'U | / ( M ulp )
  159. *>
  160. *> (35) | I - VT VT' | / ( N ulp )
  161. *>
  162. *> The "sizes" are specified by the arrays MM(1:NSIZES) and
  163. *> NN(1:NSIZES); the value of each element pair (MM(j),NN(j))
  164. *> specifies one size. The "types" are specified by a logical array
  165. *> DOTYPE( 1:NTYPES ); if DOTYPE(j) is .TRUE., then matrix type "j"
  166. *> will be generated.
  167. *> Currently, the list of possible types is:
  168. *>
  169. *> (1) The zero matrix.
  170. *> (2) The identity matrix.
  171. *> (3) A matrix of the form U D V, where U and V are orthogonal and
  172. *> D has evenly spaced entries 1, ..., ULP with random signs
  173. *> on the diagonal.
  174. *> (4) Same as (3), but multiplied by the underflow-threshold / ULP.
  175. *> (5) Same as (3), but multiplied by the overflow-threshold * ULP.
  176. *> \endverbatim
  177. *
  178. * Arguments:
  179. * ==========
  180. *
  181. *> \param[in] NSIZES
  182. *> \verbatim
  183. *> NSIZES is INTEGER
  184. *> The number of matrix sizes (M,N) contained in the vectors
  185. *> MM and NN.
  186. *> \endverbatim
  187. *>
  188. *> \param[in] MM
  189. *> \verbatim
  190. *> MM is INTEGER array, dimension (NSIZES)
  191. *> The values of the matrix row dimension M.
  192. *> \endverbatim
  193. *>
  194. *> \param[in] NN
  195. *> \verbatim
  196. *> NN is INTEGER array, dimension (NSIZES)
  197. *> The values of the matrix column dimension N.
  198. *> \endverbatim
  199. *>
  200. *> \param[in] NTYPES
  201. *> \verbatim
  202. *> NTYPES is INTEGER
  203. *> The number of elements in DOTYPE. If it is zero, DDRVBD
  204. *> does nothing. It must be at least zero. If it is MAXTYP+1
  205. *> and NSIZES is 1, then an additional type, MAXTYP+1 is
  206. *> defined, which is to use whatever matrices are in A and B.
  207. *> This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
  208. *> DOTYPE(MAXTYP+1) is .TRUE. .
  209. *> \endverbatim
  210. *>
  211. *> \param[in] DOTYPE
  212. *> \verbatim
  213. *> DOTYPE is LOGICAL array, dimension (NTYPES)
  214. *> If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
  215. *> of type j will be generated. If NTYPES is smaller than the
  216. *> maximum number of types defined (PARAMETER MAXTYP), then
  217. *> types NTYPES+1 through MAXTYP will not be generated. If
  218. *> NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
  219. *> DOTYPE(NTYPES) will be ignored.
  220. *> \endverbatim
  221. *>
  222. *> \param[in,out] ISEED
  223. *> \verbatim
  224. *> ISEED is INTEGER array, dimension (4)
  225. *> On entry, the seed of the random number generator. The array
  226. *> elements should be between 0 and 4095; if not they will be
  227. *> reduced mod 4096. Also, ISEED(4) must be odd.
  228. *> On exit, ISEED is changed and can be used in the next call to
  229. *> DDRVBD to continue the same random number sequence.
  230. *> \endverbatim
  231. *>
  232. *> \param[in] THRESH
  233. *> \verbatim
  234. *> THRESH is DOUBLE PRECISION
  235. *> The threshold value for the test ratios. A result is
  236. *> included in the output file if RESULT >= THRESH. The test
  237. *> ratios are scaled to be O(1), so THRESH should be a small
  238. *> multiple of 1, e.g., 10 or 100. To have every test ratio
  239. *> printed, use THRESH = 0.
  240. *> \endverbatim
  241. *>
  242. *> \param[out] A
  243. *> \verbatim
  244. *> A is DOUBLE PRECISION array, dimension (LDA,NMAX)
  245. *> where NMAX is the maximum value of N in NN.
  246. *> \endverbatim
  247. *>
  248. *> \param[in] LDA
  249. *> \verbatim
  250. *> LDA is INTEGER
  251. *> The leading dimension of the array A. LDA >= max(1,MMAX),
  252. *> where MMAX is the maximum value of M in MM.
  253. *> \endverbatim
  254. *>
  255. *> \param[out] U
  256. *> \verbatim
  257. *> U is DOUBLE PRECISION array, dimension (LDU,MMAX)
  258. *> \endverbatim
  259. *>
  260. *> \param[in] LDU
  261. *> \verbatim
  262. *> LDU is INTEGER
  263. *> The leading dimension of the array U. LDU >= max(1,MMAX).
  264. *> \endverbatim
  265. *>
  266. *> \param[out] VT
  267. *> \verbatim
  268. *> VT is DOUBLE PRECISION array, dimension (LDVT,NMAX)
  269. *> \endverbatim
  270. *>
  271. *> \param[in] LDVT
  272. *> \verbatim
  273. *> LDVT is INTEGER
  274. *> The leading dimension of the array VT. LDVT >= max(1,NMAX).
  275. *> \endverbatim
  276. *>
  277. *> \param[out] ASAV
  278. *> \verbatim
  279. *> ASAV is DOUBLE PRECISION array, dimension (LDA,NMAX)
  280. *> \endverbatim
  281. *>
  282. *> \param[out] USAV
  283. *> \verbatim
  284. *> USAV is DOUBLE PRECISION array, dimension (LDU,MMAX)
  285. *> \endverbatim
  286. *>
  287. *> \param[out] VTSAV
  288. *> \verbatim
  289. *> VTSAV is DOUBLE PRECISION array, dimension (LDVT,NMAX)
  290. *> \endverbatim
  291. *>
  292. *> \param[out] S
  293. *> \verbatim
  294. *> S is DOUBLE PRECISION array, dimension
  295. *> (max(min(MM,NN)))
  296. *> \endverbatim
  297. *>
  298. *> \param[out] SSAV
  299. *> \verbatim
  300. *> SSAV is DOUBLE PRECISION array, dimension
  301. *> (max(min(MM,NN)))
  302. *> \endverbatim
  303. *>
  304. *> \param[out] E
  305. *> \verbatim
  306. *> E is DOUBLE PRECISION array, dimension
  307. *> (max(min(MM,NN)))
  308. *> \endverbatim
  309. *>
  310. *> \param[out] WORK
  311. *> \verbatim
  312. *> WORK is DOUBLE PRECISION array, dimension (LWORK)
  313. *> \endverbatim
  314. *>
  315. *> \param[in] LWORK
  316. *> \verbatim
  317. *> LWORK is INTEGER
  318. *> The number of entries in WORK. This must be at least
  319. *> max(3*MN+MX,5*MN-4)+2*MN**2 for all pairs
  320. *> pairs (MN,MX)=( min(MM(j),NN(j), max(MM(j),NN(j)) )
  321. *> \endverbatim
  322. *>
  323. *> \param[out] IWORK
  324. *> \verbatim
  325. *> IWORK is INTEGER array, dimension at least 8*min(M,N)
  326. *> \endverbatim
  327. *>
  328. *> \param[in] NOUT
  329. *> \verbatim
  330. *> NOUT is INTEGER
  331. *> The FORTRAN unit number for printing out error messages
  332. *> (e.g., if a routine returns IINFO not equal to 0.)
  333. *> \endverbatim
  334. *>
  335. *> \param[out] INFO
  336. *> \verbatim
  337. *> INFO is INTEGER
  338. *> If 0, then everything ran OK.
  339. *> -1: NSIZES < 0
  340. *> -2: Some MM(j) < 0
  341. *> -3: Some NN(j) < 0
  342. *> -4: NTYPES < 0
  343. *> -7: THRESH < 0
  344. *> -10: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
  345. *> -12: LDU < 1 or LDU < MMAX.
  346. *> -14: LDVT < 1 or LDVT < NMAX, where NMAX is max( NN(j) ).
  347. *> -21: LWORK too small.
  348. *> If DLATMS, or DGESVD returns an error code, the
  349. *> absolute value of it is returned.
  350. *> \endverbatim
  351. *
  352. * Authors:
  353. * ========
  354. *
  355. *> \author Univ. of Tennessee
  356. *> \author Univ. of California Berkeley
  357. *> \author Univ. of Colorado Denver
  358. *> \author NAG Ltd.
  359. *
  360. *> \date June 2016
  361. *
  362. *> \ingroup double_eig
  363. *
  364. * =====================================================================
  365. SUBROUTINE DDRVBD( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH,
  366. $ A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S,
  367. $ SSAV, E, WORK, LWORK, IWORK, NOUT, INFO )
  368. *
  369. IMPLICIT NONE
  370. *
  371. * -- LAPACK test routine (version 3.7.0) --
  372. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  373. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  374. * June 2016
  375. *
  376. * .. Scalar Arguments ..
  377. INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUT, NSIZES,
  378. $ NTYPES
  379. DOUBLE PRECISION THRESH
  380. * ..
  381. * .. Array Arguments ..
  382. LOGICAL DOTYPE( * )
  383. INTEGER ISEED( 4 ), IWORK( * ), MM( * ), NN( * )
  384. DOUBLE PRECISION A( LDA, * ), ASAV( LDA, * ), E( * ), S( * ),
  385. $ SSAV( * ), U( LDU, * ), USAV( LDU, * ),
  386. $ VT( LDVT, * ), VTSAV( LDVT, * ), WORK( * )
  387. * ..
  388. *
  389. * =====================================================================
  390. *
  391. * .. Parameters ..
  392. DOUBLE PRECISION ZERO, ONE, TWO, HALF
  393. PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
  394. $ HALF = 0.5D0 )
  395. INTEGER MAXTYP
  396. PARAMETER ( MAXTYP = 5 )
  397. * ..
  398. * .. Local Scalars ..
  399. LOGICAL BADMM, BADNN
  400. CHARACTER JOBQ, JOBU, JOBVT, RANGE
  401. CHARACTER*3 PATH
  402. INTEGER I, IINFO, IJQ, IJU, IJVT, IL,IU, IWS, IWTMP,
  403. $ ITEMP, J, JSIZE, JTYPE, LSWORK, M, MINWRK,
  404. $ MMAX, MNMAX, MNMIN, MTYPES, N, NFAIL,
  405. $ NMAX, NS, NSI, NSV, NTEST
  406. DOUBLE PRECISION ANORM, DIF, DIV, OVFL, RTUNFL, ULP,
  407. $ ULPINV, UNFL, VL, VU
  408. * ..
  409. * .. Local Scalars for DGESVDQ ..
  410. INTEGER LIWORK, LRWORK, NUMRANK
  411. * ..
  412. * .. Local Arrays for DGESVDQ ..
  413. DOUBLE PRECISION RWORK( 2 )
  414. * ..
  415. * .. Local Arrays ..
  416. CHARACTER CJOB( 4 ), CJOBR( 3 ), CJOBV( 2 )
  417. INTEGER IOLDSD( 4 ), ISEED2( 4 )
  418. DOUBLE PRECISION RESULT( 39 )
  419. * ..
  420. * .. External Functions ..
  421. DOUBLE PRECISION DLAMCH, DLARND
  422. EXTERNAL DLAMCH, DLARND
  423. * ..
  424. * .. External Subroutines ..
  425. EXTERNAL ALASVM, DBDT01, DGEJSV, DGESDD, DGESVD,
  426. $ DGESVDQ, DGESVDX, DGESVJ, DLABAD, DLACPY,
  427. $ DLASET, DLATMS, DORT01, DORT03, XERBLA
  428. * ..
  429. * .. Intrinsic Functions ..
  430. INTRINSIC ABS, DBLE, INT, MAX, MIN
  431. * ..
  432. * .. Scalars in Common ..
  433. LOGICAL LERR, OK
  434. CHARACTER*32 SRNAMT
  435. INTEGER INFOT, NUNIT
  436. * ..
  437. * .. Common blocks ..
  438. COMMON / INFOC / INFOT, NUNIT, OK, LERR
  439. COMMON / SRNAMC / SRNAMT
  440. * ..
  441. * .. Data statements ..
  442. DATA CJOB / 'N', 'O', 'S', 'A' /
  443. DATA CJOBR / 'A', 'V', 'I' /
  444. DATA CJOBV / 'N', 'V' /
  445. * ..
  446. * .. Executable Statements ..
  447. *
  448. * Check for errors
  449. *
  450. INFO = 0
  451. BADMM = .FALSE.
  452. BADNN = .FALSE.
  453. MMAX = 1
  454. NMAX = 1
  455. MNMAX = 1
  456. MINWRK = 1
  457. DO 10 J = 1, NSIZES
  458. MMAX = MAX( MMAX, MM( J ) )
  459. IF( MM( J ).LT.0 )
  460. $ BADMM = .TRUE.
  461. NMAX = MAX( NMAX, NN( J ) )
  462. IF( NN( J ).LT.0 )
  463. $ BADNN = .TRUE.
  464. MNMAX = MAX( MNMAX, MIN( MM( J ), NN( J ) ) )
  465. MINWRK = MAX( MINWRK, MAX( 3*MIN( MM( J ),
  466. $ NN( J ) )+MAX( MM( J ), NN( J ) ), 5*MIN( MM( J ),
  467. $ NN( J )-4 ) )+2*MIN( MM( J ), NN( J ) )**2 )
  468. 10 CONTINUE
  469. *
  470. * Check for errors
  471. *
  472. IF( NSIZES.LT.0 ) THEN
  473. INFO = -1
  474. ELSE IF( BADMM ) THEN
  475. INFO = -2
  476. ELSE IF( BADNN ) THEN
  477. INFO = -3
  478. ELSE IF( NTYPES.LT.0 ) THEN
  479. INFO = -4
  480. ELSE IF( LDA.LT.MAX( 1, MMAX ) ) THEN
  481. INFO = -10
  482. ELSE IF( LDU.LT.MAX( 1, MMAX ) ) THEN
  483. INFO = -12
  484. ELSE IF( LDVT.LT.MAX( 1, NMAX ) ) THEN
  485. INFO = -14
  486. ELSE IF( MINWRK.GT.LWORK ) THEN
  487. INFO = -21
  488. END IF
  489. *
  490. IF( INFO.NE.0 ) THEN
  491. CALL XERBLA( 'DDRVBD', -INFO )
  492. RETURN
  493. END IF
  494. *
  495. * Initialize constants
  496. *
  497. PATH( 1: 1 ) = 'Double precision'
  498. PATH( 2: 3 ) = 'BD'
  499. NFAIL = 0
  500. NTEST = 0
  501. UNFL = DLAMCH( 'Safe minimum' )
  502. OVFL = ONE / UNFL
  503. CALL DLABAD( UNFL, OVFL )
  504. ULP = DLAMCH( 'Precision' )
  505. RTUNFL = SQRT( UNFL )
  506. ULPINV = ONE / ULP
  507. INFOT = 0
  508. *
  509. * Loop over sizes, types
  510. *
  511. DO 240 JSIZE = 1, NSIZES
  512. M = MM( JSIZE )
  513. N = NN( JSIZE )
  514. MNMIN = MIN( M, N )
  515. *
  516. IF( NSIZES.NE.1 ) THEN
  517. MTYPES = MIN( MAXTYP, NTYPES )
  518. ELSE
  519. MTYPES = MIN( MAXTYP+1, NTYPES )
  520. END IF
  521. *
  522. DO 230 JTYPE = 1, MTYPES
  523. IF( .NOT.DOTYPE( JTYPE ) )
  524. $ GO TO 230
  525. *
  526. DO 20 J = 1, 4
  527. IOLDSD( J ) = ISEED( J )
  528. 20 CONTINUE
  529. *
  530. * Compute "A"
  531. *
  532. IF( MTYPES.GT.MAXTYP )
  533. $ GO TO 30
  534. *
  535. IF( JTYPE.EQ.1 ) THEN
  536. *
  537. * Zero matrix
  538. *
  539. CALL DLASET( 'Full', M, N, ZERO, ZERO, A, LDA )
  540. *
  541. ELSE IF( JTYPE.EQ.2 ) THEN
  542. *
  543. * Identity matrix
  544. *
  545. CALL DLASET( 'Full', M, N, ZERO, ONE, A, LDA )
  546. *
  547. ELSE
  548. *
  549. * (Scaled) random matrix
  550. *
  551. IF( JTYPE.EQ.3 )
  552. $ ANORM = ONE
  553. IF( JTYPE.EQ.4 )
  554. $ ANORM = UNFL / ULP
  555. IF( JTYPE.EQ.5 )
  556. $ ANORM = OVFL*ULP
  557. CALL DLATMS( M, N, 'U', ISEED, 'N', S, 4, DBLE( MNMIN ),
  558. $ ANORM, M-1, N-1, 'N', A, LDA, WORK, IINFO )
  559. IF( IINFO.NE.0 ) THEN
  560. WRITE( NOUT, FMT = 9996 )'Generator', IINFO, M, N,
  561. $ JTYPE, IOLDSD
  562. INFO = ABS( IINFO )
  563. RETURN
  564. END IF
  565. END IF
  566. *
  567. 30 CONTINUE
  568. CALL DLACPY( 'F', M, N, A, LDA, ASAV, LDA )
  569. *
  570. * Do for minimal and adequate (for blocking) workspace
  571. *
  572. DO 220 IWS = 1, 4
  573. *
  574. DO 40 J = 1, 32
  575. RESULT( J ) = -ONE
  576. 40 CONTINUE
  577. *
  578. * Test DGESVD: Factorize A
  579. *
  580. IWTMP = MAX( 3*MIN( M, N )+MAX( M, N ), 5*MIN( M, N ) )
  581. LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
  582. LSWORK = MIN( LSWORK, LWORK )
  583. LSWORK = MAX( LSWORK, 1 )
  584. IF( IWS.EQ.4 )
  585. $ LSWORK = LWORK
  586. *
  587. IF( IWS.GT.1 )
  588. $ CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
  589. SRNAMT = 'DGESVD'
  590. CALL DGESVD( 'A', 'A', M, N, A, LDA, SSAV, USAV, LDU,
  591. $ VTSAV, LDVT, WORK, LSWORK, IINFO )
  592. IF( IINFO.NE.0 ) THEN
  593. WRITE( NOUT, FMT = 9995 )'GESVD', IINFO, M, N, JTYPE,
  594. $ LSWORK, IOLDSD
  595. INFO = ABS( IINFO )
  596. RETURN
  597. END IF
  598. *
  599. * Do tests 1--4
  600. *
  601. CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
  602. $ VTSAV, LDVT, WORK, RESULT( 1 ) )
  603. IF( M.NE.0 .AND. N.NE.0 ) THEN
  604. CALL DORT01( 'Columns', M, M, USAV, LDU, WORK, LWORK,
  605. $ RESULT( 2 ) )
  606. CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK, LWORK,
  607. $ RESULT( 3 ) )
  608. END IF
  609. RESULT( 4 ) = ZERO
  610. DO 50 I = 1, MNMIN - 1
  611. IF( SSAV( I ).LT.SSAV( I+1 ) )
  612. $ RESULT( 4 ) = ULPINV
  613. IF( SSAV( I ).LT.ZERO )
  614. $ RESULT( 4 ) = ULPINV
  615. 50 CONTINUE
  616. IF( MNMIN.GE.1 ) THEN
  617. IF( SSAV( MNMIN ).LT.ZERO )
  618. $ RESULT( 4 ) = ULPINV
  619. END IF
  620. *
  621. * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
  622. *
  623. RESULT( 5 ) = ZERO
  624. RESULT( 6 ) = ZERO
  625. RESULT( 7 ) = ZERO
  626. DO 80 IJU = 0, 3
  627. DO 70 IJVT = 0, 3
  628. IF( ( IJU.EQ.3 .AND. IJVT.EQ.3 ) .OR.
  629. $ ( IJU.EQ.1 .AND. IJVT.EQ.1 ) )GO TO 70
  630. JOBU = CJOB( IJU+1 )
  631. JOBVT = CJOB( IJVT+1 )
  632. CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
  633. SRNAMT = 'DGESVD'
  634. CALL DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
  635. $ VT, LDVT, WORK, LSWORK, IINFO )
  636. *
  637. * Compare U
  638. *
  639. DIF = ZERO
  640. IF( M.GT.0 .AND. N.GT.0 ) THEN
  641. IF( IJU.EQ.1 ) THEN
  642. CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
  643. $ LDU, A, LDA, WORK, LWORK, DIF,
  644. $ IINFO )
  645. ELSE IF( IJU.EQ.2 ) THEN
  646. CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
  647. $ LDU, U, LDU, WORK, LWORK, DIF,
  648. $ IINFO )
  649. ELSE IF( IJU.EQ.3 ) THEN
  650. CALL DORT03( 'C', M, M, M, MNMIN, USAV, LDU,
  651. $ U, LDU, WORK, LWORK, DIF,
  652. $ IINFO )
  653. END IF
  654. END IF
  655. RESULT( 5 ) = MAX( RESULT( 5 ), DIF )
  656. *
  657. * Compare VT
  658. *
  659. DIF = ZERO
  660. IF( M.GT.0 .AND. N.GT.0 ) THEN
  661. IF( IJVT.EQ.1 ) THEN
  662. CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
  663. $ LDVT, A, LDA, WORK, LWORK, DIF,
  664. $ IINFO )
  665. ELSE IF( IJVT.EQ.2 ) THEN
  666. CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
  667. $ LDVT, VT, LDVT, WORK, LWORK,
  668. $ DIF, IINFO )
  669. ELSE IF( IJVT.EQ.3 ) THEN
  670. CALL DORT03( 'R', N, N, N, MNMIN, VTSAV,
  671. $ LDVT, VT, LDVT, WORK, LWORK,
  672. $ DIF, IINFO )
  673. END IF
  674. END IF
  675. RESULT( 6 ) = MAX( RESULT( 6 ), DIF )
  676. *
  677. * Compare S
  678. *
  679. DIF = ZERO
  680. DIV = MAX( MNMIN*ULP*S( 1 ), UNFL )
  681. DO 60 I = 1, MNMIN - 1
  682. IF( SSAV( I ).LT.SSAV( I+1 ) )
  683. $ DIF = ULPINV
  684. IF( SSAV( I ).LT.ZERO )
  685. $ DIF = ULPINV
  686. DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
  687. 60 CONTINUE
  688. RESULT( 7 ) = MAX( RESULT( 7 ), DIF )
  689. 70 CONTINUE
  690. 80 CONTINUE
  691. *
  692. * Test DGESDD: Factorize A
  693. *
  694. IWTMP = 5*MNMIN*MNMIN + 9*MNMIN + MAX( M, N )
  695. LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
  696. LSWORK = MIN( LSWORK, LWORK )
  697. LSWORK = MAX( LSWORK, 1 )
  698. IF( IWS.EQ.4 )
  699. $ LSWORK = LWORK
  700. *
  701. CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
  702. SRNAMT = 'DGESDD'
  703. CALL DGESDD( 'A', M, N, A, LDA, SSAV, USAV, LDU, VTSAV,
  704. $ LDVT, WORK, LSWORK, IWORK, IINFO )
  705. IF( IINFO.NE.0 ) THEN
  706. WRITE( NOUT, FMT = 9995 )'GESDD', IINFO, M, N, JTYPE,
  707. $ LSWORK, IOLDSD
  708. INFO = ABS( IINFO )
  709. RETURN
  710. END IF
  711. *
  712. * Do tests 8--11
  713. *
  714. CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
  715. $ VTSAV, LDVT, WORK, RESULT( 8 ) )
  716. IF( M.NE.0 .AND. N.NE.0 ) THEN
  717. CALL DORT01( 'Columns', M, M, USAV, LDU, WORK, LWORK,
  718. $ RESULT( 9 ) )
  719. CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK, LWORK,
  720. $ RESULT( 10 ) )
  721. END IF
  722. RESULT( 11 ) = ZERO
  723. DO 90 I = 1, MNMIN - 1
  724. IF( SSAV( I ).LT.SSAV( I+1 ) )
  725. $ RESULT( 11 ) = ULPINV
  726. IF( SSAV( I ).LT.ZERO )
  727. $ RESULT( 11 ) = ULPINV
  728. 90 CONTINUE
  729. IF( MNMIN.GE.1 ) THEN
  730. IF( SSAV( MNMIN ).LT.ZERO )
  731. $ RESULT( 11 ) = ULPINV
  732. END IF
  733. *
  734. * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
  735. *
  736. RESULT( 12 ) = ZERO
  737. RESULT( 13 ) = ZERO
  738. RESULT( 14 ) = ZERO
  739. DO 110 IJQ = 0, 2
  740. JOBQ = CJOB( IJQ+1 )
  741. CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
  742. SRNAMT = 'DGESDD'
  743. CALL DGESDD( JOBQ, M, N, A, LDA, S, U, LDU, VT, LDVT,
  744. $ WORK, LSWORK, IWORK, IINFO )
  745. *
  746. * Compare U
  747. *
  748. DIF = ZERO
  749. IF( M.GT.0 .AND. N.GT.0 ) THEN
  750. IF( IJQ.EQ.1 ) THEN
  751. IF( M.GE.N ) THEN
  752. CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
  753. $ LDU, A, LDA, WORK, LWORK, DIF,
  754. $ INFO )
  755. ELSE
  756. CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
  757. $ LDU, U, LDU, WORK, LWORK, DIF,
  758. $ INFO )
  759. END IF
  760. ELSE IF( IJQ.EQ.2 ) THEN
  761. CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV, LDU,
  762. $ U, LDU, WORK, LWORK, DIF, INFO )
  763. END IF
  764. END IF
  765. RESULT( 12 ) = MAX( RESULT( 12 ), DIF )
  766. *
  767. * Compare VT
  768. *
  769. DIF = ZERO
  770. IF( M.GT.0 .AND. N.GT.0 ) THEN
  771. IF( IJQ.EQ.1 ) THEN
  772. IF( M.GE.N ) THEN
  773. CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
  774. $ LDVT, VT, LDVT, WORK, LWORK,
  775. $ DIF, INFO )
  776. ELSE
  777. CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
  778. $ LDVT, A, LDA, WORK, LWORK, DIF,
  779. $ INFO )
  780. END IF
  781. ELSE IF( IJQ.EQ.2 ) THEN
  782. CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
  783. $ LDVT, VT, LDVT, WORK, LWORK, DIF,
  784. $ INFO )
  785. END IF
  786. END IF
  787. RESULT( 13 ) = MAX( RESULT( 13 ), DIF )
  788. *
  789. * Compare S
  790. *
  791. DIF = ZERO
  792. DIV = MAX( MNMIN*ULP*S( 1 ), UNFL )
  793. DO 100 I = 1, MNMIN - 1
  794. IF( SSAV( I ).LT.SSAV( I+1 ) )
  795. $ DIF = ULPINV
  796. IF( SSAV( I ).LT.ZERO )
  797. $ DIF = ULPINV
  798. DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
  799. 100 CONTINUE
  800. RESULT( 14 ) = MAX( RESULT( 14 ), DIF )
  801. 110 CONTINUE
  802. *
  803. * Test DGESVDQ
  804. * Note: DGESVDQ only works for M >= N
  805. *
  806. RESULT( 36 ) = ZERO
  807. RESULT( 37 ) = ZERO
  808. RESULT( 38 ) = ZERO
  809. RESULT( 39 ) = ZERO
  810. *
  811. IF( M.GE.N ) THEN
  812. IWTMP = 5*MNMIN*MNMIN + 9*MNMIN + MAX( M, N )
  813. LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
  814. LSWORK = MIN( LSWORK, LWORK )
  815. LSWORK = MAX( LSWORK, 1 )
  816. IF( IWS.EQ.4 )
  817. $ LSWORK = LWORK
  818. *
  819. CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
  820. SRNAMT = 'DGESVDQ'
  821. *
  822. LRWORK = 2
  823. LIWORK = MAX( N, 1 )
  824. CALL DGESVDQ( 'H', 'N', 'N', 'A', 'A',
  825. $ M, N, A, LDA, SSAV, USAV, LDU,
  826. $ VTSAV, LDVT, NUMRANK, IWORK, LIWORK,
  827. $ WORK, LWORK, RWORK, LRWORK, IINFO )
  828. *
  829. IF( IINFO.NE.0 ) THEN
  830. WRITE( NOUT, FMT = 9995 )'DGESVDQ', IINFO, M, N,
  831. $ JTYPE, LSWORK, IOLDSD
  832. INFO = ABS( IINFO )
  833. RETURN
  834. END IF
  835. *
  836. * Do tests 36--39
  837. *
  838. CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
  839. $ VTSAV, LDVT, WORK, RESULT( 36 ) )
  840. IF( M.NE.0 .AND. N.NE.0 ) THEN
  841. CALL DORT01( 'Columns', M, M, USAV, LDU, WORK,
  842. $ LWORK, RESULT( 37 ) )
  843. CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK,
  844. $ LWORK, RESULT( 38 ) )
  845. END IF
  846. RESULT( 39 ) = ZERO
  847. DO 199 I = 1, MNMIN - 1
  848. IF( SSAV( I ).LT.SSAV( I+1 ) )
  849. $ RESULT( 39 ) = ULPINV
  850. IF( SSAV( I ).LT.ZERO )
  851. $ RESULT( 39 ) = ULPINV
  852. 199 CONTINUE
  853. IF( MNMIN.GE.1 ) THEN
  854. IF( SSAV( MNMIN ).LT.ZERO )
  855. $ RESULT( 39 ) = ULPINV
  856. END IF
  857. END IF
  858. *
  859. * Test DGESVJ
  860. * Note: DGESVJ only works for M >= N
  861. *
  862. RESULT( 15 ) = ZERO
  863. RESULT( 16 ) = ZERO
  864. RESULT( 17 ) = ZERO
  865. RESULT( 18 ) = ZERO
  866. *
  867. IF( M.GE.N ) THEN
  868. IWTMP = 5*MNMIN*MNMIN + 9*MNMIN + MAX( M, N )
  869. LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
  870. LSWORK = MIN( LSWORK, LWORK )
  871. LSWORK = MAX( LSWORK, 1 )
  872. IF( IWS.EQ.4 )
  873. $ LSWORK = LWORK
  874. *
  875. CALL DLACPY( 'F', M, N, ASAV, LDA, USAV, LDA )
  876. SRNAMT = 'DGESVJ'
  877. CALL DGESVJ( 'G', 'U', 'V', M, N, USAV, LDA, SSAV,
  878. & 0, A, LDVT, WORK, LWORK, INFO )
  879. *
  880. * DGESVJ returns V not VT
  881. *
  882. DO J=1,N
  883. DO I=1,N
  884. VTSAV(J,I) = A(I,J)
  885. END DO
  886. END DO
  887. *
  888. IF( IINFO.NE.0 ) THEN
  889. WRITE( NOUT, FMT = 9995 )'GESVJ', IINFO, M, N,
  890. $ JTYPE, LSWORK, IOLDSD
  891. INFO = ABS( IINFO )
  892. RETURN
  893. END IF
  894. *
  895. * Do tests 15--18
  896. *
  897. CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
  898. $ VTSAV, LDVT, WORK, RESULT( 15 ) )
  899. IF( M.NE.0 .AND. N.NE.0 ) THEN
  900. CALL DORT01( 'Columns', M, M, USAV, LDU, WORK,
  901. $ LWORK, RESULT( 16 ) )
  902. CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK,
  903. $ LWORK, RESULT( 17 ) )
  904. END IF
  905. RESULT( 18 ) = ZERO
  906. DO 120 I = 1, MNMIN - 1
  907. IF( SSAV( I ).LT.SSAV( I+1 ) )
  908. $ RESULT( 18 ) = ULPINV
  909. IF( SSAV( I ).LT.ZERO )
  910. $ RESULT( 18 ) = ULPINV
  911. 120 CONTINUE
  912. IF( MNMIN.GE.1 ) THEN
  913. IF( SSAV( MNMIN ).LT.ZERO )
  914. $ RESULT( 18 ) = ULPINV
  915. END IF
  916. END IF
  917. *
  918. * Test DGEJSV
  919. * Note: DGEJSV only works for M >= N
  920. *
  921. RESULT( 19 ) = ZERO
  922. RESULT( 20 ) = ZERO
  923. RESULT( 21 ) = ZERO
  924. RESULT( 22 ) = ZERO
  925. IF( M.GE.N ) THEN
  926. IWTMP = 5*MNMIN*MNMIN + 9*MNMIN + MAX( M, N )
  927. LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
  928. LSWORK = MIN( LSWORK, LWORK )
  929. LSWORK = MAX( LSWORK, 1 )
  930. IF( IWS.EQ.4 )
  931. $ LSWORK = LWORK
  932. *
  933. CALL DLACPY( 'F', M, N, ASAV, LDA, VTSAV, LDA )
  934. SRNAMT = 'DGEJSV'
  935. CALL DGEJSV( 'G', 'U', 'V', 'R', 'N', 'N',
  936. & M, N, VTSAV, LDA, SSAV, USAV, LDU, A, LDVT,
  937. & WORK, LWORK, IWORK, INFO )
  938. *
  939. * DGEJSV returns V not VT
  940. *
  941. DO 140 J=1,N
  942. DO 130 I=1,N
  943. VTSAV(J,I) = A(I,J)
  944. 130 END DO
  945. 140 END DO
  946. *
  947. IF( IINFO.NE.0 ) THEN
  948. WRITE( NOUT, FMT = 9995 )'GEJSV', IINFO, M, N,
  949. $ JTYPE, LSWORK, IOLDSD
  950. INFO = ABS( IINFO )
  951. RETURN
  952. END IF
  953. *
  954. * Do tests 19--22
  955. *
  956. CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
  957. $ VTSAV, LDVT, WORK, RESULT( 19 ) )
  958. IF( M.NE.0 .AND. N.NE.0 ) THEN
  959. CALL DORT01( 'Columns', M, M, USAV, LDU, WORK,
  960. $ LWORK, RESULT( 20 ) )
  961. CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK,
  962. $ LWORK, RESULT( 21 ) )
  963. END IF
  964. RESULT( 22 ) = ZERO
  965. DO 150 I = 1, MNMIN - 1
  966. IF( SSAV( I ).LT.SSAV( I+1 ) )
  967. $ RESULT( 22 ) = ULPINV
  968. IF( SSAV( I ).LT.ZERO )
  969. $ RESULT( 22 ) = ULPINV
  970. 150 CONTINUE
  971. IF( MNMIN.GE.1 ) THEN
  972. IF( SSAV( MNMIN ).LT.ZERO )
  973. $ RESULT( 22 ) = ULPINV
  974. END IF
  975. END IF
  976. *
  977. * Test DGESVDX
  978. *
  979. CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
  980. CALL DGESVDX( 'V', 'V', 'A', M, N, A, LDA,
  981. $ VL, VU, IL, IU, NS, SSAV, USAV, LDU,
  982. $ VTSAV, LDVT, WORK, LWORK, IWORK,
  983. $ IINFO )
  984. IF( IINFO.NE.0 ) THEN
  985. WRITE( NOUT, FMT = 9995 )'GESVDX', IINFO, M, N,
  986. $ JTYPE, LSWORK, IOLDSD
  987. INFO = ABS( IINFO )
  988. RETURN
  989. END IF
  990. *
  991. * Do tests 23--29
  992. *
  993. RESULT( 23 ) = ZERO
  994. RESULT( 24 ) = ZERO
  995. RESULT( 25 ) = ZERO
  996. CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
  997. $ VTSAV, LDVT, WORK, RESULT( 23 ) )
  998. IF( M.NE.0 .AND. N.NE.0 ) THEN
  999. CALL DORT01( 'Columns', M, M, USAV, LDU, WORK, LWORK,
  1000. $ RESULT( 24 ) )
  1001. CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK, LWORK,
  1002. $ RESULT( 25 ) )
  1003. END IF
  1004. RESULT( 26 ) = ZERO
  1005. DO 160 I = 1, MNMIN - 1
  1006. IF( SSAV( I ).LT.SSAV( I+1 ) )
  1007. $ RESULT( 26 ) = ULPINV
  1008. IF( SSAV( I ).LT.ZERO )
  1009. $ RESULT( 26 ) = ULPINV
  1010. 160 CONTINUE
  1011. IF( MNMIN.GE.1 ) THEN
  1012. IF( SSAV( MNMIN ).LT.ZERO )
  1013. $ RESULT( 26 ) = ULPINV
  1014. END IF
  1015. *
  1016. * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
  1017. *
  1018. RESULT( 27 ) = ZERO
  1019. RESULT( 28 ) = ZERO
  1020. RESULT( 29 ) = ZERO
  1021. DO 180 IJU = 0, 1
  1022. DO 170 IJVT = 0, 1
  1023. IF( ( IJU.EQ.0 .AND. IJVT.EQ.0 ) .OR.
  1024. $ ( IJU.EQ.1 .AND. IJVT.EQ.1 ) )GO TO 170
  1025. JOBU = CJOBV( IJU+1 )
  1026. JOBVT = CJOBV( IJVT+1 )
  1027. RANGE = CJOBR( 1 )
  1028. CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
  1029. CALL DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA,
  1030. $ VL, VU, IL, IU, NS, S, U, LDU,
  1031. $ VT, LDVT, WORK, LWORK, IWORK,
  1032. $ IINFO )
  1033. *
  1034. * Compare U
  1035. *
  1036. DIF = ZERO
  1037. IF( M.GT.0 .AND. N.GT.0 ) THEN
  1038. IF( IJU.EQ.1 ) THEN
  1039. CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
  1040. $ LDU, U, LDU, WORK, LWORK, DIF,
  1041. $ IINFO )
  1042. END IF
  1043. END IF
  1044. RESULT( 27 ) = MAX( RESULT( 27 ), DIF )
  1045. *
  1046. * Compare VT
  1047. *
  1048. DIF = ZERO
  1049. IF( M.GT.0 .AND. N.GT.0 ) THEN
  1050. IF( IJVT.EQ.1 ) THEN
  1051. CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
  1052. $ LDVT, VT, LDVT, WORK, LWORK,
  1053. $ DIF, IINFO )
  1054. END IF
  1055. END IF
  1056. RESULT( 28 ) = MAX( RESULT( 28 ), DIF )
  1057. *
  1058. * Compare S
  1059. *
  1060. DIF = ZERO
  1061. DIV = MAX( MNMIN*ULP*S( 1 ), UNFL )
  1062. DO 190 I = 1, MNMIN - 1
  1063. IF( SSAV( I ).LT.SSAV( I+1 ) )
  1064. $ DIF = ULPINV
  1065. IF( SSAV( I ).LT.ZERO )
  1066. $ DIF = ULPINV
  1067. DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
  1068. 190 CONTINUE
  1069. RESULT( 29 ) = MAX( RESULT( 29 ), DIF )
  1070. 170 CONTINUE
  1071. 180 CONTINUE
  1072. *
  1073. * Do tests 30--32: DGESVDX( 'V', 'V', 'I' )
  1074. *
  1075. DO 200 I = 1, 4
  1076. ISEED2( I ) = ISEED( I )
  1077. 200 CONTINUE
  1078. IF( MNMIN.LE.1 ) THEN
  1079. IL = 1
  1080. IU = MAX( 1, MNMIN )
  1081. ELSE
  1082. IL = 1 + INT( ( MNMIN-1 )*DLARND( 1, ISEED2 ) )
  1083. IU = 1 + INT( ( MNMIN-1 )*DLARND( 1, ISEED2 ) )
  1084. IF( IU.LT.IL ) THEN
  1085. ITEMP = IU
  1086. IU = IL
  1087. IL = ITEMP
  1088. END IF
  1089. END IF
  1090. CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
  1091. CALL DGESVDX( 'V', 'V', 'I', M, N, A, LDA,
  1092. $ VL, VU, IL, IU, NSI, S, U, LDU,
  1093. $ VT, LDVT, WORK, LWORK, IWORK,
  1094. $ IINFO )
  1095. IF( IINFO.NE.0 ) THEN
  1096. WRITE( NOUT, FMT = 9995 )'GESVDX', IINFO, M, N,
  1097. $ JTYPE, LSWORK, IOLDSD
  1098. INFO = ABS( IINFO )
  1099. RETURN
  1100. END IF
  1101. *
  1102. RESULT( 30 ) = ZERO
  1103. RESULT( 31 ) = ZERO
  1104. RESULT( 32 ) = ZERO
  1105. CALL DBDT05( M, N, ASAV, LDA, S, NSI, U, LDU,
  1106. $ VT, LDVT, WORK, RESULT( 30 ) )
  1107. CALL DORT01( 'Columns', M, NSI, U, LDU, WORK, LWORK,
  1108. $ RESULT( 31 ) )
  1109. CALL DORT01( 'Rows', NSI, N, VT, LDVT, WORK, LWORK,
  1110. $ RESULT( 32 ) )
  1111. *
  1112. * Do tests 33--35: DGESVDX( 'V', 'V', 'V' )
  1113. *
  1114. IF( MNMIN.GT.0 .AND. NSI.GT.1 ) THEN
  1115. IF( IL.NE.1 ) THEN
  1116. VU = SSAV( IL ) +
  1117. $ MAX( HALF*ABS( SSAV( IL )-SSAV( IL-1 ) ),
  1118. $ ULP*ANORM, TWO*RTUNFL )
  1119. ELSE
  1120. VU = SSAV( 1 ) +
  1121. $ MAX( HALF*ABS( SSAV( NS )-SSAV( 1 ) ),
  1122. $ ULP*ANORM, TWO*RTUNFL )
  1123. END IF
  1124. IF( IU.NE.NS ) THEN
  1125. VL = SSAV( IU ) - MAX( ULP*ANORM, TWO*RTUNFL,
  1126. $ HALF*ABS( SSAV( IU+1 )-SSAV( IU ) ) )
  1127. ELSE
  1128. VL = SSAV( NS ) - MAX( ULP*ANORM, TWO*RTUNFL,
  1129. $ HALF*ABS( SSAV( NS )-SSAV( 1 ) ) )
  1130. END IF
  1131. VL = MAX( VL,ZERO )
  1132. VU = MAX( VU,ZERO )
  1133. IF( VL.GE.VU ) VU = MAX( VU*2, VU+VL+HALF )
  1134. ELSE
  1135. VL = ZERO
  1136. VU = ONE
  1137. END IF
  1138. CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
  1139. CALL DGESVDX( 'V', 'V', 'V', M, N, A, LDA,
  1140. $ VL, VU, IL, IU, NSV, S, U, LDU,
  1141. $ VT, LDVT, WORK, LWORK, IWORK,
  1142. $ IINFO )
  1143. IF( IINFO.NE.0 ) THEN
  1144. WRITE( NOUT, FMT = 9995 )'GESVDX', IINFO, M, N,
  1145. $ JTYPE, LSWORK, IOLDSD
  1146. INFO = ABS( IINFO )
  1147. RETURN
  1148. END IF
  1149. *
  1150. RESULT( 33 ) = ZERO
  1151. RESULT( 34 ) = ZERO
  1152. RESULT( 35 ) = ZERO
  1153. CALL DBDT05( M, N, ASAV, LDA, S, NSV, U, LDU,
  1154. $ VT, LDVT, WORK, RESULT( 33 ) )
  1155. CALL DORT01( 'Columns', M, NSV, U, LDU, WORK, LWORK,
  1156. $ RESULT( 34 ) )
  1157. CALL DORT01( 'Rows', NSV, N, VT, LDVT, WORK, LWORK,
  1158. $ RESULT( 35 ) )
  1159. *
  1160. * End of Loop -- Check for RESULT(j) > THRESH
  1161. *
  1162. DO 210 J = 1, 39
  1163. IF( RESULT( J ).GE.THRESH ) THEN
  1164. IF( NFAIL.EQ.0 ) THEN
  1165. WRITE( NOUT, FMT = 9999 )
  1166. WRITE( NOUT, FMT = 9998 )
  1167. END IF
  1168. WRITE( NOUT, FMT = 9997 )M, N, JTYPE, IWS, IOLDSD,
  1169. $ J, RESULT( J )
  1170. NFAIL = NFAIL + 1
  1171. END IF
  1172. 210 CONTINUE
  1173. NTEST = NTEST + 39
  1174. 220 CONTINUE
  1175. 230 CONTINUE
  1176. 240 CONTINUE
  1177. *
  1178. * Summary
  1179. *
  1180. CALL ALASVM( PATH, NOUT, NFAIL, NTEST, 0 )
  1181. *
  1182. 9999 FORMAT( ' SVD -- Real Singular Value Decomposition Driver ',
  1183. $ / ' Matrix types (see DDRVBD for details):',
  1184. $ / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
  1185. $ / ' 3 = Evenly spaced singular values near 1',
  1186. $ / ' 4 = Evenly spaced singular values near underflow',
  1187. $ / ' 5 = Evenly spaced singular values near overflow', / /
  1188. $ ' Tests performed: ( A is dense, U and V are orthogonal,',
  1189. $ / 19X, ' S is an array, and Upartial, VTpartial, and',
  1190. $ / 19X, ' Spartial are partially computed U, VT and S),', / )
  1191. 9998 FORMAT( ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
  1192. $ / ' 2 = | I - U**T U | / ( M ulp ) ',
  1193. $ / ' 3 = | I - VT VT**T | / ( N ulp ) ',
  1194. $ / ' 4 = 0 if S contains min(M,N) nonnegative values in',
  1195. $ ' decreasing order, else 1/ulp',
  1196. $ / ' 5 = | U - Upartial | / ( M ulp )',
  1197. $ / ' 6 = | VT - VTpartial | / ( N ulp )',
  1198. $ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
  1199. $ / ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
  1200. $ / ' 9 = | I - U**T U | / ( M ulp ) ',
  1201. $ / '10 = | I - VT VT**T | / ( N ulp ) ',
  1202. $ / '11 = 0 if S contains min(M,N) nonnegative values in',
  1203. $ ' decreasing order, else 1/ulp',
  1204. $ / '12 = | U - Upartial | / ( M ulp )',
  1205. $ / '13 = | VT - VTpartial | / ( N ulp )',
  1206. $ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )',
  1207. $ / '15 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
  1208. $ / '16 = | I - U**T U | / ( M ulp ) ',
  1209. $ / '17 = | I - VT VT**T | / ( N ulp ) ',
  1210. $ / '18 = 0 if S contains min(M,N) nonnegative values in',
  1211. $ ' decreasing order, else 1/ulp',
  1212. $ / '19 = | U - Upartial | / ( M ulp )',
  1213. $ / '20 = | VT - VTpartial | / ( N ulp )',
  1214. $ / '21 = | S - Spartial | / ( min(M,N) ulp |S| )',
  1215. $ / '22 = 0 if S contains min(M,N) nonnegative values in',
  1216. $ ' decreasing order, else 1/ulp',
  1217. $ / '23 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ),',
  1218. $ ' DGESVDX(V,V,A) ',
  1219. $ / '24 = | I - U**T U | / ( M ulp ) ',
  1220. $ / '25 = | I - VT VT**T | / ( N ulp ) ',
  1221. $ / '26 = 0 if S contains min(M,N) nonnegative values in',
  1222. $ ' decreasing order, else 1/ulp',
  1223. $ / '27 = | U - Upartial | / ( M ulp )',
  1224. $ / '28 = | VT - VTpartial | / ( N ulp )',
  1225. $ / '29 = | S - Spartial | / ( min(M,N) ulp |S| )',
  1226. $ / '30 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp ),',
  1227. $ ' DGESVDX(V,V,I) ',
  1228. $ / '31 = | I - U**T U | / ( M ulp ) ',
  1229. $ / '32 = | I - VT VT**T | / ( N ulp ) ',
  1230. $ / '33 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp ),',
  1231. $ ' DGESVDX(V,V,V) ',
  1232. $ / '34 = | I - U**T U | / ( M ulp ) ',
  1233. $ / '35 = | I - VT VT**T | / ( N ulp ) ',
  1234. $ ' DGESVDQ(H,N,N,A,A',
  1235. $ / '36 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
  1236. $ / '37 = | I - U**T U | / ( M ulp ) ',
  1237. $ / '38 = | I - VT VT**T | / ( N ulp ) ',
  1238. $ / '39 = 0 if S contains min(M,N) nonnegative values in',
  1239. $ ' decreasing order, else 1/ulp',
  1240. $ / / )
  1241. 9997 FORMAT( ' M=', I5, ', N=', I5, ', type ', I1, ', IWS=', I1,
  1242. $ ', seed=', 4( I4, ',' ), ' test(', I2, ')=', G11.4 )
  1243. 9996 FORMAT( ' DDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
  1244. $ I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
  1245. $ I5, ')' )
  1246. 9995 FORMAT( ' DDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
  1247. $ I6, ', N=', I6, ', JTYPE=', I6, ', LSWORK=', I6, / 9X,
  1248. $ 'ISEED=(', 3( I5, ',' ), I5, ')' )
  1249. *
  1250. RETURN
  1251. *
  1252. * End of DDRVBD
  1253. *
  1254. END