      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
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
