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.

clatmt.c 59 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977
  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 CLATMT */
  245. /* =========== DOCUMENTATION =========== */
  246. /* Online html documentation available at */
  247. /* http://www.netlib.org/lapack/explore-html/ */
  248. /* Definition: */
  249. /* =========== */
  250. /* SUBROUTINE CLATMT( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, */
  251. /* RANK, KL, KU, PACK, A, LDA, WORK, INFO ) */
  252. /* REAL COND, DMAX */
  253. /* INTEGER INFO, KL, KU, LDA, M, MODE, N, RANK */
  254. /* CHARACTER DIST, PACK, SYM */
  255. /* COMPLEX A( LDA, * ), WORK( * ) */
  256. /* REAL D( * ) */
  257. /* INTEGER ISEED( 4 ) */
  258. /* > \par Purpose: */
  259. /* ============= */
  260. /* > */
  261. /* > \verbatim */
  262. /* > */
  263. /* > CLATMT generates random matrices with specified singular values */
  264. /* > (or hermitian with specified eigenvalues) */
  265. /* > for testing LAPACK programs. */
  266. /* > */
  267. /* > CLATMT 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 CLATMT */
  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:RANK)=1.0/COND */
  395. /* > MODE = 2 sets D(1:RANK-1)=1 and D(RANK)=1.0/COND */
  396. /* > MODE = 3 sets D(I)=COND**(-(I-1)/(RANK-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] RANK */
  433. /* > \verbatim */
  434. /* > RANK is INTEGER */
  435. /* > The rank of matrix to be generated for modes 1,2,3 only. */
  436. /* > D( RANK+1:N ) = 0. */
  437. /* > Not modified. */
  438. /* > \endverbatim */
  439. /* > */
  440. /* > \param[in] KL */
  441. /* > \verbatim */
  442. /* > KL is INTEGER */
  443. /* > This specifies the lower bandwidth of the matrix. For */
  444. /* > example, KL=0 implies upper triangular, KL=1 implies upper */
  445. /* > Hessenberg, and KL being at least M-1 means that the matrix */
  446. /* > has full lower bandwidth. KL must equal KU if the matrix */
  447. /* > is symmetric or hermitian. */
  448. /* > Not modified. */
  449. /* > \endverbatim */
  450. /* > */
  451. /* > \param[in] KU */
  452. /* > \verbatim */
  453. /* > KU is INTEGER */
  454. /* > This specifies the upper bandwidth of the matrix. For */
  455. /* > example, KU=0 implies lower triangular, KU=1 implies lower */
  456. /* > Hessenberg, and KU being at least N-1 means that the matrix */
  457. /* > has full upper bandwidth. KL must equal KU if the matrix */
  458. /* > is symmetric or hermitian. */
  459. /* > Not modified. */
  460. /* > \endverbatim */
  461. /* > */
  462. /* > \param[in] PACK */
  463. /* > \verbatim */
  464. /* > PACK is CHARACTER*1 */
  465. /* > This specifies packing of matrix as follows: */
  466. /* > 'N' => no packing */
  467. /* > 'U' => zero out all subdiagonal entries (if symmetric */
  468. /* > or hermitian) */
  469. /* > 'L' => zero out all superdiagonal entries (if symmetric */
  470. /* > or hermitian) */
  471. /* > 'C' => store the upper triangle columnwise (only if the */
  472. /* > matrix is symmetric, hermitian, or upper triangular) */
  473. /* > 'R' => store the lower triangle columnwise (only if the */
  474. /* > matrix is symmetric, hermitian, or lower triangular) */
  475. /* > 'B' => store the lower triangle in band storage scheme */
  476. /* > (only if the matrix is symmetric, hermitian, or */
  477. /* > lower triangular) */
  478. /* > 'Q' => store the upper triangle in band storage scheme */
  479. /* > (only if the matrix is symmetric, hermitian, or */
  480. /* > upper triangular) */
  481. /* > 'Z' => store the entire matrix in band storage scheme */
  482. /* > (pivoting can be provided for by using this */
  483. /* > option to store A in the trailing rows of */
  484. /* > the allocated storage) */
  485. /* > */
  486. /* > Using these options, the various LAPACK packed and banded */
  487. /* > storage schemes can be obtained: */
  488. /* > GB - use 'Z' */
  489. /* > PB, SB, HB, or TB - use 'B' or 'Q' */
  490. /* > PP, SP, HB, or TP - use 'C' or 'R' */
  491. /* > */
  492. /* > If two calls to CLATMT differ only in the PACK parameter, */
  493. /* > they will generate mathematically equivalent matrices. */
  494. /* > Not modified. */
  495. /* > \endverbatim */
  496. /* > */
  497. /* > \param[in,out] A */
  498. /* > \verbatim */
  499. /* > A is COMPLEX array, dimension ( LDA, N ) */
  500. /* > On exit A is the desired test matrix. A is first generated */
  501. /* > in full (unpacked) form, and then packed, if so specified */
  502. /* > by PACK. Thus, the first M elements of the first N */
  503. /* > columns will always be modified. If PACK specifies a */
  504. /* > packed or banded storage scheme, all LDA elements of the */
  505. /* > first N columns will be modified; the elements of the */
  506. /* > array which do not correspond to elements of the generated */
  507. /* > matrix are set to zero. */
  508. /* > Modified. */
  509. /* > \endverbatim */
  510. /* > */
  511. /* > \param[in] LDA */
  512. /* > \verbatim */
  513. /* > LDA is INTEGER */
  514. /* > LDA specifies the first dimension of A as declared in the */
  515. /* > calling program. If PACK='N', 'U', 'L', 'C', or 'R', then */
  516. /* > LDA must be at least M. If PACK='B' or 'Q', then LDA must */
  517. /* > be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)). */
  518. /* > If PACK='Z', LDA must be large enough to hold the packed */
  519. /* > array: MIN( KU, N-1) + MIN( KL, M-1) + 1. */
  520. /* > Not modified. */
  521. /* > \endverbatim */
  522. /* > */
  523. /* > \param[out] WORK */
  524. /* > \verbatim */
  525. /* > WORK is COMPLEX array, dimension ( 3*MAX( N, M ) ) */
  526. /* > Workspace. */
  527. /* > Modified. */
  528. /* > \endverbatim */
  529. /* > */
  530. /* > \param[out] INFO */
  531. /* > \verbatim */
  532. /* > INFO is INTEGER */
  533. /* > Error code. On exit, INFO will be set to one of the */
  534. /* > following values: */
  535. /* > 0 => normal return */
  536. /* > -1 => M negative or unequal to N and SYM='S', 'H', or 'P' */
  537. /* > -2 => N negative */
  538. /* > -3 => DIST illegal string */
  539. /* > -5 => SYM illegal string */
  540. /* > -7 => MODE not in range -6 to 6 */
  541. /* > -8 => COND less than 1.0, and MODE neither -6, 0 nor 6 */
  542. /* > -10 => KL negative */
  543. /* > -11 => KU negative, or SYM is not 'N' and KU is not equal to */
  544. /* > KL */
  545. /* > -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N'; */
  546. /* > or PACK='C' or 'Q' and SYM='N' and KL is not zero; */
  547. /* > or PACK='R' or 'B' and SYM='N' and KU is not zero; */
  548. /* > or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not */
  549. /* > N. */
  550. /* > -14 => LDA is less than M, or PACK='Z' and LDA is less than */
  551. /* > MIN(KU,N-1) + MIN(KL,M-1) + 1. */
  552. /* > 1 => Error return from SLATM7 */
  553. /* > 2 => Cannot scale to DMAX (f2cmax. sing. value is 0) */
  554. /* > 3 => Error return from CLAGGE, CLAGHE or CLAGSY */
  555. /* > \endverbatim */
  556. /* Authors: */
  557. /* ======== */
  558. /* > \author Univ. of Tennessee */
  559. /* > \author Univ. of California Berkeley */
  560. /* > \author Univ. of Colorado Denver */
  561. /* > \author NAG Ltd. */
  562. /* > \date December 2016 */
  563. /* > \ingroup complex_matgen */
  564. /* ===================================================================== */
  565. /* Subroutine */ void clatmt_(integer *m, integer *n, char *dist, integer *
  566. iseed, char *sym, real *d__, integer *mode, real *cond, real *dmax__,
  567. integer *rank, integer *kl, integer *ku, char *pack, complex *a,
  568. integer *lda, complex *work, integer *info)
  569. {
  570. /* System generated locals */
  571. integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5, i__6;
  572. real r__1, r__2, r__3;
  573. complex q__1, q__2, q__3;
  574. logical L__1;
  575. /* Local variables */
  576. integer ilda, icol;
  577. real temp;
  578. logical csym;
  579. integer irow, isym;
  580. complex c__;
  581. integer i__, j, k;
  582. complex s;
  583. real alpha, angle, realc;
  584. integer ipack, ioffg;
  585. extern logical lsame_(char *, char *);
  586. integer iinfo;
  587. extern /* Subroutine */ void sscal_(integer *, real *, real *, integer *);
  588. complex ctemp;
  589. integer idist, mnmin;
  590. complex extra;
  591. integer iskew;
  592. complex dummy;
  593. extern /* Subroutine */ void slatm7_(integer *, real *, integer *, integer
  594. *, integer *, real *, integer *, integer *, integer *);
  595. integer ic, jc, nc;
  596. extern /* Subroutine */ void clagge_(integer *, integer *, integer *,
  597. integer *, real *, complex *, integer *, integer *, complex *,
  598. integer *), claghe_(integer *, integer *, real *, complex *,
  599. integer *, integer *, complex *, integer *);
  600. integer il;
  601. complex ct;
  602. integer iendch, ir, jr, ipackg, mr;
  603. //extern /* Complex */ VOID clarnd_(complex *, integer *, integer *);
  604. extern complex clarnd_(integer *, integer *);
  605. integer minlda;
  606. complex st;
  607. extern /* Subroutine */ void claset_(char *, integer *, integer *, complex
  608. *, complex *, complex *, integer *), clartg_(complex *,
  609. complex *, real *, complex *, complex *);
  610. extern int xerbla_(char *, integer *, ftnlen);
  611. extern void clagsy_(integer *, integer *, real *, complex *,
  612. integer *, integer *, complex *, integer *);
  613. extern real slarnd_(integer *, integer *);
  614. extern /* Subroutine */ void clarot_(logical *, logical *, logical *,
  615. integer *, complex *, complex *, complex *, integer *, complex *,
  616. complex *);
  617. integer ioffst, irsign;
  618. logical givens, iltemp, ilextr, topdwn;
  619. integer ir1, ir2, isympk, jch, llb, jkl, jku, uub;
  620. /* -- LAPACK computational routine (version 3.7.0) -- */
  621. /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
  622. /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
  623. /* December 2016 */
  624. /* ===================================================================== */
  625. /* 1) Decode and Test the input parameters. */
  626. /* Initialize flags & seed. */
  627. /* Parameter adjustments */
  628. --iseed;
  629. --d__;
  630. a_dim1 = *lda;
  631. a_offset = 1 + a_dim1 * 1;
  632. a -= a_offset;
  633. --work;
  634. /* Function Body */
  635. *info = 0;
  636. /* Quick return if possible */
  637. if (*m == 0 || *n == 0) {
  638. return;
  639. }
  640. /* Decode DIST */
  641. if (lsame_(dist, "U")) {
  642. idist = 1;
  643. } else if (lsame_(dist, "S")) {
  644. idist = 2;
  645. } else if (lsame_(dist, "N")) {
  646. idist = 3;
  647. } else {
  648. idist = -1;
  649. }
  650. /* Decode SYM */
  651. if (lsame_(sym, "N")) {
  652. isym = 1;
  653. irsign = 0;
  654. csym = FALSE_;
  655. } else if (lsame_(sym, "P")) {
  656. isym = 2;
  657. irsign = 0;
  658. csym = FALSE_;
  659. } else if (lsame_(sym, "S")) {
  660. isym = 2;
  661. irsign = 0;
  662. csym = TRUE_;
  663. } else if (lsame_(sym, "H")) {
  664. isym = 2;
  665. irsign = 1;
  666. csym = FALSE_;
  667. } else {
  668. isym = -1;
  669. }
  670. /* Decode PACK */
  671. isympk = 0;
  672. if (lsame_(pack, "N")) {
  673. ipack = 0;
  674. } else if (lsame_(pack, "U")) {
  675. ipack = 1;
  676. isympk = 1;
  677. } else if (lsame_(pack, "L")) {
  678. ipack = 2;
  679. isympk = 1;
  680. } else if (lsame_(pack, "C")) {
  681. ipack = 3;
  682. isympk = 2;
  683. } else if (lsame_(pack, "R")) {
  684. ipack = 4;
  685. isympk = 3;
  686. } else if (lsame_(pack, "B")) {
  687. ipack = 5;
  688. isympk = 3;
  689. } else if (lsame_(pack, "Q")) {
  690. ipack = 6;
  691. isympk = 2;
  692. } else if (lsame_(pack, "Z")) {
  693. ipack = 7;
  694. } else {
  695. ipack = -1;
  696. }
  697. /* Set certain internal parameters */
  698. mnmin = f2cmin(*m,*n);
  699. /* Computing MIN */
  700. i__1 = *kl, i__2 = *m - 1;
  701. llb = f2cmin(i__1,i__2);
  702. /* Computing MIN */
  703. i__1 = *ku, i__2 = *n - 1;
  704. uub = f2cmin(i__1,i__2);
  705. /* Computing MIN */
  706. i__1 = *m, i__2 = *n + llb;
  707. mr = f2cmin(i__1,i__2);
  708. /* Computing MIN */
  709. i__1 = *n, i__2 = *m + uub;
  710. nc = f2cmin(i__1,i__2);
  711. if (ipack == 5 || ipack == 6) {
  712. minlda = uub + 1;
  713. } else if (ipack == 7) {
  714. minlda = llb + uub + 1;
  715. } else {
  716. minlda = *m;
  717. }
  718. /* Use Givens rotation method if bandwidth small enough, */
  719. /* or if LDA is too small to store the matrix unpacked. */
  720. givens = FALSE_;
  721. if (isym == 1) {
  722. /* Computing MAX */
  723. i__1 = 1, i__2 = mr + nc;
  724. if ((real) (llb + uub) < (real) f2cmax(i__1,i__2) * .3f) {
  725. givens = TRUE_;
  726. }
  727. } else {
  728. if (llb << 1 < *m) {
  729. givens = TRUE_;
  730. }
  731. }
  732. if (*lda < *m && *lda >= minlda) {
  733. givens = TRUE_;
  734. }
  735. /* Set INFO if an error */
  736. if (*m < 0) {
  737. *info = -1;
  738. } else if (*m != *n && isym != 1) {
  739. *info = -1;
  740. } else if (*n < 0) {
  741. *info = -2;
  742. } else if (idist == -1) {
  743. *info = -3;
  744. } else if (isym == -1) {
  745. *info = -5;
  746. } else if (abs(*mode) > 6) {
  747. *info = -7;
  748. } else if (*mode != 0 && abs(*mode) != 6 && *cond < 1.f) {
  749. *info = -8;
  750. } else if (*kl < 0) {
  751. *info = -10;
  752. } else if (*ku < 0 || isym != 1 && *kl != *ku) {
  753. *info = -11;
  754. } else if (ipack == -1 || isympk == 1 && isym == 1 || isympk == 2 && isym
  755. == 1 && *kl > 0 || isympk == 3 && isym == 1 && *ku > 0 || isympk
  756. != 0 && *m != *n) {
  757. *info = -12;
  758. } else if (*lda < f2cmax(1,minlda)) {
  759. *info = -14;
  760. }
  761. if (*info != 0) {
  762. i__1 = -(*info);
  763. xerbla_("CLATMT", &i__1, 6);
  764. return;
  765. }
  766. /* Initialize random number generator */
  767. for (i__ = 1; i__ <= 4; ++i__) {
  768. iseed[i__] = (i__1 = iseed[i__], abs(i__1)) % 4096;
  769. /* L100: */
  770. }
  771. if (iseed[4] % 2 != 1) {
  772. ++iseed[4];
  773. }
  774. /* 2) Set up D if indicated. */
  775. /* Compute D according to COND and MODE */
  776. slatm7_(mode, cond, &irsign, &idist, &iseed[1], &d__[1], &mnmin, rank, &
  777. iinfo);
  778. if (iinfo != 0) {
  779. *info = 1;
  780. return;
  781. }
  782. /* Choose Top-Down if D is (apparently) increasing, */
  783. /* Bottom-Up if D is (apparently) decreasing. */
  784. if (abs(d__[1]) <= (r__1 = d__[*rank], abs(r__1))) {
  785. topdwn = TRUE_;
  786. } else {
  787. topdwn = FALSE_;
  788. }
  789. if (*mode != 0 && abs(*mode) != 6) {
  790. /* Scale by DMAX */
  791. temp = abs(d__[1]);
  792. i__1 = *rank;
  793. for (i__ = 2; i__ <= i__1; ++i__) {
  794. /* Computing MAX */
  795. r__2 = temp, r__3 = (r__1 = d__[i__], abs(r__1));
  796. temp = f2cmax(r__2,r__3);
  797. /* L110: */
  798. }
  799. if (temp > 0.f) {
  800. alpha = *dmax__ / temp;
  801. } else {
  802. *info = 2;
  803. return;
  804. }
  805. sscal_(rank, &alpha, &d__[1], &c__1);
  806. }
  807. claset_("Full", lda, n, &c_b1, &c_b1, &a[a_offset], lda);
  808. /* 3) Generate Banded Matrix using Givens rotations. */
  809. /* Also the special case of UUB=LLB=0 */
  810. /* Compute Addressing constants to cover all */
  811. /* storage formats. Whether GE, HE, SY, GB, HB, or SB, */
  812. /* upper or lower triangle or both, */
  813. /* the (i,j)-th element is in */
  814. /* A( i - ISKEW*j + IOFFST, j ) */
  815. if (ipack > 4) {
  816. ilda = *lda - 1;
  817. iskew = 1;
  818. if (ipack > 5) {
  819. ioffst = uub + 1;
  820. } else {
  821. ioffst = 1;
  822. }
  823. } else {
  824. ilda = *lda;
  825. iskew = 0;
  826. ioffst = 0;
  827. }
  828. /* IPACKG is the format that the matrix is generated in. If this is */
  829. /* different from IPACK, then the matrix must be repacked at the */
  830. /* end. It also signals how to compute the norm, for scaling. */
  831. ipackg = 0;
  832. /* Diagonal Matrix -- We are done, unless it */
  833. /* is to be stored HP/SP/PP/TP (PACK='R' or 'C') */
  834. if (llb == 0 && uub == 0) {
  835. i__1 = mnmin;
  836. for (j = 1; j <= i__1; ++j) {
  837. i__2 = (1 - iskew) * j + ioffst + j * a_dim1;
  838. i__3 = j;
  839. q__1.r = d__[i__3], q__1.i = 0.f;
  840. a[i__2].r = q__1.r, a[i__2].i = q__1.i;
  841. /* L120: */
  842. }
  843. if (ipack <= 2 || ipack >= 5) {
  844. ipackg = ipack;
  845. }
  846. } else if (givens) {
  847. /* Check whether to use Givens rotations, */
  848. /* Householder transformations, or nothing. */
  849. if (isym == 1) {
  850. /* Non-symmetric -- A = U D V */
  851. if (ipack > 4) {
  852. ipackg = ipack;
  853. } else {
  854. ipackg = 0;
  855. }
  856. i__1 = mnmin;
  857. for (j = 1; j <= i__1; ++j) {
  858. i__2 = (1 - iskew) * j + ioffst + j * a_dim1;
  859. i__3 = j;
  860. q__1.r = d__[i__3], q__1.i = 0.f;
  861. a[i__2].r = q__1.r, a[i__2].i = q__1.i;
  862. /* L130: */
  863. }
  864. if (topdwn) {
  865. jkl = 0;
  866. i__1 = uub;
  867. for (jku = 1; jku <= i__1; ++jku) {
  868. /* Transform from bandwidth JKL, JKU-1 to JKL, JKU */
  869. /* Last row actually rotated is M */
  870. /* Last column actually rotated is MIN( M+JKU, N ) */
  871. /* Computing MIN */
  872. i__3 = *m + jku;
  873. i__2 = f2cmin(i__3,*n) + jkl - 1;
  874. for (jr = 1; jr <= i__2; ++jr) {
  875. extra.r = 0.f, extra.i = 0.f;
  876. angle = slarnd_(&c__1, &iseed[1]) *
  877. 6.2831853071795864769252867663f;
  878. r__1 = cos(angle);
  879. //clarnd_(&q__2, &c__5, &iseed[1]);
  880. q__2=clarnd_(&c__5, &iseed[1]);
  881. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  882. c__.r = q__1.r, c__.i = q__1.i;
  883. r__1 = sin(angle);
  884. //clarnd_(&q__2, &c__5, &iseed[1]);
  885. q__2=clarnd_(&c__5, &iseed[1]);
  886. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  887. s.r = q__1.r, s.i = q__1.i;
  888. /* Computing MAX */
  889. i__3 = 1, i__4 = jr - jkl;
  890. icol = f2cmax(i__3,i__4);
  891. if (jr < *m) {
  892. /* Computing MIN */
  893. i__3 = *n, i__4 = jr + jku;
  894. il = f2cmin(i__3,i__4) + 1 - icol;
  895. L__1 = jr > jkl;
  896. clarot_(&c_true, &L__1, &c_false, &il, &c__, &s, &
  897. a[jr - iskew * icol + ioffst + icol *
  898. a_dim1], &ilda, &extra, &dummy);
  899. }
  900. /* Chase "EXTRA" back up */
  901. ir = jr;
  902. ic = icol;
  903. i__3 = -jkl - jku;
  904. for (jch = jr - jkl; i__3 < 0 ? jch >= 1 : jch <= 1;
  905. jch += i__3) {
  906. if (ir < *m) {
  907. clartg_(&a[ir + 1 - iskew * (ic + 1) + ioffst
  908. + (ic + 1) * a_dim1], &extra, &realc,
  909. &s, &dummy);
  910. //clarnd_(&q__1, &c__5, &iseed[1]);
  911. q__1=clarnd_(&c__5, &iseed[1]);
  912. dummy.r = q__1.r, dummy.i = q__1.i;
  913. q__2.r = realc * dummy.r, q__2.i = realc *
  914. dummy.i;
  915. r_cnjg(&q__1, &q__2);
  916. c__.r = q__1.r, c__.i = q__1.i;
  917. q__3.r = -s.r, q__3.i = -s.i;
  918. q__2.r = q__3.r * dummy.r - q__3.i * dummy.i,
  919. q__2.i = q__3.r * dummy.i + q__3.i *
  920. dummy.r;
  921. r_cnjg(&q__1, &q__2);
  922. s.r = q__1.r, s.i = q__1.i;
  923. }
  924. /* Computing MAX */
  925. i__4 = 1, i__5 = jch - jku;
  926. irow = f2cmax(i__4,i__5);
  927. il = ir + 2 - irow;
  928. ctemp.r = 0.f, ctemp.i = 0.f;
  929. iltemp = jch > jku;
  930. clarot_(&c_false, &iltemp, &c_true, &il, &c__, &s,
  931. &a[irow - iskew * ic + ioffst + ic *
  932. a_dim1], &ilda, &ctemp, &extra);
  933. if (iltemp) {
  934. clartg_(&a[irow + 1 - iskew * (ic + 1) +
  935. ioffst + (ic + 1) * a_dim1], &ctemp, &
  936. realc, &s, &dummy);
  937. //clarnd_(&q__1, &c__5, &iseed[1]);
  938. q__1=clarnd_(&c__5, &iseed[1]);
  939. dummy.r = q__1.r, dummy.i = q__1.i;
  940. q__2.r = realc * dummy.r, q__2.i = realc *
  941. dummy.i;
  942. r_cnjg(&q__1, &q__2);
  943. c__.r = q__1.r, c__.i = q__1.i;
  944. q__3.r = -s.r, q__3.i = -s.i;
  945. q__2.r = q__3.r * dummy.r - q__3.i * dummy.i,
  946. q__2.i = q__3.r * dummy.i + q__3.i *
  947. dummy.r;
  948. r_cnjg(&q__1, &q__2);
  949. s.r = q__1.r, s.i = q__1.i;
  950. /* Computing MAX */
  951. i__4 = 1, i__5 = jch - jku - jkl;
  952. icol = f2cmax(i__4,i__5);
  953. il = ic + 2 - icol;
  954. extra.r = 0.f, extra.i = 0.f;
  955. L__1 = jch > jku + jkl;
  956. clarot_(&c_true, &L__1, &c_true, &il, &c__, &
  957. s, &a[irow - iskew * icol + ioffst +
  958. icol * a_dim1], &ilda, &extra, &ctemp)
  959. ;
  960. ic = icol;
  961. ir = irow;
  962. }
  963. /* L140: */
  964. }
  965. /* L150: */
  966. }
  967. /* L160: */
  968. }
  969. jku = uub;
  970. i__1 = llb;
  971. for (jkl = 1; jkl <= i__1; ++jkl) {
  972. /* Transform from bandwidth JKL-1, JKU to JKL, JKU */
  973. /* Computing MIN */
  974. i__3 = *n + jkl;
  975. i__2 = f2cmin(i__3,*m) + jku - 1;
  976. for (jc = 1; jc <= i__2; ++jc) {
  977. extra.r = 0.f, extra.i = 0.f;
  978. angle = slarnd_(&c__1, &iseed[1]) *
  979. 6.2831853071795864769252867663f;
  980. r__1 = cos(angle);
  981. //clarnd_(&q__2, &c__5, &iseed[1]);
  982. q__2=clarnd_(&c__5, &iseed[1]);
  983. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  984. c__.r = q__1.r, c__.i = q__1.i;
  985. r__1 = sin(angle);
  986. //clarnd_(&q__2, &c__5, &iseed[1]);
  987. q__2=clarnd_(&c__5, &iseed[1]);
  988. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  989. s.r = q__1.r, s.i = q__1.i;
  990. /* Computing MAX */
  991. i__3 = 1, i__4 = jc - jku;
  992. irow = f2cmax(i__3,i__4);
  993. if (jc < *n) {
  994. /* Computing MIN */
  995. i__3 = *m, i__4 = jc + jkl;
  996. il = f2cmin(i__3,i__4) + 1 - irow;
  997. L__1 = jc > jku;
  998. clarot_(&c_false, &L__1, &c_false, &il, &c__, &s,
  999. &a[irow - iskew * jc + ioffst + jc *
  1000. a_dim1], &ilda, &extra, &dummy);
  1001. }
  1002. /* Chase "EXTRA" back up */
  1003. ic = jc;
  1004. ir = irow;
  1005. i__3 = -jkl - jku;
  1006. for (jch = jc - jku; i__3 < 0 ? jch >= 1 : jch <= 1;
  1007. jch += i__3) {
  1008. if (ic < *n) {
  1009. clartg_(&a[ir + 1 - iskew * (ic + 1) + ioffst
  1010. + (ic + 1) * a_dim1], &extra, &realc,
  1011. &s, &dummy);
  1012. //clarnd_(&q__1, &c__5, &iseed[1]);
  1013. q__1=clarnd_(&c__5, &iseed[1]);
  1014. dummy.r = q__1.r, dummy.i = q__1.i;
  1015. q__2.r = realc * dummy.r, q__2.i = realc *
  1016. dummy.i;
  1017. r_cnjg(&q__1, &q__2);
  1018. c__.r = q__1.r, c__.i = q__1.i;
  1019. q__3.r = -s.r, q__3.i = -s.i;
  1020. q__2.r = q__3.r * dummy.r - q__3.i * dummy.i,
  1021. q__2.i = q__3.r * dummy.i + q__3.i *
  1022. dummy.r;
  1023. r_cnjg(&q__1, &q__2);
  1024. s.r = q__1.r, s.i = q__1.i;
  1025. }
  1026. /* Computing MAX */
  1027. i__4 = 1, i__5 = jch - jkl;
  1028. icol = f2cmax(i__4,i__5);
  1029. il = ic + 2 - icol;
  1030. ctemp.r = 0.f, ctemp.i = 0.f;
  1031. iltemp = jch > jkl;
  1032. clarot_(&c_true, &iltemp, &c_true, &il, &c__, &s,
  1033. &a[ir - iskew * icol + ioffst + icol *
  1034. a_dim1], &ilda, &ctemp, &extra);
  1035. if (iltemp) {
  1036. clartg_(&a[ir + 1 - iskew * (icol + 1) +
  1037. ioffst + (icol + 1) * a_dim1], &ctemp,
  1038. &realc, &s, &dummy);
  1039. //clarnd_(&q__1, &c__5, &iseed[1]);
  1040. q__1=clarnd_(&c__5, &iseed[1]);
  1041. dummy.r = q__1.r, dummy.i = q__1.i;
  1042. q__2.r = realc * dummy.r, q__2.i = realc *
  1043. dummy.i;
  1044. r_cnjg(&q__1, &q__2);
  1045. c__.r = q__1.r, c__.i = q__1.i;
  1046. q__3.r = -s.r, q__3.i = -s.i;
  1047. q__2.r = q__3.r * dummy.r - q__3.i * dummy.i,
  1048. q__2.i = q__3.r * dummy.i + q__3.i *
  1049. dummy.r;
  1050. r_cnjg(&q__1, &q__2);
  1051. s.r = q__1.r, s.i = q__1.i;
  1052. /* Computing MAX */
  1053. i__4 = 1, i__5 = jch - jkl - jku;
  1054. irow = f2cmax(i__4,i__5);
  1055. il = ir + 2 - irow;
  1056. extra.r = 0.f, extra.i = 0.f;
  1057. L__1 = jch > jkl + jku;
  1058. clarot_(&c_false, &L__1, &c_true, &il, &c__, &
  1059. s, &a[irow - iskew * icol + ioffst +
  1060. icol * a_dim1], &ilda, &extra, &ctemp)
  1061. ;
  1062. ic = icol;
  1063. ir = irow;
  1064. }
  1065. /* L170: */
  1066. }
  1067. /* L180: */
  1068. }
  1069. /* L190: */
  1070. }
  1071. } else {
  1072. /* Bottom-Up -- Start at the bottom right. */
  1073. jkl = 0;
  1074. i__1 = uub;
  1075. for (jku = 1; jku <= i__1; ++jku) {
  1076. /* Transform from bandwidth JKL, JKU-1 to JKL, JKU */
  1077. /* First row actually rotated is M */
  1078. /* First column actually rotated is MIN( M+JKU, N ) */
  1079. /* Computing MIN */
  1080. i__2 = *m, i__3 = *n + jkl;
  1081. iendch = f2cmin(i__2,i__3) - 1;
  1082. /* Computing MIN */
  1083. i__2 = *m + jku;
  1084. i__3 = 1 - jkl;
  1085. for (jc = f2cmin(i__2,*n) - 1; jc >= i__3; --jc) {
  1086. extra.r = 0.f, extra.i = 0.f;
  1087. angle = slarnd_(&c__1, &iseed[1]) *
  1088. 6.2831853071795864769252867663f;
  1089. r__1 = cos(angle);
  1090. //clarnd_(&q__2, &c__5, &iseed[1]);
  1091. q__2=clarnd_(&c__5, &iseed[1]);
  1092. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1093. c__.r = q__1.r, c__.i = q__1.i;
  1094. r__1 = sin(angle);
  1095. //clarnd_(&q__2, &c__5, &iseed[1]);
  1096. q__2=clarnd_(&c__5, &iseed[1]);
  1097. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1098. s.r = q__1.r, s.i = q__1.i;
  1099. /* Computing MAX */
  1100. i__2 = 1, i__4 = jc - jku + 1;
  1101. irow = f2cmax(i__2,i__4);
  1102. if (jc > 0) {
  1103. /* Computing MIN */
  1104. i__2 = *m, i__4 = jc + jkl + 1;
  1105. il = f2cmin(i__2,i__4) + 1 - irow;
  1106. L__1 = jc + jkl < *m;
  1107. clarot_(&c_false, &c_false, &L__1, &il, &c__, &s,
  1108. &a[irow - iskew * jc + ioffst + jc *
  1109. a_dim1], &ilda, &dummy, &extra);
  1110. }
  1111. /* Chase "EXTRA" back down */
  1112. ic = jc;
  1113. i__2 = iendch;
  1114. i__4 = jkl + jku;
  1115. for (jch = jc + jkl; i__4 < 0 ? jch >= i__2 : jch <=
  1116. i__2; jch += i__4) {
  1117. ilextr = ic > 0;
  1118. if (ilextr) {
  1119. clartg_(&a[jch - iskew * ic + ioffst + ic *
  1120. a_dim1], &extra, &realc, &s, &dummy);
  1121. //clarnd_(&q__1, &c__5, &iseed[1]);
  1122. q__1=clarnd_(&c__5, &iseed[1]);
  1123. dummy.r = q__1.r, dummy.i = q__1.i;
  1124. q__1.r = realc * dummy.r, q__1.i = realc *
  1125. dummy.i;
  1126. c__.r = q__1.r, c__.i = q__1.i;
  1127. q__1.r = s.r * dummy.r - s.i * dummy.i,
  1128. q__1.i = s.r * dummy.i + s.i *
  1129. dummy.r;
  1130. s.r = q__1.r, s.i = q__1.i;
  1131. }
  1132. ic = f2cmax(1,ic);
  1133. /* Computing MIN */
  1134. i__5 = *n - 1, i__6 = jch + jku;
  1135. icol = f2cmin(i__5,i__6);
  1136. iltemp = jch + jku < *n;
  1137. ctemp.r = 0.f, ctemp.i = 0.f;
  1138. i__5 = icol + 2 - ic;
  1139. clarot_(&c_true, &ilextr, &iltemp, &i__5, &c__, &
  1140. s, &a[jch - iskew * ic + ioffst + ic *
  1141. a_dim1], &ilda, &extra, &ctemp);
  1142. if (iltemp) {
  1143. clartg_(&a[jch - iskew * icol + ioffst + icol
  1144. * a_dim1], &ctemp, &realc, &s, &dummy)
  1145. ;
  1146. //clarnd_(&q__1, &c__5, &iseed[1]);
  1147. q__1=clarnd_(&c__5, &iseed[1]);
  1148. dummy.r = q__1.r, dummy.i = q__1.i;
  1149. q__1.r = realc * dummy.r, q__1.i = realc *
  1150. dummy.i;
  1151. c__.r = q__1.r, c__.i = q__1.i;
  1152. q__1.r = s.r * dummy.r - s.i * dummy.i,
  1153. q__1.i = s.r * dummy.i + s.i *
  1154. dummy.r;
  1155. s.r = q__1.r, s.i = q__1.i;
  1156. /* Computing MIN */
  1157. i__5 = iendch, i__6 = jch + jkl + jku;
  1158. il = f2cmin(i__5,i__6) + 2 - jch;
  1159. extra.r = 0.f, extra.i = 0.f;
  1160. L__1 = jch + jkl + jku <= iendch;
  1161. clarot_(&c_false, &c_true, &L__1, &il, &c__, &
  1162. s, &a[jch - iskew * icol + ioffst +
  1163. icol * a_dim1], &ilda, &ctemp, &extra)
  1164. ;
  1165. ic = icol;
  1166. }
  1167. /* L200: */
  1168. }
  1169. /* L210: */
  1170. }
  1171. /* L220: */
  1172. }
  1173. jku = uub;
  1174. i__1 = llb;
  1175. for (jkl = 1; jkl <= i__1; ++jkl) {
  1176. /* Transform from bandwidth JKL-1, JKU to JKL, JKU */
  1177. /* First row actually rotated is MIN( N+JKL, M ) */
  1178. /* First column actually rotated is N */
  1179. /* Computing MIN */
  1180. i__3 = *n, i__4 = *m + jku;
  1181. iendch = f2cmin(i__3,i__4) - 1;
  1182. /* Computing MIN */
  1183. i__3 = *n + jkl;
  1184. i__4 = 1 - jku;
  1185. for (jr = f2cmin(i__3,*m) - 1; jr >= i__4; --jr) {
  1186. extra.r = 0.f, extra.i = 0.f;
  1187. angle = slarnd_(&c__1, &iseed[1]) *
  1188. 6.2831853071795864769252867663f;
  1189. r__1 = cos(angle);
  1190. //clarnd_(&q__2, &c__5, &iseed[1]);
  1191. q__2=clarnd_(&c__5, &iseed[1]);
  1192. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1193. c__.r = q__1.r, c__.i = q__1.i;
  1194. r__1 = sin(angle);
  1195. //clarnd_(&q__2, &c__5, &iseed[1]);
  1196. q__2=clarnd_(&c__5, &iseed[1]);
  1197. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1198. s.r = q__1.r, s.i = q__1.i;
  1199. /* Computing MAX */
  1200. i__3 = 1, i__2 = jr - jkl + 1;
  1201. icol = f2cmax(i__3,i__2);
  1202. if (jr > 0) {
  1203. /* Computing MIN */
  1204. i__3 = *n, i__2 = jr + jku + 1;
  1205. il = f2cmin(i__3,i__2) + 1 - icol;
  1206. L__1 = jr + jku < *n;
  1207. clarot_(&c_true, &c_false, &L__1, &il, &c__, &s, &
  1208. a[jr - iskew * icol + ioffst + icol *
  1209. a_dim1], &ilda, &dummy, &extra);
  1210. }
  1211. /* Chase "EXTRA" back down */
  1212. ir = jr;
  1213. i__3 = iendch;
  1214. i__2 = jkl + jku;
  1215. for (jch = jr + jku; i__2 < 0 ? jch >= i__3 : jch <=
  1216. i__3; jch += i__2) {
  1217. ilextr = ir > 0;
  1218. if (ilextr) {
  1219. clartg_(&a[ir - iskew * jch + ioffst + jch *
  1220. a_dim1], &extra, &realc, &s, &dummy);
  1221. //clarnd_(&q__1, &c__5, &iseed[1]);
  1222. q__1=clarnd_(&c__5, &iseed[1]);
  1223. dummy.r = q__1.r, dummy.i = q__1.i;
  1224. q__1.r = realc * dummy.r, q__1.i = realc *
  1225. dummy.i;
  1226. c__.r = q__1.r, c__.i = q__1.i;
  1227. q__1.r = s.r * dummy.r - s.i * dummy.i,
  1228. q__1.i = s.r * dummy.i + s.i *
  1229. dummy.r;
  1230. s.r = q__1.r, s.i = q__1.i;
  1231. }
  1232. ir = f2cmax(1,ir);
  1233. /* Computing MIN */
  1234. i__5 = *m - 1, i__6 = jch + jkl;
  1235. irow = f2cmin(i__5,i__6);
  1236. iltemp = jch + jkl < *m;
  1237. ctemp.r = 0.f, ctemp.i = 0.f;
  1238. i__5 = irow + 2 - ir;
  1239. clarot_(&c_false, &ilextr, &iltemp, &i__5, &c__, &
  1240. s, &a[ir - iskew * jch + ioffst + jch *
  1241. a_dim1], &ilda, &extra, &ctemp);
  1242. if (iltemp) {
  1243. clartg_(&a[irow - iskew * jch + ioffst + jch *
  1244. a_dim1], &ctemp, &realc, &s, &dummy);
  1245. //clarnd_(&q__1, &c__5, &iseed[1]);
  1246. q__1=clarnd_(&c__5, &iseed[1]);
  1247. dummy.r = q__1.r, dummy.i = q__1.i;
  1248. q__1.r = realc * dummy.r, q__1.i = realc *
  1249. dummy.i;
  1250. c__.r = q__1.r, c__.i = q__1.i;
  1251. q__1.r = s.r * dummy.r - s.i * dummy.i,
  1252. q__1.i = s.r * dummy.i + s.i *
  1253. dummy.r;
  1254. s.r = q__1.r, s.i = q__1.i;
  1255. /* Computing MIN */
  1256. i__5 = iendch, i__6 = jch + jkl + jku;
  1257. il = f2cmin(i__5,i__6) + 2 - jch;
  1258. extra.r = 0.f, extra.i = 0.f;
  1259. L__1 = jch + jkl + jku <= iendch;
  1260. clarot_(&c_true, &c_true, &L__1, &il, &c__, &
  1261. s, &a[irow - iskew * jch + ioffst +
  1262. jch * a_dim1], &ilda, &ctemp, &extra);
  1263. ir = irow;
  1264. }
  1265. /* L230: */
  1266. }
  1267. /* L240: */
  1268. }
  1269. /* L250: */
  1270. }
  1271. }
  1272. } else {
  1273. /* Symmetric -- A = U D U' */
  1274. /* Hermitian -- A = U D U* */
  1275. ipackg = ipack;
  1276. ioffg = ioffst;
  1277. if (topdwn) {
  1278. /* Top-Down -- Generate Upper triangle only */
  1279. if (ipack >= 5) {
  1280. ipackg = 6;
  1281. ioffg = uub + 1;
  1282. } else {
  1283. ipackg = 1;
  1284. }
  1285. i__1 = mnmin;
  1286. for (j = 1; j <= i__1; ++j) {
  1287. i__4 = (1 - iskew) * j + ioffg + j * a_dim1;
  1288. i__2 = j;
  1289. q__1.r = d__[i__2], q__1.i = 0.f;
  1290. a[i__4].r = q__1.r, a[i__4].i = q__1.i;
  1291. /* L260: */
  1292. }
  1293. i__1 = uub;
  1294. for (k = 1; k <= i__1; ++k) {
  1295. i__4 = *n - 1;
  1296. for (jc = 1; jc <= i__4; ++jc) {
  1297. /* Computing MAX */
  1298. i__2 = 1, i__3 = jc - k;
  1299. irow = f2cmax(i__2,i__3);
  1300. /* Computing MIN */
  1301. i__2 = jc + 1, i__3 = k + 2;
  1302. il = f2cmin(i__2,i__3);
  1303. extra.r = 0.f, extra.i = 0.f;
  1304. i__2 = jc - iskew * (jc + 1) + ioffg + (jc + 1) *
  1305. a_dim1;
  1306. ctemp.r = a[i__2].r, ctemp.i = a[i__2].i;
  1307. angle = slarnd_(&c__1, &iseed[1]) *
  1308. 6.2831853071795864769252867663f;
  1309. r__1 = cos(angle);
  1310. //clarnd_(&q__2, &c__5, &iseed[1]);
  1311. q__2=clarnd_(&c__5, &iseed[1]);
  1312. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1313. c__.r = q__1.r, c__.i = q__1.i;
  1314. r__1 = sin(angle);
  1315. //clarnd_(&q__2, &c__5, &iseed[1]);
  1316. q__2=clarnd_(&c__5, &iseed[1]);
  1317. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1318. s.r = q__1.r, s.i = q__1.i;
  1319. if (csym) {
  1320. ct.r = c__.r, ct.i = c__.i;
  1321. st.r = s.r, st.i = s.i;
  1322. } else {
  1323. r_cnjg(&q__1, &ctemp);
  1324. ctemp.r = q__1.r, ctemp.i = q__1.i;
  1325. r_cnjg(&q__1, &c__);
  1326. ct.r = q__1.r, ct.i = q__1.i;
  1327. r_cnjg(&q__1, &s);
  1328. st.r = q__1.r, st.i = q__1.i;
  1329. }
  1330. L__1 = jc > k;
  1331. clarot_(&c_false, &L__1, &c_true, &il, &c__, &s, &a[
  1332. irow - iskew * jc + ioffg + jc * a_dim1], &
  1333. ilda, &extra, &ctemp);
  1334. /* Computing MIN */
  1335. i__3 = k, i__5 = *n - jc;
  1336. i__2 = f2cmin(i__3,i__5) + 1;
  1337. clarot_(&c_true, &c_true, &c_false, &i__2, &ct, &st, &
  1338. a[(1 - iskew) * jc + ioffg + jc * a_dim1], &
  1339. ilda, &ctemp, &dummy);
  1340. /* Chase EXTRA back up the matrix */
  1341. icol = jc;
  1342. i__2 = -k;
  1343. for (jch = jc - k; i__2 < 0 ? jch >= 1 : jch <= 1;
  1344. jch += i__2) {
  1345. clartg_(&a[jch + 1 - iskew * (icol + 1) + ioffg +
  1346. (icol + 1) * a_dim1], &extra, &realc, &s,
  1347. &dummy);
  1348. //clarnd_(&q__1, &c__5, &iseed[1]);
  1349. q__1=clarnd_(&c__5, &iseed[1]);
  1350. dummy.r = q__1.r, dummy.i = q__1.i;
  1351. q__2.r = realc * dummy.r, q__2.i = realc *
  1352. dummy.i;
  1353. r_cnjg(&q__1, &q__2);
  1354. c__.r = q__1.r, c__.i = q__1.i;
  1355. q__3.r = -s.r, q__3.i = -s.i;
  1356. q__2.r = q__3.r * dummy.r - q__3.i * dummy.i,
  1357. q__2.i = q__3.r * dummy.i + q__3.i *
  1358. dummy.r;
  1359. r_cnjg(&q__1, &q__2);
  1360. s.r = q__1.r, s.i = q__1.i;
  1361. i__3 = jch - iskew * (jch + 1) + ioffg + (jch + 1)
  1362. * a_dim1;
  1363. ctemp.r = a[i__3].r, ctemp.i = a[i__3].i;
  1364. if (csym) {
  1365. ct.r = c__.r, ct.i = c__.i;
  1366. st.r = s.r, st.i = s.i;
  1367. } else {
  1368. r_cnjg(&q__1, &ctemp);
  1369. ctemp.r = q__1.r, ctemp.i = q__1.i;
  1370. r_cnjg(&q__1, &c__);
  1371. ct.r = q__1.r, ct.i = q__1.i;
  1372. r_cnjg(&q__1, &s);
  1373. st.r = q__1.r, st.i = q__1.i;
  1374. }
  1375. i__3 = k + 2;
  1376. clarot_(&c_true, &c_true, &c_true, &i__3, &c__, &
  1377. s, &a[(1 - iskew) * jch + ioffg + jch *
  1378. a_dim1], &ilda, &ctemp, &extra);
  1379. /* Computing MAX */
  1380. i__3 = 1, i__5 = jch - k;
  1381. irow = f2cmax(i__3,i__5);
  1382. /* Computing MIN */
  1383. i__3 = jch + 1, i__5 = k + 2;
  1384. il = f2cmin(i__3,i__5);
  1385. extra.r = 0.f, extra.i = 0.f;
  1386. L__1 = jch > k;
  1387. clarot_(&c_false, &L__1, &c_true, &il, &ct, &st, &
  1388. a[irow - iskew * jch + ioffg + jch *
  1389. a_dim1], &ilda, &extra, &ctemp);
  1390. icol = jch;
  1391. /* L270: */
  1392. }
  1393. /* L280: */
  1394. }
  1395. /* L290: */
  1396. }
  1397. /* If we need lower triangle, copy from upper. Note that */
  1398. /* the order of copying is chosen to work for 'q' -> 'b' */
  1399. if (ipack != ipackg && ipack != 3) {
  1400. i__1 = *n;
  1401. for (jc = 1; jc <= i__1; ++jc) {
  1402. irow = ioffst - iskew * jc;
  1403. if (csym) {
  1404. /* Computing MIN */
  1405. i__2 = *n, i__3 = jc + uub;
  1406. i__4 = f2cmin(i__2,i__3);
  1407. for (jr = jc; jr <= i__4; ++jr) {
  1408. i__2 = jr + irow + jc * a_dim1;
  1409. i__3 = jc - iskew * jr + ioffg + jr * a_dim1;
  1410. a[i__2].r = a[i__3].r, a[i__2].i = a[i__3].i;
  1411. /* L300: */
  1412. }
  1413. } else {
  1414. /* Computing MIN */
  1415. i__2 = *n, i__3 = jc + uub;
  1416. i__4 = f2cmin(i__2,i__3);
  1417. for (jr = jc; jr <= i__4; ++jr) {
  1418. i__2 = jr + irow + jc * a_dim1;
  1419. r_cnjg(&q__1, &a[jc - iskew * jr + ioffg + jr
  1420. * a_dim1]);
  1421. a[i__2].r = q__1.r, a[i__2].i = q__1.i;
  1422. /* L310: */
  1423. }
  1424. }
  1425. /* L320: */
  1426. }
  1427. if (ipack == 5) {
  1428. i__1 = *n;
  1429. for (jc = *n - uub + 1; jc <= i__1; ++jc) {
  1430. i__4 = uub + 1;
  1431. for (jr = *n + 2 - jc; jr <= i__4; ++jr) {
  1432. i__2 = jr + jc * a_dim1;
  1433. a[i__2].r = 0.f, a[i__2].i = 0.f;
  1434. /* L330: */
  1435. }
  1436. /* L340: */
  1437. }
  1438. }
  1439. if (ipackg == 6) {
  1440. ipackg = ipack;
  1441. } else {
  1442. ipackg = 0;
  1443. }
  1444. }
  1445. } else {
  1446. /* Bottom-Up -- Generate Lower triangle only */
  1447. if (ipack >= 5) {
  1448. ipackg = 5;
  1449. if (ipack == 6) {
  1450. ioffg = 1;
  1451. }
  1452. } else {
  1453. ipackg = 2;
  1454. }
  1455. i__1 = mnmin;
  1456. for (j = 1; j <= i__1; ++j) {
  1457. i__4 = (1 - iskew) * j + ioffg + j * a_dim1;
  1458. i__2 = j;
  1459. q__1.r = d__[i__2], q__1.i = 0.f;
  1460. a[i__4].r = q__1.r, a[i__4].i = q__1.i;
  1461. /* L350: */
  1462. }
  1463. i__1 = uub;
  1464. for (k = 1; k <= i__1; ++k) {
  1465. for (jc = *n - 1; jc >= 1; --jc) {
  1466. /* Computing MIN */
  1467. i__4 = *n + 1 - jc, i__2 = k + 2;
  1468. il = f2cmin(i__4,i__2);
  1469. extra.r = 0.f, extra.i = 0.f;
  1470. i__4 = (1 - iskew) * jc + 1 + ioffg + jc * a_dim1;
  1471. ctemp.r = a[i__4].r, ctemp.i = a[i__4].i;
  1472. angle = slarnd_(&c__1, &iseed[1]) *
  1473. 6.2831853071795864769252867663f;
  1474. r__1 = cos(angle);
  1475. //clarnd_(&q__2, &c__5, &iseed[1]);
  1476. q__2=clarnd_(&c__5, &iseed[1]);
  1477. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1478. c__.r = q__1.r, c__.i = q__1.i;
  1479. r__1 = sin(angle);
  1480. //clarnd_(&q__2, &c__5, &iseed[1]);
  1481. q__2=clarnd_(&c__5, &iseed[1]);
  1482. q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
  1483. s.r = q__1.r, s.i = q__1.i;
  1484. if (csym) {
  1485. ct.r = c__.r, ct.i = c__.i;
  1486. st.r = s.r, st.i = s.i;
  1487. } else {
  1488. r_cnjg(&q__1, &ctemp);
  1489. ctemp.r = q__1.r, ctemp.i = q__1.i;
  1490. r_cnjg(&q__1, &c__);
  1491. ct.r = q__1.r, ct.i = q__1.i;
  1492. r_cnjg(&q__1, &s);
  1493. st.r = q__1.r, st.i = q__1.i;
  1494. }
  1495. L__1 = *n - jc > k;
  1496. clarot_(&c_false, &c_true, &L__1, &il, &c__, &s, &a[(
  1497. 1 - iskew) * jc + ioffg + jc * a_dim1], &ilda,
  1498. &ctemp, &extra);
  1499. /* Computing MAX */
  1500. i__4 = 1, i__2 = jc - k + 1;
  1501. icol = f2cmax(i__4,i__2);
  1502. i__4 = jc + 2 - icol;
  1503. clarot_(&c_true, &c_false, &c_true, &i__4, &ct, &st, &
  1504. a[jc - iskew * icol + ioffg + icol * a_dim1],
  1505. &ilda, &dummy, &ctemp);
  1506. /* Chase EXTRA back down the matrix */
  1507. icol = jc;
  1508. i__4 = *n - 1;
  1509. i__2 = k;
  1510. for (jch = jc + k; i__2 < 0 ? jch >= i__4 : jch <=
  1511. i__4; jch += i__2) {
  1512. clartg_(&a[jch - iskew * icol + ioffg + icol *
  1513. a_dim1], &extra, &realc, &s, &dummy);
  1514. //clarnd_(&q__1, &c__5, &iseed[1]);
  1515. q__1=clarnd_(&c__5, &iseed[1]);
  1516. dummy.r = q__1.r, dummy.i = q__1.i;
  1517. q__1.r = realc * dummy.r, q__1.i = realc *
  1518. dummy.i;
  1519. c__.r = q__1.r, c__.i = q__1.i;
  1520. q__1.r = s.r * dummy.r - s.i * dummy.i, q__1.i =
  1521. s.r * dummy.i + s.i * dummy.r;
  1522. s.r = q__1.r, s.i = q__1.i;
  1523. i__3 = (1 - iskew) * jch + 1 + ioffg + jch *
  1524. a_dim1;
  1525. ctemp.r = a[i__3].r, ctemp.i = a[i__3].i;
  1526. if (csym) {
  1527. ct.r = c__.r, ct.i = c__.i;
  1528. st.r = s.r, st.i = s.i;
  1529. } else {
  1530. r_cnjg(&q__1, &ctemp);
  1531. ctemp.r = q__1.r, ctemp.i = q__1.i;
  1532. r_cnjg(&q__1, &c__);
  1533. ct.r = q__1.r, ct.i = q__1.i;
  1534. r_cnjg(&q__1, &s);
  1535. st.r = q__1.r, st.i = q__1.i;
  1536. }
  1537. i__3 = k + 2;
  1538. clarot_(&c_true, &c_true, &c_true, &i__3, &c__, &
  1539. s, &a[jch - iskew * icol + ioffg + icol *
  1540. a_dim1], &ilda, &extra, &ctemp);
  1541. /* Computing MIN */
  1542. i__3 = *n + 1 - jch, i__5 = k + 2;
  1543. il = f2cmin(i__3,i__5);
  1544. extra.r = 0.f, extra.i = 0.f;
  1545. L__1 = *n - jch > k;
  1546. clarot_(&c_false, &c_true, &L__1, &il, &ct, &st, &
  1547. a[(1 - iskew) * jch + ioffg + jch *
  1548. a_dim1], &ilda, &ctemp, &extra);
  1549. icol = jch;
  1550. /* L360: */
  1551. }
  1552. /* L370: */
  1553. }
  1554. /* L380: */
  1555. }
  1556. /* If we need upper triangle, copy from lower. Note that */
  1557. /* the order of copying is chosen to work for 'b' -> 'q' */
  1558. if (ipack != ipackg && ipack != 4) {
  1559. for (jc = *n; jc >= 1; --jc) {
  1560. irow = ioffst - iskew * jc;
  1561. if (csym) {
  1562. /* Computing MAX */
  1563. i__2 = 1, i__4 = jc - uub;
  1564. i__1 = f2cmax(i__2,i__4);
  1565. for (jr = jc; jr >= i__1; --jr) {
  1566. i__2 = jr + irow + jc * a_dim1;
  1567. i__4 = jc - iskew * jr + ioffg + jr * a_dim1;
  1568. a[i__2].r = a[i__4].r, a[i__2].i = a[i__4].i;
  1569. /* L390: */
  1570. }
  1571. } else {
  1572. /* Computing MAX */
  1573. i__2 = 1, i__4 = jc - uub;
  1574. i__1 = f2cmax(i__2,i__4);
  1575. for (jr = jc; jr >= i__1; --jr) {
  1576. i__2 = jr + irow + jc * a_dim1;
  1577. r_cnjg(&q__1, &a[jc - iskew * jr + ioffg + jr
  1578. * a_dim1]);
  1579. a[i__2].r = q__1.r, a[i__2].i = q__1.i;
  1580. /* L400: */
  1581. }
  1582. }
  1583. /* L410: */
  1584. }
  1585. if (ipack == 6) {
  1586. i__1 = uub;
  1587. for (jc = 1; jc <= i__1; ++jc) {
  1588. i__2 = uub + 1 - jc;
  1589. for (jr = 1; jr <= i__2; ++jr) {
  1590. i__4 = jr + jc * a_dim1;
  1591. a[i__4].r = 0.f, a[i__4].i = 0.f;
  1592. /* L420: */
  1593. }
  1594. /* L430: */
  1595. }
  1596. }
  1597. if (ipackg == 5) {
  1598. ipackg = ipack;
  1599. } else {
  1600. ipackg = 0;
  1601. }
  1602. }
  1603. }
  1604. /* Ensure that the diagonal is real if Hermitian */
  1605. if (! csym) {
  1606. i__1 = *n;
  1607. for (jc = 1; jc <= i__1; ++jc) {
  1608. irow = ioffst + (1 - iskew) * jc;
  1609. i__2 = irow + jc * a_dim1;
  1610. i__4 = irow + jc * a_dim1;
  1611. r__1 = a[i__4].r;
  1612. q__1.r = r__1, q__1.i = 0.f;
  1613. a[i__2].r = q__1.r, a[i__2].i = q__1.i;
  1614. /* L440: */
  1615. }
  1616. }
  1617. }
  1618. } else {
  1619. /* 4) Generate Banded Matrix by first */
  1620. /* Rotating by random Unitary matrices, */
  1621. /* then reducing the bandwidth using Householder */
  1622. /* transformations. */
  1623. /* Note: we should get here only if LDA .ge. N */
  1624. if (isym == 1) {
  1625. /* Non-symmetric -- A = U D V */
  1626. clagge_(&mr, &nc, &llb, &uub, &d__[1], &a[a_offset], lda, &iseed[
  1627. 1], &work[1], &iinfo);
  1628. } else {
  1629. /* Symmetric -- A = U D U' or */
  1630. /* Hermitian -- A = U D U* */
  1631. if (csym) {
  1632. clagsy_(m, &llb, &d__[1], &a[a_offset], lda, &iseed[1], &work[
  1633. 1], &iinfo);
  1634. } else {
  1635. claghe_(m, &llb, &d__[1], &a[a_offset], lda, &iseed[1], &work[
  1636. 1], &iinfo);
  1637. }
  1638. }
  1639. if (iinfo != 0) {
  1640. *info = 3;
  1641. return;
  1642. }
  1643. }
  1644. /* 5) Pack the matrix */
  1645. if (ipack != ipackg) {
  1646. if (ipack == 1) {
  1647. /* 'U' -- Upper triangular, not packed */
  1648. i__1 = *m;
  1649. for (j = 1; j <= i__1; ++j) {
  1650. i__2 = *m;
  1651. for (i__ = j + 1; i__ <= i__2; ++i__) {
  1652. i__4 = i__ + j * a_dim1;
  1653. a[i__4].r = 0.f, a[i__4].i = 0.f;
  1654. /* L450: */
  1655. }
  1656. /* L460: */
  1657. }
  1658. } else if (ipack == 2) {
  1659. /* 'L' -- Lower triangular, not packed */
  1660. i__1 = *m;
  1661. for (j = 2; j <= i__1; ++j) {
  1662. i__2 = j - 1;
  1663. for (i__ = 1; i__ <= i__2; ++i__) {
  1664. i__4 = i__ + j * a_dim1;
  1665. a[i__4].r = 0.f, a[i__4].i = 0.f;
  1666. /* L470: */
  1667. }
  1668. /* L480: */
  1669. }
  1670. } else if (ipack == 3) {
  1671. /* 'C' -- Upper triangle packed Columnwise. */
  1672. icol = 1;
  1673. irow = 0;
  1674. i__1 = *m;
  1675. for (j = 1; j <= i__1; ++j) {
  1676. i__2 = j;
  1677. for (i__ = 1; i__ <= i__2; ++i__) {
  1678. ++irow;
  1679. if (irow > *lda) {
  1680. irow = 1;
  1681. ++icol;
  1682. }
  1683. i__4 = irow + icol * a_dim1;
  1684. i__3 = i__ + j * a_dim1;
  1685. a[i__4].r = a[i__3].r, a[i__4].i = a[i__3].i;
  1686. /* L490: */
  1687. }
  1688. /* L500: */
  1689. }
  1690. } else if (ipack == 4) {
  1691. /* 'R' -- Lower triangle packed Columnwise. */
  1692. icol = 1;
  1693. irow = 0;
  1694. i__1 = *m;
  1695. for (j = 1; j <= i__1; ++j) {
  1696. i__2 = *m;
  1697. for (i__ = j; i__ <= i__2; ++i__) {
  1698. ++irow;
  1699. if (irow > *lda) {
  1700. irow = 1;
  1701. ++icol;
  1702. }
  1703. i__4 = irow + icol * a_dim1;
  1704. i__3 = i__ + j * a_dim1;
  1705. a[i__4].r = a[i__3].r, a[i__4].i = a[i__3].i;
  1706. /* L510: */
  1707. }
  1708. /* L520: */
  1709. }
  1710. } else if (ipack >= 5) {
  1711. /* 'B' -- The lower triangle is packed as a band matrix. */
  1712. /* 'Q' -- The upper triangle is packed as a band matrix. */
  1713. /* 'Z' -- The whole matrix is packed as a band matrix. */
  1714. if (ipack == 5) {
  1715. uub = 0;
  1716. }
  1717. if (ipack == 6) {
  1718. llb = 0;
  1719. }
  1720. i__1 = uub;
  1721. for (j = 1; j <= i__1; ++j) {
  1722. /* Computing MIN */
  1723. i__2 = j + llb;
  1724. for (i__ = f2cmin(i__2,*m); i__ >= 1; --i__) {
  1725. i__2 = i__ - j + uub + 1 + j * a_dim1;
  1726. i__4 = i__ + j * a_dim1;
  1727. a[i__2].r = a[i__4].r, a[i__2].i = a[i__4].i;
  1728. /* L530: */
  1729. }
  1730. /* L540: */
  1731. }
  1732. i__1 = *n;
  1733. for (j = uub + 2; j <= i__1; ++j) {
  1734. /* Computing MIN */
  1735. i__4 = j + llb;
  1736. i__2 = f2cmin(i__4,*m);
  1737. for (i__ = j - uub; i__ <= i__2; ++i__) {
  1738. i__4 = i__ - j + uub + 1 + j * a_dim1;
  1739. i__3 = i__ + j * a_dim1;
  1740. a[i__4].r = a[i__3].r, a[i__4].i = a[i__3].i;
  1741. /* L550: */
  1742. }
  1743. /* L560: */
  1744. }
  1745. }
  1746. /* If packed, zero out extraneous elements. */
  1747. /* Symmetric/Triangular Packed -- */
  1748. /* zero out everything after A(IROW,ICOL) */
  1749. if (ipack == 3 || ipack == 4) {
  1750. i__1 = *m;
  1751. for (jc = icol; jc <= i__1; ++jc) {
  1752. i__2 = *lda;
  1753. for (jr = irow + 1; jr <= i__2; ++jr) {
  1754. i__4 = jr + jc * a_dim1;
  1755. a[i__4].r = 0.f, a[i__4].i = 0.f;
  1756. /* L570: */
  1757. }
  1758. irow = 0;
  1759. /* L580: */
  1760. }
  1761. } else if (ipack >= 5) {
  1762. /* Packed Band -- */
  1763. /* 1st row is now in A( UUB+2-j, j), zero above it */
  1764. /* m-th row is now in A( M+UUB-j,j), zero below it */
  1765. /* last non-zero diagonal is now in A( UUB+LLB+1,j ), */
  1766. /* zero below it, too. */
  1767. ir1 = uub + llb + 2;
  1768. ir2 = uub + *m + 2;
  1769. i__1 = *n;
  1770. for (jc = 1; jc <= i__1; ++jc) {
  1771. i__2 = uub + 1 - jc;
  1772. for (jr = 1; jr <= i__2; ++jr) {
  1773. i__4 = jr + jc * a_dim1;
  1774. a[i__4].r = 0.f, a[i__4].i = 0.f;
  1775. /* L590: */
  1776. }
  1777. /* Computing MAX */
  1778. /* Computing MIN */
  1779. i__3 = ir1, i__5 = ir2 - jc;
  1780. i__2 = 1, i__4 = f2cmin(i__3,i__5);
  1781. i__6 = *lda;
  1782. for (jr = f2cmax(i__2,i__4); jr <= i__6; ++jr) {
  1783. i__2 = jr + jc * a_dim1;
  1784. a[i__2].r = 0.f, a[i__2].i = 0.f;
  1785. /* L600: */
  1786. }
  1787. /* L610: */
  1788. }
  1789. }
  1790. }
  1791. return;
  1792. /* End of CLATMT */
  1793. } /* clatmt_ */