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

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