!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      subroutine odr_leastsquare(NP,FCN,BETA,N,X,M,Y,NQ,
     &weitx,weity,iderivative,shortx,shorty,fvalue,INFO)
      implicit none
!if derivatives are provided, set iderivative to 1, otherwise set it to 0.
!for ordinary least square regression, set INFO to 0.
!for explicit orthorgonal distance regression, set INFO to 1.
!the content of INFO is destroyed on return
C  ODRPACK ARGUMENT DEFINITIONS
C      ==> FCN      NAME OF THE USER SUPPLIED FUNCTION SUBROUTINE
C      ==> N        NUMBER OF OBSERVATIONS 
C      ==> M        COLUMNS OF DATA IN THE EXPLANATORY VARIABLE
C      ==> NP       NUMBER OF PARAMETERS
C      ==> NQ       NUMBER OF RESPONSES PER OBSERVATION
C     <==> BETA     FUNCTION PARAMETERS
C      ==> Y        RESPONSE VARIABLE
C      ==> X        EXPLANATORY VARIABLE
C      ==> LWORK    DIMENSION OF VECTOR WORK
C      ==> LIWORK   DIMENSION OF VECTOR IWORK
C     <==  INFO     STOPPING CONDITION 
       
C  VARIABLE DECLARATIONS
      INTEGER INFO,M,N,NP,NQ,iderivative,LWORK,LIWORK
      double precision weity(N,NQ),weitx(N,M),shorty(N,NQ),
     &shortx(N,M),fvalue,BETA(NP),X(N,M),Y(N,NQ)
      EXTERNAL FCN
      LWORK=18+11*NP+NP**2+M+M**2+4*N*NQ+6*N*M+2*N*NQ*NP+
     &2*N*NQ*M+NQ**2+5*NQ+NQ*(NP+M)+N*1*NQ
      LIWORK=20+NP+NQ*(NP+M)
      call odr_interface(NP,FCN,BETA,N,X(1:N,1:M),M,Y(1:N,1:NQ),NQ,
     &LWORK,LIWORK,weitx(1:N,1:M),weity(1:N,1:NQ),iderivative,
     &shortx(1:N,1:M),shorty(1:N,1:NQ),fvalue,INFO)
      return
      end subroutine odr_leastsquare
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
      subroutine odr_interface(NP,FCN,BETA,N,X,M,Y,NQ,LWORK,
     &LIWORK,weitx,weity,iderivative,shortx,shorty,fvalue,INFO)
      implicit none
!if derivatives are provided, set iderivative to 1, otherwise set it to 0.
!for ordinary least square regression, set INFO to 0.
!for explicit orthorgonal distance regression, set INFO to 1.
!the content of INFO is destroyed on return
C  ODRPACK ARGUMENT DEFINITIONS
C      ==> FCN      NAME OF THE USER SUPPLIED FUNCTION SUBROUTINE
C      ==> N        NUMBER OF OBSERVATIONS 
C      ==> M        COLUMNS OF DATA IN THE EXPLANATORY VARIABLE
C      ==> NP       NUMBER OF PARAMETERS
C      ==> NQ       NUMBER OF RESPONSES PER OBSERVATION
C     <==> BETA     FUNCTION PARAMETERS
C      ==> Y        RESPONSE VARIABLE
C      ==> LDY      LEADING DIMENSION OF ARRAY Y
C      ==> X        EXPLANATORY VARIABLE
C      ==> LDX      LEADING DIMENSION OF ARRAY X
C      ==> WE       "EPSILON" WEIGHTS
C      ==> LDWE     LEADING DIMENSION OF ARRAY WE
C      ==> LD2WE    SECOND DIMENSION OF ARRAY WE
C      ==> WD       "DELTA" WEIGHTS
C      ==> LDWD     LEADING DIMENSION OF ARRAY WD
C      ==> LD2WD    SECOND DIMENSION OF ARRAY WD
C      ==> IFIXB    INDICATORS FOR "FIXING" PARAMETERS (BETA)
C      ==> IFIXX    INDICATORS FOR "FIXING" EXPLANATORY VARIABLE (X)
C      ==> LDIFX    LEADING DIMENSION OF ARRAY IFIXX
C      ==> JOB      TASK TO BE PERFORMED 
C      ==> NDIGIT   GOOD DIGITS IN SUBROUTINE FUNCTION RESULTS
C      ==> TAUFAC   TRUST REGION INITIALIZATION FACTOR
C      ==> SSTOL    SUM OF SQUARES CONVERGENCE CRITERION
C      ==> PARTOL   PARAMETER CONVERGENCE CRITERION
C      ==> MAXIT    MAXIMUM NUMBER OF ITERATIONS
C      ==> IPRINT   PRINT CONTROL 
C      ==> LUNERR   LOGICAL UNIT FOR ERROR REPORTS 
C      ==> LUNRPT   LOGICAL UNIT FOR COMPUTATION REPORTS 
C      ==> STPB     STEP SIZES FOR FINITE DIFFERENCE DERIVATIVES WRT BETA
C      ==> STPD     STEP SIZES FOR FINITE DIFFERENCE DERIVATIVES WRT DELTA
C      ==> LDSTPD   LEADING DIMENSION OF ARRAY STPD
C      ==> SCLB     SCALE VALUES FOR PARAMETERS BETA
C      ==> SCLD     SCALE VALUES FOR ERRORS DELTA IN EXPLANATORY VARIABLE 
C      ==> LDSCLD   LEADING DIMENSION OF ARRAY SCLD
C     <==> WORK     DOUBLE PRECISION WORK VECTOR
C      ==> LWORK    DIMENSION OF VECTOR WORK
C     <==  IWORK    INTEGER WORK VECTOR
C      ==> LIWORK   DIMENSION OF VECTOR IWORK
C     <==  INFO     STOPPING CONDITION 

