@@ -252,271 +252,6 @@ typedef struct Namelist Namelist;
252252//#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);}
253253#define mymaxloc_ (w ,s ,e ,n ) dmaxloc_(w,*(s),*(e),n)
254254
255- /* procedure parameter types for -A and -C++ */
256-
257- #define F2C_proc_par_types 1
258- #ifdef __cplusplus
259- typedef logical (* L_fp )(...);
260- #else
261- typedef logical (* L_fp )();
262- #endif
263- #if 0
264- static float spow_ui (float x , integer n ) {
265- float pow = 1.0 ; unsigned long int u ;
266- if (n != 0 ) {
267- if (n < 0 ) n = - n , x = 1 /x ;
268- for (u = n ; ; ) {
269- if (u & 01 ) pow *= x ;
270- if (u >>= 1 ) x *= x ;
271- else break ;
272- }
273- }
274- return pow ;
275- }
276- static double dpow_ui (double x , integer n ) {
277- double pow = 1.0 ; unsigned long int u ;
278- if (n != 0 ) {
279- if (n < 0 ) n = - n , x = 1 /x ;
280- for (u = n ; ; ) {
281- if (u & 01 ) pow *= x ;
282- if (u >>= 1 ) x *= x ;
283- else break ;
284- }
285- }
286- return pow ;
287- }
288- #ifdef _MSC_VER
289- static _Fcomplex cpow_ui (complex x , integer n ) {
290- complex pow = {1.0 ,0.0 }; unsigned long int u ;
291- if (n != 0 ) {
292- if (n < 0 ) n = - n , x .r = 1 /x .r , x .i = 1 /x .i ;
293- for (u = n ; ; ) {
294- if (u & 01 ) pow .r *= x .r , pow .i *= x .i ;
295- if (u >>= 1 ) x .r *= x .r , x .i *= x .i ;
296- else break ;
297- }
298- }
299- _Fcomplex p = {pow .r , pow .i };
300- return p ;
301- }
302- #else
303- static _Complex float cpow_ui (_Complex float x , integer n ) {
304- _Complex float pow = 1.0 ; unsigned long int u ;
305- if (n != 0 ) {
306- if (n < 0 ) n = - n , x = 1 /x ;
307- for (u = n ; ; ) {
308- if (u & 01 ) pow *= x ;
309- if (u >>= 1 ) x *= x ;
310- else break ;
311- }
312- }
313- return pow ;
314- }
315- #endif
316- #ifdef _MSC_VER
317- static _Dcomplex zpow_ui (_Dcomplex x , integer n ) {
318- _Dcomplex pow = {1.0 ,0.0 }; unsigned long int u ;
319- if (n != 0 ) {
320- if (n < 0 ) n = - n , x ._Val [0 ] = 1 /x ._Val [0 ], x ._Val [1 ] = 1 /x ._Val [1 ];
321- for (u = n ; ; ) {
322- if (u & 01 ) pow ._Val [0 ] *= x ._Val [0 ], pow ._Val [1 ] *= x ._Val [1 ];
323- if (u >>= 1 ) x ._Val [0 ] *= x ._Val [0 ], x ._Val [1 ] *= x ._Val [1 ];
324- else break ;
325- }
326- }
327- _Dcomplex p = {pow ._Val [0 ], pow ._Val [1 ]};
328- return p ;
329- }
330- #else
331- static _Complex double zpow_ui (_Complex double x , integer n ) {
332- _Complex double pow = 1.0 ; unsigned long int u ;
333- if (n != 0 ) {
334- if (n < 0 ) n = - n , x = 1 /x ;
335- for (u = n ; ; ) {
336- if (u & 01 ) pow *= x ;
337- if (u >>= 1 ) x *= x ;
338- else break ;
339- }
340- }
341- return pow ;
342- }
343- #endif
344- static integer pow_ii (integer x , integer n ) {
345- integer pow ; unsigned long int u ;
346- if (n <= 0 ) {
347- if (n == 0 || x == 1 ) pow = 1 ;
348- else if (x != -1 ) pow = x == 0 ? 1 /x : 0 ;
349- else n = - n ;
350- }
351- if ((n > 0 ) || !(n == 0 || x == 1 || x != -1 )) {
352- u = n ;
353- for (pow = 1 ; ; ) {
354- if (u & 01 ) pow *= x ;
355- if (u >>= 1 ) x *= x ;
356- else break ;
357- }
358- }
359- return pow ;
360- }
361- static integer dmaxloc_ (double * w , integer s , integer e , integer * n )
362- {
363- double m ; integer i , mi ;
364- for (m = w [s - 1 ], mi = s , i = s + 1 ; i <=e ; i ++ )
365- if (w [i - 1 ]> m ) mi = i ,m = w [i - 1 ];
366- return mi - s + 1 ;
367- }
368- static integer smaxloc_ (float * w , integer s , integer e , integer * n )
369- {
370- float m ; integer i , mi ;
371- for (m = w [s - 1 ], mi = s , i = s + 1 ; i <=e ; i ++ )
372- if (w [i - 1 ]> m ) mi = i ,m = w [i - 1 ];
373- return mi - s + 1 ;
374- }
375- static inline void cdotc_ (complex * z , integer * n_ , complex * x , integer * incx_ , complex * y , integer * incy_ ) {
376- integer n = * n_ , incx = * incx_ , incy = * incy_ , i ;
377- #ifdef _MSC_VER
378- _Fcomplex zdotc = {0.0 , 0.0 };
379- if (incx == 1 && incy == 1 ) {
380- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
381- zdotc ._Val [0 ] += conjf (Cf (& x [i ]))._Val [0 ] * Cf (& y [i ])._Val [0 ];
382- zdotc ._Val [1 ] += conjf (Cf (& x [i ]))._Val [1 ] * Cf (& y [i ])._Val [1 ];
383- }
384- } else {
385- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
386- zdotc ._Val [0 ] += conjf (Cf (& x [i * incx ]))._Val [0 ] * Cf (& y [i * incy ])._Val [0 ];
387- zdotc ._Val [1 ] += conjf (Cf (& x [i * incx ]))._Val [1 ] * Cf (& y [i * incy ])._Val [1 ];
388- }
389- }
390- pCf (z ) = zdotc ;
391- }
392- #else
393- _Complex float zdotc = 0.0 ;
394- if (incx == 1 && incy == 1 ) {
395- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
396- zdotc += conjf (Cf (& x [i ])) * Cf (& y [i ]);
397- }
398- } else {
399- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
400- zdotc += conjf (Cf (& x [i * incx ])) * Cf (& y [i * incy ]);
401- }
402- }
403- pCf (z ) = zdotc ;
404- }
405- #endif
406- static inline void zdotc_ (doublecomplex * z , integer * n_ , doublecomplex * x , integer * incx_ , doublecomplex * y , integer * incy_ ) {
407- integer n = * n_ , incx = * incx_ , incy = * incy_ , i ;
408- #ifdef _MSC_VER
409- _Dcomplex zdotc = {0.0 , 0.0 };
410- if (incx == 1 && incy == 1 ) {
411- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
412- zdotc ._Val [0 ] += conj (Cd (& x [i ]))._Val [0 ] * Cd (& y [i ])._Val [0 ];
413- zdotc ._Val [1 ] += conj (Cd (& x [i ]))._Val [1 ] * Cd (& y [i ])._Val [1 ];
414- }
415- } else {
416- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
417- zdotc ._Val [0 ] += conj (Cd (& x [i * incx ]))._Val [0 ] * Cd (& y [i * incy ])._Val [0 ];
418- zdotc ._Val [1 ] += conj (Cd (& x [i * incx ]))._Val [1 ] * Cd (& y [i * incy ])._Val [1 ];
419- }
420- }
421- pCd (z ) = zdotc ;
422- }
423- #else
424- _Complex double zdotc = 0.0 ;
425- if (incx == 1 && incy == 1 ) {
426- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
427- zdotc += conj (Cd (& x [i ])) * Cd (& y [i ]);
428- }
429- } else {
430- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
431- zdotc += conj (Cd (& x [i * incx ])) * Cd (& y [i * incy ]);
432- }
433- }
434- pCd (z ) = zdotc ;
435- }
436- #endif
437- static inline void cdotu_ (complex * z , integer * n_ , complex * x , integer * incx_ , complex * y , integer * incy_ ) {
438- integer n = * n_ , incx = * incx_ , incy = * incy_ , i ;
439- #ifdef _MSC_VER
440- _Fcomplex zdotc = {0.0 , 0.0 };
441- if (incx == 1 && incy == 1 ) {
442- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
443- zdotc ._Val [0 ] += Cf (& x [i ])._Val [0 ] * Cf (& y [i ])._Val [0 ];
444- zdotc ._Val [1 ] += Cf (& x [i ])._Val [1 ] * Cf (& y [i ])._Val [1 ];
445- }
446- } else {
447- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
448- zdotc ._Val [0 ] += Cf (& x [i * incx ])._Val [0 ] * Cf (& y [i * incy ])._Val [0 ];
449- zdotc ._Val [1 ] += Cf (& x [i * incx ])._Val [1 ] * Cf (& y [i * incy ])._Val [1 ];
450- }
451- }
452- pCf (z ) = zdotc ;
453- }
454- #else
455- _Complex float zdotc = 0.0 ;
456- if (incx == 1 && incy == 1 ) {
457- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
458- zdotc += Cf (& x [i ]) * Cf (& y [i ]);
459- }
460- } else {
461- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
462- zdotc += Cf (& x [i * incx ]) * Cf (& y [i * incy ]);
463- }
464- }
465- pCf (z ) = zdotc ;
466- }
467- #endif
468- static inline void zdotu_ (doublecomplex * z , integer * n_ , doublecomplex * x , integer * incx_ , doublecomplex * y , integer * incy_ ) {
469- integer n = * n_ , incx = * incx_ , incy = * incy_ , i ;
470- #ifdef _MSC_VER
471- _Dcomplex zdotc = {0.0 , 0.0 };
472- if (incx == 1 && incy == 1 ) {
473- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
474- zdotc ._Val [0 ] += Cd (& x [i ])._Val [0 ] * Cd (& y [i ])._Val [0 ];
475- zdotc ._Val [1 ] += Cd (& x [i ])._Val [1 ] * Cd (& y [i ])._Val [1 ];
476- }
477- } else {
478- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
479- zdotc ._Val [0 ] += Cd (& x [i * incx ])._Val [0 ] * Cd (& y [i * incy ])._Val [0 ];
480- zdotc ._Val [1 ] += Cd (& x [i * incx ])._Val [1 ] * Cd (& y [i * incy ])._Val [1 ];
481- }
482- }
483- pCd (z ) = zdotc ;
484- }
485- #else
486- _Complex double zdotc = 0.0 ;
487- if (incx == 1 && incy == 1 ) {
488- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
489- zdotc += Cd (& x [i ]) * Cd (& y [i ]);
490- }
491- } else {
492- for (i = 0 ;i < n ;i ++ ) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
493- zdotc += Cd (& x [i * incx ]) * Cd (& y [i * incy ]);
494- }
495- }
496- pCd (z ) = zdotc ;
497- }
498- #endif
499- #endif
500- /* -- translated by f2c (version 20000121).
501- You must link the resulting object file with the libraries:
502- -lf2c -lm (in that order)
503- */
504-
505-
506-
507- /* -- translated by f2c (version 20200916).
508- You must link the resulting object file with libf2c:
509- on Microsoft Windows system, link with libf2c.lib;
510- on Linux or Unix systems, link with .../path/to/libf2c.a -lm
511- or, if you install libf2c.a in a standard place, with -lf2c -lm
512- -- in that order, at the end of the command line, as in
513- cc *.o -lf2c -lm
514- Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
515-
516- http://www.netlib.org/f2c/libf2c.zip
517- */
518-
519-
520255
521256/* Table of constant values */
522257
@@ -575,7 +310,7 @@ static integer c__1 = 1;
575310 extern real slamch_ (char * );
576311 integer nblock ;
577312 extern real clanhs_ (char * , integer * , complex * , integer * , real * );
578- extern /* Subroutine */ void claset_ (char * , integer * , integer * , complex
313+ extern /* Subroutine */ int claset_ (char * , integer * , integer * , complex
579314 * , complex * , complex * , integer * ), clartg_ (complex * ,
580315 complex * , real * , complex * , complex * );
581316 real safmin ;
@@ -849,6 +584,8 @@ static integer c__1 = 1;
849584 k = istop ;
850585 while (k >= istart2 ) {
851586 if (c_abs (& b [k + k * b_dim1 ]) < btol ) {
587+ i__2 = k + k * b_dim1 ;
588+ b [i__2 ].r = 0.f , b [i__2 ].i = 0.f ;
852589/* A diagonal element of B is negligible, move it */
853590/* to the top and deflate it */
854591 i__2 = istart2 + 1 ;
0 commit comments