249 lines
8.5 KiB
FortranFixed
249 lines
8.5 KiB
FortranFixed
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
|
|
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
|