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.

claror.c 19 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659
  1. #include <math.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <stdio.h>
  5. #include <complex.h>
  6. #ifdef complex
  7. #undef complex
  8. #endif
  9. #ifdef I
  10. #undef I
  11. #endif
  12. #if defined(_WIN64)
  13. typedef long long BLASLONG;
  14. typedef unsigned long long BLASULONG;
  15. #else
  16. typedef long BLASLONG;
  17. typedef unsigned long BLASULONG;
  18. #endif
  19. #ifdef LAPACK_ILP64
  20. typedef BLASLONG blasint;
  21. #if defined(_WIN64)
  22. #define blasabs(x) llabs(x)
  23. #else
  24. #define blasabs(x) labs(x)
  25. #endif
  26. #else
  27. typedef int blasint;
  28. #define blasabs(x) abs(x)
  29. #endif
  30. typedef blasint integer;
  31. typedef unsigned int uinteger;
  32. typedef char *address;
  33. typedef short int shortint;
  34. typedef float real;
  35. typedef double doublereal;
  36. typedef struct { real r, i; } complex;
  37. typedef struct { doublereal r, i; } doublecomplex;
  38. #ifdef _MSC_VER
  39. static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
  40. static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
  41. static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
  42. static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
  43. #else
  44. static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
  45. static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
  46. static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
  47. static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
  48. #endif
  49. #define pCf(z) (*_pCf(z))
  50. #define pCd(z) (*_pCd(z))
  51. typedef blasint logical;
  52. typedef char logical1;
  53. typedef char integer1;
  54. #define TRUE_ (1)
  55. #define FALSE_ (0)
  56. /* Extern is for use with -E */
  57. #ifndef Extern
  58. #define Extern extern
  59. #endif
  60. /* I/O stuff */
  61. typedef int flag;
  62. typedef int ftnlen;
  63. typedef int ftnint;
  64. /*external read, write*/
  65. typedef struct
  66. { flag cierr;
  67. ftnint ciunit;
  68. flag ciend;
  69. char *cifmt;
  70. ftnint cirec;
  71. } cilist;
  72. /*internal read, write*/
  73. typedef struct
  74. { flag icierr;
  75. char *iciunit;
  76. flag iciend;
  77. char *icifmt;
  78. ftnint icirlen;
  79. ftnint icirnum;
  80. } icilist;
  81. /*open*/
  82. typedef struct
  83. { flag oerr;
  84. ftnint ounit;
  85. char *ofnm;
  86. ftnlen ofnmlen;
  87. char *osta;
  88. char *oacc;
  89. char *ofm;
  90. ftnint orl;
  91. char *oblnk;
  92. } olist;
  93. /*close*/
  94. typedef struct
  95. { flag cerr;
  96. ftnint cunit;
  97. char *csta;
  98. } cllist;
  99. /*rewind, backspace, endfile*/
  100. typedef struct
  101. { flag aerr;
  102. ftnint aunit;
  103. } alist;
  104. /* inquire */
  105. typedef struct
  106. { flag inerr;
  107. ftnint inunit;
  108. char *infile;
  109. ftnlen infilen;
  110. ftnint *inex; /*parameters in standard's order*/
  111. ftnint *inopen;
  112. ftnint *innum;
  113. ftnint *innamed;
  114. char *inname;
  115. ftnlen innamlen;
  116. char *inacc;
  117. ftnlen inacclen;
  118. char *inseq;
  119. ftnlen inseqlen;
  120. char *indir;
  121. ftnlen indirlen;
  122. char *infmt;
  123. ftnlen infmtlen;
  124. char *inform;
  125. ftnint informlen;
  126. char *inunf;
  127. ftnlen inunflen;
  128. ftnint *inrecl;
  129. ftnint *innrec;
  130. char *inblank;
  131. ftnlen inblanklen;
  132. } inlist;
  133. #define VOID void
  134. union Multitype { /* for multiple entry points */
  135. integer1 g;
  136. shortint h;
  137. integer i;
  138. /* longint j; */
  139. real r;
  140. doublereal d;
  141. complex c;
  142. doublecomplex z;
  143. };
  144. typedef union Multitype Multitype;
  145. struct Vardesc { /* for Namelist */
  146. char *name;
  147. char *addr;
  148. ftnlen *dims;
  149. int type;
  150. };
  151. typedef struct Vardesc Vardesc;
  152. struct Namelist {
  153. char *name;
  154. Vardesc **vars;
  155. int nvars;
  156. };
  157. typedef struct Namelist Namelist;
  158. #define abs(x) ((x) >= 0 ? (x) : -(x))
  159. #define dabs(x) (fabs(x))
  160. #define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
  161. #define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
  162. #define dmin(a,b) (f2cmin(a,b))
  163. #define dmax(a,b) (f2cmax(a,b))
  164. #define bit_test(a,b) ((a) >> (b) & 1)
  165. #define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
  166. #define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
  167. #define abort_() { sig_die("Fortran abort routine called", 1); }
  168. #define c_abs(z) (cabsf(Cf(z)))
  169. #define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
  170. #ifdef _MSC_VER
  171. #define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
  172. #define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/df(b)._Val[1]);}
  173. #else
  174. #define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
  175. #define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
  176. #endif
  177. #define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
  178. #define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
  179. #define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
  180. //#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
  181. #define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
  182. #define d_abs(x) (fabs(*(x)))
  183. #define d_acos(x) (acos(*(x)))
  184. #define d_asin(x) (asin(*(x)))
  185. #define d_atan(x) (atan(*(x)))
  186. #define d_atn2(x, y) (atan2(*(x),*(y)))
  187. #define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
  188. #define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
  189. #define d_cos(x) (cos(*(x)))
  190. #define d_cosh(x) (cosh(*(x)))
  191. #define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
  192. #define d_exp(x) (exp(*(x)))
  193. #define d_imag(z) (cimag(Cd(z)))
  194. #define r_imag(z) (cimagf(Cf(z)))
  195. #define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
  196. #define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
  197. #define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
  198. #define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
  199. #define d_log(x) (log(*(x)))
  200. #define d_mod(x, y) (fmod(*(x), *(y)))
  201. #define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
  202. #define d_nint(x) u_nint(*(x))
  203. #define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
  204. #define d_sign(a,b) u_sign(*(a),*(b))
  205. #define r_sign(a,b) u_sign(*(a),*(b))
  206. #define d_sin(x) (sin(*(x)))
  207. #define d_sinh(x) (sinh(*(x)))
  208. #define d_sqrt(x) (sqrt(*(x)))
  209. #define d_tan(x) (tan(*(x)))
  210. #define d_tanh(x) (tanh(*(x)))
  211. #define i_abs(x) abs(*(x))
  212. #define i_dnnt(x) ((integer)u_nint(*(x)))
  213. #define i_len(s, n) (n)
  214. #define i_nint(x) ((integer)u_nint(*(x)))
  215. #define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
  216. #define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
  217. #define pow_si(B,E) spow_ui(*(B),*(E))
  218. #define pow_ri(B,E) spow_ui(*(B),*(E))
  219. #define pow_di(B,E) dpow_ui(*(B),*(E))
  220. #define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
  221. #define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
  222. #define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
  223. #define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
  224. #define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
  225. #define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
  226. #define sig_die(s, kill) { exit(1); }
  227. #define s_stop(s, n) {exit(0);}
  228. #define z_abs(z) (cabs(Cd(z)))
  229. #define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
  230. #define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
  231. #define myexit_() break;
  232. #define mycycle() continue;
  233. #define myceiling(w) {ceil(w)}
  234. #define myhuge(w) {HUGE_VAL}
  235. //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);}
  236. #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)}
  237. /* procedure parameter types for -A and -C++ */
  238. /* Table of constant values */
  239. static complex c_b1 = {0.f,0.f};
  240. static complex c_b2 = {1.f,0.f};
  241. static integer c__3 = 3;
  242. static integer c__1 = 1;
  243. /* > \brief \b CLAROR */
  244. /* =========== DOCUMENTATION =========== */
  245. /* Online html documentation available at */
  246. /* http://www.netlib.org/lapack/explore-html/ */
  247. /* Definition: */
  248. /* =========== */
  249. /* SUBROUTINE CLAROR( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO ) */
  250. /* CHARACTER INIT, SIDE */
  251. /* INTEGER INFO, LDA, M, N */
  252. /* INTEGER ISEED( 4 ) */
  253. /* COMPLEX A( LDA, * ), X( * ) */
  254. /* > \par Purpose: */
  255. /* ============= */
  256. /* > */
  257. /* > \verbatim */
  258. /* > */
  259. /* > CLAROR pre- or post-multiplies an M by N matrix A by a random */
  260. /* > unitary matrix U, overwriting A. A may optionally be */
  261. /* > initialized to the identity matrix before multiplying by U. */
  262. /* > U is generated using the method of G.W. Stewart */
  263. /* > ( SIAM J. Numer. Anal. 17, 1980, pp. 403-409 ). */
  264. /* > (BLAS-2 version) */
  265. /* > \endverbatim */
  266. /* Arguments: */
  267. /* ========== */
  268. /* > \param[in] SIDE */
  269. /* > \verbatim */
  270. /* > SIDE is CHARACTER*1 */
  271. /* > SIDE specifies whether A is multiplied on the left or right */
  272. /* > by U. */
  273. /* > SIDE = 'L' Multiply A on the left (premultiply) by U */
  274. /* > SIDE = 'R' Multiply A on the right (postmultiply) by UC> SIDE = 'C' Multiply A on the lef
  275. t by U and the right by UC> SIDE = 'T' Multiply A on the left by U and the right by U' */
  276. /* > Not modified. */
  277. /* > \endverbatim */
  278. /* > */
  279. /* > \param[in] INIT */
  280. /* > \verbatim */
  281. /* > INIT is CHARACTER*1 */
  282. /* > INIT specifies whether or not A should be initialized to */
  283. /* > the identity matrix. */
  284. /* > INIT = 'I' Initialize A to (a section of) the */
  285. /* > identity matrix before applying U. */
  286. /* > INIT = 'N' No initialization. Apply U to the */
  287. /* > input matrix A. */
  288. /* > */
  289. /* > INIT = 'I' may be used to generate square (i.e., unitary) */
  290. /* > or rectangular orthogonal matrices (orthogonality being */
  291. /* > in the sense of CDOTC): */
  292. /* > */
  293. /* > For square matrices, M=N, and SIDE many be either 'L' or */
  294. /* > 'R'; the rows will be orthogonal to each other, as will the */
  295. /* > columns. */
  296. /* > For rectangular matrices where M < N, SIDE = 'R' will */
  297. /* > produce a dense matrix whose rows will be orthogonal and */
  298. /* > whose columns will not, while SIDE = 'L' will produce a */
  299. /* > matrix whose rows will be orthogonal, and whose first M */
  300. /* > columns will be orthogonal, the remaining columns being */
  301. /* > zero. */
  302. /* > For matrices where M > N, just use the previous */
  303. /* > explanation, interchanging 'L' and 'R' and "rows" and */
  304. /* > "columns". */
  305. /* > */
  306. /* > Not modified. */
  307. /* > \endverbatim */
  308. /* > */
  309. /* > \param[in] M */
  310. /* > \verbatim */
  311. /* > M is INTEGER */
  312. /* > Number of rows of A. Not modified. */
  313. /* > \endverbatim */
  314. /* > */
  315. /* > \param[in] N */
  316. /* > \verbatim */
  317. /* > N is INTEGER */
  318. /* > Number of columns of A. Not modified. */
  319. /* > \endverbatim */
  320. /* > */
  321. /* > \param[in,out] A */
  322. /* > \verbatim */
  323. /* > A is COMPLEX array, dimension ( LDA, N ) */
  324. /* > Input and output array. Overwritten by U A ( if SIDE = 'L' ) */
  325. /* > or by A U ( if SIDE = 'R' ) */
  326. /* > or by U A U* ( if SIDE = 'C') */
  327. /* > or by U A U' ( if SIDE = 'T') on exit. */
  328. /* > \endverbatim */
  329. /* > */
  330. /* > \param[in] LDA */
  331. /* > \verbatim */
  332. /* > LDA is INTEGER */
  333. /* > Leading dimension of A. Must be at least MAX ( 1, M ). */
  334. /* > Not modified. */
  335. /* > \endverbatim */
  336. /* > */
  337. /* > \param[in,out] ISEED */
  338. /* > \verbatim */
  339. /* > ISEED is INTEGER array, dimension ( 4 ) */
  340. /* > On entry ISEED specifies the seed of the random number */
  341. /* > generator. The array elements should be between 0 and 4095; */
  342. /* > if not they will be reduced mod 4096. Also, ISEED(4) must */
  343. /* > be odd. The random number generator uses a linear */
  344. /* > congruential sequence limited to small integers, and so */
  345. /* > should produce machine independent random numbers. The */
  346. /* > values of ISEED are changed on exit, and can be used in the */
  347. /* > next call to CLAROR to continue the same random number */
  348. /* > sequence. */
  349. /* > Modified. */
  350. /* > \endverbatim */
  351. /* > */
  352. /* > \param[out] X */
  353. /* > \verbatim */
  354. /* > X is COMPLEX array, dimension ( 3*MAX( M, N ) ) */
  355. /* > Workspace. Of length: */
  356. /* > 2*M + N if SIDE = 'L', */
  357. /* > 2*N + M if SIDE = 'R', */
  358. /* > 3*N if SIDE = 'C' or 'T'. */
  359. /* > Modified. */
  360. /* > \endverbatim */
  361. /* > */
  362. /* > \param[out] INFO */
  363. /* > \verbatim */
  364. /* > INFO is INTEGER */
  365. /* > An error flag. It is set to: */
  366. /* > 0 if no error. */
  367. /* > 1 if CLARND returned a bad random number (installation */
  368. /* > problem) */
  369. /* > -1 if SIDE is not L, R, C, or T. */
  370. /* > -3 if M is negative. */
  371. /* > -4 if N is negative or if SIDE is C or T and N is not equal */
  372. /* > to M. */
  373. /* > -6 if LDA is less than M. */
  374. /* > \endverbatim */
  375. /* Authors: */
  376. /* ======== */
  377. /* > \author Univ. of Tennessee */
  378. /* > \author Univ. of California Berkeley */
  379. /* > \author Univ. of Colorado Denver */
  380. /* > \author NAG Ltd. */
  381. /* > \date December 2016 */
  382. /* > \ingroup complex_matgen */
  383. /* ===================================================================== */
  384. /* Subroutine */ void claror_(char *side, char *init, integer *m, integer *n,
  385. complex *a, integer *lda, integer *iseed, complex *x, integer *info)
  386. {
  387. /* System generated locals */
  388. integer a_dim1, a_offset, i__1, i__2, i__3;
  389. complex q__1, q__2;
  390. /* Local variables */
  391. integer kbeg, jcol;
  392. real xabs;
  393. integer irow, j;
  394. extern /* Subroutine */ void cgerc_(integer *, integer *, complex *,
  395. complex *, integer *, complex *, integer *, complex *, integer *),
  396. cscal_(integer *, complex *, complex *, integer *);
  397. extern logical lsame_(char *, char *);
  398. extern /* Subroutine */ void cgemv_(char *, integer *, integer *, complex *
  399. , complex *, integer *, complex *, integer *, complex *, complex *
  400. , integer *);
  401. complex csign;
  402. integer ixfrm, itype, nxfrm;
  403. real xnorm;
  404. extern real scnrm2_(integer *, complex *, integer *);
  405. extern /* Subroutine */ void clacgv_(integer *, complex *, integer *);
  406. //extern /* Complex */ VOID clarnd_(complex *, integer *, integer *);
  407. extern complex clarnd_(integer *, integer *);
  408. extern /* Subroutine */ void claset_(char *, integer *, integer *, complex
  409. *, complex *, complex *, integer *);
  410. extern int xerbla_(char *, integer *, ftnlen);
  411. real factor;
  412. complex xnorms;
  413. /* -- LAPACK auxiliary routine (version 3.7.0) -- */
  414. /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
  415. /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
  416. /* December 2016 */
  417. /* ===================================================================== */
  418. /* Parameter adjustments */
  419. a_dim1 = *lda;
  420. a_offset = 1 + a_dim1 * 1;
  421. a -= a_offset;
  422. --iseed;
  423. --x;
  424. /* Function Body */
  425. *info = 0;
  426. if (*n == 0 || *m == 0) {
  427. return;
  428. }
  429. itype = 0;
  430. if (lsame_(side, "L")) {
  431. itype = 1;
  432. } else if (lsame_(side, "R")) {
  433. itype = 2;
  434. } else if (lsame_(side, "C")) {
  435. itype = 3;
  436. } else if (lsame_(side, "T")) {
  437. itype = 4;
  438. }
  439. /* Check for argument errors. */
  440. if (itype == 0) {
  441. *info = -1;
  442. } else if (*m < 0) {
  443. *info = -3;
  444. } else if (*n < 0 || itype == 3 && *n != *m) {
  445. *info = -4;
  446. } else if (*lda < *m) {
  447. *info = -6;
  448. }
  449. if (*info != 0) {
  450. i__1 = -(*info);
  451. xerbla_("CLAROR", &i__1, 6);
  452. return;
  453. }
  454. if (itype == 1) {
  455. nxfrm = *m;
  456. } else {
  457. nxfrm = *n;
  458. }
  459. /* Initialize A to the identity matrix if desired */
  460. if (lsame_(init, "I")) {
  461. claset_("Full", m, n, &c_b1, &c_b2, &a[a_offset], lda);
  462. }
  463. /* If no rotation possible, still multiply by */
  464. /* a random complex number from the circle |x| = 1 */
  465. /* 2) Compute Rotation by computing Householder */
  466. /* Transformations H(2), H(3), ..., H(n). Note that the */
  467. /* order in which they are computed is irrelevant. */
  468. i__1 = nxfrm;
  469. for (j = 1; j <= i__1; ++j) {
  470. i__2 = j;
  471. x[i__2].r = 0.f, x[i__2].i = 0.f;
  472. /* L40: */
  473. }
  474. i__1 = nxfrm;
  475. for (ixfrm = 2; ixfrm <= i__1; ++ixfrm) {
  476. kbeg = nxfrm - ixfrm + 1;
  477. /* Generate independent normal( 0, 1 ) random numbers */
  478. i__2 = nxfrm;
  479. for (j = kbeg; j <= i__2; ++j) {
  480. i__3 = j;
  481. //clarnd_(&q__1, &c__3, &iseed[1]);
  482. q__1=clarnd_(&c__3, &iseed[1]);
  483. x[i__3].r = q__1.r, x[i__3].i = q__1.i;
  484. /* L50: */
  485. }
  486. /* Generate a Householder transformation from the random vector X */
  487. xnorm = scnrm2_(&ixfrm, &x[kbeg], &c__1);
  488. xabs = c_abs(&x[kbeg]);
  489. if (xabs != 0.f) {
  490. i__2 = kbeg;
  491. q__1.r = x[i__2].r / xabs, q__1.i = x[i__2].i / xabs;
  492. csign.r = q__1.r, csign.i = q__1.i;
  493. } else {
  494. csign.r = 1.f, csign.i = 0.f;
  495. }
  496. q__1.r = xnorm * csign.r, q__1.i = xnorm * csign.i;
  497. xnorms.r = q__1.r, xnorms.i = q__1.i;
  498. i__2 = nxfrm + kbeg;
  499. q__1.r = -csign.r, q__1.i = -csign.i;
  500. x[i__2].r = q__1.r, x[i__2].i = q__1.i;
  501. factor = xnorm * (xnorm + xabs);
  502. if (abs(factor) < 1e-20f) {
  503. *info = 1;
  504. i__2 = -(*info);
  505. xerbla_("CLAROR", &i__2, 6);
  506. return;
  507. } else {
  508. factor = 1.f / factor;
  509. }
  510. i__2 = kbeg;
  511. i__3 = kbeg;
  512. q__1.r = x[i__3].r + xnorms.r, q__1.i = x[i__3].i + xnorms.i;
  513. x[i__2].r = q__1.r, x[i__2].i = q__1.i;
  514. /* Apply Householder transformation to A */
  515. if (itype == 1 || itype == 3 || itype == 4) {
  516. /* Apply H(k) on the left of A */
  517. cgemv_("C", &ixfrm, n, &c_b2, &a[kbeg + a_dim1], lda, &x[kbeg], &
  518. c__1, &c_b1, &x[(nxfrm << 1) + 1], &c__1);
  519. q__2.r = factor, q__2.i = 0.f;
  520. q__1.r = -q__2.r, q__1.i = -q__2.i;
  521. cgerc_(&ixfrm, n, &q__1, &x[kbeg], &c__1, &x[(nxfrm << 1) + 1], &
  522. c__1, &a[kbeg + a_dim1], lda);
  523. }
  524. if (itype >= 2 && itype <= 4) {
  525. /* Apply H(k)* (or H(k)') on the right of A */
  526. if (itype == 4) {
  527. clacgv_(&ixfrm, &x[kbeg], &c__1);
  528. }
  529. cgemv_("N", m, &ixfrm, &c_b2, &a[kbeg * a_dim1 + 1], lda, &x[kbeg]
  530. , &c__1, &c_b1, &x[(nxfrm << 1) + 1], &c__1);
  531. q__2.r = factor, q__2.i = 0.f;
  532. q__1.r = -q__2.r, q__1.i = -q__2.i;
  533. cgerc_(m, &ixfrm, &q__1, &x[(nxfrm << 1) + 1], &c__1, &x[kbeg], &
  534. c__1, &a[kbeg * a_dim1 + 1], lda);
  535. }
  536. /* L60: */
  537. }
  538. //clarnd_(&q__1, &c__3, &iseed[1]);
  539. q__1=clarnd_(&c__3, &iseed[1]);
  540. x[1].r = q__1.r, x[1].i = q__1.i;
  541. xabs = c_abs(&x[1]);
  542. if (xabs != 0.f) {
  543. q__1.r = x[1].r / xabs, q__1.i = x[1].i / xabs;
  544. csign.r = q__1.r, csign.i = q__1.i;
  545. } else {
  546. csign.r = 1.f, csign.i = 0.f;
  547. }
  548. i__1 = nxfrm << 1;
  549. x[i__1].r = csign.r, x[i__1].i = csign.i;
  550. /* Scale the matrix A by D. */
  551. if (itype == 1 || itype == 3 || itype == 4) {
  552. i__1 = *m;
  553. for (irow = 1; irow <= i__1; ++irow) {
  554. r_cnjg(&q__1, &x[nxfrm + irow]);
  555. cscal_(n, &q__1, &a[irow + a_dim1], lda);
  556. /* L70: */
  557. }
  558. }
  559. if (itype == 2 || itype == 3) {
  560. i__1 = *n;
  561. for (jcol = 1; jcol <= i__1; ++jcol) {
  562. cscal_(m, &x[nxfrm + jcol], &a[jcol * a_dim1 + 1], &c__1);
  563. /* L80: */
  564. }
  565. }
  566. if (itype == 4) {
  567. i__1 = *n;
  568. for (jcol = 1; jcol <= i__1; ++jcol) {
  569. r_cnjg(&q__1, &x[nxfrm + jcol]);
  570. cscal_(m, &q__1, &a[jcol * a_dim1 + 1], &c__1);
  571. /* L90: */
  572. }
  573. }
  574. return;
  575. /* End of CLAROR */
  576. } /* claror_ */