C C .................................................................. C C SUBROUTINE DARAT C C PURPOSE C CALCULATE BEST RATIONAL APPROXIMATION OF A DISCRETE C FUNCTION IN THE LEAST SQUARES SENSE C C USAGE C CALL DARAT(DATI,N,WORK,P,IP,IQ,IER) C C DESCRIPTION OF PARAMETERS C DATI - TWODIMENSIONAL ARRAY WITH 3 COLUMNS AND N ROWS C THE FIRST COLUMN MUST CONTAIN THE GIVEN ARGUMENTS, C THE SECOND COLUMN THE GIVEN FUNCTION VALUES AND C THE THIRD COLUMN THE GIVEN WEIGHTS IF ANY. C IF NO WEIGHTS ARE TO BE USED THEN THE THIRD C COLUMN MAY BE DROPPED , EXCEPT THE FIRST ELEMENT C WHICH MUST CONTAIN A NONPOSITIVE VALUE C DATI MUST BE OF DOUBLE PRECISION C N - NUMBER OF NODES OF THE GIVEN DISCRETE FUNCTION C WORK - WORKING STORAGE WHICH IS OF DIMENSION C (IP+IQ)*(IP+IQ+1)+4*N+1 AT LEAST. C ON RETURN THE VALUES OF THE NUMERATOR ARE CONTAINED C IN WORK(N+1) UP TO WORK(2*N), WHILE THE VALUES OF C THE DENOMINATOR ARE STORED IN WORK(2*N+1) UP TO C WORK(3*N) C WORK MUST BE OF DOUBLE PRECISION C P - RESULTANT COEFFICIENT VECTOR OF DENOMINATOR AND C NUMERATOR. THE DENOMINATOR IS STORED IN FIRST IQ C LOCATIONS, THE NUMERATOR IN THE FOLLOWING IP C LOCATIONS. C COEFFICIENTS ARE ORDERED FROM LOW TO HIGH. C P MUST BE OF DOUBLE PRECISION C IP - DIMENSION OF THE NUMERATOR (INPUT VALUE) C IQ - DIMENSION OF THE DENOMINATOR (INPUT VALUE) C IER - RESULTANT ERROR PARAMETER C IER =-1 MEANS FORMAL ERRORS C IER = 0 MEANS NO ERRORS C IER = 1,2 MEANS POOR CONVERGENCE OF ITERATION C IER IS ALSO USED AS INPUT VALUE C A NONZERO INPUT VALUE INDICATES AVAILABILITY OF AN C INITIAL APPROXIMATION STORED IN P C C REMARKS C THE COEFFICIENT VECTORS OF THE DENOMINATOR AND NUMERATOR C OF THE RATIONAL APPROXIMATION ARE BOTH STORED IN P C STARTING WITH LOW POWERS (DENOMINATOR FIRST). C IP+IQ MUST NOT EXCEED N, ALL THREE VALUES MUST BE POSITIVE. C SINCE CHEBYSHEV POLYNOMIALS ARE USED AS FUNDAMENTAL C FUNCTIONS, THE ARGUMENTS SHOULD BE REDUCED TO THE INTERVAL C (-1,1). THIS CAN ALWAYS BE ACCOMPLISHED BY MEANS OF A LINEAR C TRANSFORMATION OF THE ORIGINALLY GIVEN ARGUMENTS. C IF A FIT IN OTHER FUNCTIONS IS REQUIRED, DCNP AND DCNPS MUST C BE REPLACED BY SUBROUTINES WHICH ARE OF ANALOGOUS DESIGN. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C DAPLL, DAPFS, DFRAT, DCNPS, DCNP C DCNP IS REQUIRED WITHIN DFRAT C C METHOD C THE ITERATIVE SCHEME USED FOR CALCULATION OF THE C APPROXIMATION IS REPEATED SOLUTION OF THE NORMAL EQUATIONS C WHICH ARE OBTAINED BY LINEARIZATION. C A REFINED TECHNIQUE OF THIS LINEAR LEAST SQUARES APPROACH C IS USED WHICH GUARANTEES THAT THE DENOMINATOR IS FREE OF C ZEROES WITHIN THE APPROXIMATION INTERVAL. C FOR REFERENCE SEE C D.BRAESS, UEBER DAEMPFUNG BEI MINIMALISIERUNGSVERFAHREN, C COMPUTING(1966), VOL.1, ED.3, PP.264-272. C D.W.MARQUARDT, AN ALGORITHM FOR LEAST-SQUARES ESTIMATION C OF NONLINEAR PARAMETERS, C JSIAM(1963), VOL.11, ED.2, PP.431-441. C C .................................................................. C SUBROUTINE DARAT(DATI,N,WORK,P,IP,IQ,IER) C C EXTERNAL DFRAT C C DIMENSIONED LOCAL VARIABLE DIMENSION IERV(3) C C DIMENSIONED DUMMY VARIABLES DIMENSION DATI(1),WORK(1),P(1) DOUBLE PRECISION DATI,WORK,P,T,OSUM,DIAG,RELAX,SUM,SSOE,SAVE C C INITIALIZE TESTVALUES LIMIT=20 ETA=1.E-29 EPS=1.E-14 C C CHECK FOR FORMAL ERRORS IF(N)4,4,1 1 IF(IP)4,4,2 2 IF(IQ)4,4,3 3 IPQ=IP+IQ IF(N-IPQ)4,5,5 C C ERROR RETURN IN CASE OF FORMAL ERRORS 4 IER=-1 RETURN C C INITIALIZE ITERATION PROCESS 5 KOUNT=0 IERV(2)=IP IERV(3)=IQ NDP=N+N+1 NNE=NDP+NDP IX=IPQ-1 IQP1=IQ+1 IRHS=NNE+IPQ*IX/2 IEND=IRHS+IX C C TEST FOR AVAILABILITY OF AN INITIAL APPROXIMATION IF(IER)8,6,8 C C INITIALIZE NUMERATOR AND DENOMINATOR 6 DO 7 I=2,IPQ 7 P(I)=0.D0 P(1)=1.D0 C C CALCULATE VALUES OF NUMERATOR AND DENOMINATOR FOR INITIAL C APPROXIMATION 8 DO 9 J=1,N T=DATI(J) I=J+N CALL DCNPS(WORK(I),T,P(IQP1),IP) K=I+N 9 CALL DCNPS(WORK(K),T,P,IQ) C C SET UP NORMAL EQUATIONS (MAIN LOOP OF ITERATION) 10 CALL DAPLL(DFRAT,N,IX,WORK,WORK(IEND+1),DATI,IERV) C C CHECK FOR ZERO DENOMINATOR IF(IERV(1))4,11,4 11 INCR=0 RELAX=2.D0 C C RESTORE MATRIX IN WORKING STORAGE 12 J=IEND DO 13 I=NNE,IEND J=J+1 13 WORK(I)=WORK(J) IF(KOUNT)14,14,15 C C SAVE SQUARE SUM OF ERRORS 14 OSUM=WORK(IEND) DIAG=OSUM*EPS K=IQ C C ADD CONSTANT TO DIAGONAL IF(WORK(NNE))17,17,19 15 IF(INCR)19,19,16 16 K=IPQ 17 J=NNE-1 DO 18 I=1,K WORK(J)=WORK(J)+DIAG 18 J=J+I C C SOLVE NORMAL EQUATIONS 19 CALL DAPFS(WORK(NNE),IX,IRES,1,EPS,ETA,IER) C C CHECK FOR FAILURE OF EQUATION SOLVER IF(IRES)4,4,20 C C TEST FOR DEFECTIVE NORMALEQUATIONS 20 IF(IRES-IX)21,24,24 21 IF(INCR)22,22,23 22 DIAG=DIAG*0.125D0 23 DIAG=DIAG+DIAG INCR=INCR+1 C C START WITH OVER RELAXATION RELAX=8.D0 IF(INCR-LIMIT)12,45,45 C C CALCULATE VALUES OF CHANGE OF NUMERATOR AND DENOMINATOR 24 L=NDP J=NNE+IRES*(IRES-1)/2-1 K=J+IQ WORK(J)=0.D0 IRQ=IQ IRP=IRES-IQ+1 IF(IRP)25,26,26 25 IRQ=IRES+1 26 DO 29 I=1,N T=DATI(I) WORK(I)=0.D0 CALL DCNPS(WORK(I),T,WORK(K),IRP) M=L+N CALL DCNPS(WORK(M),T,WORK(J),IRQ) IF(WORK(M)*WORK(L))27,29,29 27 SUM=WORK(L)/WORK(M) IF(RELAX+SUM)29,29,28 28 RELAX=-SUM 29 L=L+1 C C MODIFY RELAXATION FACTOR IF NECESSARY SSOE=OSUM ITER=LIMIT 30 SUM=0.D0 RELAX=RELAX*0.5D0 DO 32 I=1,N M=I+N K=M+N L=K+N SAVE=DATI(M)-(WORK(M)+RELAX*WORK(I))/(WORK(K)+RELAX*WORK(L)) SAVE=SAVE*SAVE IF(DATI(NDP))32,32,31 31 SAVE=SAVE*DATI(K) 32 SUM=SUM+SAVE IF(ITER)45,33,33 33 ITER=ITER-1 IF(SUM-OSUM)34,37,35 34 OSUM=SUM GOTO 30 C C TEST FOR IMPROVEMENT 35 IF(OSUM-SSOE)36,30,30 36 RELAX=RELAX+RELAX 37 T=0. SAVE=0.D0 K=IRES+1 DO 38 I=2,K J=J+1 T=T+DABS(P(I)) P(I)=P(I)+RELAX*WORK(J) 38 SAVE=SAVE+DABS(P(I)) C C UPDATE CURRENT VALUES OF NUMERATOR AND DENOMINATOR DO 39 I=1,N J=I+N K=J+N L=K+N WORK(J)=WORK(J)+RELAX*WORK(I) 39 WORK(K)=WORK(K)+RELAX*WORK(L) C C TEST FOR CONVERGENCE IF(INCR)40,40,42 40 IF(SSOE-OSUM-RELAX*OSUM*DBLE(EPS))46,46,41 41 IF(DABS(T-SAVE)-RELAX*SAVE*DBLE(EPS))46,46,42 42 IF(OSUM-SAVE*DBLE(ETA))46,46,43 43 KOUNT=KOUNT+1 IF(KOUNT-LIMIT)10,44,44 C C ERROR RETURN IN CASE OF POOR CONVERGENCE 44 IER=2 RETURN 45 IER=1 RETURN C C NORMAL RETURN 46 IER=0 RETURN END