/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:22 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "spolz.h" #include #include /* PARAMETER translations */ #define C43 (-.4375e0) #define C75 .75e0 #define C95 .95e0 #define HALF .5e0 #define ONE 1.e0 #define ZERO 0.e0 /* end of PARAMETER translations */ void /*FUNCTION*/ spolz( float a[], long ndeg, float z[], float *h, long *ierr) { #define H(I_,J_) (*(h+(I_)*(ndeg)+(J_))) LOGICAL32 more, notlas; long int _l0, en, enm2, i, its, j, k, l, ll, low, m, mm, mp2, n, na; float c, f, g, p, q, r, s, t, w, x, y, zz; static float b2, base, machep; static LOGICAL32 first = TRUE; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const A = &a[0] - 1; float *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2001-05-25 SPOLZ Krogh Minor change for making .f90 version. *>> 1996-04-27 SPOLZ Krogh Changes to use .C. and C%%. *>> 1996-03-30 SPOLZ Krogh Added external statement, MIN0 => MIN. *>> 1995-01-25 SPOLZ Krogh Automate C conversion. *>> 1995-11-17 SPOLZ Krogh SFTRAN converted to Fortran 77 *>> 1994-10-19 SPOLZ Krogh Changes to use M77CON *>> 1992-05-11 SPOLZ CLL *>> 1988-11-16 CLL More editing of spec stmts. *>> 1988-06-07 SPOLZ CLL Reordered spec stmts for ANSI standard. *>> 1987-09-16 SPOLZ Lawson Initial code. *--S replaces "?": ?POLZ, * ------------------------------------------------------------------ * * Given coefficients A(1),...,A(NDEG+1) this subr computes the * NDEG roots of the polynomial * A(1)*X**NDEG + ... + A(NDEG+1) * storing the roots as complex numbers in the array Z( ). * Require NDEG .ge. 1 and A(1) .ne. 0. * * ------------------------------------------------------------------ *-- Begin mask code changes * * Argument Definitions * -------------------- * * A( ) (In) Contains the coefficients of a polynomial, high * order coefficient first with A(1).ne.0. The contents * of this array will not be modified by the subroutine. * * NDEG (In) Degree of the polynomial. * * Z( ) (Out) Contains the polynomial roots stored as complex * numbers. The real and imaginary parts of the Jth roots * will be stored in Z(2*J-1) and Z(2*J) respectively. * * H( ) (Scratch) Array of work space. * * IERR (Out) Error flag. Set by the subroutine to 0 on normal * termination. Set to -1 if A(1)=0. Set to -2 if NDEG * .le. 0. Set to J > 0 if the iteration count limit * has been exceeded and roots 1 through J have not been * determined. * * C.L.Lawson & S.Y.Chan, JPL, June 3,1986. * 1992-05-11 CLL IERR was not being set when N = 0 or 1. Fixed this. * Added type stmts for all variables. *-- End mask code changes * ------------------------------------------------------------------ * */ /*++ Default NO_COMPLEX = .C. | (.N. == 'D') *++ Default COMPLEX = ~NO_COMPLEX *++ CODE for COMPLEX is inactive * COMPLEX Z(*) *++ CODE for NO_COMPLEX is active */ /*++ END */ /* ------------------------------------------------------------------ * */ if (first) { /* Set MACHEP = machine dependent parameter specifying the * relative precision of floating point arithmetic. * BASE = machine dependent parameter specifying the base * of the machine floating point representation. * */ first = FALSE; machep = FLT_EPSILON; base = FLT_RADIX; b2 = base*base; } *ierr = 0; if (ndeg <= 0) { *ierr = -2; ermsg( "SPOLZ", *ierr, 0, "NDEG .LE. 0", '.' ); return; } if (A[1] == ZERO) { *ierr = -1; ermsg( "SPOLZ", *ierr, 0, "A(1) .EQ. ZERO", '.' ); return; } n = ndeg; *ierr = 0; /* Build first row of companion matrix. * */ for (i = 2; i <= (ndeg + 1); i++) { H(i - 2,0) = -(A[i]/A[1]); } /* Extract any exact zero roots and set N = degree of * remaining polynomial. * */ for (j = ndeg; j >= 1; j--) { if (H(j - 1,0) != ZERO) goto L_30; /*++ Code for COMPLEX is inactive * Z(J) = ZERO *++ Code for NO_COMPLEX is active */ Z[2*j - 1] = ZERO; Z[2*j] = ZERO; /*++ End */ n -= 1; } L_30: ; /* Special for N = 0 or 1. * */ if (n == 0) return; if (n == 1) { Z[1] = H(0,0); return; } /* Build rows 2 thru N of the companion matrix. * */ for (i = 2; i <= n; i++) { for (j = 1; j <= n; j++) { H(j - 1,i - 1) = ZERO; } H(i - 2,i - 1) = ONE; } /* ***************** BALANCE THE MATRIX *********************** * * This is an adaption of the EISPACK subroutine BALANC to * the special case of a companion matrix. The EISPACK BALANCE * is a translation of the ALGOL procedure BALANCE, NUM. MATH. * 13, 293-304(1969) by Parlett and Reinsch. * HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). * * ********** ITERATIVE LOOP FOR NORM REDUCTION ********** */ L_60: ; more = FALSE; for (i = 1; i <= n; i++) { /* Compute R = sum of magnitudes in row I skipping diagonal. * C = sum of magnitudes in col I skipping diagonal. */ if (i == 1) { r = fabsf( H(1,0) ); for (j = 3; j <= n; j++) { r += fabsf( H(j - 1,0) ); } c = fabsf( H(0,1) ); } else { r = fabsf( H(i - 2,i - 1) ); c = fabsf( H(i - 1,0) ); if (i != n) { c += fabsf( H(i - 1,i) ); } } /* Determine column scale factor, F. * */ g = r/base; f = ONE; s = c + r; L_80: if (c < g) { f *= base; c *= b2; goto L_80; } g = r*base; L_90: if (c >= g) { f /= base; c /= b2; goto L_90; } /* Will the factor F have a significant effect ? * */ if ((c + r)/f < C95*s) { /* Yes, so do the scaling. * */ g = ONE/f; more = TRUE; /* Scale Row I * */ if (i == 1) { for (j = 1; j <= n; j++) { H(j - 1,0) *= g; } } else { H(i - 2,i - 1) *= g; } /* Scale Column I * */ H(i - 1,0) *= f; if (i != n) H(i - 1,i) *= f; } } if (more) goto L_60; /* ***************** QR EIGENVALUE ALGORITHM *********************** * * This is the EISPACK subroutine HQR that uses the QR * algorithm to compute all eigenvalues of an upper * Hessenberg matrix. Original ALGOL code was due to Martin, * Peters, and Wilkinson, Numer. Math., 14, 219-231(1970). * */ low = 1; en = n; t = ZERO; /* ********** SEARCH FOR NEXT EIGENVALUES ********** */ L_160: if (en < low) goto L_1001; its = 0; na = en - 1; enm2 = na - 1; /* ********** LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT * FOR L=EN STEP -1 UNTIL LOW DO -- ********** */ L_170: for (ll = low; ll <= en; ll++) { l = en + low - ll; if (l == low) goto L_200; if (fabsf( H(l - 2,l - 1) ) <= machep*(fabsf( H(l - 2,l - 2) ) + fabsf( H(l - 1,l - 1) ))) goto L_200; } /* ********** FORM SHIFT ********** */ L_200: x = H(en - 1,en - 1); if (l == en) goto L_370; y = H(na - 1,na - 1); w = H(na - 1,en - 1)*H(en - 1,na - 1); if (l == na) goto L_380; if (its == 30) goto L_1000; if (its != 10 && its != 20) goto L_230; /* ********** FORM EXCEPTIONAL SHIFT ********** */ t += x; for (i = low; i <= en; i++) { H(i - 1,i - 1) -= x; } s = fabsf( H(na - 1,en - 1) ) + fabsf( H(enm2 - 1,na - 1) ); x = C75*s; y = x; w = C43*s*s; L_230: its += 1; /* ********** LOOK FOR TWO CONSECUTIVE SMALL * SUB-DIAGONAL ELEMENTS. * FOR M=EN-2 STEP -1 UNTIL L DO -- ********** */ for (mm = l; mm <= enm2; mm++) { m = enm2 + l - mm; zz = H(m - 1,m - 1); r = x - zz; s = y - zz; p = (r*s - w)/H(m - 1,m) + H(m,m - 1); q = H(m,m) - zz - r - s; r = H(m,m + 1); s = fabsf( p ) + fabsf( q ) + fabsf( r ); p /= s; q /= s; r /= s; if (m == l) goto L_250; if (fabsf( H(m - 2,m - 1) )*(fabsf( q ) + fabsf( r )) <= machep* fabsf( p )*(fabsf( H(m - 2,m - 2) ) + fabsf( zz ) + fabsf( H(m,m) ))) goto L_250; } L_250: mp2 = m + 2; for (i = mp2; i <= en; i++) { H(i - 3,i - 1) = ZERO; if (i == mp2) goto L_260; H(i - 4,i - 1) = ZERO; L_260: ; } /* ********** DOUBLE QR STEP INVOLVING ROWS L TO EN AND * COLUMNS M TO EN ********** */ for (k = m; k <= na; k++) { notlas = k != na; if (k == m) goto L_270; p = H(k - 2,k - 1); q = H(k - 2,k); r = ZERO; if (notlas) r = H(k - 2,k + 1); x = fabsf( p ) + fabsf( q ) + fabsf( r ); if (x == ZERO) goto L_360; p /= x; q /= x; r /= x; L_270: s = signf( sqrtf( p*p + q*q + r*r ), p ); if (k == m) goto L_280; H(k - 2,k - 1) = -s*x; goto L_290; L_280: if (l != m) H(k - 2,k - 1) = -H(k - 2,k - 1); L_290: p += s; x = p/s; y = q/s; zz = r/s; q /= p; r /= p; /* ********** ROW MODIFICATION ********** */ for (j = k; j <= en; j++) { p = H(j - 1,k - 1) + q*H(j - 1,k); if (!notlas) goto L_300; p += r*H(j - 1,k + 1); H(j - 1,k + 1) += -p*zz; L_300: H(j - 1,k) += -p*y; H(j - 1,k - 1) += -p*x; } j = min( en, k + 3 ); /* ********** COLUMN MODIFICATION ********** */ for (i = l; i <= j; i++) { p = x*H(k - 1,i - 1) + y*H(k,i - 1); if (!notlas) goto L_320; p += zz*H(k + 1,i - 1); H(k + 1,i - 1) += -p*r; L_320: H(k,i - 1) += -p*q; H(k - 1,i - 1) -= p; } L_360: ; } goto L_170; /* ********** ONE ROOT FOUND ********** *++ Code for COMPLEX is inactive * 370 Z(EN) = CMPLX(X+T,ZERO) *++ Code for NO_COMPLEX is active */ L_370: Z[2*en - 1] = x + t; Z[2*en] = ZERO; /*++ End */ en = na; goto L_160; /* ********** TWO ROOTS FOUND ********** */ L_380: p = (y - x)*HALF; q = p*p + w; zz = sqrtf( fabsf( q ) ); x += t; if (q < ZERO) goto L_420; /* ********** PAIR OF REALS ********** */ zz = p + signf( zz, p ); /*++ Code for COMPLEX is inactive * Z(NA) = CMPLX(X+ZZ,ZERO) * Z(EN) = Z(NA) * IF (ZZ .NE. ZERO) Z(EN) = CMPLX(X-W/ZZ,ZERO) *++ Code for NO_COMPLEX is active */ Z[2*na - 1] = x + zz; Z[2*na] = ZERO; Z[2*en - 1] = Z[2*na - 1]; Z[2*en] = Z[2*na]; if (zz != ZERO) { Z[2*en - 1] = x - w/zz; Z[2*en] = ZERO; } /*++ End */ goto L_430; /* ********** COMPLEX PAIR ********** *++ Code for COMPLEX is inactive * 420 Z(NA) = CMPLX(X+P,ZZ) * Z(EN) = CMPLX(X+P,-ZZ) *++ Code for NO_COMPLEX is active */ L_420: Z[2*na - 1] = x + p; Z[2*na] = zz; Z[2*en - 1] = x + p; Z[2*en] = -zz; /*++ End */ L_430: en = enm2; goto L_160; /* ********** SET ERROR -- NO CONVERGENCE TO AN * EIGENVALUE AFTER 30 ITERATIONS ********** */ L_1000: *ierr = en; L_1001: ; if (*ierr != 0) { ermsg( "SPOLZ", *ierr, 0, "Convergence failure", '.' ); } return; #undef H } /* end of function */