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.

cdrvbd.f 49 kB

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