/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:23 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "zpolz.h"
#include <float.h>
#include <stdlib.h>
/* PARAMETER translations */
#define C95 .95e0
#define ONE 1.e0
#define ZERO 0.e0
/* end of PARAMETER translations */
void /*FUNCTION*/ zpolz(
double a[][2],
long ndeg,
double z[][2],
double *h,
long *ierr)
{
#define H(I_,J_,K_) (*(h+(I_)*(ndeg)*(ndeg)+(J_)*(ndeg)+(K_)))
LOGICAL32 more;
long int _l0, i, j, n;
double c, f, g, r, s, temp[2];
static double b2, base;
static LOGICAL32 first = TRUE;
/* OFFSET Vectors w/subscript range: 1 to dimension */
double *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 ZPOLZ Krogh Changes to use .C. and C%%.
*>> 1996-03-30 ZPOLZ Krogh Added external statement, removed "!..."
*>> 1995-01-18 ZPOLZ Krogh More M77CON for conversion to C.
*>> 1995-11-17 ZPOLZ Krogh Added M77CON statements for conversion to C
*>> 1995-11-17 ZPOLZ 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 ZPOLZ Lawson Initial code.
*--Z Replaces "?": ?POLZ, ?QUO
*--D (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( "ZPOLZ", *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( "ZPOLZ", *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 */
zquo( &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 = fabs( H(0,1,0) ) + fabs( H(1,1,0) );
for (j = 3; j <= n; j++)
{
r += fabs( H(0,j - 1,0) ) + fabs( H(1,j - 1,0) );
}
c = fabs( H(0,0,1) ) + fabs( H(1,0,1) );
}
else
{
r = fabs( H(0,i - 2,i - 1) ) + fabs( H(1,i - 2,i - 1) );
c = fabs( H(0,i - 1,0) ) + fabs( H(1,i - 1,0) );
if (i != n)
{
c += fabs( H(0,i - 1,i) ) + fabs( 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;
dcomqr( ndeg, n, 1, n, &H(0,0,0), &H(1,0,0), (double*)z, ierr );
if (*ierr != 0)
{
ermsg( "ZPOLZ", *ierr, 0, "Convergence failure", '.' );
}
return;
#undef H
} /* end of function */