C  PARAMETER DECLARATIONS AND SPECIFICATIONS
      INTEGER    LDIFX,LDSCLD,LDSTPD,LDWD,LDWE,LDX,LDY,LD2WD,LD2WE,
     +           LIWORK,LWORK
C  VARIABLE DECLARATIONS 
      INTEGER          I,INFO,IPRINT,J,JOB,L,LUNERR,LUNRPT,M,MAXIT,N,
     +                 NDIGIT,NP,NQ
      INTEGER          IFIXB(NP),IFIXX(N,M),IWORK(LIWORK)
      DOUBLE PRECISION PARTOL,SSTOL,TAUFAC
      DOUBLE PRECISION BETA(NP),SCLB(NP),SCLD(1,M),
     +                 STPB(NP),STPD(1,M),
     +                 WD(N,1,M),WE(N,1,NQ),
     +                 WORK(LWORK),X(N,M),Y(N,NQ)
!------------For using information in WORK----------------------------
      LOGICAL
     +ISODR
      INTEGER
     +   DELTAI,EPSI,XPLUSI,FNI,SDI,VCVI,
     +   RVARI,WSSI,WSSDEI,WSSEPI,RCONDI,ETAI,
     +   OLMAVI,TAUI,ALPHAI,ACTRSI,PNORMI,RNORSI,PRERSI,
     +   PARTLI,SSTOLI,TAUFCI,EPSMAI,
     +   BETA0I,BETACI,BETASI,BETANI,SI,SSI,SSFI,QRAUXI,UI,
     +   FSI,FJACBI,WE1I,DIFFI,
     +   DELTSI,DELTNI,TI,TTI,OMEGAI,FJACDI,
     +   WRK1I,WRK2I,WRK3I,WRK4I,WRK5I,WRK6I,WRK7I,
     +   LWKMN
c
      integer i1,i2,i3,i4,i5,iderivative
      double precision weity(N,NQ),weitx(N,M),shorty(N,NQ),
     &shortx(N,M),fvalue
      EXTERNAL FCN
c
C  SPECIFY DEFAULT VALUES FOR DODRC ARGUMENTS

      LDY=N
      LDX=N
      LDWE=N
      LD2WE=1
      LDWD=N
      LD2WD=1
      LDIFX=N
      LDSTPD=1
      LDSCLD=1

      WE(1,1,1)  = -1.0D0
      WD(1,1,1)  = -1.0D0
      IFIXB(1)   = -1
!      IFIXX(1,1) = -1

      if(INFO.eq.0)then
!explicit ordinary least square fitting
       ISODR=.false.
       if(iderivative.eq.0)then
!no derivatives provided, using central finite difference
         JOB=13
       else
