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