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.

clatms.c 58 kB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969
  1. #include <math.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <stdio.h>
  5. #include <complex.h>
  6. #ifdef complex
  7. #undef complex
  8. #endif
  9. #ifdef I
  10. #undef I
  11. #endif
  12. #if defined(_WIN64)
  13. typedef long long BLASLONG;
  14. typedef unsigned long long BLASULONG;
  15. #else
  16. typedef long BLASLONG;
  17. typedef unsigned long BLASULONG;
  18. #endif
  19. #ifdef LAPACK_ILP64
  20. typedef BLASLONG blasint;
  21. #if defined(_WIN64)
  22. #define blasabs(x) llabs(x)
  23. #else
  24. #define blasabs(x) labs(x)
  25. #endif
  26. #else
  27. typedef int blasint;
  28. #define blasabs(x) abs(x)
  29. #endif
  30. typedef blasint integer;
  31. typedef unsigned int uinteger;
  32. typedef char *address;
  33. typedef short int shortint;
  34. typedef float real;
  35. typedef double doublereal;
  36. typedef struct { real r, i; } complex;
  37. typedef struct { doublereal r, i; } doublecomplex;
  38. #ifdef _MSC_VER
  39. static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
  40. static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
  41. static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
  42. static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
  43. #else
  44. static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
  45. static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
  46. static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
  47. static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
  48. #endif
  49. #define pCf(z) (*_pCf(z))
  50. #define pCd(z) (*_pCd(z))
  51. typedef blasint logical;
  52. typedef char logical1;
  53. typedef char integer1;
  54. #define TRUE_ (1)
  55. #define FALSE_ (0)
  56. /* Extern is for use with -E */
  57. #ifndef Extern
  58. #define Extern extern
  59. #endif
  60. /* I/O stuff */
  61. typedef int flag;
  62. typedef int ftnlen;
  63. typedef int ftnint;
  64. /*external read, write*/
  65. typedef struct
  66. { flag cierr;
  67. ftnint ciunit;
  68. flag ciend;
  69. char *cifmt;
  70. ftnint cirec;
  71. } cilist;
  72. /*internal read, write*/
  73. typedef struct
  74. { flag icierr;
  75. char *iciunit;
  76. flag iciend;
  77. char *icifmt;
  78. ftnint icirlen;
  79. ftnint icirnum;
  80. } icilist;
  81. /*open*/
  82. typedef struct
  83. { flag oerr;
  84. ftnint ounit;
  85. char *ofnm;
  86. ftnlen ofnmlen;
  87. char *osta;
  88. char *oacc;
  89. char *ofm;
  90. ftnint orl;
  91. char *oblnk;
  92. } olist;
  93. /*close*/
  94. typedef struct
  95. { flag cerr;
  96. ftnint cunit;
  97. char *csta;
  98. } cllist;
  99. /*rewind, backspace, endfile*/
  100. typedef struct
  101. { flag aerr;
  102. ftnint aunit;
  103. } alist;
  104. /* inquire */
  105. typedef struct
  106. { flag inerr;
  107. ftnint inunit;
  108. char *infile;
  109. ftnlen infilen;
  110. ftnint *inex; /*parameters in standard's order*/
  111. ftnint *inopen;
  112. ftnint *innum;
  113. ftnint *innamed;
  114. char *inname;
  115. ftnlen innamlen;
  116. char *inacc;
  117. ftnlen inacclen;
  118. char *inseq;
  119. ftnlen inseqlen;
  120. char *indir;
  121. ftnlen indirlen;
  122. char *infmt;
  123. ftnlen infmtlen;
  124. char *inform;
  125. ftnint informlen;
  126. char *inunf;
  127. ftnlen inunflen;
  128. ftnint *inrecl;
  129. ftnint *innrec;
  130. char *inblank;
  131. ftnlen inblanklen;
  132. } inlist;
  133. #define VOID void
  134. union Multitype { /* for multiple entry points */
  135. integer1 g;
  136. shortint h;
  137. integer i;
  138. /* longint j; */
  139. real r;
  140. doublereal d;
  141. complex c;
  142. doublecomplex z;
  143. };
  144. typedef union Multitype Multitype;
  145. struct Vardesc { /* for Namelist */
  146. char *name;
  147. char *addr;
  148. ftnlen *dims;
  149. int type;
  150. };
  151. typedef struct Vardesc Vardesc;
  152. struct Namelist {
  153. char *name;
  154. Vardesc **vars;
  155. int nvars;
  156. };
  157. typedef struct Namelist Namelist;
  158. #define abs(x) ((x) >= 0 ? (x) : -(x))
  159. #define dabs(x) (fabs(x))
  160. #define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
  161. #define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
  162. #define dmin(a,b) (f2cmin(a,b))
  163. #define dmax(a,b) (f2cmax(a,b))
  164. #define bit_test(a,b) ((a) >> (b) & 1)
  165. #define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
  166. #define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
  167. #define abort_() { sig_die("Fortran abort routine called", 1); }
  168. #define c_abs(z) (cabsf(Cf(z)))
  169. #define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
  170. #ifdef _MSC_VER
  171. #define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
  172. #define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/df(b)._Val[1]);}
  173. #else
  174. #define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
  175. #define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
  176. #endif
  177. #define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
  178. #define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
  179. #define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
  180. //#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
  181. #define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
  182. #define d_abs(x) (fabs(*(x)))
  183. #define d_acos(x) (acos(*(x)))
  184. #define d_asin(x) (asin(*(x)))
  185. #define d_atan(x) (atan(*(x)))
  186. #define d_atn2(x, y) (atan2(*(x),*(y)))
  187. #define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
  188. #define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
  189. #define d_cos(x) (cos(*(x)))
  190. #define d_cosh(x) (cosh(*(x)))
  191. #define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
  192. #define d_exp(x) (exp(*(x)))
  193. #define d_imag(z) (cimag(Cd(z)))
  194. #define r_imag(z) (cimagf(Cf(z)))
  195. #define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
  196. #define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
  197. #define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
  198. #define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
  199. #define d_log(x) (log(*(x)))
  200. #define d_mod(x, y) (fmod(*(x), *(y)))
  201. #define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
  202. #define d_nint(x) u_nint(*(x))
  203. #define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
  204. #define d_sign(a,b) u_sign(*(a),*(b))
  205. #define r_sign(a,b) u_sign(*(a),*(b))
  206. #define d_sin(x) (sin(*(x)))
  207. #define d_sinh(x) (sinh(*(x)))
  208. #define d_sqrt(x) (sqrt(*(x)))
  209. #define d_tan(x) (tan(*(x)))
  210. #define d_tanh(x) (tanh(*(x)))
  211. #define i_abs(x) abs(*(x))
  212. #define i_dnnt(x) ((integer)u_nint(*(x)))
  213. #define i_len(s, n) (n)
  214. #define i_nint(x) ((integer)u_nint(*(x)))
  215. #define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
  216. #define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
  217. #define pow_si(B,E) spow_ui(*(B),*(E))
  218. #define pow_ri(B,E) spow_ui(*(B),*(E))
  219. #define pow_di(B,E) dpow_ui(*(B),*(E))
  220. #define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
  221. #define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
  222. #define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
  223. #define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
  224. #define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
  225. #define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
  226. #define sig_die(s, kill) { exit(1); }
  227. #define s_stop(s, n) {exit(0);}
  228. #define z_abs(z) (cabs(Cd(z)))
  229. #define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
  230. #define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
  231. #define myexit_() break;
  232. #define mycycle() continue;
  233. #define myceiling(w) {ceil(w)}
  234. #define myhuge(w) {HUGE_VAL}
  235. //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);}
  236. #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)}
  237. /* procedure parameter types for -A and -C++ */
  238. /* Table of constant values */
  239. static complex c_b1 = {0.f,0.f};
  240. static integer c__1 = 1;
  241. static integer c__5 = 5;
  242. static logical c_true = TRUE_;
  243. static logical c_false = FALSE_;
  244. /* > \brief \b CLATMS */
  245. /* =========== DOCUMENTATION =========== */
  246. /* Online html documentation available at */
  247. /* http://www.netlib.org/lapack/explore-html/ */
  248. /* Definition: */
  249. /* =========== */
  250. /* SUBROUTINE CLATMS( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, */
  251. /* KL, KU, PACK, A, LDA, WORK, INFO ) */
  252. /* CHARACTER DIST, PACK, SYM */
  253. /* INTEGER INFO, KL, KU, LDA, M, MODE, N */
  254. /* REAL COND, DMAX */
  255. /* INTEGER ISEED( 4 ) */
  256. /* REAL D( * ) */
  257. /* COMPLEX A( LDA, * ), WORK( * ) */
  258. /* > \par Purpose: */
  259. /* ============= */
  260. /* > */
  261. /* > \verbatim */
  262. /* > */
  263. /* > CLATMS generates random matrices with specified singular values */
  264. /* > (or hermitian with specified eigenvalues) */
  265. /* > for testing LAPACK programs. */
  266. /* > */
  267. /* > CLATMS operates by applying the following sequence of */
  268. /* > operations: */
  269. /* > */
  270. /* > Set the diagonal to D, where D may be input or */
  271. /* > computed according to MODE, COND, DMAX, and SYM */
  272. /* > as described below. */
  273. /* > */
  274. /* > Generate a matrix with the appropriate band structure, by one */
  275. /* > of two methods: */
  276. /* > */
  277. /* > Method A: */
  278. /* > Generate a dense M x N matrix by multiplying D on the left */
  279. /* > and the right by random unitary matrices, then: */
  280. /* > */
  281. /* > Reduce the bandwidth according to KL and KU, using */
  282. /* > Householder transformations. */
  283. /* > */
  284. /* > Method B: */
  285. /* > Convert the bandwidth-0 (i.e., diagonal) matrix to a */
  286. /* > bandwidth-1 matrix using Givens rotations, "chasing" */
  287. /* > out-of-band elements back, much as in QR; then convert */
  288. /* > the bandwidth-1 to a bandwidth-2 matrix, etc. Note */
  289. /* > that for reasonably small bandwidths (relative to M and */
  290. /* > N) this requires less storage, as a dense matrix is not */
  291. /* > generated. Also, for hermitian or symmetric matrices, */
  292. /* > only one triangle is generated. */
  293. /* > */
  294. /* > Method A is chosen if the bandwidth is a large fraction of the */
  295. /* > order of the matrix, and LDA is at least M (so a dense */
  296. /* > matrix can be stored.) Method B is chosen if the bandwidth */
  297. /* > is small (< 1/2 N for hermitian or symmetric, < .3 N+M for */
  298. /* > non-symmetric), or LDA is less than M and not less than the */
  299. /* > bandwidth. */
  300. /* > */
  301. /* > Pack the matrix if desired. Options specified by PACK are: */
  302. /* > no packing */
  303. /* > zero out upper half (if hermitian) */
  304. /* > zero out lower half (if hermitian) */
  305. /* > store the upper half columnwise (if hermitian or upper */
  306. /* > triangular) */
  307. /* > store the lower half columnwise (if hermitian or lower */
  308. /* > triangular) */
  309. /* > store the lower triangle in banded format (if hermitian or */
  310. /* > lower triangular) */
  311. /* > store the upper triangle in banded format (if hermitian or */
  312. /* > upper triangular) */
  313. /* > store the entire matrix in banded format */
  314. /* > If Method B is chosen, and band format is specified, then the */
  315. /* > matrix will be generated in the band format, so no repacking */
  316. /* > will be necessary. */
  317. /* > \endverbatim */
  318. /* Arguments: */
  319. /* ========== */
  320. /* > \param[in] M */
  321. /* > \verbatim */
  322. /* > M is INTEGER */
  323. /* > The number of rows of A. Not modified. */
  324. /* > \endverbatim */
  325. /* > */
  326. /* > \param[in] N */
  327. /* > \verbatim */
  328. /* > N is INTEGER */
  329. /* > The number of columns of A. N must equal M if the matrix */
  330. /* > is symmetric or hermitian (i.e., if SYM is not 'N') */
  331. /* > Not modified. */
  332. /* > \endverbatim */
  333. /* > */
  334. /* > \param[in] DIST */
  335. /* > \verbatim */
  336. /* > DIST is CHARACTER*1 */
  337. /* > On entry, DIST specifies the type of distribution to be used */
  338. /* > to generate the random eigen-/singular values. */
  339. /* > 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform ) */
  340. /* > 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric ) */
  341. /* > 'N' => NORMAL( 0, 1 ) ( 'N' for normal ) */
  342. /* > Not modified. */
  343. /* > \endverbatim */
  344. /* > */
  345. /* > \param[in,out] ISEED */
  346. /* > \verbatim */
  347. /* > ISEED is INTEGER array, dimension ( 4 ) */
  348. /* > On entry ISEED specifies the seed of the random number */
  349. /* > generator. They should lie between 0 and 4095 inclusive, */
  350. /* > and ISEED(4) should be odd. The random number generator */
  351. /* > uses a linear congruential sequence limited to small */
  352. /* > integers, and so should produce machine independent */
  353. /* > random numbers. The values of ISEED are changed on */
  354. /* > exit, and can be used in the next call to CLATMS */
  355. /* > to continue the same random number sequence. */
  356. /* > Changed on exit. */
  357. /* > \endverbatim */
  358. /* > */
  359. /* > \param[in] SYM */
  360. /* > \verbatim */
  361. /* > SYM is CHARACTER*1 */
  362. /* > If SYM='H', the generated matrix is hermitian, with */
  363. /* > eigenvalues specified by D, COND, MODE, and DMAX; they */
  364. /* > may be positive, negative, or zero. */
  365. /* > If SYM='P', the generated matrix is hermitian, with */
  366. /* > eigenvalues (= singular values) specified by D, COND, */
  367. /* > MODE, and DMAX; they will not be negative. */
  368. /* > If SYM='N', the generated matrix is nonsymmetric, with */
  369. /* > singular values specified by D, COND, MODE, and DMAX; */
  370. /* > they will not be negative. */
  371. /* > If SYM='S', the generated matrix is (complex) symmetric, */
  372. /* > with singular values specified by D, COND, MODE, and */
  373. /* > DMAX; they will not be negative. */
  374. /* > Not modified. */
  375. /* > \endverbatim */
  376. /* > */
  377. /* > \param[in,out] D */
  378. /* > \verbatim */
  379. /* > D is REAL array, dimension ( MIN( M, N ) ) */
  380. /* > This array is used to specify the singular values or */
  381. /* > eigenvalues of A (see SYM, above.) If MODE=0, then D is */
  382. /* > assumed to contain the singular/eigenvalues, otherwise */
  383. /* > they will be computed according to MODE, COND, and DMAX, */
  384. /* > and placed in D. */
  385. /* > Modified if MODE is nonzero. */
  386. /* > \endverbatim */
  387. /* > */
  388. /* > \param[in] MODE */
  389. /* > \verbatim */
  390. /* > MODE is INTEGER */
  391. /* > On entry this describes how the singular/eigenvalues are to */
  392. /* > be specified: */
  393. /* > MODE = 0 means use D as input */
  394. /* > MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND */
  395. /* > MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND */
  396. /* > MODE = 3 sets D(I)=COND**(-(I-1)/(N-1)) */
  397. /* > MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND) */
  398. /* > MODE = 5 sets D to random numbers in the range */
  399. /* > ( 1/COND , 1 ) such that their logarithms */
  400. /* > are uniformly distributed. */
  401. /* > MODE = 6 set D to random numbers from same distribution */
  402. /* > as the rest of the matrix. */
  403. /* > MODE < 0 has the same meaning as ABS(MODE), except that */
  404. /* > the order of the elements of D is reversed. */
  405. /* > Thus if MODE is positive, D has entries ranging from */
  406. /* > 1 to 1/COND, if negative, from 1/COND to 1, */
  407. /* > If SYM='H', and MODE is neither 0, 6, nor -6, then */
  408. /* > the elements of D will also be multiplied by a random */
  409. /* > sign (i.e., +1 or -1.) */
  410. /* > Not modified. */
  411. /* > \endverbatim */
  412. /* > */
  413. /* > \param[in] COND */
  414. /* > \verbatim */
  415. /* > COND is REAL */
  416. /* > On entry, this is used as described under MODE above. */
  417. /* > If used, it must be >= 1. Not modified. */
  418. /* > \endverbatim */
  419. /* > */
  420. /* > \param[in] DMAX */
  421. /* > \verbatim */
  422. /* > DMAX is REAL */
  423. /* > If MODE is neither -6, 0 nor 6, the contents of D, as */
  424. /* > computed according to MODE and COND, will be scaled by */
  425. /* > DMAX / f2cmax(abs(D(i))); thus, the maximum absolute eigen- or */
  426. /* > singular value (which is to say the norm) will be abs(DMAX). */
  427. /* > Note that DMAX need not be positive: if DMAX is negative */
  428. /* > (or zero), D will be scaled by a negative number (or zero). */
  429. /* > Not modified. */
  430. /* > \endverbatim */
  431. /* > */
  432. /* > \param[in] KL */
  433. /* > \verbatim */
  434. /* > KL is INTEGER */
  435. /* > This specifies the lower bandwidth of the matrix. For */
  436. /* > example, KL=0 implies upper triangular, KL=1 implies upper */
  437. /* > Hessenberg, and KL being at least M-1 means that the matrix */
  438. /* > has full lower bandwidth. KL must equal KU if the matrix */
  439. /* > is symmetric or hermitian. */
  440. /* > Not modified. */
  441. /* > \endverbatim */
  442. /* > */
  443. /* > \param[in] KU */
  444. /* > \verbatim */
  445. /* > KU is INTEGER */
  446. /* > This specifies the upper bandwidth of the matrix. For */
  447. /* > example, KU=0 implies lower triangular, KU=1 implies lower */
  448. /* > Hessenberg, and KU being at least N-1 means that the matrix */
  449. /* > has full upper bandwidth. KL must equal KU if the matrix */
  450. /* > is symmetric or hermitian. */
  451. /* > Not modified. */
  452. /* > \endverbatim */
  453. /* > */
  454. /* > \param[in] PACK */
  455. /* > \verbatim */
  456. /* > PACK is CHARACTER*1 */
  457. /* > This specifies packing of matrix as follows: */
  458. /* > 'N' => no packing */
  459. /* > 'U' => zero out all subdiagonal entries (if symmetric */
  460. /* > or hermitian) */
  461. /* > 'L' => zero out all superdiagonal entries (if symmetric */
  462. /* > or hermitian) */
  463. /* > 'C' => store the upper triangle columnwise (only if the */
  464. /* > matrix is symmetric, hermitian, or upper triangular) */
  465. /* > 'R' => store the lower triangle columnwise (only if the */
  466. /* > matrix is symmetric, hermitian, or lower triangular) */
  467. /* > 'B' => store the lower triangle in band storage scheme */
  468. /* > (only if the matrix is symmetric, hermitian, or */
  469. /* > lower triangular) */
  470. /* > 'Q' => store the upper triangle in band storage scheme */
  471. /* > (only if the matrix is symmetric, hermitian, or */
  472. /* > upper triangular) */
  473. /* > 'Z' => store the entire matrix in band storage scheme */
  474. /* > (pivoting can be provided for by using this */
  475. /* > option to store A in the trailing rows of */
  476. /* > the allocated storage) */
  477. /* > */
  478. /* > Using these options, the various LAPACK packed and banded */
  479. /* > storage schemes can be obtained: */
  480. /* > GB - use 'Z' */
  481. /* > PB, SB, HB, or TB - use 'B' or 'Q' */
  482. /* > PP, SP, HB, or TP - use 'C' or 'R' */
  483. /* > */
  484. /* > If two calls to CLATMS differ only in the PACK parameter, */
  485. /* > they will generate mathematically equivalent matrices. */
  486. /* > Not modified. */
  487. /* > \endverbatim */
  488. /* > */
  489. /* > \param[in,out] A */
  490. /* > \verbatim */
  491. /* > A is COMPLEX array, dimension ( LDA, N ) */
  492. /* > On exit A is the desired test matrix. A is first generated */
  493. /* > in full (unpacked) form, and then packed, if so specified */
  494. /* > by PACK. Thus, the first M elements of the first N */
  495. /* > columns will always be modified. If PACK specifies a */
  496. /* > packed or banded storage scheme, all LDA elements of the */
  497. /* > first N columns will be modified; the elements of the */
  498. /* > array which do not correspond to elements of the generated */
  499. /* > matrix are set to zero. */
  500. /* > Modified. */
  501. /* > \endverbatim */
  502. /* > */
  503. /* > \param[in] LDA */
  504. /* > \verbatim */
  505. /* > LDA is INTEGER */
  506. /* > LDA specifies the first dimension of A as declared in the */
  507. /* > calling program. If PACK='N', 'U', 'L', 'C', or 'R', then */
  508. /* > LDA must be at least M. If PACK='B' or 'Q', then LDA must */
  509. /* > be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)). */
  510. /* > If PACK='Z', LDA must be large enough to hold the packed */
  511. /* > array: MIN( KU, N-1) + MIN( KL, M-1) + 1. */
  512. /* > Not modified. */
  513. /* > \endverbatim */
  514. /* > */
  515. /* > \param[out] WORK */
  516. /* > \verbatim */
  517. /* > WORK is COMPLEX array, dimension ( 3*MAX( N, M ) ) */
  518. /* > Workspace. */
  519. /* > Modified. */
  520. /* > \endverbatim */
  521. /* > */
  522. /* > \param[out] INFO */
  523. /* > \verbatim */
  524. /* > INFO is INTEGER */
  525. /* > Error code. On exit, INFO will be set to one of the */
  526. /* > following values: */
  527. /* > 0 => normal return */
  528. /* > -1 => M negative or unequal to N and SYM='S', 'H', or 'P' */
  529. /* > -2 => N negative */
  530. /* > -3 => DIST illegal string */
  531. /* > -5 => SYM illegal string */
  532. /* > -7 => MODE not in range -6 to 6 */
  533. /* > -8 => COND less than 1.0, and MODE neither -6, 0 nor 6 */
  534. /* > -10 => KL negative */
  535. /* > -11 => KU negative, or SYM is not 'N' and KU is not equal to */
  536. /* > KL */
  537. /* > -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N'; */
  538. /* > or PACK='C' or 'Q' and SYM='N' and KL is not zero; */
  539. /* > or PACK='R' or 'B' and SYM='N' and KU is not zero; */
  540. /* > or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not */
  541. /* > N. */
  542. /* > -14 => LDA is less than M, or PACK='Z' and LDA is less than */
  543. /* > MIN(KU,N-1) + MIN(KL,M-1) + 1. */
  544. /* > 1 => Error return from SLATM1 */
  545. /* > 2 => Cannot scale to DMAX (f2cmax. sing. value is 0) */
  546. /* > 3 => Error return from CLAGGE, CLAGHE or CLAGSY */
  547. /* > \endverbatim */
  548. /* Authors: */
  549. /* ======== */
  550. /* > \author Univ. of Tennessee */
  551. /* > \author Univ. of California Berkeley */
  552. /* > \author Univ. of Colorado Denver */
  553. /* > \author NAG Ltd. */
  554. /* > \date December 2016 */
  555. /* > \ingroup complex_matgen */
  556. /* ===================================================================== */
  557. /* Subroutine */ void clatms_(integer *m, integer *n, char *dist, integer *
  558. iseed, char *sym, real *d__, integer *mode, real *cond, real *dmax__,
  559. integer *kl, integer *ku, char *pack, complex *a, integer *lda,
  560. complex *work, integer *info)
  561. {
  562. /* System generated locals */
  563. integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5, i__6;
  564. real r__1, r__2, r__3;
  565. complex q__1, q__2, q__3;
  566. logical L__1;
  567. /* Local variables */
  568. integer ilda, icol;
  569. real temp;
  570. logical csym;
  571. integer irow, isym;
  572. complex c__;
  573. integer i__, j, k;
  574. complex s;
  575. real alpha, angle;
  576. integer ipack;
  577. real realc;
  578. integer ioffg;
  579. extern logical lsame_(char *, char *);
  580. integer iinfo;
  581. extern /* Subroutine */ void sscal_(integer *, real *, real *, integer *);
  582. complex ctemp;
  583. integer idist, mnmin, iskew;
  584. complex extra, dummy;
  585. extern /* Subroutine */ void slatm1_(integer *, real *, integer *, integer
  586. *, integer *, real *, integer *, integer *);
  587. integer ic, jc, nc;
  588. extern /* Subroutine */ void clagge_(integer *, integer *, integer *,
  589. integer *, real *, complex *, integer *, integer *, complex *,
  590. integer *), claghe_(integer *, integer *, real *, complex *,
  591. integer *, integer *, complex *, integer *);
  592. integer il;
  593. complex ct;
  594. integer iendch, ir, jr, ipackg, mr;
  595. //extern /* Complex */ VOID clarnd_(complex *, integer *, integer *);
  596. extern complex clarnd_(integer *, integer *);
  597. integer minlda;
  598. complex st;
  599. extern /* Subroutine */ void claset_(char *, integer *, integer *, complex
  600. *, complex *, complex *, integer *), clartg_(complex *,
  601. complex *, real *, complex *, complex *);
  602. extern int xerbla_(char *, integer *, ftnlen);
  603. extern void clagsy_(integer *, integer *, real *, complex *,
  604. integer *, integer *, complex *, integer *);
  605. extern real slarnd_(integer *, integer *);
  606. extern /* Subroutine */ void clarot_(logical *, logical *, logical *,
  607. integer *, complex *, complex *, complex *, integer *, complex *,
  608. complex *);
  609. logical iltemp, givens;
  610. integer ioffst, irsign;
  611. logical ilextr, topdwn;
  612. integer ir1, ir2, isympk, jch, llb, jkl, jku, uub;
  613. /* -- LAPACK computational routine (version 3.7.0) -- */
  614. /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
  615. /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
  616. /* December 2016 */
  617. /* ===================================================================== */
  618. /* 1) Decode and Test the input parameters. */
  619. /* Initialize flags & seed. */
  620. /* Parameter adjustments */
  621. --iseed;
  622. --d__;
  623. a_dim1 = *lda;
  624. a_offset = 1 + a_dim1 * 1;
  625. a -= a_offset;
  626. --work;
  627. /* Function Body */
  628. *info = 0;
  629. /* Quick return if possible */
  630. if (*m == 0 || *n == 0) {
  631. return;
  632. }
  633. /* Decode DIST */
  634. if (lsame_(dist, "U")) {
  635. idist = 1;
  636. } else if (lsame_(dist, "S")) {
  637. idist = 2;
  638. } else if (lsame_(dist, "N")) {
  639. idist = 3;
  640. } else {
  641. idist = -1;
  642. }
  643. /* Decode SYM */
  644. if (lsame_(sym, "N")) {
  645. isym = 1;
  646. irsign = 0;
  647. csym = FALSE_;
  648. } else if (lsame_(sym, "P")) {
  649. isym = 2;
  650. irsign = 0;
  651. csym = FALSE_;
  652. } else if (lsame_(sym, "S")) {
  653. isym = 2;
  654. irsign = 0;
  655. csym = TRUE_;
  656. } else if (lsame_(sym, "H")) {
  657. isym = 2;
  658. irsign = 1;
  659. csym = FALSE_;
  660. } else {
  661. isym = -1;
  662. }
  663. /* Decode PACK */
  664. isympk = 0;
  665. if (lsame_(pack, "N")) {
  666. ipack = 0;
  667. } else if (lsame_(pack, "U")) {
  668. ipack = 1;
  669. isympk = 1;
  670. } else if (lsame_(pack, "L")) {
  671. ipack = 2;
  672. isympk = 1;
  673. } else if (lsame_(pack, "C")) {
  674. ipack = 3;
  675. isympk = 2;
  676. } else if (lsame_(pack, "R")) {
  677. ipack = 4;
  678. isympk = 3;
  679. } else if (lsame_(pack, "B")) {
  680. ipack = 5;
  681. isympk = 3;
  682. } else if (lsame_(pack, "Q")) {
  683. ipack = 6;
  684. isympk = 2;
  685. } else if (lsame_(pack, "Z")) {
  686. ipack = 7;
  687. } else {
  688. ipack = -1;
  689. }
  690. /* Set certain internal parameters */
  691. mnmin = f2cmin(*m,*n);
  692. /* Computing MIN */
  693. i__1 = *kl, i__2 = *m - 1;
  694. llb = f2cmin(i__1,i__2);
  695. /* Computing MIN */
  696. i__1 = *ku, i__2 = *n - 1;
  697. uub = f2cmin(i__1,i__2);
  698. /* Computing MIN */
  699. i__1 = *m, i__2 = *n + llb;
  700. mr = f2cmin(i__1,i__2);
  701. /* Computing MIN */
  702. i__1 = *n, i__2 = *m + uub;
  703. nc = f2cmin(i__1,i__2);
  704. if (ipack == 5 || ipack == 6) {
  705. minlda = uub + 1;
  706. } else if (ipack == 7) {
  707. minlda = llb + uub + 1;
  708. } else {
  709. minlda = *m;
  710. }
  711. /* Use Givens rotation method if bandwidth small enough, */
  712. /* or if LDA is too small to store the matrix unpacked. */
  713. givens = FALSE_;
  714. if (isym == 1) {
  715. /* Computing MAX */
  716. i__1 = 1, i__2 = mr + nc;
  717. if ((real) (llb + uub) < (real) f2cmax(i__1,i__2) * .3f) {
  718. givens = TRUE_;
  719. }
  720. } else {
  721. if (llb << 1 < *m) {
  722. givens = TRUE_;
  723. }
  724. }
  725. if (*lda < *m && *lda >= minlda) {
  726. givens = TRUE_;
  727. }
  728. /* Set INFO if an error */
  729. if (*m < 0) {
  730. *info = -1;
  731. } else if (*m != *n && isym != 1) {
  732. *info = -1;
  733. } else if (*n < 0) {
  734. *info = -2;
  735. } else if (idist == -1) {
  736. *info = -3;
  737. } else if (isym == -1) {
  738. *info = -5;
  739. } else if (abs(*mode) > 6) {
  740. *info = -7;
  741. } else if (*mode != 0 && abs(*mode) != 6 && *cond < 1.f) {
  742. *info = -8;
  743. } else if (*kl < 0) {
  744. *info = -10;
  745. } else if (*ku < 0 || isym != 1 && *kl != *ku) {
  746. *info = -11;
  747. } else if (ipack == -1 || isympk == 1 && isym == 1 || isympk == 2 && isym
  748. == 1 && *kl > 0 || isympk == 3 && isym == 1 && *ku > 0 || isympk
  749. != 0 && *m != *n) {
  750. *info = -12;
  751. } else if (*lda < f2cmax(1,minlda)) {
  752. *info = -14;
  753. }
  754. if (*info != 0) {
  755. i__1 = -(*info);
  756. xerbla_("CLATMS", &i__1, 6);
  757. return;
  758. }
  759. /* Initialize random number generator */
  760. for (i__ = 1; i__ <= 4; ++i__) {
  761. iseed[i__] = (i__1 = iseed[i__], abs(i__1)) % 4096;
  762. /* L10: */
  763. }
  764. if (iseed[4] % 2 != 1) {
  765. ++iseed[4];
  766. }
  767. /* 2) Set up D if indicated. */
  768. /* Compute D according to COND and MODE */
  769. slatm1_(mode, cond, &irsign, &idist, &iseed[1], &d__[1], &mnmin, &iinfo);
  770. if (iinfo != 0) {
  771. *info = 1;
  772. return;
  773. }
  774. /* Choose Top-Down if D is (apparently) increasing, */
  775. /* Bottom-Up if D is (apparently) decreasing. */
  776. if (abs(d__[1]) <= (r__1 = d__[mnmin], abs(r__1))) {
  777. topdwn = TRUE_;
  778. } else {
  779. topdwn = FALSE_;
  780. }
  781. if (*mode != 0 && abs(*mode) != 6) {
  782. /* Scale by DMAX */
  783. temp = abs(d__[1]);
  784. i__1 = mnmin;
  785. for (i__ = 2; i__ <= i__1; ++i__) {
  786. /* Computing MAX */
  787. r__2 = temp, r__3 = (r__1 = d__[i__], abs(r__1));
  788. temp = f2cmax(r__2,r__3);
  789. /* L20: */
  790. }
  791. if (temp > 0.f) {
  792. alpha = *dmax__ / temp;
  793. } else {
  794. *info = 2;
  795. return;
  796. }
  797. sscal_(&mnmin, &alpha, &d__[1], &c__1);
  798. }
  799. claset_("Full", lda, n, &c_b1, &c_b1, &a[a_offset], lda);
  800. /* 3) Generate Banded Matrix using Givens rotations. */
  801. /* Also the special case of UUB=LLB=0 */
  802. /* Compute Addressing constants to cover all */
  803. /* storage formats. Whether GE, HE, SY, GB, HB, or SB, */
  804. /* upper or lower triangle or both, */
  805. /* the (i,j)-th element is in */
  806. /* A( i - ISKEW*j + IOFFST, j ) */
  807. if (ipack > 4) {
  808. ilda = *lda - 1;
  809. iskew = 1;
  810. if (ipack > 5) {
  811. ioffst = uub + 1;
  812. } else {
  813. ioffst = 1;
  814. }
  815. } else {
  816. ilda = *lda;
  817. iskew = 0;
  818. ioffst = 0;
  819. }
  820. /* IPACKG is the format that the matrix is generated in. If this is */
  821. /* different from IPACK, then the matrix must be repacked at the */
  822. /* end. It also signals how to compute the norm, for scaling. */
  823. ipackg = 0;
  824. /* Diagonal Matrix -- We are done, unless it */
  825. /* is to be stored HP/SP/PP/TP (PACK='R' or 'C') */
  826. if (llb == 0 && uub == 0) {
  827. i__1 = mnmin;
  828. for (j = 1; j <= i__1; ++j) {
  829. i__2 = (1 - iskew) * j + ioffst + j * a_dim1;
  830. i__3 = j;
  831. q__1.r = d__[i__3], q__1.i = 0.f;
  832. a[i__2].r = q__1.r, a[i__2].i = q__1.i;
  833. /* L30: */
  834. }
  835. if (ipack <= 2 || ipack >= 5) {
  836. ipackg = ipack;
  837. }
  838. } else if (givens) {
  839. /* Check whether to use Givens rotations, */
  840. /* Householder transformations, or nothing. */
  841. if (isym == 1) {
  842. /* Non-symmetric -- A = U D V */
  843. if (ipack > 4) {
  844. ipackg = ipack;
  845. } else {
  846. ipackg = 0;
  847. }
  848. i__1 = mnmin;
  849. for (j = 1; j <= i__1; ++j) {
  850. i__2 = (1 - iskew) * j + ioffst + j * a_dim1;
  851. i__3 = j;
  852. q__1.r = d__[i__3], q__1.i = 0.f;
  853. a[i__2].r = q__1.r, a[i__2].i = q__1.i;
  854. /* L40: */
  855. }
  856. if (topdwn) {
  857. jkl = 0;
  858. i__1 = uub;
  859. for (jku = 1; jku <= i__1; ++jku) {
  860. /* Transform from bandwidth JKL, JKU-1 to JKL, JKU */
  861. /* Last row actually rotated is M */
  862. /* Last column actually rotated is MIN( M+JKU, N ) */
  863. /* Computing MIN */
  864. i__3 = *m + jku;
  865. i__2 = f2cmin(i__3,*n) + jkl - 1;
  866. for (jr = 1; jr <= i__2; ++jr) {
  867. extra.r = 0.f, extra.i = 0.f;
  868. angle = slarnd_(&c__1, &iseed[1]) *
  869. 6.2831853071795864769252867663f;
  870. r__1 = cos(angle);
  871. //clarnd_(&q__2, &c__5, &iseed[1]);
  872. q__2=clarnd_(&c__5, &iseed[1]);
  873. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  874. c__.r = q__1.r, c__.i = q__1.i;
  875. r__1 = sin(angle);
  876. //clarnd_(&q__2, &c__5, &iseed[1]);
  877. q__2=clarnd_(&c__5, &iseed[1]);
  878. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  879. s.r = q__1.r, s.i = q__1.i;
  880. /* Computing MAX */
  881. i__3 = 1, i__4 = jr - jkl;
  882. icol = f2cmax(i__3,i__4);
  883. if (jr < *m) {
  884. /* Computing MIN */
  885. i__3 = *n, i__4 = jr + jku;
  886. il = f2cmin(i__3,i__4) + 1 - icol;
  887. L__1 = jr > jkl;
  888. clarot_(&c_true, &L__1, &c_false, &il, &c__, &s, &
  889. a[jr - iskew * icol + ioffst + icol *
  890. a_dim1], &ilda, &extra, &dummy);
  891. }
  892. /* Chase "EXTRA" back up */
  893. ir = jr;
  894. ic = icol;
  895. i__3 = -jkl - jku;
  896. for (jch = jr - jkl; i__3 < 0 ? jch >= 1 : jch <= 1;
  897. jch += i__3) {
  898. if (ir < *m) {
  899. clartg_(&a[ir + 1 - iskew * (ic + 1) + ioffst
  900. + (ic + 1) * a_dim1], &extra, &realc,
  901. &s, &dummy);
  902. //clarnd_(&q__1, &c__5, &iseed[1]);
  903. q__1=clarnd_(&c__5, &iseed[1]);
  904. dummy.r = q__1.r, dummy.i = q__1.i;
  905. q__2.r = realc * dummy.r, q__2.i = realc *
  906. dummy.i;
  907. r_cnjg(&q__1, &q__2);
  908. c__.r = q__1.r, c__.i = q__1.i;
  909. q__3.r = -s.r, q__3.i = -s.i;
  910. q__2.r = q__3.r * dummy.r - q__3.i * dummy.i,
  911. q__2.i = q__3.r * dummy.i + q__3.i *
  912. dummy.r;
  913. r_cnjg(&q__1, &q__2);
  914. s.r = q__1.r, s.i = q__1.i;
  915. }
  916. /* Computing MAX */
  917. i__4 = 1, i__5 = jch - jku;
  918. irow = f2cmax(i__4,i__5);
  919. il = ir + 2 - irow;
  920. ctemp.r = 0.f, ctemp.i = 0.f;
  921. iltemp = jch > jku;
  922. clarot_(&c_false, &iltemp, &c_true, &il, &c__, &s,
  923. &a[irow - iskew * ic + ioffst + ic *
  924. a_dim1], &ilda, &ctemp, &extra);
  925. if (iltemp) {
  926. clartg_(&a[irow + 1 - iskew * (ic + 1) +
  927. ioffst + (ic + 1) * a_dim1], &ctemp, &
  928. realc, &s, &dummy);
  929. //clarnd_(&q__1, &c__5, &iseed[1]);
  930. q__1=clarnd_(&c__5, &iseed[1]);
  931. dummy.r = q__1.r, dummy.i = q__1.i;
  932. q__2.r = realc * dummy.r, q__2.i = realc *
  933. dummy.i;
  934. r_cnjg(&q__1, &q__2);
  935. c__.r = q__1.r, c__.i = q__1.i;
  936. q__3.r = -s.r, q__3.i = -s.i;
  937. q__2.r = q__3.r * dummy.r - q__3.i * dummy.i,
  938. q__2.i = q__3.r * dummy.i + q__3.i *
  939. dummy.r;
  940. r_cnjg(&q__1, &q__2);
  941. s.r = q__1.r, s.i = q__1.i;
  942. /* Computing MAX */
  943. i__4 = 1, i__5 = jch - jku - jkl;
  944. icol = f2cmax(i__4,i__5);
  945. il = ic + 2 - icol;
  946. extra.r = 0.f, extra.i = 0.f;
  947. L__1 = jch > jku + jkl;
  948. clarot_(&c_true, &L__1, &c_true, &il, &c__, &
  949. s, &a[irow - iskew * icol + ioffst +
  950. icol * a_dim1], &ilda, &extra, &ctemp)
  951. ;
  952. ic = icol;
  953. ir = irow;
  954. }
  955. /* L50: */
  956. }
  957. /* L60: */
  958. }
  959. /* L70: */
  960. }
  961. jku = uub;
  962. i__1 = llb;
  963. for (jkl = 1; jkl <= i__1; ++jkl) {
  964. /* Transform from bandwidth JKL-1, JKU to JKL, JKU */
  965. /* Computing MIN */
  966. i__3 = *n + jkl;
  967. i__2 = f2cmin(i__3,*m) + jku - 1;
  968. for (jc = 1; jc <= i__2; ++jc) {
  969. extra.r = 0.f, extra.i = 0.f;
  970. angle = slarnd_(&c__1, &iseed[1]) *
  971. 6.2831853071795864769252867663f;
  972. r__1 = cos(angle);
  973. //clarnd_(&q__2, &c__5, &iseed[1]);
  974. q__2=clarnd_(&c__5, &iseed[1]);
  975. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  976. c__.r = q__1.r, c__.i = q__1.i;
  977. r__1 = sin(angle);
  978. //clarnd_(&q__2, &c__5, &iseed[1]);
  979. q__2=clarnd_(&c__5, &iseed[1]);
  980. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  981. s.r = q__1.r, s.i = q__1.i;
  982. /* Computing MAX */
  983. i__3 = 1, i__4 = jc - jku;
  984. irow = f2cmax(i__3,i__4);
  985. if (jc < *n) {
  986. /* Computing MIN */
  987. i__3 = *m, i__4 = jc + jkl;
  988. il = f2cmin(i__3,i__4) + 1 - irow;
  989. L__1 = jc > jku;
  990. clarot_(&c_false, &L__1, &c_false, &il, &c__, &s,
  991. &a[irow - iskew * jc + ioffst + jc *
  992. a_dim1], &ilda, &extra, &dummy);
  993. }
  994. /* Chase "EXTRA" back up */
  995. ic = jc;
  996. ir = irow;
  997. i__3 = -jkl - jku;
  998. for (jch = jc - jku; i__3 < 0 ? jch >= 1 : jch <= 1;
  999. jch += i__3) {
  1000. if (ic < *n) {
  1001. clartg_(&a[ir + 1 - iskew * (ic + 1) + ioffst
  1002. + (ic + 1) * a_dim1], &extra, &realc,
  1003. &s, &dummy);
  1004. //clarnd_(&q__1, &c__5, &iseed[1]);
  1005. q__1=clarnd_(&c__5, &iseed[1]);
  1006. dummy.r = q__1.r, dummy.i = q__1.i;
  1007. q__2.r = realc * dummy.r, q__2.i = realc *
  1008. dummy.i;
  1009. r_cnjg(&q__1, &q__2);
  1010. c__.r = q__1.r, c__.i = q__1.i;
  1011. q__3.r = -s.r, q__3.i = -s.i;
  1012. q__2.r = q__3.r * dummy.r - q__3.i * dummy.i,
  1013. q__2.i = q__3.r * dummy.i + q__3.i *
  1014. dummy.r;
  1015. r_cnjg(&q__1, &q__2);
  1016. s.r = q__1.r, s.i = q__1.i;
  1017. }
  1018. /* Computing MAX */
  1019. i__4 = 1, i__5 = jch - jkl;
  1020. icol = f2cmax(i__4,i__5);
  1021. il = ic + 2 - icol;
  1022. ctemp.r = 0.f, ctemp.i = 0.f;
  1023. iltemp = jch > jkl;
  1024. clarot_(&c_true, &iltemp, &c_true, &il, &c__, &s,
  1025. &a[ir - iskew * icol + ioffst + icol *
  1026. a_dim1], &ilda, &ctemp, &extra);
  1027. if (iltemp) {
  1028. clartg_(&a[ir + 1 - iskew * (icol + 1) +
  1029. ioffst + (icol + 1) * a_dim1], &ctemp,
  1030. &realc, &s, &dummy);
  1031. //clarnd_(&q__1, &c__5, &iseed[1]);
  1032. q__1=clarnd_(&c__5, &iseed[1]);
  1033. dummy.r = q__1.r, dummy.i = q__1.i;
  1034. q__2.r = realc * dummy.r, q__2.i = realc *
  1035. dummy.i;
  1036. r_cnjg(&q__1, &q__2);
  1037. c__.r = q__1.r, c__.i = q__1.i;
  1038. q__3.r = -s.r, q__3.i = -s.i;
  1039. q__2.r = q__3.r * dummy.r - q__3.i * dummy.i,
  1040. q__2.i = q__3.r * dummy.i + q__3.i *
  1041. dummy.r;
  1042. r_cnjg(&q__1, &q__2);
  1043. s.r = q__1.r, s.i = q__1.i;
  1044. /* Computing MAX */
  1045. i__4 = 1, i__5 = jch - jkl - jku;
  1046. irow = f2cmax(i__4,i__5);
  1047. il = ir + 2 - irow;
  1048. extra.r = 0.f, extra.i = 0.f;
  1049. L__1 = jch > jkl + jku;
  1050. clarot_(&c_false, &L__1, &c_true, &il, &c__, &
  1051. s, &a[irow - iskew * icol + ioffst +
  1052. icol * a_dim1], &ilda, &extra, &ctemp)
  1053. ;
  1054. ic = icol;
  1055. ir = irow;
  1056. }
  1057. /* L80: */
  1058. }
  1059. /* L90: */
  1060. }
  1061. /* L100: */
  1062. }
  1063. } else {
  1064. /* Bottom-Up -- Start at the bottom right. */
  1065. jkl = 0;
  1066. i__1 = uub;
  1067. for (jku = 1; jku <= i__1; ++jku) {
  1068. /* Transform from bandwidth JKL, JKU-1 to JKL, JKU */
  1069. /* First row actually rotated is M */
  1070. /* First column actually rotated is MIN( M+JKU, N ) */
  1071. /* Computing MIN */
  1072. i__2 = *m, i__3 = *n + jkl;
  1073. iendch = f2cmin(i__2,i__3) - 1;
  1074. /* Computing MIN */
  1075. i__2 = *m + jku;
  1076. i__3 = 1 - jkl;
  1077. for (jc = f2cmin(i__2,*n) - 1; jc >= i__3; --jc) {
  1078. extra.r = 0.f, extra.i = 0.f;
  1079. angle = slarnd_(&c__1, &iseed[1]) *
  1080. 6.2831853071795864769252867663f;
  1081. r__1 = cos(angle);
  1082. //clarnd_(&q__2, &c__5, &iseed[1]);
  1083. q__2=clarnd_(&c__5, &iseed[1]);
  1084. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1085. c__.r = q__1.r, c__.i = q__1.i;
  1086. r__1 = sin(angle);
  1087. //clarnd_(&q__2, &c__5, &iseed[1]);
  1088. q__2=clarnd_(&c__5, &iseed[1]);
  1089. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1090. s.r = q__1.r, s.i = q__1.i;
  1091. /* Computing MAX */
  1092. i__2 = 1, i__4 = jc - jku + 1;
  1093. irow = f2cmax(i__2,i__4);
  1094. if (jc > 0) {
  1095. /* Computing MIN */
  1096. i__2 = *m, i__4 = jc + jkl + 1;
  1097. il = f2cmin(i__2,i__4) + 1 - irow;
  1098. L__1 = jc + jkl < *m;
  1099. clarot_(&c_false, &c_false, &L__1, &il, &c__, &s,
  1100. &a[irow - iskew * jc + ioffst + jc *
  1101. a_dim1], &ilda, &dummy, &extra);
  1102. }
  1103. /* Chase "EXTRA" back down */
  1104. ic = jc;
  1105. i__2 = iendch;
  1106. i__4 = jkl + jku;
  1107. for (jch = jc + jkl; i__4 < 0 ? jch >= i__2 : jch <=
  1108. i__2; jch += i__4) {
  1109. ilextr = ic > 0;
  1110. if (ilextr) {
  1111. clartg_(&a[jch - iskew * ic + ioffst + ic *
  1112. a_dim1], &extra, &realc, &s, &dummy);
  1113. //clarnd_(&q__1, &c__5, &iseed[1]);
  1114. q__1=clarnd_(&c__5, &iseed[1]);
  1115. dummy.r = q__1.r, dummy.i = q__1.i;
  1116. q__1.r = realc * dummy.r, q__1.i = realc *
  1117. dummy.i;
  1118. c__.r = q__1.r, c__.i = q__1.i;
  1119. q__1.r = s.r * dummy.r - s.i * dummy.i,
  1120. q__1.i = s.r * dummy.i + s.i *
  1121. dummy.r;
  1122. s.r = q__1.r, s.i = q__1.i;
  1123. }
  1124. ic = f2cmax(1,ic);
  1125. /* Computing MIN */
  1126. i__5 = *n - 1, i__6 = jch + jku;
  1127. icol = f2cmin(i__5,i__6);
  1128. iltemp = jch + jku < *n;
  1129. ctemp.r = 0.f, ctemp.i = 0.f;
  1130. i__5 = icol + 2 - ic;
  1131. clarot_(&c_true, &ilextr, &iltemp, &i__5, &c__, &
  1132. s, &a[jch - iskew * ic + ioffst + ic *
  1133. a_dim1], &ilda, &extra, &ctemp);
  1134. if (iltemp) {
  1135. clartg_(&a[jch - iskew * icol + ioffst + icol
  1136. * a_dim1], &ctemp, &realc, &s, &dummy)
  1137. ;
  1138. //clarnd_(&q__1, &c__5, &iseed[1]);
  1139. q__1=clarnd_(&c__5, &iseed[1]);
  1140. dummy.r = q__1.r, dummy.i = q__1.i;
  1141. q__1.r = realc * dummy.r, q__1.i = realc *
  1142. dummy.i;
  1143. c__.r = q__1.r, c__.i = q__1.i;
  1144. q__1.r = s.r * dummy.r - s.i * dummy.i,
  1145. q__1.i = s.r * dummy.i + s.i *
  1146. dummy.r;
  1147. s.r = q__1.r, s.i = q__1.i;
  1148. /* Computing MIN */
  1149. i__5 = iendch, i__6 = jch + jkl + jku;
  1150. il = f2cmin(i__5,i__6) + 2 - jch;
  1151. extra.r = 0.f, extra.i = 0.f;
  1152. L__1 = jch + jkl + jku <= iendch;
  1153. clarot_(&c_false, &c_true, &L__1, &il, &c__, &
  1154. s, &a[jch - iskew * icol + ioffst +
  1155. icol * a_dim1], &ilda, &ctemp, &extra)
  1156. ;
  1157. ic = icol;
  1158. }
  1159. /* L110: */
  1160. }
  1161. /* L120: */
  1162. }
  1163. /* L130: */
  1164. }
  1165. jku = uub;
  1166. i__1 = llb;
  1167. for (jkl = 1; jkl <= i__1; ++jkl) {
  1168. /* Transform from bandwidth JKL-1, JKU to JKL, JKU */
  1169. /* First row actually rotated is MIN( N+JKL, M ) */
  1170. /* First column actually rotated is N */
  1171. /* Computing MIN */
  1172. i__3 = *n, i__4 = *m + jku;
  1173. iendch = f2cmin(i__3,i__4) - 1;
  1174. /* Computing MIN */
  1175. i__3 = *n + jkl;
  1176. i__4 = 1 - jku;
  1177. for (jr = f2cmin(i__3,*m) - 1; jr >= i__4; --jr) {
  1178. extra.r = 0.f, extra.i = 0.f;
  1179. angle = slarnd_(&c__1, &iseed[1]) *
  1180. 6.2831853071795864769252867663f;
  1181. r__1 = cos(angle);
  1182. //clarnd_(&q__2, &c__5, &iseed[1]);
  1183. q__2=clarnd_(&c__5, &iseed[1]);
  1184. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1185. c__.r = q__1.r, c__.i = q__1.i;
  1186. r__1 = sin(angle);
  1187. //clarnd_(&q__2, &c__5, &iseed[1]);
  1188. q__2=clarnd_(&c__5, &iseed[1]);
  1189. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1190. s.r = q__1.r, s.i = q__1.i;
  1191. /* Computing MAX */
  1192. i__3 = 1, i__2 = jr - jkl + 1;
  1193. icol = f2cmax(i__3,i__2);
  1194. if (jr > 0) {
  1195. /* Computing MIN */
  1196. i__3 = *n, i__2 = jr + jku + 1;
  1197. il = f2cmin(i__3,i__2) + 1 - icol;
  1198. L__1 = jr + jku < *n;
  1199. clarot_(&c_true, &c_false, &L__1, &il, &c__, &s, &
  1200. a[jr - iskew * icol + ioffst + icol *
  1201. a_dim1], &ilda, &dummy, &extra);
  1202. }
  1203. /* Chase "EXTRA" back down */
  1204. ir = jr;
  1205. i__3 = iendch;
  1206. i__2 = jkl + jku;
  1207. for (jch = jr + jku; i__2 < 0 ? jch >= i__3 : jch <=
  1208. i__3; jch += i__2) {
  1209. ilextr = ir > 0;
  1210. if (ilextr) {
  1211. clartg_(&a[ir - iskew * jch + ioffst + jch *
  1212. a_dim1], &extra, &realc, &s, &dummy);
  1213. //clarnd_(&q__1, &c__5, &iseed[1]);
  1214. q__1=clarnd_(&c__5, &iseed[1]);
  1215. dummy.r = q__1.r, dummy.i = q__1.i;
  1216. q__1.r = realc * dummy.r, q__1.i = realc *
  1217. dummy.i;
  1218. c__.r = q__1.r, c__.i = q__1.i;
  1219. q__1.r = s.r * dummy.r - s.i * dummy.i,
  1220. q__1.i = s.r * dummy.i + s.i *
  1221. dummy.r;
  1222. s.r = q__1.r, s.i = q__1.i;
  1223. }
  1224. ir = f2cmax(1,ir);
  1225. /* Computing MIN */
  1226. i__5 = *m - 1, i__6 = jch + jkl;
  1227. irow = f2cmin(i__5,i__6);
  1228. iltemp = jch + jkl < *m;
  1229. ctemp.r = 0.f, ctemp.i = 0.f;
  1230. i__5 = irow + 2 - ir;
  1231. clarot_(&c_false, &ilextr, &iltemp, &i__5, &c__, &
  1232. s, &a[ir - iskew * jch + ioffst + jch *
  1233. a_dim1], &ilda, &extra, &ctemp);
  1234. if (iltemp) {
  1235. clartg_(&a[irow - iskew * jch + ioffst + jch *
  1236. a_dim1], &ctemp, &realc, &s, &dummy);
  1237. //clarnd_(&q__1, &c__5, &iseed[1]);
  1238. q__1=clarnd_(&c__5, &iseed[1]);
  1239. dummy.r = q__1.r, dummy.i = q__1.i;
  1240. q__1.r = realc * dummy.r, q__1.i = realc *
  1241. dummy.i;
  1242. c__.r = q__1.r, c__.i = q__1.i;
  1243. q__1.r = s.r * dummy.r - s.i * dummy.i,
  1244. q__1.i = s.r * dummy.i + s.i *
  1245. dummy.r;
  1246. s.r = q__1.r, s.i = q__1.i;
  1247. /* Computing MIN */
  1248. i__5 = iendch, i__6 = jch + jkl + jku;
  1249. il = f2cmin(i__5,i__6) + 2 - jch;
  1250. extra.r = 0.f, extra.i = 0.f;
  1251. L__1 = jch + jkl + jku <= iendch;
  1252. clarot_(&c_true, &c_true, &L__1, &il, &c__, &
  1253. s, &a[irow - iskew * jch + ioffst +
  1254. jch * a_dim1], &ilda, &ctemp, &extra);
  1255. ir = irow;
  1256. }
  1257. /* L140: */
  1258. }
  1259. /* L150: */
  1260. }
  1261. /* L160: */
  1262. }
  1263. }
  1264. } else {
  1265. /* Symmetric -- A = U D U' */
  1266. /* Hermitian -- A = U D U* */
  1267. ipackg = ipack;
  1268. ioffg = ioffst;
  1269. if (topdwn) {
  1270. /* Top-Down -- Generate Upper triangle only */
  1271. if (ipack >= 5) {
  1272. ipackg = 6;
  1273. ioffg = uub + 1;
  1274. } else {
  1275. ipackg = 1;
  1276. }
  1277. i__1 = mnmin;
  1278. for (j = 1; j <= i__1; ++j) {
  1279. i__4 = (1 - iskew) * j + ioffg + j * a_dim1;
  1280. i__2 = j;
  1281. q__1.r = d__[i__2], q__1.i = 0.f;
  1282. a[i__4].r = q__1.r, a[i__4].i = q__1.i;
  1283. /* L170: */
  1284. }
  1285. i__1 = uub;
  1286. for (k = 1; k <= i__1; ++k) {
  1287. i__4 = *n - 1;
  1288. for (jc = 1; jc <= i__4; ++jc) {
  1289. /* Computing MAX */
  1290. i__2 = 1, i__3 = jc - k;
  1291. irow = f2cmax(i__2,i__3);
  1292. /* Computing MIN */
  1293. i__2 = jc + 1, i__3 = k + 2;
  1294. il = f2cmin(i__2,i__3);
  1295. extra.r = 0.f, extra.i = 0.f;
  1296. i__2 = jc - iskew * (jc + 1) + ioffg + (jc + 1) *
  1297. a_dim1;
  1298. ctemp.r = a[i__2].r, ctemp.i = a[i__2].i;
  1299. angle = slarnd_(&c__1, &iseed[1]) *
  1300. 6.2831853071795864769252867663f;
  1301. r__1 = cos(angle);
  1302. //clarnd_(&q__2, &c__5, &iseed[1]);
  1303. q__2=clarnd_(&c__5, &iseed[1]);
  1304. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1305. c__.r = q__1.r, c__.i = q__1.i;
  1306. r__1 = sin(angle);
  1307. //clarnd_(&q__2, &c__5, &iseed[1]);
  1308. q__2=clarnd_(&c__5, &iseed[1]);
  1309. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1310. s.r = q__1.r, s.i = q__1.i;
  1311. if (csym) {
  1312. ct.r = c__.r, ct.i = c__.i;
  1313. st.r = s.r, st.i = s.i;
  1314. } else {
  1315. r_cnjg(&q__1, &ctemp);
  1316. ctemp.r = q__1.r, ctemp.i = q__1.i;
  1317. r_cnjg(&q__1, &c__);
  1318. ct.r = q__1.r, ct.i = q__1.i;
  1319. r_cnjg(&q__1, &s);
  1320. st.r = q__1.r, st.i = q__1.i;
  1321. }
  1322. L__1 = jc > k;
  1323. clarot_(&c_false, &L__1, &c_true, &il, &c__, &s, &a[
  1324. irow - iskew * jc + ioffg + jc * a_dim1], &
  1325. ilda, &extra, &ctemp);
  1326. /* Computing MIN */
  1327. i__3 = k, i__5 = *n - jc;
  1328. i__2 = f2cmin(i__3,i__5) + 1;
  1329. clarot_(&c_true, &c_true, &c_false, &i__2, &ct, &st, &
  1330. a[(1 - iskew) * jc + ioffg + jc * a_dim1], &
  1331. ilda, &ctemp, &dummy);
  1332. /* Chase EXTRA back up the matrix */
  1333. icol = jc;
  1334. i__2 = -k;
  1335. for (jch = jc - k; i__2 < 0 ? jch >= 1 : jch <= 1;
  1336. jch += i__2) {
  1337. clartg_(&a[jch + 1 - iskew * (icol + 1) + ioffg +
  1338. (icol + 1) * a_dim1], &extra, &realc, &s,
  1339. &dummy);
  1340. //clarnd_(&q__1, &c__5, &iseed[1]);
  1341. q__1=clarnd_(&c__5, &iseed[1]);
  1342. dummy.r = q__1.r, dummy.i = q__1.i;
  1343. q__2.r = realc * dummy.r, q__2.i = realc *
  1344. dummy.i;
  1345. r_cnjg(&q__1, &q__2);
  1346. c__.r = q__1.r, c__.i = q__1.i;
  1347. q__3.r = -s.r, q__3.i = -s.i;
  1348. q__2.r = q__3.r * dummy.r - q__3.i * dummy.i,
  1349. q__2.i = q__3.r * dummy.i + q__3.i *
  1350. dummy.r;
  1351. r_cnjg(&q__1, &q__2);
  1352. s.r = q__1.r, s.i = q__1.i;
  1353. i__3 = jch - iskew * (jch + 1) + ioffg + (jch + 1)
  1354. * a_dim1;
  1355. ctemp.r = a[i__3].r, ctemp.i = a[i__3].i;
  1356. if (csym) {
  1357. ct.r = c__.r, ct.i = c__.i;
  1358. st.r = s.r, st.i = s.i;
  1359. } else {
  1360. r_cnjg(&q__1, &ctemp);
  1361. ctemp.r = q__1.r, ctemp.i = q__1.i;
  1362. r_cnjg(&q__1, &c__);
  1363. ct.r = q__1.r, ct.i = q__1.i;
  1364. r_cnjg(&q__1, &s);
  1365. st.r = q__1.r, st.i = q__1.i;
  1366. }
  1367. i__3 = k + 2;
  1368. clarot_(&c_true, &c_true, &c_true, &i__3, &c__, &
  1369. s, &a[(1 - iskew) * jch + ioffg + jch *
  1370. a_dim1], &ilda, &ctemp, &extra);
  1371. /* Computing MAX */
  1372. i__3 = 1, i__5 = jch - k;
  1373. irow = f2cmax(i__3,i__5);
  1374. /* Computing MIN */
  1375. i__3 = jch + 1, i__5 = k + 2;
  1376. il = f2cmin(i__3,i__5);
  1377. extra.r = 0.f, extra.i = 0.f;
  1378. L__1 = jch > k;
  1379. clarot_(&c_false, &L__1, &c_true, &il, &ct, &st, &
  1380. a[irow - iskew * jch + ioffg + jch *
  1381. a_dim1], &ilda, &extra, &ctemp);
  1382. icol = jch;
  1383. /* L180: */
  1384. }
  1385. /* L190: */
  1386. }
  1387. /* L200: */
  1388. }
  1389. /* If we need lower triangle, copy from upper. Note that */
  1390. /* the order of copying is chosen to work for 'q' -> 'b' */
  1391. if (ipack != ipackg && ipack != 3) {
  1392. i__1 = *n;
  1393. for (jc = 1; jc <= i__1; ++jc) {
  1394. irow = ioffst - iskew * jc;
  1395. if (csym) {
  1396. /* Computing MIN */
  1397. i__2 = *n, i__3 = jc + uub;
  1398. i__4 = f2cmin(i__2,i__3);
  1399. for (jr = jc; jr <= i__4; ++jr) {
  1400. i__2 = jr + irow + jc * a_dim1;
  1401. i__3 = jc - iskew * jr + ioffg + jr * a_dim1;
  1402. a[i__2].r = a[i__3].r, a[i__2].i = a[i__3].i;
  1403. /* L210: */
  1404. }
  1405. } else {
  1406. /* Computing MIN */
  1407. i__2 = *n, i__3 = jc + uub;
  1408. i__4 = f2cmin(i__2,i__3);
  1409. for (jr = jc; jr <= i__4; ++jr) {
  1410. i__2 = jr + irow + jc * a_dim1;
  1411. r_cnjg(&q__1, &a[jc - iskew * jr + ioffg + jr
  1412. * a_dim1]);
  1413. a[i__2].r = q__1.r, a[i__2].i = q__1.i;
  1414. /* L220: */
  1415. }
  1416. }
  1417. /* L230: */
  1418. }
  1419. if (ipack == 5) {
  1420. i__1 = *n;
  1421. for (jc = *n - uub + 1; jc <= i__1; ++jc) {
  1422. i__4 = uub + 1;
  1423. for (jr = *n + 2 - jc; jr <= i__4; ++jr) {
  1424. i__2 = jr + jc * a_dim1;
  1425. a[i__2].r = 0.f, a[i__2].i = 0.f;
  1426. /* L240: */
  1427. }
  1428. /* L250: */
  1429. }
  1430. }
  1431. if (ipackg == 6) {
  1432. ipackg = ipack;
  1433. } else {
  1434. ipackg = 0;
  1435. }
  1436. }
  1437. } else {
  1438. /* Bottom-Up -- Generate Lower triangle only */
  1439. if (ipack >= 5) {
  1440. ipackg = 5;
  1441. if (ipack == 6) {
  1442. ioffg = 1;
  1443. }
  1444. } else {
  1445. ipackg = 2;
  1446. }
  1447. i__1 = mnmin;
  1448. for (j = 1; j <= i__1; ++j) {
  1449. i__4 = (1 - iskew) * j + ioffg + j * a_dim1;
  1450. i__2 = j;
  1451. q__1.r = d__[i__2], q__1.i = 0.f;
  1452. a[i__4].r = q__1.r, a[i__4].i = q__1.i;
  1453. /* L260: */
  1454. }
  1455. i__1 = uub;
  1456. for (k = 1; k <= i__1; ++k) {
  1457. for (jc = *n - 1; jc >= 1; --jc) {
  1458. /* Computing MIN */
  1459. i__4 = *n + 1 - jc, i__2 = k + 2;
  1460. il = f2cmin(i__4,i__2);
  1461. extra.r = 0.f, extra.i = 0.f;
  1462. i__4 = (1 - iskew) * jc + 1 + ioffg + jc * a_dim1;
  1463. ctemp.r = a[i__4].r, ctemp.i = a[i__4].i;
  1464. angle = slarnd_(&c__1, &iseed[1]) *
  1465. 6.2831853071795864769252867663f;
  1466. r__1 = cos(angle);
  1467. //clarnd_(&q__2, &c__5, &iseed[1]);
  1468. q__2=clarnd_(&c__5, &iseed[1]);
  1469. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1470. c__.r = q__1.r, c__.i = q__1.i;
  1471. r__1 = sin(angle);
  1472. //clarnd_(&q__2, &c__5, &iseed[1]);
  1473. q__2=clarnd_(&c__5, &iseed[1]);
  1474. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1475. s.r = q__1.r, s.i = q__1.i;
  1476. if (csym) {
  1477. ct.r = c__.r, ct.i = c__.i;
  1478. st.r = s.r, st.i = s.i;
  1479. } else {
  1480. r_cnjg(&q__1, &ctemp);
  1481. ctemp.r = q__1.r, ctemp.i = q__1.i;
  1482. r_cnjg(&q__1, &c__);
  1483. ct.r = q__1.r, ct.i = q__1.i;
  1484. r_cnjg(&q__1, &s);
  1485. st.r = q__1.r, st.i = q__1.i;
  1486. }
  1487. L__1 = *n - jc > k;
  1488. clarot_(&c_false, &c_true, &L__1, &il, &c__, &s, &a[(
  1489. 1 - iskew) * jc + ioffg + jc * a_dim1], &ilda,
  1490. &ctemp, &extra);
  1491. /* Computing MAX */
  1492. i__4 = 1, i__2 = jc - k + 1;
  1493. icol = f2cmax(i__4,i__2);
  1494. i__4 = jc + 2 - icol;
  1495. clarot_(&c_true, &c_false, &c_true, &i__4, &ct, &st, &
  1496. a[jc - iskew * icol + ioffg + icol * a_dim1],
  1497. &ilda, &dummy, &ctemp);
  1498. /* Chase EXTRA back down the matrix */
  1499. icol = jc;
  1500. i__4 = *n - 1;
  1501. i__2 = k;
  1502. for (jch = jc + k; i__2 < 0 ? jch >= i__4 : jch <=
  1503. i__4; jch += i__2) {
  1504. clartg_(&a[jch - iskew * icol + ioffg + icol *
  1505. a_dim1], &extra, &realc, &s, &dummy);
  1506. //clarnd_(&q__1, &c__5, &iseed[1]);
  1507. q__1=clarnd_(&c__5, &iseed[1]);
  1508. dummy.r = q__1.r, dummy.i = q__1.i;
  1509. q__1.r = realc * dummy.r, q__1.i = realc *
  1510. dummy.i;
  1511. c__.r = q__1.r, c__.i = q__1.i;
  1512. q__1.r = s.r * dummy.r - s.i * dummy.i, q__1.i =
  1513. s.r * dummy.i + s.i * dummy.r;
  1514. s.r = q__1.r, s.i = q__1.i;
  1515. i__3 = (1 - iskew) * jch + 1 + ioffg + jch *
  1516. a_dim1;
  1517. ctemp.r = a[i__3].r, ctemp.i = a[i__3].i;
  1518. if (csym) {
  1519. ct.r = c__.r, ct.i = c__.i;
  1520. st.r = s.r, st.i = s.i;
  1521. } else {
  1522. r_cnjg(&q__1, &ctemp);
  1523. ctemp.r = q__1.r, ctemp.i = q__1.i;
  1524. r_cnjg(&q__1, &c__);
  1525. ct.r = q__1.r, ct.i = q__1.i;
  1526. r_cnjg(&q__1, &s);
  1527. st.r = q__1.r, st.i = q__1.i;
  1528. }
  1529. i__3 = k + 2;
  1530. clarot_(&c_true, &c_true, &c_true, &i__3, &c__, &
  1531. s, &a[jch - iskew * icol + ioffg + icol *
  1532. a_dim1], &ilda, &extra, &ctemp);
  1533. /* Computing MIN */
  1534. i__3 = *n + 1 - jch, i__5 = k + 2;
  1535. il = f2cmin(i__3,i__5);
  1536. extra.r = 0.f, extra.i = 0.f;
  1537. L__1 = *n - jch > k;
  1538. clarot_(&c_false, &c_true, &L__1, &il, &ct, &st, &
  1539. a[(1 - iskew) * jch + ioffg + jch *
  1540. a_dim1], &ilda, &ctemp, &extra);
  1541. icol = jch;
  1542. /* L270: */
  1543. }
  1544. /* L280: */
  1545. }
  1546. /* L290: */
  1547. }
  1548. /* If we need upper triangle, copy from lower. Note that */
  1549. /* the order of copying is chosen to work for 'b' -> 'q' */
  1550. if (ipack != ipackg && ipack != 4) {
  1551. for (jc = *n; jc >= 1; --jc) {
  1552. irow = ioffst - iskew * jc;
  1553. if (csym) {
  1554. /* Computing MAX */
  1555. i__2 = 1, i__4 = jc - uub;
  1556. i__1 = f2cmax(i__2,i__4);
  1557. for (jr = jc; jr >= i__1; --jr) {
  1558. i__2 = jr + irow + jc * a_dim1;
  1559. i__4 = jc - iskew * jr + ioffg + jr * a_dim1;
  1560. a[i__2].r = a[i__4].r, a[i__2].i = a[i__4].i;
  1561. /* L300: */
  1562. }
  1563. } else {
  1564. /* Computing MAX */
  1565. i__2 = 1, i__4 = jc - uub;
  1566. i__1 = f2cmax(i__2,i__4);
  1567. for (jr = jc; jr >= i__1; --jr) {
  1568. i__2 = jr + irow + jc * a_dim1;
  1569. r_cnjg(&q__1, &a[jc - iskew * jr + ioffg + jr
  1570. * a_dim1]);
  1571. a[i__2].r = q__1.r, a[i__2].i = q__1.i;
  1572. /* L310: */
  1573. }
  1574. }
  1575. /* L320: */
  1576. }
  1577. if (ipack == 6) {
  1578. i__1 = uub;
  1579. for (jc = 1; jc <= i__1; ++jc) {
  1580. i__2 = uub + 1 - jc;
  1581. for (jr = 1; jr <= i__2; ++jr) {
  1582. i__4 = jr + jc * a_dim1;
  1583. a[i__4].r = 0.f, a[i__4].i = 0.f;
  1584. /* L330: */
  1585. }
  1586. /* L340: */
  1587. }
  1588. }
  1589. if (ipackg == 5) {
  1590. ipackg = ipack;
  1591. } else {
  1592. ipackg = 0;
  1593. }
  1594. }
  1595. }
  1596. /* Ensure that the diagonal is real if Hermitian */
  1597. if (! csym) {
  1598. i__1 = *n;
  1599. for (jc = 1; jc <= i__1; ++jc) {
  1600. irow = ioffst + (1 - iskew) * jc;
  1601. i__2 = irow + jc * a_dim1;
  1602. i__4 = irow + jc * a_dim1;
  1603. r__1 = a[i__4].r;
  1604. q__1.r = r__1, q__1.i = 0.f;
  1605. a[i__2].r = q__1.r, a[i__2].i = q__1.i;
  1606. /* L350: */
  1607. }
  1608. }
  1609. }
  1610. } else {
  1611. /* 4) Generate Banded Matrix by first */
  1612. /* Rotating by random Unitary matrices, */
  1613. /* then reducing the bandwidth using Householder */
  1614. /* transformations. */
  1615. /* Note: we should get here only if LDA .ge. N */
  1616. if (isym == 1) {
  1617. /* Non-symmetric -- A = U D V */
  1618. clagge_(&mr, &nc, &llb, &uub, &d__[1], &a[a_offset], lda, &iseed[
  1619. 1], &work[1], &iinfo);
  1620. } else {
  1621. /* Symmetric -- A = U D U' or */
  1622. /* Hermitian -- A = U D U* */
  1623. if (csym) {
  1624. clagsy_(m, &llb, &d__[1], &a[a_offset], lda, &iseed[1], &work[
  1625. 1], &iinfo);
  1626. } else {
  1627. claghe_(m, &llb, &d__[1], &a[a_offset], lda, &iseed[1], &work[
  1628. 1], &iinfo);
  1629. }
  1630. }
  1631. if (iinfo != 0) {
  1632. *info = 3;
  1633. return;
  1634. }
  1635. }
  1636. /* 5) Pack the matrix */
  1637. if (ipack != ipackg) {
  1638. if (ipack == 1) {
  1639. /* 'U' -- Upper triangular, not packed */
  1640. i__1 = *m;
  1641. for (j = 1; j <= i__1; ++j) {
  1642. i__2 = *m;
  1643. for (i__ = j + 1; i__ <= i__2; ++i__) {
  1644. i__4 = i__ + j * a_dim1;
  1645. a[i__4].r = 0.f, a[i__4].i = 0.f;
  1646. /* L360: */
  1647. }
  1648. /* L370: */
  1649. }
  1650. } else if (ipack == 2) {
  1651. /* 'L' -- Lower triangular, not packed */
  1652. i__1 = *m;
  1653. for (j = 2; j <= i__1; ++j) {
  1654. i__2 = j - 1;
  1655. for (i__ = 1; i__ <= i__2; ++i__) {
  1656. i__4 = i__ + j * a_dim1;
  1657. a[i__4].r = 0.f, a[i__4].i = 0.f;
  1658. /* L380: */
  1659. }
  1660. /* L390: */
  1661. }
  1662. } else if (ipack == 3) {
  1663. /* 'C' -- Upper triangle packed Columnwise. */
  1664. icol = 1;
  1665. irow = 0;
  1666. i__1 = *m;
  1667. for (j = 1; j <= i__1; ++j) {
  1668. i__2 = j;
  1669. for (i__ = 1; i__ <= i__2; ++i__) {
  1670. ++irow;
  1671. if (irow > *lda) {
  1672. irow = 1;
  1673. ++icol;
  1674. }
  1675. i__4 = irow + icol * a_dim1;
  1676. i__3 = i__ + j * a_dim1;
  1677. a[i__4].r = a[i__3].r, a[i__4].i = a[i__3].i;
  1678. /* L400: */
  1679. }
  1680. /* L410: */
  1681. }
  1682. } else if (ipack == 4) {
  1683. /* 'R' -- Lower triangle packed Columnwise. */
  1684. icol = 1;
  1685. irow = 0;
  1686. i__1 = *m;
  1687. for (j = 1; j <= i__1; ++j) {
  1688. i__2 = *m;
  1689. for (i__ = j; i__ <= i__2; ++i__) {
  1690. ++irow;
  1691. if (irow > *lda) {
  1692. irow = 1;
  1693. ++icol;
  1694. }
  1695. i__4 = irow + icol * a_dim1;
  1696. i__3 = i__ + j * a_dim1;
  1697. a[i__4].r = a[i__3].r, a[i__4].i = a[i__3].i;
  1698. /* L420: */
  1699. }
  1700. /* L430: */
  1701. }
  1702. } else if (ipack >= 5) {
  1703. /* 'B' -- The lower triangle is packed as a band matrix. */
  1704. /* 'Q' -- The upper triangle is packed as a band matrix. */
  1705. /* 'Z' -- The whole matrix is packed as a band matrix. */
  1706. if (ipack == 5) {
  1707. uub = 0;
  1708. }
  1709. if (ipack == 6) {
  1710. llb = 0;
  1711. }
  1712. i__1 = uub;
  1713. for (j = 1; j <= i__1; ++j) {
  1714. /* Computing MIN */
  1715. i__2 = j + llb;
  1716. for (i__ = f2cmin(i__2,*m); i__ >= 1; --i__) {
  1717. i__2 = i__ - j + uub + 1 + j * a_dim1;
  1718. i__4 = i__ + j * a_dim1;
  1719. a[i__2].r = a[i__4].r, a[i__2].i = a[i__4].i;
  1720. /* L440: */
  1721. }
  1722. /* L450: */
  1723. }
  1724. i__1 = *n;
  1725. for (j = uub + 2; j <= i__1; ++j) {
  1726. /* Computing MIN */
  1727. i__4 = j + llb;
  1728. i__2 = f2cmin(i__4,*m);
  1729. for (i__ = j - uub; i__ <= i__2; ++i__) {
  1730. i__4 = i__ - j + uub + 1 + j * a_dim1;
  1731. i__3 = i__ + j * a_dim1;
  1732. a[i__4].r = a[i__3].r, a[i__4].i = a[i__3].i;
  1733. /* L460: */
  1734. }
  1735. /* L470: */
  1736. }
  1737. }
  1738. /* If packed, zero out extraneous elements. */
  1739. /* Symmetric/Triangular Packed -- */
  1740. /* zero out everything after A(IROW,ICOL) */
  1741. if (ipack == 3 || ipack == 4) {
  1742. i__1 = *m;
  1743. for (jc = icol; jc <= i__1; ++jc) {
  1744. i__2 = *lda;
  1745. for (jr = irow + 1; jr <= i__2; ++jr) {
  1746. i__4 = jr + jc * a_dim1;
  1747. a[i__4].r = 0.f, a[i__4].i = 0.f;
  1748. /* L480: */
  1749. }
  1750. irow = 0;
  1751. /* L490: */
  1752. }
  1753. } else if (ipack >= 5) {
  1754. /* Packed Band -- */
  1755. /* 1st row is now in A( UUB+2-j, j), zero above it */
  1756. /* m-th row is now in A( M+UUB-j,j), zero below it */
  1757. /* last non-zero diagonal is now in A( UUB+LLB+1,j ), */
  1758. /* zero below it, too. */
  1759. ir1 = uub + llb + 2;
  1760. ir2 = uub + *m + 2;
  1761. i__1 = *n;
  1762. for (jc = 1; jc <= i__1; ++jc) {
  1763. i__2 = uub + 1 - jc;
  1764. for (jr = 1; jr <= i__2; ++jr) {
  1765. i__4 = jr + jc * a_dim1;
  1766. a[i__4].r = 0.f, a[i__4].i = 0.f;
  1767. /* L500: */
  1768. }
  1769. /* Computing MAX */
  1770. /* Computing MIN */
  1771. i__3 = ir1, i__5 = ir2 - jc;
  1772. i__2 = 1, i__4 = f2cmin(i__3,i__5);
  1773. i__6 = *lda;
  1774. for (jr = f2cmax(i__2,i__4); jr <= i__6; ++jr) {
  1775. i__2 = jr + jc * a_dim1;
  1776. a[i__2].r = 0.f, a[i__2].i = 0.f;
  1777. /* L510: */
  1778. }
  1779. /* L520: */
  1780. }
  1781. }
  1782. }
  1783. return;
  1784. /* End of CLATMS */
  1785. } /* clatms_ */