!don't check derivatives
         JOB=43
!check derivatives
!         JOB=23
        endif
      endif
      if(INFO.eq.1)then
!explicit orthogonal distance regression
        ISODR=.true.
        if(iderivative.eq.0)then
!no derivatives provided, using central finite difference
          JOB=10
        else
!don't check derivatives
          JOB=40
!check derivatives
!         JOB=20
        endif
      endif
      if(INFO.eq.-1)then
!implicit orthogonal distance regression
        ISODR=.true.
        if(iderivative.eq.0)then
!no derivatives provided, using central finite difference
          JOB=11
        else
!don't check derivatives
          JOB=31
!check derivatives
!          JOB=21
        endif
      endif
      NDIGIT     = -1
      TAUFAC     = -1.0D0
      SSTOL      = -1.0D0
      PARTOL     = -1.0D0
      MAXIT      = -1
      IPRINT     = -1
      IPRINT=0
      LUNERR     = -1
      LUNRPT     = -1
      STPB(1)    = -1.0D0
      STPD(1,1)  = -1.0D0
      SCLB(1)    = -1.0D0
      SCLD(1,1)  = -1.0D0

      MAXIT      = 200000
C  SET UP ODRPACK REPORT FILES
      LUNERR  =   107
      LUNRPT  =   108
      LWKMN=LWORK
c
      do I=1,N
        do i1=1,M
          WD(I,1,i1)=weitx(I,i1)
        enddo
        do i1=1,NQ
          WE(I,1,i1)=weity(I,i1)
        enddo
      enddo
C  READ PROBLEM DATA, AND SET NONDEFAULT VALUE FOR ARGUMENT IFIXX
      DO 10 I=1,N
         DO 15 J=1, M
           IFIXX(I,J) = 1
15       CONTINUE
10    CONTINUE
60    CALL DODRC(FCN,
     +           N,M,NP,NQ,
     +           BETA,
     +           Y,LDY,X,LDX,
     +           WE,LDWE,LD2WE,WD,LDWD,LD2WD,
     +           IFIXB,IFIXX,LDIFX,
     +           JOB,NDIGIT,TAUFAC,
     +           SSTOL,PARTOL,MAXIT,
     +           IPRINT,LUNERR,LUNRPT,
     +           STPB,STPD,LDSTPD,
     +           SCLB,SCLD,LDSCLD,
     +           WORK,LWORK,IWORK,LIWORK,
     +           INFO)
      i1=mod(INFO,10)
      i2=(mod(INFO,100)-i1)/10
      i3=(mod(INFO,1000)-mod(INFO,100))/100
      i4=(mod(INFO,10000)-mod(INFO,1000))/1000
      i5=(INFO-mod(INFO,10000))/10000
      CALL DWINF
     +   (N,M,NP,NQ,LDWE,LD2WE,ISODR,
     +   DELTAI,EPSI,XPLUSI,FNI,SDI,VCVI,
     +   RVARI,WSSI,WSSDEI,WSSEPI,RCONDI,ETAI,
     +   OLMAVI,TAUI,ALPHAI,ACTRSI,PNORMI,RNORSI,PRERSI,
     +   PARTLI,SSTOLI,TAUFCI,EPSMAI,
     +   BETA0I,BETACI,BETASI,BETANI,SI,SSI,SSFI,QRAUXI,UI,
     +   FSI,FJACBI,WE1I,DIFFI,
     +   DELTSI,DELTNI,TI,TTI,OMEGAI,FJACDI,
     +   WRK1I,WRK2I,WRK3I,WRK4I,WRK5I,WRK6I,WRK7I,
     +   LWKMN)
      fvalue=0.0d0
      do I=1,N
        do J=1,M
          shortx(I,J)=WORK(XPLUSI-1+I+(J-1)*N)
          fvalue=fvalue+weitx(I,J)*WORK(DELTAI-1+I+(J-1)*N)
     +      *WORK(DELTAI-1+I+(J-1)*N)
        enddo
        do J=1,NQ
          shorty(I,J)=WORK(FNI-1+I+(J-1)*N)
          fvalue=fvalue+weity(I,J)*WORK(EPSI-1+I+(J-1)*N)
     +*WORK(EPSI-1+I+(J-1)*N)
        enddo
      enddo
      return
      END
