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.

slatms.c 48 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641
  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 integer c__1 = 1;
  240. static real c_b22 = 0.f;
  241. static logical c_true = TRUE_;
  242. static logical c_false = FALSE_;
  243. /* > \brief \b SLATMS */
  244. /* =========== DOCUMENTATION =========== */
  245. /* Online html documentation available at */
  246. /* http://www.netlib.org/lapack/explore-html/ */
  247. /* Definition: */
  248. /* =========== */
  249. /* SUBROUTINE SLATMS( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, */
  250. /* KL, KU, PACK, A, LDA, WORK, INFO ) */
  251. /* CHARACTER DIST, PACK, SYM */
  252. /* INTEGER INFO, KL, KU, LDA, M, MODE, N */
  253. /* REAL COND, DMAX */
  254. /* INTEGER ISEED( 4 ) */
  255. /* REAL A( LDA, * ), D( * ), WORK( * ) */
  256. /* > \par Purpose: */
  257. /* ============= */
  258. /* > */
  259. /* > \verbatim */
  260. /* > */
  261. /* > SLATMS generates random matrices with specified singular values */
  262. /* > (or symmetric/hermitian with specified eigenvalues) */
  263. /* > for testing LAPACK programs. */
  264. /* > */
  265. /* > SLATMS operates by applying the following sequence of */
  266. /* > operations: */
  267. /* > */
  268. /* > Set the diagonal to D, where D may be input or */
  269. /* > computed according to MODE, COND, DMAX, and SYM */
  270. /* > as described below. */
  271. /* > */
  272. /* > Generate a matrix with the appropriate band structure, by one */
  273. /* > of two methods: */
  274. /* > */
  275. /* > Method A: */
  276. /* > Generate a dense M x N matrix by multiplying D on the left */
  277. /* > and the right by random unitary matrices, then: */
  278. /* > */
  279. /* > Reduce the bandwidth according to KL and KU, using */
  280. /* > Householder transformations. */
  281. /* > */
  282. /* > Method B: */
  283. /* > Convert the bandwidth-0 (i.e., diagonal) matrix to a */
  284. /* > bandwidth-1 matrix using Givens rotations, "chasing" */
  285. /* > out-of-band elements back, much as in QR; then */
  286. /* > convert the bandwidth-1 to a bandwidth-2 matrix, etc. */
  287. /* > Note that for reasonably small bandwidths (relative to */
  288. /* > M and N) this requires less storage, as a dense matrix */
  289. /* > is not generated. Also, for symmetric matrices, only */
  290. /* > one triangle is generated. */
  291. /* > */
  292. /* > Method A is chosen if the bandwidth is a large fraction of the */
  293. /* > order of the matrix, and LDA is at least M (so a dense */
  294. /* > matrix can be stored.) Method B is chosen if the bandwidth */
  295. /* > is small (< 1/2 N for symmetric, < .3 N+M for */
  296. /* > non-symmetric), or LDA is less than M and not less than the */
  297. /* > bandwidth. */
  298. /* > */
  299. /* > Pack the matrix if desired. Options specified by PACK are: */
  300. /* > no packing */
  301. /* > zero out upper half (if symmetric) */
  302. /* > zero out lower half (if symmetric) */
  303. /* > store the upper half columnwise (if symmetric or upper */
  304. /* > triangular) */
  305. /* > store the lower half columnwise (if symmetric or lower */
  306. /* > triangular) */
  307. /* > store the lower triangle in banded format (if symmetric */
  308. /* > or lower triangular) */
  309. /* > store the upper triangle in banded format (if symmetric */
  310. /* > or upper triangular) */
  311. /* > store the entire matrix in banded format */
  312. /* > If Method B is chosen, and band format is specified, then the */
  313. /* > matrix will be generated in the band format, so no repacking */
  314. /* > will be necessary. */
  315. /* > \endverbatim */
  316. /* Arguments: */
  317. /* ========== */
  318. /* > \param[in] M */
  319. /* > \verbatim */
  320. /* > M is INTEGER */
  321. /* > The number of rows of A. Not modified. */
  322. /* > \endverbatim */
  323. /* > */
  324. /* > \param[in] N */
  325. /* > \verbatim */
  326. /* > N is INTEGER */
  327. /* > The number of columns of A. Not modified. */
  328. /* > \endverbatim */
  329. /* > */
  330. /* > \param[in] DIST */
  331. /* > \verbatim */
  332. /* > DIST is CHARACTER*1 */
  333. /* > On entry, DIST specifies the type of distribution to be used */
  334. /* > to generate the random eigen-/singular values. */
  335. /* > 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform ) */
  336. /* > 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric ) */
  337. /* > 'N' => NORMAL( 0, 1 ) ( 'N' for normal ) */
  338. /* > Not modified. */
  339. /* > \endverbatim */
  340. /* > */
  341. /* > \param[in,out] ISEED */
  342. /* > \verbatim */
  343. /* > ISEED is INTEGER array, dimension ( 4 ) */
  344. /* > On entry ISEED specifies the seed of the random number */
  345. /* > generator. They should lie between 0 and 4095 inclusive, */
  346. /* > and ISEED(4) should be odd. The random number generator */
  347. /* > uses a linear congruential sequence limited to small */
  348. /* > integers, and so should produce machine independent */
  349. /* > random numbers. The values of ISEED are changed on */
  350. /* > exit, and can be used in the next call to SLATMS */
  351. /* > to continue the same random number sequence. */
  352. /* > Changed on exit. */
  353. /* > \endverbatim */
  354. /* > */
  355. /* > \param[in] SYM */
  356. /* > \verbatim */
  357. /* > SYM is CHARACTER*1 */
  358. /* > If SYM='S' or 'H', the generated matrix is symmetric, with */
  359. /* > eigenvalues specified by D, COND, MODE, and DMAX; they */
  360. /* > may be positive, negative, or zero. */
  361. /* > If SYM='P', the generated matrix is symmetric, with */
  362. /* > eigenvalues (= singular values) specified by D, COND, */
  363. /* > MODE, and DMAX; they will not be negative. */
  364. /* > If SYM='N', the generated matrix is nonsymmetric, with */
  365. /* > singular values specified by D, COND, MODE, and DMAX; */
  366. /* > they will not be negative. */
  367. /* > Not modified. */
  368. /* > \endverbatim */
  369. /* > */
  370. /* > \param[in,out] D */
  371. /* > \verbatim */
  372. /* > D is REAL array, dimension ( MIN( M , N ) ) */
  373. /* > This array is used to specify the singular values or */
  374. /* > eigenvalues of A (see SYM, above.) If MODE=0, then D is */
  375. /* > assumed to contain the singular/eigenvalues, otherwise */
  376. /* > they will be computed according to MODE, COND, and DMAX, */
  377. /* > and placed in D. */
  378. /* > Modified if MODE is nonzero. */
  379. /* > \endverbatim */
  380. /* > */
  381. /* > \param[in] MODE */
  382. /* > \verbatim */
  383. /* > MODE is INTEGER */
  384. /* > On entry this describes how the singular/eigenvalues are to */
  385. /* > be specified: */
  386. /* > MODE = 0 means use D as input */
  387. /* > MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND */
  388. /* > MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND */
  389. /* > MODE = 3 sets D(I)=COND**(-(I-1)/(N-1)) */
  390. /* > MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND) */
  391. /* > MODE = 5 sets D to random numbers in the range */
  392. /* > ( 1/COND , 1 ) such that their logarithms */
  393. /* > are uniformly distributed. */
  394. /* > MODE = 6 set D to random numbers from same distribution */
  395. /* > as the rest of the matrix. */
  396. /* > MODE < 0 has the same meaning as ABS(MODE), except that */
  397. /* > the order of the elements of D is reversed. */
  398. /* > Thus if MODE is positive, D has entries ranging from */
  399. /* > 1 to 1/COND, if negative, from 1/COND to 1, */
  400. /* > If SYM='S' or 'H', and MODE is neither 0, 6, nor -6, then */
  401. /* > the elements of D will also be multiplied by a random */
  402. /* > sign (i.e., +1 or -1.) */
  403. /* > Not modified. */
  404. /* > \endverbatim */
  405. /* > */
  406. /* > \param[in] COND */
  407. /* > \verbatim */
  408. /* > COND is REAL */
  409. /* > On entry, this is used as described under MODE above. */
  410. /* > If used, it must be >= 1. Not modified. */
  411. /* > \endverbatim */
  412. /* > */
  413. /* > \param[in] DMAX */
  414. /* > \verbatim */
  415. /* > DMAX is REAL */
  416. /* > If MODE is neither -6, 0 nor 6, the contents of D, as */
  417. /* > computed according to MODE and COND, will be scaled by */
  418. /* > DMAX / f2cmax(abs(D(i))); thus, the maximum absolute eigen- or */
  419. /* > singular value (which is to say the norm) will be abs(DMAX). */
  420. /* > Note that DMAX need not be positive: if DMAX is negative */
  421. /* > (or zero), D will be scaled by a negative number (or zero). */
  422. /* > Not modified. */
  423. /* > \endverbatim */
  424. /* > */
  425. /* > \param[in] KL */
  426. /* > \verbatim */
  427. /* > KL is INTEGER */
  428. /* > This specifies the lower bandwidth of the matrix. For */
  429. /* > example, KL=0 implies upper triangular, KL=1 implies upper */
  430. /* > Hessenberg, and KL being at least M-1 means that the matrix */
  431. /* > has full lower bandwidth. KL must equal KU if the matrix */
  432. /* > is symmetric. */
  433. /* > Not modified. */
  434. /* > \endverbatim */
  435. /* > */
  436. /* > \param[in] KU */
  437. /* > \verbatim */
  438. /* > KU is INTEGER */
  439. /* > This specifies the upper bandwidth of the matrix. For */
  440. /* > example, KU=0 implies lower triangular, KU=1 implies lower */
  441. /* > Hessenberg, and KU being at least N-1 means that the matrix */
  442. /* > has full upper bandwidth. KL must equal KU if the matrix */
  443. /* > is symmetric. */
  444. /* > Not modified. */
  445. /* > \endverbatim */
  446. /* > */
  447. /* > \param[in] PACK */
  448. /* > \verbatim */
  449. /* > PACK is CHARACTER*1 */
  450. /* > This specifies packing of matrix as follows: */
  451. /* > 'N' => no packing */
  452. /* > 'U' => zero out all subdiagonal entries (if symmetric) */
  453. /* > 'L' => zero out all superdiagonal entries (if symmetric) */
  454. /* > 'C' => store the upper triangle columnwise */
  455. /* > (only if the matrix is symmetric or upper triangular) */
  456. /* > 'R' => store the lower triangle columnwise */
  457. /* > (only if the matrix is symmetric or lower triangular) */
  458. /* > 'B' => store the lower triangle in band storage scheme */
  459. /* > (only if matrix symmetric or lower triangular) */
  460. /* > 'Q' => store the upper triangle in band storage scheme */
  461. /* > (only if matrix symmetric or upper triangular) */
  462. /* > 'Z' => store the entire matrix in band storage scheme */
  463. /* > (pivoting can be provided for by using this */
  464. /* > option to store A in the trailing rows of */
  465. /* > the allocated storage) */
  466. /* > */
  467. /* > Using these options, the various LAPACK packed and banded */
  468. /* > storage schemes can be obtained: */
  469. /* > GB - use 'Z' */
  470. /* > PB, SB or TB - use 'B' or 'Q' */
  471. /* > PP, SP or TP - use 'C' or 'R' */
  472. /* > */
  473. /* > If two calls to SLATMS differ only in the PACK parameter, */
  474. /* > they will generate mathematically equivalent matrices. */
  475. /* > Not modified. */
  476. /* > \endverbatim */
  477. /* > */
  478. /* > \param[in,out] A */
  479. /* > \verbatim */
  480. /* > A is REAL array, dimension ( LDA, N ) */
  481. /* > On exit A is the desired test matrix. A is first generated */
  482. /* > in full (unpacked) form, and then packed, if so specified */
  483. /* > by PACK. Thus, the first M elements of the first N */
  484. /* > columns will always be modified. If PACK specifies a */
  485. /* > packed or banded storage scheme, all LDA elements of the */
  486. /* > first N columns will be modified; the elements of the */
  487. /* > array which do not correspond to elements of the generated */
  488. /* > matrix are set to zero. */
  489. /* > Modified. */
  490. /* > \endverbatim */
  491. /* > */
  492. /* > \param[in] LDA */
  493. /* > \verbatim */
  494. /* > LDA is INTEGER */
  495. /* > LDA specifies the first dimension of A as declared in the */
  496. /* > calling program. If PACK='N', 'U', 'L', 'C', or 'R', then */
  497. /* > LDA must be at least M. If PACK='B' or 'Q', then LDA must */
  498. /* > be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)). */
  499. /* > If PACK='Z', LDA must be large enough to hold the packed */
  500. /* > array: MIN( KU, N-1) + MIN( KL, M-1) + 1. */
  501. /* > Not modified. */
  502. /* > \endverbatim */
  503. /* > */
  504. /* > \param[out] WORK */
  505. /* > \verbatim */
  506. /* > WORK is REAL array, dimension ( 3*MAX( N , M ) ) */
  507. /* > Workspace. */
  508. /* > Modified. */
  509. /* > \endverbatim */
  510. /* > */
  511. /* > \param[out] INFO */
  512. /* > \verbatim */
  513. /* > INFO is INTEGER */
  514. /* > Error code. On exit, INFO will be set to one of the */
  515. /* > following values: */
  516. /* > 0 => normal return */
  517. /* > -1 => M negative or unequal to N and SYM='S', 'H', or 'P' */
  518. /* > -2 => N negative */
  519. /* > -3 => DIST illegal string */
  520. /* > -5 => SYM illegal string */
  521. /* > -7 => MODE not in range -6 to 6 */
  522. /* > -8 => COND less than 1.0, and MODE neither -6, 0 nor 6 */
  523. /* > -10 => KL negative */
  524. /* > -11 => KU negative, or SYM='S' or 'H' and KU not equal to KL */
  525. /* > -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N'; */
  526. /* > or PACK='C' or 'Q' and SYM='N' and KL is not zero; */
  527. /* > or PACK='R' or 'B' and SYM='N' and KU is not zero; */
  528. /* > or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not */
  529. /* > N. */
  530. /* > -14 => LDA is less than M, or PACK='Z' and LDA is less than */
  531. /* > MIN(KU,N-1) + MIN(KL,M-1) + 1. */
  532. /* > 1 => Error return from SLATM1 */
  533. /* > 2 => Cannot scale to DMAX (f2cmax. sing. value is 0) */
  534. /* > 3 => Error return from SLAGGE or SLAGSY */
  535. /* > \endverbatim */
  536. /* Authors: */
  537. /* ======== */
  538. /* > \author Univ. of Tennessee */
  539. /* > \author Univ. of California Berkeley */
  540. /* > \author Univ. of Colorado Denver */
  541. /* > \author NAG Ltd. */
  542. /* > \date December 2016 */
  543. /* > \ingroup real_matgen */
  544. /* ===================================================================== */
  545. /* Subroutine */ void slatms_(integer *m, integer *n, char *dist, integer *
  546. iseed, char *sym, real *d__, integer *mode, real *cond, real *dmax__,
  547. integer *kl, integer *ku, char *pack, real *a, integer *lda, real *
  548. work, integer *info)
  549. {
  550. /* System generated locals */
  551. integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5, i__6;
  552. real r__1, r__2, r__3;
  553. logical L__1;
  554. /* Local variables */
  555. integer ilda, icol;
  556. real temp;
  557. integer irow, isym;
  558. real c__;
  559. integer i__, j, k;
  560. real s, alpha, angle;
  561. integer ipack, ioffg;
  562. extern logical lsame_(char *, char *);
  563. integer iinfo;
  564. extern /* Subroutine */ void sscal_(integer *, real *, real *, integer *);
  565. integer idist, mnmin, iskew;
  566. real extra, dummy;
  567. extern /* Subroutine */ void scopy_(integer *, real *, integer *, real *,
  568. integer *), slatm1_(integer *, real *, integer *, integer *,
  569. integer *, real *, integer *, integer *);
  570. integer ic, jc, nc, il, iendch, ir, jr, ipackg, mr;
  571. extern /* Subroutine */ void slagge_(integer *, integer *, integer *,
  572. integer *, real *, real *, integer *, integer *, real *, integer *
  573. );
  574. integer minlda;
  575. extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
  576. extern real slarnd_(integer *, integer *);
  577. logical iltemp, givens;
  578. integer ioffst, irsign;
  579. extern /* Subroutine */ void slartg_(real *, real *, real *, real *, real *
  580. ), slaset_(char *, integer *, integer *, real *, real *, real *,
  581. integer *), slagsy_(integer *, integer *, real *, real *,
  582. integer *, integer *, real *, integer *), slarot_(logical *,
  583. logical *, logical *, integer *, real *, real *, real *, integer *
  584. , real *, real *);
  585. logical ilextr, topdwn;
  586. integer ir1, ir2, isympk, jch, llb, jkl, jku, uub;
  587. /* -- LAPACK computational routine (version 3.7.0) -- */
  588. /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
  589. /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
  590. /* December 2016 */
  591. /* ===================================================================== */
  592. /* 1) Decode and Test the input parameters. */
  593. /* Initialize flags & seed. */
  594. /* Parameter adjustments */
  595. --iseed;
  596. --d__;
  597. a_dim1 = *lda;
  598. a_offset = 1 + a_dim1 * 1;
  599. a -= a_offset;
  600. --work;
  601. /* Function Body */
  602. *info = 0;
  603. /* Quick return if possible */
  604. if (*m == 0 || *n == 0) {
  605. return;
  606. }
  607. /* Decode DIST */
  608. if (lsame_(dist, "U")) {
  609. idist = 1;
  610. } else if (lsame_(dist, "S")) {
  611. idist = 2;
  612. } else if (lsame_(dist, "N")) {
  613. idist = 3;
  614. } else {
  615. idist = -1;
  616. }
  617. /* Decode SYM */
  618. if (lsame_(sym, "N")) {
  619. isym = 1;
  620. irsign = 0;
  621. } else if (lsame_(sym, "P")) {
  622. isym = 2;
  623. irsign = 0;
  624. } else if (lsame_(sym, "S")) {
  625. isym = 2;
  626. irsign = 1;
  627. } else if (lsame_(sym, "H")) {
  628. isym = 2;
  629. irsign = 1;
  630. } else {
  631. isym = -1;
  632. }
  633. /* Decode PACK */
  634. isympk = 0;
  635. if (lsame_(pack, "N")) {
  636. ipack = 0;
  637. } else if (lsame_(pack, "U")) {
  638. ipack = 1;
  639. isympk = 1;
  640. } else if (lsame_(pack, "L")) {
  641. ipack = 2;
  642. isympk = 1;
  643. } else if (lsame_(pack, "C")) {
  644. ipack = 3;
  645. isympk = 2;
  646. } else if (lsame_(pack, "R")) {
  647. ipack = 4;
  648. isympk = 3;
  649. } else if (lsame_(pack, "B")) {
  650. ipack = 5;
  651. isympk = 3;
  652. } else if (lsame_(pack, "Q")) {
  653. ipack = 6;
  654. isympk = 2;
  655. } else if (lsame_(pack, "Z")) {
  656. ipack = 7;
  657. } else {
  658. ipack = -1;
  659. }
  660. /* Set certain internal parameters */
  661. mnmin = f2cmin(*m,*n);
  662. /* Computing MIN */
  663. i__1 = *kl, i__2 = *m - 1;
  664. llb = f2cmin(i__1,i__2);
  665. /* Computing MIN */
  666. i__1 = *ku, i__2 = *n - 1;
  667. uub = f2cmin(i__1,i__2);
  668. /* Computing MIN */
  669. i__1 = *m, i__2 = *n + llb;
  670. mr = f2cmin(i__1,i__2);
  671. /* Computing MIN */
  672. i__1 = *n, i__2 = *m + uub;
  673. nc = f2cmin(i__1,i__2);
  674. if (ipack == 5 || ipack == 6) {
  675. minlda = uub + 1;
  676. } else if (ipack == 7) {
  677. minlda = llb + uub + 1;
  678. } else {
  679. minlda = *m;
  680. }
  681. /* Use Givens rotation method if bandwidth small enough, */
  682. /* or if LDA is too small to store the matrix unpacked. */
  683. givens = FALSE_;
  684. if (isym == 1) {
  685. /* Computing MAX */
  686. i__1 = 1, i__2 = mr + nc;
  687. if ((real) (llb + uub) < (real) f2cmax(i__1,i__2) * .3f) {
  688. givens = TRUE_;
  689. }
  690. } else {
  691. if (llb << 1 < *m) {
  692. givens = TRUE_;
  693. }
  694. }
  695. if (*lda < *m && *lda >= minlda) {
  696. givens = TRUE_;
  697. }
  698. /* Set INFO if an error */
  699. if (*m < 0) {
  700. *info = -1;
  701. } else if (*m != *n && isym != 1) {
  702. *info = -1;
  703. } else if (*n < 0) {
  704. *info = -2;
  705. } else if (idist == -1) {
  706. *info = -3;
  707. } else if (isym == -1) {
  708. *info = -5;
  709. } else if (abs(*mode) > 6) {
  710. *info = -7;
  711. } else if (*mode != 0 && abs(*mode) != 6 && *cond < 1.f) {
  712. *info = -8;
  713. } else if (*kl < 0) {
  714. *info = -10;
  715. } else if (*ku < 0 || isym != 1 && *kl != *ku) {
  716. *info = -11;
  717. } else if (ipack == -1 || isympk == 1 && isym == 1 || isympk == 2 && isym
  718. == 1 && *kl > 0 || isympk == 3 && isym == 1 && *ku > 0 || isympk
  719. != 0 && *m != *n) {
  720. *info = -12;
  721. } else if (*lda < f2cmax(1,minlda)) {
  722. *info = -14;
  723. }
  724. if (*info != 0) {
  725. i__1 = -(*info);
  726. xerbla_("SLATMS", &i__1, 6);
  727. return;
  728. }
  729. /* Initialize random number generator */
  730. for (i__ = 1; i__ <= 4; ++i__) {
  731. iseed[i__] = (i__1 = iseed[i__], abs(i__1)) % 4096;
  732. /* L10: */
  733. }
  734. if (iseed[4] % 2 != 1) {
  735. ++iseed[4];
  736. }
  737. /* 2) Set up D if indicated. */
  738. /* Compute D according to COND and MODE */
  739. slatm1_(mode, cond, &irsign, &idist, &iseed[1], &d__[1], &mnmin, &iinfo);
  740. if (iinfo != 0) {
  741. *info = 1;
  742. return;
  743. }
  744. /* Choose Top-Down if D is (apparently) increasing, */
  745. /* Bottom-Up if D is (apparently) decreasing. */
  746. if (abs(d__[1]) <= (r__1 = d__[mnmin], abs(r__1))) {
  747. topdwn = TRUE_;
  748. } else {
  749. topdwn = FALSE_;
  750. }
  751. if (*mode != 0 && abs(*mode) != 6) {
  752. /* Scale by DMAX */
  753. temp = abs(d__[1]);
  754. i__1 = mnmin;
  755. for (i__ = 2; i__ <= i__1; ++i__) {
  756. /* Computing MAX */
  757. r__2 = temp, r__3 = (r__1 = d__[i__], abs(r__1));
  758. temp = f2cmax(r__2,r__3);
  759. /* L20: */
  760. }
  761. if (temp > 0.f) {
  762. alpha = *dmax__ / temp;
  763. } else {
  764. *info = 2;
  765. return;
  766. }
  767. sscal_(&mnmin, &alpha, &d__[1], &c__1);
  768. }
  769. /* 3) Generate Banded Matrix using Givens rotations. */
  770. /* Also the special case of UUB=LLB=0 */
  771. /* Compute Addressing constants to cover all */
  772. /* storage formats. Whether GE, SY, GB, or SB, */
  773. /* upper or lower triangle or both, */
  774. /* the (i,j)-th element is in */
  775. /* A( i - ISKEW*j + IOFFST, j ) */
  776. if (ipack > 4) {
  777. ilda = *lda - 1;
  778. iskew = 1;
  779. if (ipack > 5) {
  780. ioffst = uub + 1;
  781. } else {
  782. ioffst = 1;
  783. }
  784. } else {
  785. ilda = *lda;
  786. iskew = 0;
  787. ioffst = 0;
  788. }
  789. /* IPACKG is the format that the matrix is generated in. If this is */
  790. /* different from IPACK, then the matrix must be repacked at the */
  791. /* end. It also signals how to compute the norm, for scaling. */
  792. ipackg = 0;
  793. slaset_("Full", lda, n, &c_b22, &c_b22, &a[a_offset], lda);
  794. /* Diagonal Matrix -- We are done, unless it */
  795. /* is to be stored SP/PP/TP (PACK='R' or 'C') */
  796. if (llb == 0 && uub == 0) {
  797. i__1 = ilda + 1;
  798. scopy_(&mnmin, &d__[1], &c__1, &a[1 - iskew + ioffst + a_dim1], &i__1)
  799. ;
  800. if (ipack <= 2 || ipack >= 5) {
  801. ipackg = ipack;
  802. }
  803. } else if (givens) {
  804. /* Check whether to use Givens rotations, */
  805. /* Householder transformations, or nothing. */
  806. if (isym == 1) {
  807. /* Non-symmetric -- A = U D V */
  808. if (ipack > 4) {
  809. ipackg = ipack;
  810. } else {
  811. ipackg = 0;
  812. }
  813. i__1 = ilda + 1;
  814. scopy_(&mnmin, &d__[1], &c__1, &a[1 - iskew + ioffst + a_dim1], &
  815. i__1);
  816. if (topdwn) {
  817. jkl = 0;
  818. i__1 = uub;
  819. for (jku = 1; jku <= i__1; ++jku) {
  820. /* Transform from bandwidth JKL, JKU-1 to JKL, JKU */
  821. /* Last row actually rotated is M */
  822. /* Last column actually rotated is MIN( M+JKU, N ) */
  823. /* Computing MIN */
  824. i__3 = *m + jku;
  825. i__2 = f2cmin(i__3,*n) + jkl - 1;
  826. for (jr = 1; jr <= i__2; ++jr) {
  827. extra = 0.f;
  828. angle = slarnd_(&c__1, &iseed[1]) *
  829. 6.2831853071795864769252867663f;
  830. c__ = cos(angle);
  831. s = sin(angle);
  832. /* Computing MAX */
  833. i__3 = 1, i__4 = jr - jkl;
  834. icol = f2cmax(i__3,i__4);
  835. if (jr < *m) {
  836. /* Computing MIN */
  837. i__3 = *n, i__4 = jr + jku;
  838. il = f2cmin(i__3,i__4) + 1 - icol;
  839. L__1 = jr > jkl;
  840. slarot_(&c_true, &L__1, &c_false, &il, &c__, &s, &
  841. a[jr - iskew * icol + ioffst + icol *
  842. a_dim1], &ilda, &extra, &dummy);
  843. }
  844. /* Chase "EXTRA" back up */
  845. ir = jr;
  846. ic = icol;
  847. i__3 = -jkl - jku;
  848. for (jch = jr - jkl; i__3 < 0 ? jch >= 1 : jch <= 1;
  849. jch += i__3) {
  850. if (ir < *m) {
  851. slartg_(&a[ir + 1 - iskew * (ic + 1) + ioffst
  852. + (ic + 1) * a_dim1], &extra, &c__, &
  853. s, &dummy);
  854. }
  855. /* Computing MAX */
  856. i__4 = 1, i__5 = jch - jku;
  857. irow = f2cmax(i__4,i__5);
  858. il = ir + 2 - irow;
  859. temp = 0.f;
  860. iltemp = jch > jku;
  861. r__1 = -s;
  862. slarot_(&c_false, &iltemp, &c_true, &il, &c__, &
  863. r__1, &a[irow - iskew * ic + ioffst + ic *
  864. a_dim1], &ilda, &temp, &extra);
  865. if (iltemp) {
  866. slartg_(&a[irow + 1 - iskew * (ic + 1) +
  867. ioffst + (ic + 1) * a_dim1], &temp, &
  868. c__, &s, &dummy);
  869. /* Computing MAX */
  870. i__4 = 1, i__5 = jch - jku - jkl;
  871. icol = f2cmax(i__4,i__5);
  872. il = ic + 2 - icol;
  873. extra = 0.f;
  874. L__1 = jch > jku + jkl;
  875. r__1 = -s;
  876. slarot_(&c_true, &L__1, &c_true, &il, &c__, &
  877. r__1, &a[irow - iskew * icol + ioffst
  878. + icol * a_dim1], &ilda, &extra, &
  879. temp);
  880. ic = icol;
  881. ir = irow;
  882. }
  883. /* L30: */
  884. }
  885. /* L40: */
  886. }
  887. /* L50: */
  888. }
  889. jku = uub;
  890. i__1 = llb;
  891. for (jkl = 1; jkl <= i__1; ++jkl) {
  892. /* Transform from bandwidth JKL-1, JKU to JKL, JKU */
  893. /* Computing MIN */
  894. i__3 = *n + jkl;
  895. i__2 = f2cmin(i__3,*m) + jku - 1;
  896. for (jc = 1; jc <= i__2; ++jc) {
  897. extra = 0.f;
  898. angle = slarnd_(&c__1, &iseed[1]) *
  899. 6.2831853071795864769252867663f;
  900. c__ = cos(angle);
  901. s = sin(angle);
  902. /* Computing MAX */
  903. i__3 = 1, i__4 = jc - jku;
  904. irow = f2cmax(i__3,i__4);
  905. if (jc < *n) {
  906. /* Computing MIN */
  907. i__3 = *m, i__4 = jc + jkl;
  908. il = f2cmin(i__3,i__4) + 1 - irow;
  909. L__1 = jc > jku;
  910. slarot_(&c_false, &L__1, &c_false, &il, &c__, &s,
  911. &a[irow - iskew * jc + ioffst + jc *
  912. a_dim1], &ilda, &extra, &dummy);
  913. }
  914. /* Chase "EXTRA" back up */
  915. ic = jc;
  916. ir = irow;
  917. i__3 = -jkl - jku;
  918. for (jch = jc - jku; i__3 < 0 ? jch >= 1 : jch <= 1;
  919. jch += i__3) {
  920. if (ic < *n) {
  921. slartg_(&a[ir + 1 - iskew * (ic + 1) + ioffst
  922. + (ic + 1) * a_dim1], &extra, &c__, &
  923. s, &dummy);
  924. }
  925. /* Computing MAX */
  926. i__4 = 1, i__5 = jch - jkl;
  927. icol = f2cmax(i__4,i__5);
  928. il = ic + 2 - icol;
  929. temp = 0.f;
  930. iltemp = jch > jkl;
  931. r__1 = -s;
  932. slarot_(&c_true, &iltemp, &c_true, &il, &c__, &
  933. r__1, &a[ir - iskew * icol + ioffst +
  934. icol * a_dim1], &ilda, &temp, &extra);
  935. if (iltemp) {
  936. slartg_(&a[ir + 1 - iskew * (icol + 1) +
  937. ioffst + (icol + 1) * a_dim1], &temp,
  938. &c__, &s, &dummy);
  939. /* Computing MAX */
  940. i__4 = 1, i__5 = jch - jkl - jku;
  941. irow = f2cmax(i__4,i__5);
  942. il = ir + 2 - irow;
  943. extra = 0.f;
  944. L__1 = jch > jkl + jku;
  945. r__1 = -s;
  946. slarot_(&c_false, &L__1, &c_true, &il, &c__, &
  947. r__1, &a[irow - iskew * icol + ioffst
  948. + icol * a_dim1], &ilda, &extra, &
  949. temp);
  950. ic = icol;
  951. ir = irow;
  952. }
  953. /* L60: */
  954. }
  955. /* L70: */
  956. }
  957. /* L80: */
  958. }
  959. } else {
  960. /* Bottom-Up -- Start at the bottom right. */
  961. jkl = 0;
  962. i__1 = uub;
  963. for (jku = 1; jku <= i__1; ++jku) {
  964. /* Transform from bandwidth JKL, JKU-1 to JKL, JKU */
  965. /* First row actually rotated is M */
  966. /* First column actually rotated is MIN( M+JKU, N ) */
  967. /* Computing MIN */
  968. i__2 = *m, i__3 = *n + jkl;
  969. iendch = f2cmin(i__2,i__3) - 1;
  970. /* Computing MIN */
  971. i__2 = *m + jku;
  972. i__3 = 1 - jkl;
  973. for (jc = f2cmin(i__2,*n) - 1; jc >= i__3; --jc) {
  974. extra = 0.f;
  975. angle = slarnd_(&c__1, &iseed[1]) *
  976. 6.2831853071795864769252867663f;
  977. c__ = cos(angle);
  978. s = sin(angle);
  979. /* Computing MAX */
  980. i__2 = 1, i__4 = jc - jku + 1;
  981. irow = f2cmax(i__2,i__4);
  982. if (jc > 0) {
  983. /* Computing MIN */
  984. i__2 = *m, i__4 = jc + jkl + 1;
  985. il = f2cmin(i__2,i__4) + 1 - irow;
  986. L__1 = jc + jkl < *m;
  987. slarot_(&c_false, &c_false, &L__1, &il, &c__, &s,
  988. &a[irow - iskew * jc + ioffst + jc *
  989. a_dim1], &ilda, &dummy, &extra);
  990. }
  991. /* Chase "EXTRA" back down */
  992. ic = jc;
  993. i__2 = iendch;
  994. i__4 = jkl + jku;
  995. for (jch = jc + jkl; i__4 < 0 ? jch >= i__2 : jch <=
  996. i__2; jch += i__4) {
  997. ilextr = ic > 0;
  998. if (ilextr) {
  999. slartg_(&a[jch - iskew * ic + ioffst + ic *
  1000. a_dim1], &extra, &c__, &s, &dummy);
  1001. }
  1002. ic = f2cmax(1,ic);
  1003. /* Computing MIN */
  1004. i__5 = *n - 1, i__6 = jch + jku;
  1005. icol = f2cmin(i__5,i__6);
  1006. iltemp = jch + jku < *n;
  1007. temp = 0.f;
  1008. i__5 = icol + 2 - ic;
  1009. slarot_(&c_true, &ilextr, &iltemp, &i__5, &c__, &
  1010. s, &a[jch - iskew * ic + ioffst + ic *
  1011. a_dim1], &ilda, &extra, &temp);
  1012. if (iltemp) {
  1013. slartg_(&a[jch - iskew * icol + ioffst + icol
  1014. * a_dim1], &temp, &c__, &s, &dummy);
  1015. /* Computing MIN */
  1016. i__5 = iendch, i__6 = jch + jkl + jku;
  1017. il = f2cmin(i__5,i__6) + 2 - jch;
  1018. extra = 0.f;
  1019. L__1 = jch + jkl + jku <= iendch;
  1020. slarot_(&c_false, &c_true, &L__1, &il, &c__, &
  1021. s, &a[jch - iskew * icol + ioffst +
  1022. icol * a_dim1], &ilda, &temp, &extra);
  1023. ic = icol;
  1024. }
  1025. /* L90: */
  1026. }
  1027. /* L100: */
  1028. }
  1029. /* L110: */
  1030. }
  1031. jku = uub;
  1032. i__1 = llb;
  1033. for (jkl = 1; jkl <= i__1; ++jkl) {
  1034. /* Transform from bandwidth JKL-1, JKU to JKL, JKU */
  1035. /* First row actually rotated is MIN( N+JKL, M ) */
  1036. /* First column actually rotated is N */
  1037. /* Computing MIN */
  1038. i__3 = *n, i__4 = *m + jku;
  1039. iendch = f2cmin(i__3,i__4) - 1;
  1040. /* Computing MIN */
  1041. i__3 = *n + jkl;
  1042. i__4 = 1 - jku;
  1043. for (jr = f2cmin(i__3,*m) - 1; jr >= i__4; --jr) {
  1044. extra = 0.f;
  1045. angle = slarnd_(&c__1, &iseed[1]) *
  1046. 6.2831853071795864769252867663f;
  1047. c__ = cos(angle);
  1048. s = sin(angle);
  1049. /* Computing MAX */
  1050. i__3 = 1, i__2 = jr - jkl + 1;
  1051. icol = f2cmax(i__3,i__2);
  1052. if (jr > 0) {
  1053. /* Computing MIN */
  1054. i__3 = *n, i__2 = jr + jku + 1;
  1055. il = f2cmin(i__3,i__2) + 1 - icol;
  1056. L__1 = jr + jku < *n;
  1057. slarot_(&c_true, &c_false, &L__1, &il, &c__, &s, &
  1058. a[jr - iskew * icol + ioffst + icol *
  1059. a_dim1], &ilda, &dummy, &extra);
  1060. }
  1061. /* Chase "EXTRA" back down */
  1062. ir = jr;
  1063. i__3 = iendch;
  1064. i__2 = jkl + jku;
  1065. for (jch = jr + jku; i__2 < 0 ? jch >= i__3 : jch <=
  1066. i__3; jch += i__2) {
  1067. ilextr = ir > 0;
  1068. if (ilextr) {
  1069. slartg_(&a[ir - iskew * jch + ioffst + jch *
  1070. a_dim1], &extra, &c__, &s, &dummy);
  1071. }
  1072. ir = f2cmax(1,ir);
  1073. /* Computing MIN */
  1074. i__5 = *m - 1, i__6 = jch + jkl;
  1075. irow = f2cmin(i__5,i__6);
  1076. iltemp = jch + jkl < *m;
  1077. temp = 0.f;
  1078. i__5 = irow + 2 - ir;
  1079. slarot_(&c_false, &ilextr, &iltemp, &i__5, &c__, &
  1080. s, &a[ir - iskew * jch + ioffst + jch *
  1081. a_dim1], &ilda, &extra, &temp);
  1082. if (iltemp) {
  1083. slartg_(&a[irow - iskew * jch + ioffst + jch *
  1084. a_dim1], &temp, &c__, &s, &dummy);
  1085. /* Computing MIN */
  1086. i__5 = iendch, i__6 = jch + jkl + jku;
  1087. il = f2cmin(i__5,i__6) + 2 - jch;
  1088. extra = 0.f;
  1089. L__1 = jch + jkl + jku <= iendch;
  1090. slarot_(&c_true, &c_true, &L__1, &il, &c__, &
  1091. s, &a[irow - iskew * jch + ioffst +
  1092. jch * a_dim1], &ilda, &temp, &extra);
  1093. ir = irow;
  1094. }
  1095. /* L120: */
  1096. }
  1097. /* L130: */
  1098. }
  1099. /* L140: */
  1100. }
  1101. }
  1102. } else {
  1103. /* Symmetric -- A = U D U' */
  1104. ipackg = ipack;
  1105. ioffg = ioffst;
  1106. if (topdwn) {
  1107. /* Top-Down -- Generate Upper triangle only */
  1108. if (ipack >= 5) {
  1109. ipackg = 6;
  1110. ioffg = uub + 1;
  1111. } else {
  1112. ipackg = 1;
  1113. }
  1114. i__1 = ilda + 1;
  1115. scopy_(&mnmin, &d__[1], &c__1, &a[1 - iskew + ioffg + a_dim1],
  1116. &i__1);
  1117. i__1 = uub;
  1118. for (k = 1; k <= i__1; ++k) {
  1119. i__4 = *n - 1;
  1120. for (jc = 1; jc <= i__4; ++jc) {
  1121. /* Computing MAX */
  1122. i__2 = 1, i__3 = jc - k;
  1123. irow = f2cmax(i__2,i__3);
  1124. /* Computing MIN */
  1125. i__2 = jc + 1, i__3 = k + 2;
  1126. il = f2cmin(i__2,i__3);
  1127. extra = 0.f;
  1128. temp = a[jc - iskew * (jc + 1) + ioffg + (jc + 1) *
  1129. a_dim1];
  1130. angle = slarnd_(&c__1, &iseed[1]) *
  1131. 6.2831853071795864769252867663f;
  1132. c__ = cos(angle);
  1133. s = sin(angle);
  1134. L__1 = jc > k;
  1135. slarot_(&c_false, &L__1, &c_true, &il, &c__, &s, &a[
  1136. irow - iskew * jc + ioffg + jc * a_dim1], &
  1137. ilda, &extra, &temp);
  1138. /* Computing MIN */
  1139. i__3 = k, i__5 = *n - jc;
  1140. i__2 = f2cmin(i__3,i__5) + 1;
  1141. slarot_(&c_true, &c_true, &c_false, &i__2, &c__, &s, &
  1142. a[(1 - iskew) * jc + ioffg + jc * a_dim1], &
  1143. ilda, &temp, &dummy);
  1144. /* Chase EXTRA back up the matrix */
  1145. icol = jc;
  1146. i__2 = -k;
  1147. for (jch = jc - k; i__2 < 0 ? jch >= 1 : jch <= 1;
  1148. jch += i__2) {
  1149. slartg_(&a[jch + 1 - iskew * (icol + 1) + ioffg +
  1150. (icol + 1) * a_dim1], &extra, &c__, &s, &
  1151. dummy);
  1152. temp = a[jch - iskew * (jch + 1) + ioffg + (jch +
  1153. 1) * a_dim1];
  1154. i__3 = k + 2;
  1155. r__1 = -s;
  1156. slarot_(&c_true, &c_true, &c_true, &i__3, &c__, &
  1157. r__1, &a[(1 - iskew) * jch + ioffg + jch *
  1158. a_dim1], &ilda, &temp, &extra);
  1159. /* Computing MAX */
  1160. i__3 = 1, i__5 = jch - k;
  1161. irow = f2cmax(i__3,i__5);
  1162. /* Computing MIN */
  1163. i__3 = jch + 1, i__5 = k + 2;
  1164. il = f2cmin(i__3,i__5);
  1165. extra = 0.f;
  1166. L__1 = jch > k;
  1167. r__1 = -s;
  1168. slarot_(&c_false, &L__1, &c_true, &il, &c__, &
  1169. r__1, &a[irow - iskew * jch + ioffg + jch
  1170. * a_dim1], &ilda, &extra, &temp);
  1171. icol = jch;
  1172. /* L150: */
  1173. }
  1174. /* L160: */
  1175. }
  1176. /* L170: */
  1177. }
  1178. /* If we need lower triangle, copy from upper. Note that */
  1179. /* the order of copying is chosen to work for 'q' -> 'b' */
  1180. if (ipack != ipackg && ipack != 3) {
  1181. i__1 = *n;
  1182. for (jc = 1; jc <= i__1; ++jc) {
  1183. irow = ioffst - iskew * jc;
  1184. /* Computing MIN */
  1185. i__2 = *n, i__3 = jc + uub;
  1186. i__4 = f2cmin(i__2,i__3);
  1187. for (jr = jc; jr <= i__4; ++jr) {
  1188. a[jr + irow + jc * a_dim1] = a[jc - iskew * jr +
  1189. ioffg + jr * a_dim1];
  1190. /* L180: */
  1191. }
  1192. /* L190: */
  1193. }
  1194. if (ipack == 5) {
  1195. i__1 = *n;
  1196. for (jc = *n - uub + 1; jc <= i__1; ++jc) {
  1197. i__4 = uub + 1;
  1198. for (jr = *n + 2 - jc; jr <= i__4; ++jr) {
  1199. a[jr + jc * a_dim1] = 0.f;
  1200. /* L200: */
  1201. }
  1202. /* L210: */
  1203. }
  1204. }
  1205. if (ipackg == 6) {
  1206. ipackg = ipack;
  1207. } else {
  1208. ipackg = 0;
  1209. }
  1210. }
  1211. } else {
  1212. /* Bottom-Up -- Generate Lower triangle only */
  1213. if (ipack >= 5) {
  1214. ipackg = 5;
  1215. if (ipack == 6) {
  1216. ioffg = 1;
  1217. }
  1218. } else {
  1219. ipackg = 2;
  1220. }
  1221. i__1 = ilda + 1;
  1222. scopy_(&mnmin, &d__[1], &c__1, &a[1 - iskew + ioffg + a_dim1],
  1223. &i__1);
  1224. i__1 = uub;
  1225. for (k = 1; k <= i__1; ++k) {
  1226. for (jc = *n - 1; jc >= 1; --jc) {
  1227. /* Computing MIN */
  1228. i__4 = *n + 1 - jc, i__2 = k + 2;
  1229. il = f2cmin(i__4,i__2);
  1230. extra = 0.f;
  1231. temp = a[(1 - iskew) * jc + 1 + ioffg + jc * a_dim1];
  1232. angle = slarnd_(&c__1, &iseed[1]) *
  1233. 6.2831853071795864769252867663f;
  1234. c__ = cos(angle);
  1235. s = -sin(angle);
  1236. L__1 = *n - jc > k;
  1237. slarot_(&c_false, &c_true, &L__1, &il, &c__, &s, &a[(
  1238. 1 - iskew) * jc + ioffg + jc * a_dim1], &ilda,
  1239. &temp, &extra);
  1240. /* Computing MAX */
  1241. i__4 = 1, i__2 = jc - k + 1;
  1242. icol = f2cmax(i__4,i__2);
  1243. i__4 = jc + 2 - icol;
  1244. slarot_(&c_true, &c_false, &c_true, &i__4, &c__, &s, &
  1245. a[jc - iskew * icol + ioffg + icol * a_dim1],
  1246. &ilda, &dummy, &temp);
  1247. /* Chase EXTRA back down the matrix */
  1248. icol = jc;
  1249. i__4 = *n - 1;
  1250. i__2 = k;
  1251. for (jch = jc + k; i__2 < 0 ? jch >= i__4 : jch <=
  1252. i__4; jch += i__2) {
  1253. slartg_(&a[jch - iskew * icol + ioffg + icol *
  1254. a_dim1], &extra, &c__, &s, &dummy);
  1255. temp = a[(1 - iskew) * jch + 1 + ioffg + jch *
  1256. a_dim1];
  1257. i__3 = k + 2;
  1258. slarot_(&c_true, &c_true, &c_true, &i__3, &c__, &
  1259. s, &a[jch - iskew * icol + ioffg + icol *
  1260. a_dim1], &ilda, &extra, &temp);
  1261. /* Computing MIN */
  1262. i__3 = *n + 1 - jch, i__5 = k + 2;
  1263. il = f2cmin(i__3,i__5);
  1264. extra = 0.f;
  1265. L__1 = *n - jch > k;
  1266. slarot_(&c_false, &c_true, &L__1, &il, &c__, &s, &
  1267. a[(1 - iskew) * jch + ioffg + jch *
  1268. a_dim1], &ilda, &temp, &extra);
  1269. icol = jch;
  1270. /* L220: */
  1271. }
  1272. /* L230: */
  1273. }
  1274. /* L240: */
  1275. }
  1276. /* If we need upper triangle, copy from lower. Note that */
  1277. /* the order of copying is chosen to work for 'b' -> 'q' */
  1278. if (ipack != ipackg && ipack != 4) {
  1279. for (jc = *n; jc >= 1; --jc) {
  1280. irow = ioffst - iskew * jc;
  1281. /* Computing MAX */
  1282. i__2 = 1, i__4 = jc - uub;
  1283. i__1 = f2cmax(i__2,i__4);
  1284. for (jr = jc; jr >= i__1; --jr) {
  1285. a[jr + irow + jc * a_dim1] = a[jc - iskew * jr +
  1286. ioffg + jr * a_dim1];
  1287. /* L250: */
  1288. }
  1289. /* L260: */
  1290. }
  1291. if (ipack == 6) {
  1292. i__1 = uub;
  1293. for (jc = 1; jc <= i__1; ++jc) {
  1294. i__2 = uub + 1 - jc;
  1295. for (jr = 1; jr <= i__2; ++jr) {
  1296. a[jr + jc * a_dim1] = 0.f;
  1297. /* L270: */
  1298. }
  1299. /* L280: */
  1300. }
  1301. }
  1302. if (ipackg == 5) {
  1303. ipackg = ipack;
  1304. } else {
  1305. ipackg = 0;
  1306. }
  1307. }
  1308. }
  1309. }
  1310. } else {
  1311. /* 4) Generate Banded Matrix by first */
  1312. /* Rotating by random Unitary matrices, */
  1313. /* then reducing the bandwidth using Householder */
  1314. /* transformations. */
  1315. /* Note: we should get here only if LDA .ge. N */
  1316. if (isym == 1) {
  1317. /* Non-symmetric -- A = U D V */
  1318. slagge_(&mr, &nc, &llb, &uub, &d__[1], &a[a_offset], lda, &iseed[
  1319. 1], &work[1], &iinfo);
  1320. } else {
  1321. /* Symmetric -- A = U D U' */
  1322. slagsy_(m, &llb, &d__[1], &a[a_offset], lda, &iseed[1], &work[1],
  1323. &iinfo);
  1324. }
  1325. if (iinfo != 0) {
  1326. *info = 3;
  1327. return;
  1328. }
  1329. }
  1330. /* 5) Pack the matrix */
  1331. if (ipack != ipackg) {
  1332. if (ipack == 1) {
  1333. /* 'U' -- Upper triangular, not packed */
  1334. i__1 = *m;
  1335. for (j = 1; j <= i__1; ++j) {
  1336. i__2 = *m;
  1337. for (i__ = j + 1; i__ <= i__2; ++i__) {
  1338. a[i__ + j * a_dim1] = 0.f;
  1339. /* L290: */
  1340. }
  1341. /* L300: */
  1342. }
  1343. } else if (ipack == 2) {
  1344. /* 'L' -- Lower triangular, not packed */
  1345. i__1 = *m;
  1346. for (j = 2; j <= i__1; ++j) {
  1347. i__2 = j - 1;
  1348. for (i__ = 1; i__ <= i__2; ++i__) {
  1349. a[i__ + j * a_dim1] = 0.f;
  1350. /* L310: */
  1351. }
  1352. /* L320: */
  1353. }
  1354. } else if (ipack == 3) {
  1355. /* 'C' -- Upper triangle packed Columnwise. */
  1356. icol = 1;
  1357. irow = 0;
  1358. i__1 = *m;
  1359. for (j = 1; j <= i__1; ++j) {
  1360. i__2 = j;
  1361. for (i__ = 1; i__ <= i__2; ++i__) {
  1362. ++irow;
  1363. if (irow > *lda) {
  1364. irow = 1;
  1365. ++icol;
  1366. }
  1367. a[irow + icol * a_dim1] = a[i__ + j * a_dim1];
  1368. /* L330: */
  1369. }
  1370. /* L340: */
  1371. }
  1372. } else if (ipack == 4) {
  1373. /* 'R' -- Lower triangle packed Columnwise. */
  1374. icol = 1;
  1375. irow = 0;
  1376. i__1 = *m;
  1377. for (j = 1; j <= i__1; ++j) {
  1378. i__2 = *m;
  1379. for (i__ = j; i__ <= i__2; ++i__) {
  1380. ++irow;
  1381. if (irow > *lda) {
  1382. irow = 1;
  1383. ++icol;
  1384. }
  1385. a[irow + icol * a_dim1] = a[i__ + j * a_dim1];
  1386. /* L350: */
  1387. }
  1388. /* L360: */
  1389. }
  1390. } else if (ipack >= 5) {
  1391. /* 'B' -- The lower triangle is packed as a band matrix. */
  1392. /* 'Q' -- The upper triangle is packed as a band matrix. */
  1393. /* 'Z' -- The whole matrix is packed as a band matrix. */
  1394. if (ipack == 5) {
  1395. uub = 0;
  1396. }
  1397. if (ipack == 6) {
  1398. llb = 0;
  1399. }
  1400. i__1 = uub;
  1401. for (j = 1; j <= i__1; ++j) {
  1402. /* Computing MIN */
  1403. i__2 = j + llb;
  1404. for (i__ = f2cmin(i__2,*m); i__ >= 1; --i__) {
  1405. a[i__ - j + uub + 1 + j * a_dim1] = a[i__ + j * a_dim1];
  1406. /* L370: */
  1407. }
  1408. /* L380: */
  1409. }
  1410. i__1 = *n;
  1411. for (j = uub + 2; j <= i__1; ++j) {
  1412. /* Computing MIN */
  1413. i__4 = j + llb;
  1414. i__2 = f2cmin(i__4,*m);
  1415. for (i__ = j - uub; i__ <= i__2; ++i__) {
  1416. a[i__ - j + uub + 1 + j * a_dim1] = a[i__ + j * a_dim1];
  1417. /* L390: */
  1418. }
  1419. /* L400: */
  1420. }
  1421. }
  1422. /* If packed, zero out extraneous elements. */
  1423. /* Symmetric/Triangular Packed -- */
  1424. /* zero out everything after A(IROW,ICOL) */
  1425. if (ipack == 3 || ipack == 4) {
  1426. i__1 = *m;
  1427. for (jc = icol; jc <= i__1; ++jc) {
  1428. i__2 = *lda;
  1429. for (jr = irow + 1; jr <= i__2; ++jr) {
  1430. a[jr + jc * a_dim1] = 0.f;
  1431. /* L410: */
  1432. }
  1433. irow = 0;
  1434. /* L420: */
  1435. }
  1436. } else if (ipack >= 5) {
  1437. /* Packed Band -- */
  1438. /* 1st row is now in A( UUB+2-j, j), zero above it */
  1439. /* m-th row is now in A( M+UUB-j,j), zero below it */
  1440. /* last non-zero diagonal is now in A( UUB+LLB+1,j ), */
  1441. /* zero below it, too. */
  1442. ir1 = uub + llb + 2;
  1443. ir2 = uub + *m + 2;
  1444. i__1 = *n;
  1445. for (jc = 1; jc <= i__1; ++jc) {
  1446. i__2 = uub + 1 - jc;
  1447. for (jr = 1; jr <= i__2; ++jr) {
  1448. a[jr + jc * a_dim1] = 0.f;
  1449. /* L430: */
  1450. }
  1451. /* Computing MAX */
  1452. /* Computing MIN */
  1453. i__3 = ir1, i__5 = ir2 - jc;
  1454. i__2 = 1, i__4 = f2cmin(i__3,i__5);
  1455. i__6 = *lda;
  1456. for (jr = f2cmax(i__2,i__4); jr <= i__6; ++jr) {
  1457. a[jr + jc * a_dim1] = 0.f;
  1458. /* L440: */
  1459. }
  1460. /* L450: */
  1461. }
  1462. }
  1463. }
  1464. return;
  1465. /* End of SLATMS */
  1466. } /* slatms_ */