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.

slansf.c 44 kB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603
  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. static char junk[] = "\n@(#)LIBF77 VERSION 19990503\n";
  230. #define z_abs(z) (cabs(Cd(z)))
  231. #define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
  232. #define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
  233. #define myexit_() break;
  234. #define mycycle() continue;
  235. #define myceiling(w) {ceil(w)}
  236. #define myhuge(w) {HUGE_VAL}
  237. //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);}
  238. #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)}
  239. /* procedure parameter types for -A and -C++ */
  240. #define F2C_proc_par_types 1
  241. #ifdef __cplusplus
  242. typedef logical (*L_fp)(...);
  243. #else
  244. typedef logical (*L_fp)();
  245. #endif
  246. static float spow_ui(float x, integer n) {
  247. float pow=1.0; unsigned long int u;
  248. if(n != 0) {
  249. if(n < 0) n = -n, x = 1/x;
  250. for(u = n; ; ) {
  251. if(u & 01) pow *= x;
  252. if(u >>= 1) x *= x;
  253. else break;
  254. }
  255. }
  256. return pow;
  257. }
  258. static double dpow_ui(double x, integer n) {
  259. double pow=1.0; unsigned long int u;
  260. if(n != 0) {
  261. if(n < 0) n = -n, x = 1/x;
  262. for(u = n; ; ) {
  263. if(u & 01) pow *= x;
  264. if(u >>= 1) x *= x;
  265. else break;
  266. }
  267. }
  268. return pow;
  269. }
  270. #ifdef _MSC_VER
  271. static _Fcomplex cpow_ui(complex x, integer n) {
  272. complex pow={1.0,0.0}; unsigned long int u;
  273. if(n != 0) {
  274. if(n < 0) n = -n, x.r = 1/x.r, x.i=1/x.i;
  275. for(u = n; ; ) {
  276. if(u & 01) pow.r *= x.r, pow.i *= x.i;
  277. if(u >>= 1) x.r *= x.r, x.i *= x.i;
  278. else break;
  279. }
  280. }
  281. _Fcomplex p={pow.r, pow.i};
  282. return p;
  283. }
  284. #else
  285. static _Complex float cpow_ui(_Complex float x, integer n) {
  286. _Complex float pow=1.0; unsigned long int u;
  287. if(n != 0) {
  288. if(n < 0) n = -n, x = 1/x;
  289. for(u = n; ; ) {
  290. if(u & 01) pow *= x;
  291. if(u >>= 1) x *= x;
  292. else break;
  293. }
  294. }
  295. return pow;
  296. }
  297. #endif
  298. #ifdef _MSC_VER
  299. static _Dcomplex zpow_ui(_Dcomplex x, integer n) {
  300. _Dcomplex pow={1.0,0.0}; unsigned long int u;
  301. if(n != 0) {
  302. if(n < 0) n = -n, x._Val[0] = 1/x._Val[0], x._Val[1] =1/x._Val[1];
  303. for(u = n; ; ) {
  304. if(u & 01) pow._Val[0] *= x._Val[0], pow._Val[1] *= x._Val[1];
  305. if(u >>= 1) x._Val[0] *= x._Val[0], x._Val[1] *= x._Val[1];
  306. else break;
  307. }
  308. }
  309. _Dcomplex p = {pow._Val[0], pow._Val[1]};
  310. return p;
  311. }
  312. #else
  313. static _Complex double zpow_ui(_Complex double x, integer n) {
  314. _Complex double pow=1.0; unsigned long int u;
  315. if(n != 0) {
  316. if(n < 0) n = -n, x = 1/x;
  317. for(u = n; ; ) {
  318. if(u & 01) pow *= x;
  319. if(u >>= 1) x *= x;
  320. else break;
  321. }
  322. }
  323. return pow;
  324. }
  325. #endif
  326. static integer pow_ii(integer x, integer n) {
  327. integer pow; unsigned long int u;
  328. if (n <= 0) {
  329. if (n == 0 || x == 1) pow = 1;
  330. else if (x != -1) pow = x == 0 ? 1/x : 0;
  331. else n = -n;
  332. }
  333. if ((n > 0) || !(n == 0 || x == 1 || x != -1)) {
  334. u = n;
  335. for(pow = 1; ; ) {
  336. if(u & 01) pow *= x;
  337. if(u >>= 1) x *= x;
  338. else break;
  339. }
  340. }
  341. return pow;
  342. }
  343. static integer dmaxloc_(double *w, integer s, integer e, integer *n)
  344. {
  345. double m; integer i, mi;
  346. for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
  347. if (w[i-1]>m) mi=i ,m=w[i-1];
  348. return mi-s+1;
  349. }
  350. static integer smaxloc_(float *w, integer s, integer e, integer *n)
  351. {
  352. float m; integer i, mi;
  353. for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
  354. if (w[i-1]>m) mi=i ,m=w[i-1];
  355. return mi-s+1;
  356. }
  357. static inline void cdotc_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
  358. integer n = *n_, incx = *incx_, incy = *incy_, i;
  359. #ifdef _MSC_VER
  360. _Fcomplex zdotc = {0.0, 0.0};
  361. if (incx == 1 && incy == 1) {
  362. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  363. zdotc._Val[0] += conjf(Cf(&x[i]))._Val[0] * Cf(&y[i])._Val[0];
  364. zdotc._Val[1] += conjf(Cf(&x[i]))._Val[1] * Cf(&y[i])._Val[1];
  365. }
  366. } else {
  367. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  368. zdotc._Val[0] += conjf(Cf(&x[i*incx]))._Val[0] * Cf(&y[i*incy])._Val[0];
  369. zdotc._Val[1] += conjf(Cf(&x[i*incx]))._Val[1] * Cf(&y[i*incy])._Val[1];
  370. }
  371. }
  372. pCf(z) = zdotc;
  373. }
  374. #else
  375. _Complex float zdotc = 0.0;
  376. if (incx == 1 && incy == 1) {
  377. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  378. zdotc += conjf(Cf(&x[i])) * Cf(&y[i]);
  379. }
  380. } else {
  381. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  382. zdotc += conjf(Cf(&x[i*incx])) * Cf(&y[i*incy]);
  383. }
  384. }
  385. pCf(z) = zdotc;
  386. }
  387. #endif
  388. static inline void zdotc_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
  389. integer n = *n_, incx = *incx_, incy = *incy_, i;
  390. #ifdef _MSC_VER
  391. _Dcomplex zdotc = {0.0, 0.0};
  392. if (incx == 1 && incy == 1) {
  393. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  394. zdotc._Val[0] += conj(Cd(&x[i]))._Val[0] * Cd(&y[i])._Val[0];
  395. zdotc._Val[1] += conj(Cd(&x[i]))._Val[1] * Cd(&y[i])._Val[1];
  396. }
  397. } else {
  398. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  399. zdotc._Val[0] += conj(Cd(&x[i*incx]))._Val[0] * Cd(&y[i*incy])._Val[0];
  400. zdotc._Val[1] += conj(Cd(&x[i*incx]))._Val[1] * Cd(&y[i*incy])._Val[1];
  401. }
  402. }
  403. pCd(z) = zdotc;
  404. }
  405. #else
  406. _Complex double zdotc = 0.0;
  407. if (incx == 1 && incy == 1) {
  408. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  409. zdotc += conj(Cd(&x[i])) * Cd(&y[i]);
  410. }
  411. } else {
  412. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  413. zdotc += conj(Cd(&x[i*incx])) * Cd(&y[i*incy]);
  414. }
  415. }
  416. pCd(z) = zdotc;
  417. }
  418. #endif
  419. static inline void cdotu_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
  420. integer n = *n_, incx = *incx_, incy = *incy_, i;
  421. #ifdef _MSC_VER
  422. _Fcomplex zdotc = {0.0, 0.0};
  423. if (incx == 1 && incy == 1) {
  424. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  425. zdotc._Val[0] += Cf(&x[i])._Val[0] * Cf(&y[i])._Val[0];
  426. zdotc._Val[1] += Cf(&x[i])._Val[1] * Cf(&y[i])._Val[1];
  427. }
  428. } else {
  429. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  430. zdotc._Val[0] += Cf(&x[i*incx])._Val[0] * Cf(&y[i*incy])._Val[0];
  431. zdotc._Val[1] += Cf(&x[i*incx])._Val[1] * Cf(&y[i*incy])._Val[1];
  432. }
  433. }
  434. pCf(z) = zdotc;
  435. }
  436. #else
  437. _Complex float zdotc = 0.0;
  438. if (incx == 1 && incy == 1) {
  439. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  440. zdotc += Cf(&x[i]) * Cf(&y[i]);
  441. }
  442. } else {
  443. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  444. zdotc += Cf(&x[i*incx]) * Cf(&y[i*incy]);
  445. }
  446. }
  447. pCf(z) = zdotc;
  448. }
  449. #endif
  450. static inline void zdotu_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
  451. integer n = *n_, incx = *incx_, incy = *incy_, i;
  452. #ifdef _MSC_VER
  453. _Dcomplex zdotc = {0.0, 0.0};
  454. if (incx == 1 && incy == 1) {
  455. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  456. zdotc._Val[0] += Cd(&x[i])._Val[0] * Cd(&y[i])._Val[0];
  457. zdotc._Val[1] += Cd(&x[i])._Val[1] * Cd(&y[i])._Val[1];
  458. }
  459. } else {
  460. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  461. zdotc._Val[0] += Cd(&x[i*incx])._Val[0] * Cd(&y[i*incy])._Val[0];
  462. zdotc._Val[1] += Cd(&x[i*incx])._Val[1] * Cd(&y[i*incy])._Val[1];
  463. }
  464. }
  465. pCd(z) = zdotc;
  466. }
  467. #else
  468. _Complex double zdotc = 0.0;
  469. if (incx == 1 && incy == 1) {
  470. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  471. zdotc += Cd(&x[i]) * Cd(&y[i]);
  472. }
  473. } else {
  474. for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
  475. zdotc += Cd(&x[i*incx]) * Cd(&y[i*incy]);
  476. }
  477. }
  478. pCd(z) = zdotc;
  479. }
  480. #endif
  481. /* -- translated by f2c (version 20000121).
  482. You must link the resulting object file with the libraries:
  483. -lf2c -lm (in that order)
  484. */
  485. /* Table of constant values */
  486. static integer c__1 = 1;
  487. /* > \brief \b SLANSF */
  488. /* =========== DOCUMENTATION =========== */
  489. /* Online html documentation available at */
  490. /* http://www.netlib.org/lapack/explore-html/ */
  491. /* > \htmlonly */
  492. /* > Download SLANSF + dependencies */
  493. /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slansf.
  494. f"> */
  495. /* > [TGZ]</a> */
  496. /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slansf.
  497. f"> */
  498. /* > [ZIP]</a> */
  499. /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slansf.
  500. f"> */
  501. /* > [TXT]</a> */
  502. /* > \endhtmlonly */
  503. /* Definition: */
  504. /* =========== */
  505. /* REAL FUNCTION SLANSF( NORM, TRANSR, UPLO, N, A, WORK ) */
  506. /* CHARACTER NORM, TRANSR, UPLO */
  507. /* INTEGER N */
  508. /* REAL A( 0: * ), WORK( 0: * ) */
  509. /* > \par Purpose: */
  510. /* ============= */
  511. /* > */
  512. /* > \verbatim */
  513. /* > */
  514. /* > SLANSF returns the value of the one norm, or the Frobenius norm, or */
  515. /* > the infinity norm, or the element of largest absolute value of a */
  516. /* > real symmetric matrix A in RFP format. */
  517. /* > \endverbatim */
  518. /* > */
  519. /* > \return SLANSF */
  520. /* > \verbatim */
  521. /* > */
  522. /* > SLANSF = ( f2cmax(abs(A(i,j))), NORM = 'M' or 'm' */
  523. /* > ( */
  524. /* > ( norm1(A), NORM = '1', 'O' or 'o' */
  525. /* > ( */
  526. /* > ( normI(A), NORM = 'I' or 'i' */
  527. /* > ( */
  528. /* > ( normF(A), NORM = 'F', 'f', 'E' or 'e' */
  529. /* > */
  530. /* > where norm1 denotes the one norm of a matrix (maximum column sum), */
  531. /* > normI denotes the infinity norm of a matrix (maximum row sum) and */
  532. /* > normF denotes the Frobenius norm of a matrix (square root of sum of */
  533. /* > squares). Note that f2cmax(abs(A(i,j))) is not a matrix norm. */
  534. /* > \endverbatim */
  535. /* Arguments: */
  536. /* ========== */
  537. /* > \param[in] NORM */
  538. /* > \verbatim */
  539. /* > NORM is CHARACTER*1 */
  540. /* > Specifies the value to be returned in SLANSF as described */
  541. /* > above. */
  542. /* > \endverbatim */
  543. /* > */
  544. /* > \param[in] TRANSR */
  545. /* > \verbatim */
  546. /* > TRANSR is CHARACTER*1 */
  547. /* > Specifies whether the RFP format of A is normal or */
  548. /* > transposed format. */
  549. /* > = 'N': RFP format is Normal; */
  550. /* > = 'T': RFP format is Transpose. */
  551. /* > \endverbatim */
  552. /* > */
  553. /* > \param[in] UPLO */
  554. /* > \verbatim */
  555. /* > UPLO is CHARACTER*1 */
  556. /* > On entry, UPLO specifies whether the RFP matrix A came from */
  557. /* > an upper or lower triangular matrix as follows: */
  558. /* > = 'U': RFP A came from an upper triangular matrix; */
  559. /* > = 'L': RFP A came from a lower triangular matrix. */
  560. /* > \endverbatim */
  561. /* > */
  562. /* > \param[in] N */
  563. /* > \verbatim */
  564. /* > N is INTEGER */
  565. /* > The order of the matrix A. N >= 0. When N = 0, SLANSF is */
  566. /* > set to zero. */
  567. /* > \endverbatim */
  568. /* > */
  569. /* > \param[in] A */
  570. /* > \verbatim */
  571. /* > A is REAL array, dimension ( N*(N+1)/2 ); */
  572. /* > On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L') */
  573. /* > part of the symmetric matrix A stored in RFP format. See the */
  574. /* > "Notes" below for more details. */
  575. /* > Unchanged on exit. */
  576. /* > \endverbatim */
  577. /* > */
  578. /* > \param[out] WORK */
  579. /* > \verbatim */
  580. /* > WORK is REAL array, dimension (MAX(1,LWORK)), */
  581. /* > where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, */
  582. /* > WORK is not referenced. */
  583. /* > \endverbatim */
  584. /* Authors: */
  585. /* ======== */
  586. /* > \author Univ. of Tennessee */
  587. /* > \author Univ. of California Berkeley */
  588. /* > \author Univ. of Colorado Denver */
  589. /* > \author NAG Ltd. */
  590. /* > \date December 2016 */
  591. /* > \ingroup realOTHERcomputational */
  592. /* > \par Further Details: */
  593. /* ===================== */
  594. /* > */
  595. /* > \verbatim */
  596. /* > */
  597. /* > We first consider Rectangular Full Packed (RFP) Format when N is */
  598. /* > even. We give an example where N = 6. */
  599. /* > */
  600. /* > AP is Upper AP is Lower */
  601. /* > */
  602. /* > 00 01 02 03 04 05 00 */
  603. /* > 11 12 13 14 15 10 11 */
  604. /* > 22 23 24 25 20 21 22 */
  605. /* > 33 34 35 30 31 32 33 */
  606. /* > 44 45 40 41 42 43 44 */
  607. /* > 55 50 51 52 53 54 55 */
  608. /* > */
  609. /* > */
  610. /* > Let TRANSR = 'N'. RFP holds AP as follows: */
  611. /* > For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last */
  612. /* > three columns of AP upper. The lower triangle A(4:6,0:2) consists of */
  613. /* > the transpose of the first three columns of AP upper. */
  614. /* > For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first */
  615. /* > three columns of AP lower. The upper triangle A(0:2,0:2) consists of */
  616. /* > the transpose of the last three columns of AP lower. */
  617. /* > This covers the case N even and TRANSR = 'N'. */
  618. /* > */
  619. /* > RFP A RFP A */
  620. /* > */
  621. /* > 03 04 05 33 43 53 */
  622. /* > 13 14 15 00 44 54 */
  623. /* > 23 24 25 10 11 55 */
  624. /* > 33 34 35 20 21 22 */
  625. /* > 00 44 45 30 31 32 */
  626. /* > 01 11 55 40 41 42 */
  627. /* > 02 12 22 50 51 52 */
  628. /* > */
  629. /* > Now let TRANSR = 'T'. RFP A in both UPLO cases is just the */
  630. /* > transpose of RFP A above. One therefore gets: */
  631. /* > */
  632. /* > */
  633. /* > RFP A RFP A */
  634. /* > */
  635. /* > 03 13 23 33 00 01 02 33 00 10 20 30 40 50 */
  636. /* > 04 14 24 34 44 11 12 43 44 11 21 31 41 51 */
  637. /* > 05 15 25 35 45 55 22 53 54 55 22 32 42 52 */
  638. /* > */
  639. /* > */
  640. /* > We then consider Rectangular Full Packed (RFP) Format when N is */
  641. /* > odd. We give an example where N = 5. */
  642. /* > */
  643. /* > AP is Upper AP is Lower */
  644. /* > */
  645. /* > 00 01 02 03 04 00 */
  646. /* > 11 12 13 14 10 11 */
  647. /* > 22 23 24 20 21 22 */
  648. /* > 33 34 30 31 32 33 */
  649. /* > 44 40 41 42 43 44 */
  650. /* > */
  651. /* > */
  652. /* > Let TRANSR = 'N'. RFP holds AP as follows: */
  653. /* > For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last */
  654. /* > three columns of AP upper. The lower triangle A(3:4,0:1) consists of */
  655. /* > the transpose of the first two columns of AP upper. */
  656. /* > For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first */
  657. /* > three columns of AP lower. The upper triangle A(0:1,1:2) consists of */
  658. /* > the transpose of the last two columns of AP lower. */
  659. /* > This covers the case N odd and TRANSR = 'N'. */
  660. /* > */
  661. /* > RFP A RFP A */
  662. /* > */
  663. /* > 02 03 04 00 33 43 */
  664. /* > 12 13 14 10 11 44 */
  665. /* > 22 23 24 20 21 22 */
  666. /* > 00 33 34 30 31 32 */
  667. /* > 01 11 44 40 41 42 */
  668. /* > */
  669. /* > Now let TRANSR = 'T'. RFP A in both UPLO cases is just the */
  670. /* > transpose of RFP A above. One therefore gets: */
  671. /* > */
  672. /* > RFP A RFP A */
  673. /* > */
  674. /* > 02 12 22 00 01 00 10 20 30 40 50 */
  675. /* > 03 13 23 33 11 33 11 21 31 41 51 */
  676. /* > 04 14 24 34 44 43 44 22 32 42 52 */
  677. /* > \endverbatim */
  678. /* ===================================================================== */
  679. real slansf_(char *norm, char *transr, char *uplo, integer *n, real *a, real *
  680. work)
  681. {
  682. /* System generated locals */
  683. integer i__1, i__2;
  684. real ret_val, r__1;
  685. /* Local variables */
  686. real temp;
  687. integer i__, j, k, l;
  688. real s, scale;
  689. extern logical lsame_(char *, char *);
  690. real value;
  691. integer n1;
  692. real aa;
  693. extern logical sisnan_(real *);
  694. extern /* Subroutine */ void slassq_(integer *, real *, integer *, real *,
  695. real *);
  696. integer lda, ifm, noe, ilu;
  697. /* -- LAPACK computational routine (version 3.7.0) -- */
  698. /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
  699. /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
  700. /* December 2016 */
  701. /* ===================================================================== */
  702. if (*n == 0) {
  703. ret_val = 0.f;
  704. return ret_val;
  705. } else if (*n == 1) {
  706. ret_val = abs(a[0]);
  707. return ret_val;
  708. }
  709. /* set noe = 1 if n is odd. if n is even set noe=0 */
  710. noe = 1;
  711. if (*n % 2 == 0) {
  712. noe = 0;
  713. }
  714. /* set ifm = 0 when form='T or 't' and 1 otherwise */
  715. ifm = 1;
  716. if (lsame_(transr, "T")) {
  717. ifm = 0;
  718. }
  719. /* set ilu = 0 when uplo='U or 'u' and 1 otherwise */
  720. ilu = 1;
  721. if (lsame_(uplo, "U")) {
  722. ilu = 0;
  723. }
  724. /* set lda = (n+1)/2 when ifm = 0 */
  725. /* set lda = n when ifm = 1 and noe = 1 */
  726. /* set lda = n+1 when ifm = 1 and noe = 0 */
  727. if (ifm == 1) {
  728. if (noe == 1) {
  729. lda = *n;
  730. } else {
  731. /* noe=0 */
  732. lda = *n + 1;
  733. }
  734. } else {
  735. /* ifm=0 */
  736. lda = (*n + 1) / 2;
  737. }
  738. if (lsame_(norm, "M")) {
  739. /* Find f2cmax(abs(A(i,j))). */
  740. k = (*n + 1) / 2;
  741. value = 0.f;
  742. if (noe == 1) {
  743. /* n is odd */
  744. if (ifm == 1) {
  745. /* A is n by k */
  746. i__1 = k - 1;
  747. for (j = 0; j <= i__1; ++j) {
  748. i__2 = *n - 1;
  749. for (i__ = 0; i__ <= i__2; ++i__) {
  750. temp = (r__1 = a[i__ + j * lda], abs(r__1));
  751. if (value < temp || sisnan_(&temp)) {
  752. value = temp;
  753. }
  754. }
  755. }
  756. } else {
  757. /* xpose case; A is k by n */
  758. i__1 = *n - 1;
  759. for (j = 0; j <= i__1; ++j) {
  760. i__2 = k - 1;
  761. for (i__ = 0; i__ <= i__2; ++i__) {
  762. temp = (r__1 = a[i__ + j * lda], abs(r__1));
  763. if (value < temp || sisnan_(&temp)) {
  764. value = temp;
  765. }
  766. }
  767. }
  768. }
  769. } else {
  770. /* n is even */
  771. if (ifm == 1) {
  772. /* A is n+1 by k */
  773. i__1 = k - 1;
  774. for (j = 0; j <= i__1; ++j) {
  775. i__2 = *n;
  776. for (i__ = 0; i__ <= i__2; ++i__) {
  777. temp = (r__1 = a[i__ + j * lda], abs(r__1));
  778. if (value < temp || sisnan_(&temp)) {
  779. value = temp;
  780. }
  781. }
  782. }
  783. } else {
  784. /* xpose case; A is k by n+1 */
  785. i__1 = *n;
  786. for (j = 0; j <= i__1; ++j) {
  787. i__2 = k - 1;
  788. for (i__ = 0; i__ <= i__2; ++i__) {
  789. temp = (r__1 = a[i__ + j * lda], abs(r__1));
  790. if (value < temp || sisnan_(&temp)) {
  791. value = temp;
  792. }
  793. }
  794. }
  795. }
  796. }
  797. } else if (lsame_(norm, "I") || lsame_(norm, "O") || *(unsigned char *)norm == '1') {
  798. /* Find normI(A) ( = norm1(A), since A is symmetric). */
  799. if (ifm == 1) {
  800. k = *n / 2;
  801. if (noe == 1) {
  802. /* n is odd */
  803. if (ilu == 0) {
  804. i__1 = k - 1;
  805. for (i__ = 0; i__ <= i__1; ++i__) {
  806. work[i__] = 0.f;
  807. }
  808. i__1 = k;
  809. for (j = 0; j <= i__1; ++j) {
  810. s = 0.f;
  811. i__2 = k + j - 1;
  812. for (i__ = 0; i__ <= i__2; ++i__) {
  813. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  814. /* -> A(i,j+k) */
  815. s += aa;
  816. work[i__] += aa;
  817. }
  818. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  819. /* -> A(j+k,j+k) */
  820. work[j + k] = s + aa;
  821. if (i__ == k + k) {
  822. goto L10;
  823. }
  824. ++i__;
  825. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  826. /* -> A(j,j) */
  827. work[j] += aa;
  828. s = 0.f;
  829. i__2 = k - 1;
  830. for (l = j + 1; l <= i__2; ++l) {
  831. ++i__;
  832. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  833. /* -> A(l,j) */
  834. s += aa;
  835. work[l] += aa;
  836. }
  837. work[j] += s;
  838. }
  839. L10:
  840. value = work[0];
  841. i__1 = *n - 1;
  842. for (i__ = 1; i__ <= i__1; ++i__) {
  843. temp = work[i__];
  844. if (value < temp || sisnan_(&temp)) {
  845. value = temp;
  846. }
  847. }
  848. } else {
  849. /* ilu = 1 */
  850. ++k;
  851. /* k=(n+1)/2 for n odd and ilu=1 */
  852. i__1 = *n - 1;
  853. for (i__ = k; i__ <= i__1; ++i__) {
  854. work[i__] = 0.f;
  855. }
  856. for (j = k - 1; j >= 0; --j) {
  857. s = 0.f;
  858. i__1 = j - 2;
  859. for (i__ = 0; i__ <= i__1; ++i__) {
  860. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  861. /* -> A(j+k,i+k) */
  862. s += aa;
  863. work[i__ + k] += aa;
  864. }
  865. if (j > 0) {
  866. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  867. /* -> A(j+k,j+k) */
  868. s += aa;
  869. work[i__ + k] += s;
  870. /* i=j */
  871. ++i__;
  872. }
  873. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  874. /* -> A(j,j) */
  875. work[j] = aa;
  876. s = 0.f;
  877. i__1 = *n - 1;
  878. for (l = j + 1; l <= i__1; ++l) {
  879. ++i__;
  880. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  881. /* -> A(l,j) */
  882. s += aa;
  883. work[l] += aa;
  884. }
  885. work[j] += s;
  886. }
  887. value = work[0];
  888. i__1 = *n - 1;
  889. for (i__ = 1; i__ <= i__1; ++i__) {
  890. temp = work[i__];
  891. if (value < temp || sisnan_(&temp)) {
  892. value = temp;
  893. }
  894. }
  895. }
  896. } else {
  897. /* n is even */
  898. if (ilu == 0) {
  899. i__1 = k - 1;
  900. for (i__ = 0; i__ <= i__1; ++i__) {
  901. work[i__] = 0.f;
  902. }
  903. i__1 = k - 1;
  904. for (j = 0; j <= i__1; ++j) {
  905. s = 0.f;
  906. i__2 = k + j - 1;
  907. for (i__ = 0; i__ <= i__2; ++i__) {
  908. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  909. /* -> A(i,j+k) */
  910. s += aa;
  911. work[i__] += aa;
  912. }
  913. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  914. /* -> A(j+k,j+k) */
  915. work[j + k] = s + aa;
  916. ++i__;
  917. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  918. /* -> A(j,j) */
  919. work[j] += aa;
  920. s = 0.f;
  921. i__2 = k - 1;
  922. for (l = j + 1; l <= i__2; ++l) {
  923. ++i__;
  924. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  925. /* -> A(l,j) */
  926. s += aa;
  927. work[l] += aa;
  928. }
  929. work[j] += s;
  930. }
  931. value = work[0];
  932. i__1 = *n - 1;
  933. for (i__ = 1; i__ <= i__1; ++i__) {
  934. temp = work[i__];
  935. if (value < temp || sisnan_(&temp)) {
  936. value = temp;
  937. }
  938. }
  939. } else {
  940. /* ilu = 1 */
  941. i__1 = *n - 1;
  942. for (i__ = k; i__ <= i__1; ++i__) {
  943. work[i__] = 0.f;
  944. }
  945. for (j = k - 1; j >= 0; --j) {
  946. s = 0.f;
  947. i__1 = j - 1;
  948. for (i__ = 0; i__ <= i__1; ++i__) {
  949. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  950. /* -> A(j+k,i+k) */
  951. s += aa;
  952. work[i__ + k] += aa;
  953. }
  954. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  955. /* -> A(j+k,j+k) */
  956. s += aa;
  957. work[i__ + k] += s;
  958. /* i=j */
  959. ++i__;
  960. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  961. /* -> A(j,j) */
  962. work[j] = aa;
  963. s = 0.f;
  964. i__1 = *n - 1;
  965. for (l = j + 1; l <= i__1; ++l) {
  966. ++i__;
  967. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  968. /* -> A(l,j) */
  969. s += aa;
  970. work[l] += aa;
  971. }
  972. work[j] += s;
  973. }
  974. value = work[0];
  975. i__1 = *n - 1;
  976. for (i__ = 1; i__ <= i__1; ++i__) {
  977. temp = work[i__];
  978. if (value < temp || sisnan_(&temp)) {
  979. value = temp;
  980. }
  981. }
  982. }
  983. }
  984. } else {
  985. /* ifm=0 */
  986. k = *n / 2;
  987. if (noe == 1) {
  988. /* n is odd */
  989. if (ilu == 0) {
  990. n1 = k;
  991. /* n/2 */
  992. ++k;
  993. /* k is the row size and lda */
  994. i__1 = *n - 1;
  995. for (i__ = n1; i__ <= i__1; ++i__) {
  996. work[i__] = 0.f;
  997. }
  998. i__1 = n1 - 1;
  999. for (j = 0; j <= i__1; ++j) {
  1000. s = 0.f;
  1001. i__2 = k - 1;
  1002. for (i__ = 0; i__ <= i__2; ++i__) {
  1003. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1004. /* A(j,n1+i) */
  1005. work[i__ + n1] += aa;
  1006. s += aa;
  1007. }
  1008. work[j] = s;
  1009. }
  1010. /* j=n1=k-1 is special */
  1011. s = (r__1 = a[j * lda], abs(r__1));
  1012. /* A(k-1,k-1) */
  1013. i__1 = k - 1;
  1014. for (i__ = 1; i__ <= i__1; ++i__) {
  1015. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1016. /* A(k-1,i+n1) */
  1017. work[i__ + n1] += aa;
  1018. s += aa;
  1019. }
  1020. work[j] += s;
  1021. i__1 = *n - 1;
  1022. for (j = k; j <= i__1; ++j) {
  1023. s = 0.f;
  1024. i__2 = j - k - 1;
  1025. for (i__ = 0; i__ <= i__2; ++i__) {
  1026. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1027. /* A(i,j-k) */
  1028. work[i__] += aa;
  1029. s += aa;
  1030. }
  1031. /* i=j-k */
  1032. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1033. /* A(j-k,j-k) */
  1034. s += aa;
  1035. work[j - k] += s;
  1036. ++i__;
  1037. s = (r__1 = a[i__ + j * lda], abs(r__1));
  1038. /* A(j,j) */
  1039. i__2 = *n - 1;
  1040. for (l = j + 1; l <= i__2; ++l) {
  1041. ++i__;
  1042. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1043. /* A(j,l) */
  1044. work[l] += aa;
  1045. s += aa;
  1046. }
  1047. work[j] += s;
  1048. }
  1049. value = work[0];
  1050. i__1 = *n - 1;
  1051. for (i__ = 1; i__ <= i__1; ++i__) {
  1052. temp = work[i__];
  1053. if (value < temp || sisnan_(&temp)) {
  1054. value = temp;
  1055. }
  1056. }
  1057. } else {
  1058. /* ilu=1 */
  1059. ++k;
  1060. /* k=(n+1)/2 for n odd and ilu=1 */
  1061. i__1 = *n - 1;
  1062. for (i__ = k; i__ <= i__1; ++i__) {
  1063. work[i__] = 0.f;
  1064. }
  1065. i__1 = k - 2;
  1066. for (j = 0; j <= i__1; ++j) {
  1067. /* process */
  1068. s = 0.f;
  1069. i__2 = j - 1;
  1070. for (i__ = 0; i__ <= i__2; ++i__) {
  1071. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1072. /* A(j,i) */
  1073. work[i__] += aa;
  1074. s += aa;
  1075. }
  1076. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1077. /* i=j so process of A(j,j) */
  1078. s += aa;
  1079. work[j] = s;
  1080. /* is initialised here */
  1081. ++i__;
  1082. /* i=j process A(j+k,j+k) */
  1083. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1084. s = aa;
  1085. i__2 = *n - 1;
  1086. for (l = k + j + 1; l <= i__2; ++l) {
  1087. ++i__;
  1088. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1089. /* A(l,k+j) */
  1090. s += aa;
  1091. work[l] += aa;
  1092. }
  1093. work[k + j] += s;
  1094. }
  1095. /* j=k-1 is special :process col A(k-1,0:k-1) */
  1096. s = 0.f;
  1097. i__1 = k - 2;
  1098. for (i__ = 0; i__ <= i__1; ++i__) {
  1099. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1100. /* A(k,i) */
  1101. work[i__] += aa;
  1102. s += aa;
  1103. }
  1104. /* i=k-1 */
  1105. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1106. /* A(k-1,k-1) */
  1107. s += aa;
  1108. work[i__] = s;
  1109. /* done with col j=k+1 */
  1110. i__1 = *n - 1;
  1111. for (j = k; j <= i__1; ++j) {
  1112. /* process col j of A = A(j,0:k-1) */
  1113. s = 0.f;
  1114. i__2 = k - 1;
  1115. for (i__ = 0; i__ <= i__2; ++i__) {
  1116. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1117. /* A(j,i) */
  1118. work[i__] += aa;
  1119. s += aa;
  1120. }
  1121. work[j] += s;
  1122. }
  1123. value = work[0];
  1124. i__1 = *n - 1;
  1125. for (i__ = 1; i__ <= i__1; ++i__) {
  1126. temp = work[i__];
  1127. if (value < temp || sisnan_(&temp)) {
  1128. value = temp;
  1129. }
  1130. }
  1131. }
  1132. } else {
  1133. /* n is even */
  1134. if (ilu == 0) {
  1135. i__1 = *n - 1;
  1136. for (i__ = k; i__ <= i__1; ++i__) {
  1137. work[i__] = 0.f;
  1138. }
  1139. i__1 = k - 1;
  1140. for (j = 0; j <= i__1; ++j) {
  1141. s = 0.f;
  1142. i__2 = k - 1;
  1143. for (i__ = 0; i__ <= i__2; ++i__) {
  1144. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1145. /* A(j,i+k) */
  1146. work[i__ + k] += aa;
  1147. s += aa;
  1148. }
  1149. work[j] = s;
  1150. }
  1151. /* j=k */
  1152. aa = (r__1 = a[j * lda], abs(r__1));
  1153. /* A(k,k) */
  1154. s = aa;
  1155. i__1 = k - 1;
  1156. for (i__ = 1; i__ <= i__1; ++i__) {
  1157. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1158. /* A(k,k+i) */
  1159. work[i__ + k] += aa;
  1160. s += aa;
  1161. }
  1162. work[j] += s;
  1163. i__1 = *n - 1;
  1164. for (j = k + 1; j <= i__1; ++j) {
  1165. s = 0.f;
  1166. i__2 = j - 2 - k;
  1167. for (i__ = 0; i__ <= i__2; ++i__) {
  1168. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1169. /* A(i,j-k-1) */
  1170. work[i__] += aa;
  1171. s += aa;
  1172. }
  1173. /* i=j-1-k */
  1174. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1175. /* A(j-k-1,j-k-1) */
  1176. s += aa;
  1177. work[j - k - 1] += s;
  1178. ++i__;
  1179. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1180. /* A(j,j) */
  1181. s = aa;
  1182. i__2 = *n - 1;
  1183. for (l = j + 1; l <= i__2; ++l) {
  1184. ++i__;
  1185. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1186. /* A(j,l) */
  1187. work[l] += aa;
  1188. s += aa;
  1189. }
  1190. work[j] += s;
  1191. }
  1192. /* j=n */
  1193. s = 0.f;
  1194. i__1 = k - 2;
  1195. for (i__ = 0; i__ <= i__1; ++i__) {
  1196. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1197. /* A(i,k-1) */
  1198. work[i__] += aa;
  1199. s += aa;
  1200. }
  1201. /* i=k-1 */
  1202. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1203. /* A(k-1,k-1) */
  1204. s += aa;
  1205. work[i__] += s;
  1206. value = work[0];
  1207. i__1 = *n - 1;
  1208. for (i__ = 1; i__ <= i__1; ++i__) {
  1209. temp = work[i__];
  1210. if (value < temp || sisnan_(&temp)) {
  1211. value = temp;
  1212. }
  1213. }
  1214. } else {
  1215. /* ilu=1 */
  1216. i__1 = *n - 1;
  1217. for (i__ = k; i__ <= i__1; ++i__) {
  1218. work[i__] = 0.f;
  1219. }
  1220. /* j=0 is special :process col A(k:n-1,k) */
  1221. s = abs(a[0]);
  1222. /* A(k,k) */
  1223. i__1 = k - 1;
  1224. for (i__ = 1; i__ <= i__1; ++i__) {
  1225. aa = (r__1 = a[i__], abs(r__1));
  1226. /* A(k+i,k) */
  1227. work[i__ + k] += aa;
  1228. s += aa;
  1229. }
  1230. work[k] += s;
  1231. i__1 = k - 1;
  1232. for (j = 1; j <= i__1; ++j) {
  1233. /* process */
  1234. s = 0.f;
  1235. i__2 = j - 2;
  1236. for (i__ = 0; i__ <= i__2; ++i__) {
  1237. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1238. /* A(j-1,i) */
  1239. work[i__] += aa;
  1240. s += aa;
  1241. }
  1242. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1243. /* i=j-1 so process of A(j-1,j-1) */
  1244. s += aa;
  1245. work[j - 1] = s;
  1246. /* is initialised here */
  1247. ++i__;
  1248. /* i=j process A(j+k,j+k) */
  1249. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1250. s = aa;
  1251. i__2 = *n - 1;
  1252. for (l = k + j + 1; l <= i__2; ++l) {
  1253. ++i__;
  1254. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1255. /* A(l,k+j) */
  1256. s += aa;
  1257. work[l] += aa;
  1258. }
  1259. work[k + j] += s;
  1260. }
  1261. /* j=k is special :process col A(k,0:k-1) */
  1262. s = 0.f;
  1263. i__1 = k - 2;
  1264. for (i__ = 0; i__ <= i__1; ++i__) {
  1265. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1266. /* A(k,i) */
  1267. work[i__] += aa;
  1268. s += aa;
  1269. }
  1270. /* i=k-1 */
  1271. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1272. /* A(k-1,k-1) */
  1273. s += aa;
  1274. work[i__] = s;
  1275. /* done with col j=k+1 */
  1276. i__1 = *n;
  1277. for (j = k + 1; j <= i__1; ++j) {
  1278. /* process col j-1 of A = A(j-1,0:k-1) */
  1279. s = 0.f;
  1280. i__2 = k - 1;
  1281. for (i__ = 0; i__ <= i__2; ++i__) {
  1282. aa = (r__1 = a[i__ + j * lda], abs(r__1));
  1283. /* A(j-1,i) */
  1284. work[i__] += aa;
  1285. s += aa;
  1286. }
  1287. work[j - 1] += s;
  1288. }
  1289. value = work[0];
  1290. i__1 = *n - 1;
  1291. for (i__ = 1; i__ <= i__1; ++i__) {
  1292. temp = work[i__];
  1293. if (value < temp || sisnan_(&temp)) {
  1294. value = temp;
  1295. }
  1296. }
  1297. }
  1298. }
  1299. }
  1300. } else if (lsame_(norm, "F") || lsame_(norm, "E")) {
  1301. /* Find normF(A). */
  1302. k = (*n + 1) / 2;
  1303. scale = 0.f;
  1304. s = 1.f;
  1305. if (noe == 1) {
  1306. /* n is odd */
  1307. if (ifm == 1) {
  1308. /* A is normal */
  1309. if (ilu == 0) {
  1310. /* A is upper */
  1311. i__1 = k - 3;
  1312. for (j = 0; j <= i__1; ++j) {
  1313. i__2 = k - j - 2;
  1314. slassq_(&i__2, &a[k + j + 1 + j * lda], &c__1, &scale,
  1315. &s);
  1316. /* L at A(k,0) */
  1317. }
  1318. i__1 = k - 1;
  1319. for (j = 0; j <= i__1; ++j) {
  1320. i__2 = k + j - 1;
  1321. slassq_(&i__2, &a[j * lda], &c__1, &scale, &s);
  1322. /* trap U at A(0,0) */
  1323. }
  1324. s += s;
  1325. /* double s for the off diagonal elements */
  1326. i__1 = k - 1;
  1327. i__2 = lda + 1;
  1328. slassq_(&i__1, &a[k], &i__2, &scale, &s);
  1329. /* tri L at A(k,0) */
  1330. i__1 = lda + 1;
  1331. slassq_(&k, &a[k - 1], &i__1, &scale, &s);
  1332. /* tri U at A(k-1,0) */
  1333. } else {
  1334. /* ilu=1 & A is lower */
  1335. i__1 = k - 1;
  1336. for (j = 0; j <= i__1; ++j) {
  1337. i__2 = *n - j - 1;
  1338. slassq_(&i__2, &a[j + 1 + j * lda], &c__1, &scale, &s)
  1339. ;
  1340. /* trap L at A(0,0) */
  1341. }
  1342. i__1 = k - 2;
  1343. for (j = 0; j <= i__1; ++j) {
  1344. slassq_(&j, &a[(j + 1) * lda], &c__1, &scale, &s);
  1345. /* U at A(0,1) */
  1346. }
  1347. s += s;
  1348. /* double s for the off diagonal elements */
  1349. i__1 = lda + 1;
  1350. slassq_(&k, a, &i__1, &scale, &s);
  1351. /* tri L at A(0,0) */
  1352. i__1 = k - 1;
  1353. i__2 = lda + 1;
  1354. slassq_(&i__1, &a[lda], &i__2, &scale, &s);
  1355. /* tri U at A(0,1) */
  1356. }
  1357. } else {
  1358. /* A is xpose */
  1359. if (ilu == 0) {
  1360. /* A**T is upper */
  1361. i__1 = k - 2;
  1362. for (j = 1; j <= i__1; ++j) {
  1363. slassq_(&j, &a[(k + j) * lda], &c__1, &scale, &s);
  1364. /* U at A(0,k) */
  1365. }
  1366. i__1 = k - 2;
  1367. for (j = 0; j <= i__1; ++j) {
  1368. slassq_(&k, &a[j * lda], &c__1, &scale, &s);
  1369. /* k by k-1 rect. at A(0,0) */
  1370. }
  1371. i__1 = k - 2;
  1372. for (j = 0; j <= i__1; ++j) {
  1373. i__2 = k - j - 1;
  1374. slassq_(&i__2, &a[j + 1 + (j + k - 1) * lda], &c__1, &
  1375. scale, &s);
  1376. /* L at A(0,k-1) */
  1377. }
  1378. s += s;
  1379. /* double s for the off diagonal elements */
  1380. i__1 = k - 1;
  1381. i__2 = lda + 1;
  1382. slassq_(&i__1, &a[k * lda], &i__2, &scale, &s);
  1383. /* tri U at A(0,k) */
  1384. i__1 = lda + 1;
  1385. slassq_(&k, &a[(k - 1) * lda], &i__1, &scale, &s);
  1386. /* tri L at A(0,k-1) */
  1387. } else {
  1388. /* A**T is lower */
  1389. i__1 = k - 1;
  1390. for (j = 1; j <= i__1; ++j) {
  1391. slassq_(&j, &a[j * lda], &c__1, &scale, &s);
  1392. /* U at A(0,0) */
  1393. }
  1394. i__1 = *n - 1;
  1395. for (j = k; j <= i__1; ++j) {
  1396. slassq_(&k, &a[j * lda], &c__1, &scale, &s);
  1397. /* k by k-1 rect. at A(0,k) */
  1398. }
  1399. i__1 = k - 3;
  1400. for (j = 0; j <= i__1; ++j) {
  1401. i__2 = k - j - 2;
  1402. slassq_(&i__2, &a[j + 2 + j * lda], &c__1, &scale, &s)
  1403. ;
  1404. /* L at A(1,0) */
  1405. }
  1406. s += s;
  1407. /* double s for the off diagonal elements */
  1408. i__1 = lda + 1;
  1409. slassq_(&k, a, &i__1, &scale, &s);
  1410. /* tri U at A(0,0) */
  1411. i__1 = k - 1;
  1412. i__2 = lda + 1;
  1413. slassq_(&i__1, &a[1], &i__2, &scale, &s);
  1414. /* tri L at A(1,0) */
  1415. }
  1416. }
  1417. } else {
  1418. /* n is even */
  1419. if (ifm == 1) {
  1420. /* A is normal */
  1421. if (ilu == 0) {
  1422. /* A is upper */
  1423. i__1 = k - 2;
  1424. for (j = 0; j <= i__1; ++j) {
  1425. i__2 = k - j - 1;
  1426. slassq_(&i__2, &a[k + j + 2 + j * lda], &c__1, &scale,
  1427. &s);
  1428. /* L at A(k+1,0) */
  1429. }
  1430. i__1 = k - 1;
  1431. for (j = 0; j <= i__1; ++j) {
  1432. i__2 = k + j;
  1433. slassq_(&i__2, &a[j * lda], &c__1, &scale, &s);
  1434. /* trap U at A(0,0) */
  1435. }
  1436. s += s;
  1437. /* double s for the off diagonal elements */
  1438. i__1 = lda + 1;
  1439. slassq_(&k, &a[k + 1], &i__1, &scale, &s);
  1440. /* tri L at A(k+1,0) */
  1441. i__1 = lda + 1;
  1442. slassq_(&k, &a[k], &i__1, &scale, &s);
  1443. /* tri U at A(k,0) */
  1444. } else {
  1445. /* ilu=1 & A is lower */
  1446. i__1 = k - 1;
  1447. for (j = 0; j <= i__1; ++j) {
  1448. i__2 = *n - j - 1;
  1449. slassq_(&i__2, &a[j + 2 + j * lda], &c__1, &scale, &s)
  1450. ;
  1451. /* trap L at A(1,0) */
  1452. }
  1453. i__1 = k - 1;
  1454. for (j = 1; j <= i__1; ++j) {
  1455. slassq_(&j, &a[j * lda], &c__1, &scale, &s);
  1456. /* U at A(0,0) */
  1457. }
  1458. s += s;
  1459. /* double s for the off diagonal elements */
  1460. i__1 = lda + 1;
  1461. slassq_(&k, &a[1], &i__1, &scale, &s);
  1462. /* tri L at A(1,0) */
  1463. i__1 = lda + 1;
  1464. slassq_(&k, a, &i__1, &scale, &s);
  1465. /* tri U at A(0,0) */
  1466. }
  1467. } else {
  1468. /* A is xpose */
  1469. if (ilu == 0) {
  1470. /* A**T is upper */
  1471. i__1 = k - 1;
  1472. for (j = 1; j <= i__1; ++j) {
  1473. slassq_(&j, &a[(k + 1 + j) * lda], &c__1, &scale, &s);
  1474. /* U at A(0,k+1) */
  1475. }
  1476. i__1 = k - 1;
  1477. for (j = 0; j <= i__1; ++j) {
  1478. slassq_(&k, &a[j * lda], &c__1, &scale, &s);
  1479. /* k by k rect. at A(0,0) */
  1480. }
  1481. i__1 = k - 2;
  1482. for (j = 0; j <= i__1; ++j) {
  1483. i__2 = k - j - 1;
  1484. slassq_(&i__2, &a[j + 1 + (j + k) * lda], &c__1, &
  1485. scale, &s);
  1486. /* L at A(0,k) */
  1487. }
  1488. s += s;
  1489. /* double s for the off diagonal elements */
  1490. i__1 = lda + 1;
  1491. slassq_(&k, &a[(k + 1) * lda], &i__1, &scale, &s);
  1492. /* tri U at A(0,k+1) */
  1493. i__1 = lda + 1;
  1494. slassq_(&k, &a[k * lda], &i__1, &scale, &s);
  1495. /* tri L at A(0,k) */
  1496. } else {
  1497. /* A**T is lower */
  1498. i__1 = k - 1;
  1499. for (j = 1; j <= i__1; ++j) {
  1500. slassq_(&j, &a[(j + 1) * lda], &c__1, &scale, &s);
  1501. /* U at A(0,1) */
  1502. }
  1503. i__1 = *n;
  1504. for (j = k + 1; j <= i__1; ++j) {
  1505. slassq_(&k, &a[j * lda], &c__1, &scale, &s);
  1506. /* k by k rect. at A(0,k+1) */
  1507. }
  1508. i__1 = k - 2;
  1509. for (j = 0; j <= i__1; ++j) {
  1510. i__2 = k - j - 1;
  1511. slassq_(&i__2, &a[j + 1 + j * lda], &c__1, &scale, &s)
  1512. ;
  1513. /* L at A(0,0) */
  1514. }
  1515. s += s;
  1516. /* double s for the off diagonal elements */
  1517. i__1 = lda + 1;
  1518. slassq_(&k, &a[lda], &i__1, &scale, &s);
  1519. /* tri L at A(0,1) */
  1520. i__1 = lda + 1;
  1521. slassq_(&k, a, &i__1, &scale, &s);
  1522. /* tri U at A(0,0) */
  1523. }
  1524. }
  1525. }
  1526. value = scale * sqrt(s);
  1527. }
  1528. ret_val = value;
  1529. return ret_val;
  1530. /* End of SLANSF */
  1531. } /* slansf_ */