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