Re: NEC-LIST: routine for finding complex roots wanted

riccardo,

I used this quite a while ago (13 yrs). (polrt.for)
I also included an evaluation routine. (eval.for)
I believe this was for a 'wideband combline'
filter design program (wenzel) way back in '87.

happy programming,
mike

-- 
|||  Michael A. Ciminera 
|||  Exciter & RF/MMWave Instruments Group
|||  Jet Propulsion Laboratory --- Section 333
|||  818 354 4356 (tel) --- 818 354 2825 (fax)
|||  email --- mciminer_at_jazz.jpl.nasa.gov
Riccardo Leone wrote:
> 
> Does anybody know, please, where may I find a fortran (or C) routine
> for complex roots seeking, numerically, of a generic complex function?
> (that is a MoM matrix determinant) Websites, or books...
	SUBROUTINE POLRT(XCOF,COF,M,ROOTR,ROOTI,IER)
C -----------------------------------------------------------------------------
C		POLRT subroutine, which finds complex roots 
C		of polynomial with real coefficients.    
C		m.ciminera    3/24/87
C -----------------------------------------------------------------------------
	REAL XCOF(1:37),COF(1:37),ROOTR(1:37),ROOTI(1:37)
	IFIT = 0
	N = M
	IER = 0
	IF (XCOF(N+1)) 10,25,10
  10	IF (N) 15,15,32
C		SET ERROR CODE TO 1
  15	IER = 1
  20	RETURN
C		SET ERROR CODE TO 4
  25	IER = 4
	GO TO 20
C		SET ERROR CODE TO 2
  30	IER = 2
	GO TO 20
  32	IF (N - 36) 35,35,30
  35	NX = N
	NXX = N+1
	N2 = 1
	KJ1 = N+1
	DO 40 L=1,KJ1
	MT = KJ1 - L + 1
  40	COF(MT) = XCOF(L)
C		SET INITIAL VALUES
  45	X0 = .00500101
	Y0 = 0.01000101
C		ZERO INITIAL VALUE COUNTER
	IN = 0
  50	X = X0
C		INCREMENT INITIAL VALUES AND COUNTER
	X0 = -10.0 * Y0
	Y0 = -10.0 * X
C		SET X AND Y TO CURRENT VALUE
	X = X0
	Y = Y0
	IN = IN + 1
	GO TO 59
  55	IFIT = 1
	XPR = X
	YPR = Y
C		EVALUATE POLYNOMIAL AND DERIVATIVES
  59	ICT = 0
  60	UX = 0.0
	UY = 0.0
	V  = 0.0
	YT = 0.0
	XT = 1.0
	U = COF(N+1)
	IF(U) 65,130,65
  65	DO 70 I=1,N
	L = N - I + 1
	XT2 = X * XT - Y * YT
	YT2 = X * YT + Y * XT
	U = U + COF(L) * XT2
	V = V + COF(L) * YT2
	FI = I
C	^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  1 OR I ?????
	UX = UX + FI * XT * COF(L)
	UY = UY - FI * YT * COF(L)
	XT = XT2
  70	YT = YT2
	SUMSQ = UX * UX + UY * UY
	IF (SUMSQ) 75,110,75
  75	DX = (V * UY - U * UX) / SUMSQ
	X = X + DX
	DY = -(U * UY + V * UX) / SUMSQ
	Y = Y + DY
  78	IF (ABS(DY) + ABS (DX) - 1.0E-05) 100,80,80
C		STEP ITERATION COUNTER
  80	ICT = ICT + 1
	IF (ICT - 500) 60,85,85
  85	IF (IFIT) 100,90,100
  90	IF (IN - 5) 50,95,95
C		SET ERROR CODE TO 3
  95	IER = 3
	GO TO 20
C -------------------------------------
 100	DO 105 L=1,NXX
	MT = KJ1 - L + 1
	TEMP     = XCOF(MT)
	XCOF(MT) = COF(L)
 105	COF(L)   = TEMP
	ITEMP = N
	N = NX
	NX = ITEMP
	IF (IFIT) 120,55,120
 110	IF (IFIT) 115,50,115
 115	X = XPR
	Y = YPR
 120	IFIT = 0
	IF (X) 122,125,122
 122	IF (ABS(Y) - ABS(X)*1.0E-04) 135,125,125
 125	ALPHA = X + X
	SUMSQ = X * X + Y * Y
	N = N - 2
	GO TO 140
 130	X = 0.0
	NX  = NX  - 1
	NXX = NXX - 1
 135	Y = 0.0
	SUMSQ = 0.0
	ALPHA = X
	N = N - 1
 140	L1 = 1
	L2 = 2
	COF(L2) = COF(L2) + ALPHA * COF(L1)
 145	DO 150 L=2,N
 150	COF(L+1) = COF(L+1) + ALPHA * COF(L) - SUMSQ * COF(L-1)
 155	ROOTI(N2) = Y
	ROOTR(N2) = X
	N2 = N2 + 1
	IF (SUMSQ) 160,165,160
 160	Y = -Y
	SUMSQ = 0.0
	GO TO 155
 165	IF(N) 20,20,45
	END
--------------
	PROGRAM EVAL
C -----------------------------------------------------------------------------
C		Evaluates POLRT subroutine, which finds complex roots 
C		of polynomial with real coefficients.    
C			link eval,polrt
C		m.ciminera  3/23/87 
C -----------------------------------------------------------------------------
	REAL XCOF(0:36),ROOTR(1:37),ROOTI(1:37)
	XCOF(0) = 1.61
	XCOF(1) = 8.936
	XCOF(2) = 21.072
	XCOF(3) = 25.77
	XCOF(4) = 18.04
	XCOF(5) = 6.53
	XCOF(6) = 1
	M = 6
C -----------------------------------------------------------------------------
	CALL POLRT(XCOF,COF,M,ROOTR,ROOTI,IER)
	TYPE *, ' '
100	DO 120 J=1,M
	TYPE *, J,'   ROOTR = ',ROOTR(J), '   ROOTI = ',ROOTI(J)
120	CONTINUE
	TYPE *, 'IER = ',IER
	TYPE *, ' '
	END
--------------