/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:14 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "cpolz.h" #include #include /* PARAMETER translations */ #define C95 .95e0 #define ONE 1.e0 #define ZERO 0.e0 /* end of PARAMETER translations */ void /*FUNCTION*/ cpolz( float a[][2], long ndeg, float z[][2], float *h, long *ierr) { #define H(I_,J_,K_) (*(h+(I_)*(ndeg)*(ndeg)+(J_)*(ndeg)+(K_))) LOGICAL32 more; long int _l0, i, j, n; float c, f, g, r, s, temp[2]; static float b2, base; static LOGICAL32 first = TRUE; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Temp = &temp[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. *>> 1996-04-27 CPOLZ Krogh Changes to use .C. and C%%. *>> 1996-03-30 CPOLZ Krogh Added external statement, removed "!..." *>> 1995-01-18 CPOLZ Krogh More M77CON for conversion to C. *>> 1995-11-17 CPOLZ Krogh Added M77CON statements for conversion to C *>> 1995-11-17 CPOLZ Krogh Converted SFTRAN to Fortran 77. *>> 1992-05-11 CLL IERR was not being set when N = 0 or 1. Fixed this. *>> 1989-10-20 CLL Delcared all variables. *>> 1987-02-25 CPOLZ Lawson Initial code. *--C Replaces "?": ?POLZ, ?QUO *--S (Type) Replaces "?": ?COMQR *++ Default NO_COMPLEX = .C. | (.N. == 'D') *++ Default COMPLEX = ~NO_COMPLEX * ------------------------------------------------------------------ * * In the discussion below, the notation A([*,],k} should be interpreted * as the complex number A(k) if A is declared complex, and should be * interpreted as the complex number A(1,k) + i * A(2,k) if A is not * declared to be of type complex. Similar statements apply for Z(k). * Given complex 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. * * ------------------------------------------------------------------ * * Argument Definitions * -------------------- * * A( ) (In) Contains the complex coefficients of a polynomial * high order coefficient first, with A([*,]1).ne.0. The * real and imaginary parts of the Jth coefficient must * be provided in A([*],J). 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 root * will be stored in Z([*,]J). * * 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. * ------------------------------------------------------------------ */ /*++ CODE for COMPLEX is inactive * COMPLEX A(NDEG+1),TEMP,Z(NDEG) *++ CODE for NO_COMPLEX is active */ /*++ END * */ /* ------------------------------------------------------------------ * */ if (first) { /* Set BASE = machine dependent parameter specifying the base * of the machine floating point representation. * */ first = FALSE; base = FLT_RADIX; b2 = base*base; } if (ndeg <= 0) { *ierr = -2; ermsg( "CPOLZ", *ierr, 0, "NDEG .LE. 0", '.' ); return; } /*++ CODE for COMPLEX is inactive * IF (A(1) .EQ. CMPLX(ZERO, ZERO)) THEN *++ CODE for NO_COMPLEX is active */ if (a[0][0] == ZERO && a[0][1] == ZERO) { /*++ END */ *ierr = -1; ermsg( "CPOLZ", *ierr, 0, "A(*,1) .EQ. ZERO", '.' ); return; } n = ndeg; *ierr = 0; /* Build first row of companion matrix. * */ for (i = 2; i <= (n + 1); i++) { /*++ CODE for COMPLEX is inactive * TEMP = -(A(I)/A(1)) * H(1,I-1,1) = REAL(TEMP) * H(1,I-1,2) = AIMAG(TEMP) *++ CODE for NO_COMPLEX is active */ cquo( &a[i - 1][0], &a[0][0], temp ); H(0,i - 2,0) = -Temp[1]; H(1,i - 2,0) = -Temp[2]; /*++ END */ } /* Extract any exact zero roots and set N = degree of * remaining polynomial. * */ for (j = ndeg; j >= 1; j--) { /*++ CODE for COMPLEX is inactive * IF (H(1,J,1).NE.ZERO .OR. H(1,J,2).NE.ZERO) go to 40 * Z(J) = ZERO *++ CODE for NO_COMPLEX is active */ if (H(0,j - 1,0) != ZERO || H(1,j - 1,0) != ZERO) goto L_40; z[j - 1][0] = ZERO; z[j - 1][1] = ZERO; /*++ END */ n -= 1; } L_40: ; /* Special for N = 0 or 1. * */ if (n == 0) return; if (n == 1) { /*++ CODE for COMPLEX is inactive * Z(1) = CMPLX(H(1,1,1),H(1,1,2)) *++ CODE for NO_COMPLEX is active */ z[0][0] = H(0,0,0); z[0][1] = H(1,0,0); /*++ END */ return; } /* Build rows 2 thru N of the companion matrix. * */ for (i = 2; i <= n; i++) { for (j = 1; j <= n; j++) { if (j == i - 1) { H(0,j - 1,i - 1) = ONE; H(1,j - 1,i - 1) = ZERO; } else { H(0,j - 1,i - 1) = ZERO; H(1,j - 1,i - 1) = ZERO; } } } /* ***************** BALANCE THE MATRIX *********************** * * This is an adaption of the EISPACK subroutine BALANC to * the special case of a complex 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_100: ; 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(0,1,0) ) + fabsf( H(1,1,0) ); for (j = 3; j <= n; j++) { r += fabsf( H(0,j - 1,0) ) + fabsf( H(1,j - 1,0) ); } c = fabsf( H(0,0,1) ) + fabsf( H(1,0,1) ); } else { r = fabsf( H(0,i - 2,i - 1) ) + fabsf( H(1,i - 2,i - 1) ); c = fabsf( H(0,i - 1,0) ) + fabsf( H(1,i - 1,0) ); if (i != n) { c += fabsf( H(0,i - 1,i) ) + fabsf( H(1,i - 1,i) ); } } /* Determine column scale factor, F. * */ g = r/base; f = ONE; s = c + r; L_140: if (c < g) { f *= base; c *= b2; goto L_140; } g = r*base; L_160: if (c >= g) { f /= base; c /= b2; goto L_160; } /* 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(0,j - 1,0) *= g; H(1,j - 1,0) *= g; } } else { H(0,i - 2,i - 1) *= g; H(1,i - 2,i - 1) *= g; } /* Scale Column I * */ H(0,i - 1,0) *= f; H(1,i - 1,0) *= f; if (i != n) { H(0,i - 1,i) *= f; H(1,i - 1,i) *= f; } } } if (more) goto L_100; scomqr( ndeg, n, 1, n, &H(0,0,0), &H(1,0,0), (float*)z, ierr ); if (*ierr != 0) { ermsg( "CPOLZ", *ierr, 0, "Convergence failure", '.' ); } return; #undef H } /* end of function */