subroutine funkmin_generic(ndim,beta,fvalue) implicit none include 'forgenericregres.h' integer ndim double precision beta(ndim),fvalue !(in) ndim: the dimension of the parameter vector !(in) beta: the parameters !(out) fvalue: the value of the cost function at beta !----------------------------------------------------- integer i,j,idowhat,nparams double precision dydxp(nyvars,(nxvars+ndim)) ! ! check to see if parameters are out of bounds if(idobounded.eq.1)then if(betamin(1).lt.betamax(1))then do i=1,ndim if(beta(i).lt.betamin(i).or. &beta(i).gt.betamax(i))then ! parameter out of bound fvalue=1.0d+100 return endif enddo endif endif fvalue=0.0d0 if(iregrestype.eq.1)then !orthogonal distance regression do i=1,nobs call shortestdist(nyvars,nxvars,yobs(i:i,1:nyvars), & xvars(i:i,1:nxvars),xmin(i:i,1:nxvars), & xmax(i:i,1:nxvars),ndim,beta,iknowder, & shorty(i:i,1:nyvars),shortx(i:i,1:nxvars)) do j=1,nyvars fvalue=fvalue+weity(i,j)* & (shorty(i,j)-yobs(i,j))**2 enddo do j=1,nxvars fvalue=fvalue+weitx(i,j)* & (shortx(i,j)-xvars(i,j))**2 enddo enddo endif if(iregrestype.eq.0)then idowhat=0 do i=1,nobs call surffunc(nyvars,shorty(i:i,1:nyvars),nxvars, & xvars(i:i,1:nxvars),ndim,beta, & dydxp(1:nyvars,1:(nxvars+ndim)),idowhat) do j=1,nyvars fvalue=fvalue+weity(i,j)* & (shorty(i,j)-yobs(i,j))**2 enddo enddo endif if(iregrestype.eq.2)then nparams=ndim-nobs idowhat=nparams do i=1,nobs do j=1,nxvars idowhat=idowhat+1 shortx(i,j)=beta(idowhat) enddo enddo idowhat=0 do i=1,nobs call surffunc(nyvars,shorty(i:i,1:nyvars),nxvars, & shortx(i:i,1:nxvars),nparams,beta, & dydxp(1:nyvars,1:(nxvars+ndim)),idowhat) do j=1,nyvars fvalue=fvalue+weity(i,j)* & (shorty(i,j)-yobs(i,j))**2 enddo do j=1,nxvars fvalue=fvalue+weitx(i,j)* & (shortx(i,j)-xvars(i,j))**2 enddo enddo endif if(iregrestype.eq.-1)then !implicit orthogonal distance regression idowhat=0 do i=1,nobs call surffunc(nyvars,shorty(i:i,1:nyvars),nxvars, & xvars(i:i,1:nxvars),ndim,beta, & dydxp(1:nyvars,1:(nxvars+ndim)),idowhat) do j=1,nyvars fvalue=fvalue+weity(i,j)* & shorty(i,j)**2 enddo enddo endif return end subroutine funkmin_generic !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ double precision function f1dim_generic(x) implicit none double precision x CU USES funkmin_generic INTEGER j !(((((((((((((((((((((((((((((((((((((((((((((((((((( integer NMAX,ncom parameter(NMAX=1000) double precision pcom(NMAX),xicom(NMAX) COMMON /f1com/ pcom,xicom,ncom save /f1com/ !)))))))))))))))))))))))))))))))))))))))))))))))))))) double precision xt(NMAX) !----------------------------------------------------- do 11 j=1,ncom xt(j)=pcom(j)+x*xicom(j) 11 continue call funkmin_generic(ncom,xt,f1dim_generic) return END !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& SUBROUTINE FCN_generic(N,M,NP,NQ, + LDN,LDM,LDNP, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX, + IDEVAL,F,FJACB,FJACD, + ISTOP) implicit none include 'forgenericregres.h' C SUBROUTINE ARGUMENTS C ==> N NUMBER OF OBSERVATIONS C ==> M NUMBER OF COLUMNS IN EXPLANATORY VARIABLE C ==> NP NUMBER OF PARAMETERS C ==> NQ NUMBER OF RESPONSES PER OBSERVATION C ==> LDN LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING N C ==> LDM LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING M C ==> LDNP LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING NP C ==> BETA CURRENT VALUES OF PARAMETERS C ==> XPLUSD CURRENT VALUE OF EXPLANATORY VARIABLE, I.E., X + DELTA C ==> IFIXB INDICATORS FOR "FIXING" PARAMETERS (BETA) C ==> IFIXX INDICATORS FOR "FIXING" EXPLANATORY VARIABLE (X) C ==> LDIFX LEADING DIMENSION OF ARRAY IFIXX C ==> IDEVAL INDICATOR FOR SELECTING COMPUTATION TO BE PERFORMED C <== F PREDICTED FUNCTION VALUES C <== FJACB JACOBIAN WITH RESPECT TO BETA C <== FJACD JACOBIAN WITH RESPECT TO ERRORS DELTA C <== ISTOP STOPPING CONDITION, WHERE C 0 MEANS CURRENT BETA AND X+DELTA WERE C ACCEPTABLE AND VALUES WERE COMPUTED SUCCESSFULLY C 1 MEANS CURRENT BETA AND X+DELTA ARE C NOT ACCEPTABLE; ODRPACK SHOULD SELECT VALUES C CLOSER TO MOST RECENTLY USED VALUES IF POSSIBLE C -1 MEANS CURRENT BETA AND X+DELTA ARE C NOT ACCEPTABLE; ODRPACK SHOULD STOP C INPUT ARGUMENTS, NOT TO BE CHANGED BY THIS ROUTINE: INTEGER I,IDEVAL,ISTOP,L,LDIFX,LDM,LDN,LDNP,M,N,NP,NQ DOUBLE PRECISION BETA(NP),XPLUSD(LDN,M) INTEGER IFIXB(NP),IFIXX(LDIFX,M) C OUTPUT ARGUMENTS: DOUBLE PRECISION F(LDN,NQ),FJACB(LDN,LDNP,NQ),FJACD(LDN,LDM,NQ) double precision ymod(NQ),dydxp(NQ,(M+NP)) integer k,idowhat !----------------------------------------------------- if(idobounded.eq.1)then if(betamin(1).lt.betamax(1))then do I=1,NP if(BETA(I).lt.betamin(I).or.BETA(I).gt.betamax(I))then ISTOP = 1 RETURN endif enddo endif endif ISTOP=0 IF (MOD(IDEVAL,10).GE.1) THEN idowhat=0 DO 100 I = 1,N call surffunc(NQ,ymod,M,XPLUSD(I:I,1:M), & NP,BETA,dydxp(1:NQ,1:(M+NP)),idowhat) DO 110 L = 1,NQ F(I,L)=ymod(L) 110 CONTINUE 100 CONTINUE END IF C COMPUTE DERIVATIVES WITH RESPECT TO BETA IF (MOD(IDEVAL/10,10).GE.1) THEN idowhat=2 DO 200 I = 1,N call surffunc(NQ,ymod,M,XPLUSD(I:I,1:M), & NP,BETA,dydxp(1:NQ,1:(M+NP)),idowhat) DO 210 L = 1,NQ do k=1,NP FJACB(I,k,L)=dydxp(L,k) enddo 210 CONTINUE 200 CONTINUE ENDIF c compute derivatives with respect to delta IF (MOD(IDEVAL/100,10).GE.1) THEN idowhat=1 DO 300 I = 1,N call surffunc(NQ,ymod,M,XPLUSD(I:I,1:M), & NP,BETA,dydxp(1:NQ,1:(M+NP)),idowhat) DO 310 L = 1,NQ do k=1,M FJACD(I,k,L)=dydxp(L,k) enddo 310 CONTINUE 300 CONTINUE ENDIF RETURN END !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ double precision function generic_pikaia(ndim,beta01) implicit none include 'forgenericregres.h' integer ndim,i double precision beta01(ndim),beta(ndim),fvalue do i=1,ndim ! beta01(i)=(beta(i)-betamin(i))/(betamax(i)-betamin(i)) beta(i)=betamin(i)+beta01(i)*(betamax(i)-betamin(i)) enddo call funkmin_generic(ndim,beta,fvalue) generic_pikaia=1.0d0/(fvalue+0.00001d0) return end !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ subroutine distcenter(nx,x,fequ,fvalue,idowhat) implicit none include 'leastdistance.h' !idowhat=1, evaluating the system of equations and calculating the sum of squares. !idowhat=2, calculating the distance. integer nx,idowhat double precision x(nx),fequ(nx),fvalue !---------------------------------------------------------- integer i,j,ider double precision y(my),dydxp(my,(nx+nparams)), & xcopy(nx),sum,yplush(my),yminush(my),h parameter(h=1.0d-7) !==============End of Variable Declaration================== j=0 call surffunc(my,y,nx,x,nparams,params, & dydxp(1:my,1:(nx+nparams)),j) if(idowhat.eq.1)then if(iknowder.eq.1)then call surffunc(my,y,nx,x,nparams,params, & dydxp(1:my,1:(nx+nparams)),iknowder) endif if(iknowder.eq.0)then do i=1,nx xcopy(i)=x(i) enddo do i=1,nx xcopy(i)=x(i)+h call surffunc(my,yplush,nx,xcopy,nparams,params, & dydxp(1:my,1:(nx+nparams)),iknowder) xcopy(i)=x(i)-h call surffunc(my,yminush,nx,xcopy,nparams,params, & dydxp(1:my,1:(nx+nparams)),iknowder) do j=1,my dydxp(j,i)=(yplush(j)-yminush(j))/(2.0d0*h) enddo xcopy(i)=x(i) enddo endif do i=1,nx sum=0.0d0 do j=1,my sum=sum+(y(j)-targety(j))*dydxp(j,i) enddo fequ(i)=x(i)-(targetx(i)-sum) enddo fvalue=0.0d0 do i=1,nx fvalue=fvalue+fequ(i)*fequ(i) enddo endif if(idowhat.eq.2)then fvalue=0.0d0 do i=1,my fvalue=fvalue+(y(i)-targety(i))**2 enddo do i=1,nx fvalue=fvalue+(x(i)-targetx(i))**2 enddo endif return end subroutine distcenter !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ subroutine distcentersys(nunknowns,x,fequ,fsqsum) implicit none integer nunknowns,idowhat double precision x(nunknowns), & fequ(nunknowns),fsqsum parameter(idowhat=1) call distcenter(nunknowns,x,fequ,fsqsum,idowhat) return end subroutine distcentersys !----------------------------------------------------------- subroutine fsqsum_distcenter(nunknowns,x,fsqsum) implicit none integer nunknowns,idowhat double precision x(nunknowns),fsqsum, & fequ(nunknowns) parameter(idowhat=1) call distcenter(nunknowns,x,fequ,fsqsum,idowhat) return end !----------------------------------------------------------- double precision function f1dimsqsum_distcenter(x) implicit none double precision x INTEGER j,idowhat !(((((((((((((((((((((((((((((((((((((((((((((((((((( integer NMAX,ncom parameter(NMAX=1000) double precision pcom(NMAX),xicom(NMAX) COMMON /cpf1com/ pcom,xicom,ncom save /cpf1com/ !)))))))))))))))))))))))))))))))))))))))))))))))))))) double precision xt(NMAX),fequ(NMAX) parameter(idowhat=1) do 11 j=1,ncom xt(j)=pcom(j)+x*xicom(j) 11 continue call distcenter(ncom,xt,fequ, & f1dimsqsum_distcenter,idowhat) return END !------------------------------------------------------------ subroutine s2_distcenter(nunknowns,x,s2) implicit none integer nunknowns,idowhat double precision x(nunknowns),s2,fequ(nunknowns) parameter(idowhat=2) call distcenter(nunknowns,x,fequ,s2,idowhat) return end !----------------------------------------------------------- double precision function f1dims2_distcenter(x) implicit none double precision x INTEGER j,idowhat !(((((((((((((((((((((((((((((((((((((((((((((((((((( integer NMAX,ncom parameter(NMAX=1000) double precision pcom(NMAX),xicom(NMAX) COMMON /cpf1com/ pcom,xicom,ncom save /cpf1com/ !)))))))))))))))))))))))))))))))))))))))))))))))))))) double precision xt(NMAX),fequ(NMAX) parameter(idowhat=2) do 11 j=1,ncom xt(j)=pcom(j)+x*xicom(j) 11 continue call distcenter(ncom,xt,fequ, & f1dims2_distcenter,idowhat) return END !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$