Initial commit

This commit is contained in:
2016-02-03 18:52:05 +00:00
commit d40505e161
507 changed files with 91383 additions and 0 deletions
+32
View File
@@ -0,0 +1,32 @@
# Compiled Object files
*.slo
*.lo
*.o
*.obj
# Precompiled Headers
*.gch
*.pch
# Compiled Dynamic libraries
*.so
*.dylib
*.dll
# Fortran module files
*.mod
# Compiled Static libraries
*.lai
*.la
*.a
*.lib
# Executables
*.exe
*.out
*.app
# editor and mac files
*~
.DS_Store
File diff suppressed because it is too large Load Diff
+520
View File
@@ -0,0 +1,520 @@
C[BA*)
C[KA{F 5}
C[ {Iterative Methods for Linear Systems}
C[ {Iterative Methods for Linear Systems}*)
C[FE{F 5.4}
C[ {The Gau"s-Seidel Iteration}
C[ {The Gau"s-Seidel Iteration}*)
C[LE*)
c SUBROUTINE ADSOR(A,N,IA,B,X,KADAPT,EPS,KMAX,IMETH,ISWITC,
C[IX{ADSOR}*)
c * OMEGA,WORK,RES,ITNUMB,IERR)
SUBROUTINE ADSOR(A,N,IA,B,X,IERR)
C
C*****************************************************************
C *
C This program solves an inhomogeneous linear system AX = B of *
C equations with a nonsingular system matrix A. The method of *
C Jacobi is used jointly with relaxation, where the relaxation *
C parameter OMEGA is adjusted during the iteration (adaptive *
C SOR method). *
C[BE*)
C For a suitable choice of parameters (refer to the remark *
C below), this program can perform the Gauá-Seidel method or *
C a non-adaptive SOR method. *
C *
C *
C INPUT PARAMETERS: *
C ================= *
C A : 2-dimensional array A(1:IA,1:N), containing the *
C system matrix for the linear equations *
C N : size of the linear system *
C IA : leading dimension of A, as specified in the calling *
C program *
C B : N-vector B(1:N), the right hand side of the system *
C X : N-vector X(1:N) containing the starting value for *
C iteration *
C KADAPT : Number of iterations, after which the relaxation *
C parameter is to be redefined *
C EPS : desired accuracy; the iteration is stopped when the *
C maximum norm of the relative error does not exceed *
C EPS *
C KMAX : Maximal number of iterations allowed *
C IMETH : parameter that determines the method used: *
C = 0, adaptive SOR method *
C = 1, SOR method for a given relaxation parameter *
C = 2, Gauá-Seidel method *
C ISWITC : parameter that determines the convergence criterion *
C to be used: *
C = 0, none *
C = 1, row sum criterion *
C = 2, column sum criterion *
C = 3, criterion of Schmidt and v. Mises *
C OMEGA : in case IMETH=1, the optimal relaxation parameter *
C must be part of the input; otherwise only the name *
C must be declared in the callimng program. *
C *
C *
C REMARKS: *
C ======== *
C For the adaptive SOR method (IMETH=0) we recommend to set *
C KADAPT=4 or KADAPT=5. *
C If the optimal relaxationcoefficient Wopt is known for A, then*
C one should set IMETH=1 and OMEGA = Wopt, i.e., the SOR method *
C with given optimal relaxation coefficient should be used. *
C If IMETH=2, then the program performs the Gauá-Seidel method. *
C *
C *
C AUXILIARY PARAMETERS: *
C ===================== *
C WORK : 2-dim. array WORK(1:N,1:3) *
C *
C *
C OUTPUT PARAMETERS: *
C ================== *
C A : 2-dim. array A(1:IA,1:N), the input matrix A is over-*
C written by: A(I,J)=A(I,J)/A(I,I) for I,J=1, ..., N *
C B : N-vector B(1:N), the right hand side is replaced by *
C B(I)=B(I)/A(I,I); I=1,N *
C OMEGA : - if IMETH = 0, the program returns the adaptively *
C computed relaxations parameter. *
C - if IMETH = 1, the optimal relaxation parameter *
C is returned as put in externally. *
C - if IMETH = 2, then on output OMEGA = 1. *
C X : N-vector X(1:N) that contains the solution vector *
C RES : N-vector RES(1:N) containing the residuum B - AX; *
C the residuum is available even if the desired *
C accuracy EPS could not be achieved with the given *
C maximum number of iterations. *
C ITNUMB : num,bert of iterations actually performed *
C IERR : error parameter: *
C = 0, the desired convergence criterium has not been *
C met *
C = 1, the solution X has been found *
C = 2, the desired accuracy has not been achieved after*
C KMAX iterations *
C = 3, input data incorrect *
C = 4, system matrix A is numerically singular *
C *
C----------------------------------------------------------------*
C *
C Required subroutines: GAUSEI, MNORM, CONV, RESID, MACHPD *
C *
C*****************************************************************
C *
C Author : Gisela Engeln-Mllges *
C Date : 06.09.1992 *
C Source : FORTRAN 77 *
C *
C[BA*)
C*****************************************************************
C[BE*)
C
C Declarations
C
DOUBLE PRECISION A(1:IA,1:N),B(1:N),X(1:N),WORK(1:N,1:3),
* RES(1:N),EPS,OMEGA,FMACHP,HELP,DIFFN,Q,
* RELERR,SUM,XN
C
c The following 5 lines is added by GU
EPS=1.0D-06
KADAPT=4
KMAX=2000
IMETH=2
ISWITC=0
OMEGA=1.0d0
c
C Checking the inputs EPS, KMAX, IMETH and ISWITC
C
IF(EPS .LE. 0.0D0 .OR. KMAX .LT. 1 .OR. ISWITC .LT. 0 .OR.
* ISWITC .GT. 3 .OR. IMETH .LT. 0 .OR. IMETH .GT. 2) THEN
IERR=3
RETURN
ENDIF
C
C Initialize the parameters KADAPT and OMEGA depending on the method
C
IF(IMETH .EQ. 0) THEN
OMEGA=1.0D0
ELSE IF(IMETH .EQ. 1) THEN
KADAPT=KMAX
ELSE IF(IMETH .EQ. 2) THEN
KADAPT=KMAX
OMEGA=1.0D0
ENDIF
C
C Compute the machine constant and initialize the relative error bound
C
FMACHP=1.0D0
10 FMACHP=0.5D0*FMACHP
IF(MACHPD(1.0D0+FMACHP) .EQ. 1) GOTO 10
RELERR=FMACHP*8.0D0
C
C Initialize
C
Q=1.0D0
ITNUMB=0
C
C Check whether A is singular; if so, set IERR = 4.
C
DO 20 I=1,N
SUM=DABS(A(I,1))
DO 30 K=2,N
SUM=SUM+DABS(A(I,K))
30 CONTINUE
IF(SUM .EQ. 0.0D0) THEN
IERR=4
RETURN
ELSE IF(DABS(A(I,I))/SUM .LT. RELERR) THEN
IERR=4
RETURN
ENDIF
20 CONTINUE
C
C Redefine the entries in A and B: A(I,J) := A(I,J)/A(I,I)
C and B(I) := B(I)/A(I,I) .
C
DO 40 I=1,N
HELP=1.0D0/A(I,I)
DO 50 J=1,N
A(I,J)=A(I,J)*HELP
50 CONTINUE
B(I)=B(I)*HELP
40 CONTINUE
C
C Check for convergence
C
IF(ISWITC .NE. 0) THEN
CALL CONV(ISWITC,A,N,IA,IERR)
IF(IERR .EQ. 0) RETURN
ENDIF
C
C The vector RES serves as auxiliary storage for the previous solution
C vektor. Initially RES contains the staring vector.
C
DO 60 I=1,N
RES(I)=X(I)
60 CONTINUE
C
C One iteration with the Gauá-Seidel method gives the first iterate X
C
CALL GAUSEI(A,N,IA,B,OMEGA,X)
C
C Up the iteration counter
C
ITNUMB=ITNUMB+1
C
C Compute the difference of the last two iterates
C
DO 70 I=1,N
WORK(I,1)=X(I)-RES(I)
70 CONTINUE
C
C Iteration loop for the chosen method
C
DO 80 K=1,KMAX-1
C
C Check break-off criterion
C
CALL MNORM(WORK(1,1),N,DIFFN)
CALL MNORM(X,N,XN)
IF(DIFFN .LE. EPS*XN) THEN
IERR=1
ITNUMB=K
CALL RESID(A,N,IA,B,X,RES)
RETURN
ENDIF
IF(K .EQ. KMAX-1) THEN
ITNUMB=KMAX
IERR=2
CALL RESID(A,N,IA,B,X,RES)
RETURN
ENDIF
C
C RES contains the previous iterate
C
DO 90 I=1,N
RES(I)=X(I)
90 CONTINUE
C
C One iteration step using Gauá-Seidel for a fixed OMEGA
C
CALL GAUSEI(A,N,IA,B,OMEGA,X)
C
C Compute the difference of the last two iterates
C
DO 100 I=1,N
WORK(I,2)=X(I)-RES(I)
100 CONTINUE
C
C If the number of performed iterations K is divisible by KADAPT,
C then we compute Q in order to adjust the relaxation parameter;
C Q is an estimate of the spectral radius of the iteration matrix.
C
IF(MOD(K,KADAPT) .EQ. 0) THEN
DO 110 I=1,N
IF(DABS(WORK(I,1)) .LT. FMACHP) THEN
WORK(I,3)=1.0D0
ELSE
WORK(I,3)=WORK(I,2)/WORK(I,1)
ENDIF
110 CONTINUE
CALL MNORM(WORK(1,3),N,Q)
C
C If Q > 1, then the iteration counter is upped by one and
C the next Gauá-Seidel step is executed; otherwise a new
C relaxation parameter is calculated.
C
IF(Q .LE. 1.0D0) THEN
Q=MAX(Q,OMEGA-1.0D0)
OMEGA=2.0D0/(1.0D0+DSQRT(1.0D0-((Q+OMEGA-1.0D0)
* /OMEGA)**2/Q))
ENDIF
ENDIF
C
C The difference vector of the last two iterations is replaced
C by the one of the previous two iterations for the approximate solution
C
DO 120 I=1,N
WORK(I,1)=WORK(I,2)
120 CONTINUE
80 CONTINUE
END
C
C
C[BA*)
C[LE*)
SUBROUTINE GAUSEI(A,N,IA,B,OMEGA,X)
C[IX{GAUSEI}*)
C
C*****************************************************************
C *
C This subroutine performs one iteration with the Gauá-Seidel *
C method for a given relaxation parameter. *
C[BE*)
C *
C *
C INPUT PARAMETERS: *
C ================= *
C A : 2-dim. array A(1:IA, 1:N), that contains the *
C modified system matrix A : A(I,J)=A(I,J)/A(I,I) for *
C I,J=1, ..., N *
C N : order of the system *
C IA : leading dimension of A, as specified in the calling *
C program *
C B : N-vector B(1:N) with the modified right hand side: *
C B(I)=B(I)/A(I,I); I=1, ..., N *
C OMEGA : relaxation parameter *
C X : N-vector X(1:N) containing the starting vector for *
C the iteration *
C *
C *
C OUTPUT PARAMETERS: *
C ================== *
C X : N-vector X(1:N) containing the next iteration vector *
C *
C----------------------------------------------------------------*
C *
C Required subroutines: none *
C *
C*****************************************************************
C *
C Author : Gisela Engeln-Mllges *
C Date : 06.09.1992 *
C Source : FORTRAN 77 *
C *
C[BA*)
C*****************************************************************
C[BE*)
C
DOUBLE PRECISION A(1:IA,1:N),B(1:N),X(1:N),OMEGA,S
C
DO 10 I=1,N
S=B(I)
DO 20 J=1,N
S=S-A(I,J)*X(J)
20 CONTINUE
X(I)=X(I)+OMEGA*S
10 CONTINUE
RETURN
END
C
C
C[BA*)
C[LE*)
SUBROUTINE MNORM(X,N,XNORM)
C[IX{MNORM}*)
C
C*****************************************************************
C *
C This subroutine calculates the maximum norm XNORM of an *
C N-vector X. *
C *
C----------------------------------------------------------------*
C[BE*)
C *
C Required subroutines: none *
C *
C*****************************************************************
C *
C Author : Gisela Engeln-Mllges *
C Date : 06.09.1992 *
C Source : FORTRAN 77 *
C *
C*****************************************************************
C
DOUBLE PRECISION X(1:N),XNORM
C
XNORM=DABS(X(1))
DO 10 I=2,N
XNORM=DMAX1(XNORM,DABS(X(I)))
10 CONTINUE
RETURN
END
C
C
C[BA*)
C[LE*)
SUBROUTINE CONV(ISWITC,A,N,IA,IERR)
C[IX{CONV}*)
C
C*****************************************************************
C *
C This subroutine helps check convergence. *
C[BE*)
C *
C *
C INPUT PARAMETERS: *
C ================= *
C ISWITC : Parameter that determines the convergence criterion *
C to be checked: *
C = 0, none *
C = 1, row sum criterion *
C = 2, column sum criterion *
C = 3, criterion of Schmidt and v. Mises *
C A : 2-dim. array A(1:IA, 1:N), containing the matrix for *
C which we want to check convergence of the iterates *
C from the various SOR algorithms *
C N : order of the matrix A *
C IA : leading dimension of A, as prescribed in the calling *
C program *
C *
C *
C OUTPUT PARAMETERS: *
C ================== *
C IERR : error parameter: *
C = 0, the desired convergence criterion has not been *
C met *
C = 1, the desired criterion is satified *
C *
C----------------------------------------------------------------*
C *
C Required subroutines: none *
C *
C*****************************************************************
C *
C Author : Gisela Engeln-Mllges *
C Date : 06.09.1992 *
C Source : FORTRAN 77 *
C *
C[BA*)
C*****************************************************************
C[BE*)
C
DOUBLE PRECISION A(1:IA,1:N),SUM
C
C Row sum criterion
C
IF(ISWITC .EQ. 1) THEN
DO 10 I=1,N
SUM=-1.0D0
DO 20 J=1,N
SUM=SUM+DABS(A(I,J))
20 CONTINUE
IF(SUM .LT. 1.0D0) THEN
IERR=1
ELSE
IERR=0
RETURN
ENDIF
10 CONTINUE
C
C Column sum criterion
C
ELSE IF(ISWITC .EQ. 2) THEN
DO 30 J=1,N
SUM=-1.0D0
DO 40 I=1,N
SUM=SUM+DABS(A(I,J))
40 CONTINUE
IF(SUM .LT. 1.0D0) THEN
IERR=1
ELSE
IERR=0
RETURN
ENDIF
30 CONTINUE
C
C Criterion of Schmidt and v. Mises
C
ELSE IF(ISWITC .EQ. 3) THEN
SUM=-N
DO 50 I=1,N
DO 60 J=1,N
SUM=SUM+A(I,J)*A(I,J)
60 CONTINUE
50 CONTINUE
SUM=DSQRT(SUM)
IF(SUM .LT. 1.0D0) THEN
IERR=1
ELSE
IERR=0
RETURN
ENDIF
ENDIF
END
C
C
C[BA*)
C[LE*)
SUBROUTINE RESID(A,N,IA,B,X,RES)
C[IX{RESID}*)
C
C*****************************************************************
C *
C This subroutine computes the residuum RES = B - AX, where *
C both A and B are given in modified form. *
C *
C----------------------------------------------------------------*
C[BE*)
C *
C Required subroutines: none *
C *
C*****************************************************************
C *
C Author : Gisela Engeln-Mllges *
C Date : 09.06.1992 *
C Source : FORTRAN 77 *
C *
C*****************************************************************
C
DOUBLE PRECISION A(1:IA,1:N), B(1:N), X(1:N), RES(1:N),DSUM
C
DO 10 I=1,N
DSUM=B(I)
DO 20 J=1,N
DSUM=DSUM-A(I,J)*X(J)
20 CONTINUE
RES(I)=DSUM
10 CONTINUE
RETURN
END
c
C[KA{F 0}{Auxiliary Library}{Auxiliary Library}*)
INTEGER FUNCTION MACHPD(X)
C[IX{MACHPD}*)
DOUBLE PRECISION X
MACHPD=0
IF (1.0D0 .LT. X) MACHPD=1
RETURN
END
+51
View File
@@ -0,0 +1,51 @@
subroutine eigen_sym_up(N,A,W)
implicit none
!
!compute the eigenvalues and eigenvectors of a symmetrical matrix.
!A: On entry, A is a symmetrical matrix with its upper triangle filled.
! on exit, A contains the normalized eigenvectors in its columns.
!W: contains the eigenvalues in descending order.
CHARACTER JOBZ, UPLO
INTEGER INFO, LDA, LWORK, N
DOUBLE PRECISION A(N, N), W( N ), WORK(3*N-1)
double precision p
integer i,j
JOBZ='V'
UPLO='U'
LWORK=3*N-1
LDA=N
call DSYEV(JOBZ,UPLO,N,A(1:N,1:N),LDA,W,WORK,LWORK,INFO)
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, the algorithm failed to converge; i
* off-diagonal elements of an intermediate tridiagonal
* form did not converge to zero.
if(INFO.lt.0)then
write(*,*)'The ',-INFO,
& 'th argument in DSYEV has an illegal value'
stop
endif
if(INFO.gt.0)then
write(*,*)'The algorithm failed to converge'
stop
endif
! Change the eigenvalue array from ascending to descending order and rearrange
! the eigen vectors accordingly.
!---------------------------------------------
do i=1,N/2
p=W(i)
W(i)=W(N-i+1)
W(N-i+1)=p
do j=1,N
p=A(j,i)
A(j,i)=A(j,N-i+1)
A(j,N-i+1)=p
enddo
enddo
!---------------------------------------------
return
end
+201
View File
@@ -0,0 +1,201 @@
SUBROUTINE lfit(x,y,sig,ndat,a,ma,npc,funcs,ierr)
implicit none
! ierr = 1, ok; =0 singular matrix
INTEGER ma,ia(ma),npc,ndat,MMAX,ierr
double precision chisq,a(ma),covar(npc,npc),
& sig(ndat),x(ndat),y(ndat)
EXTERNAL funcs
PARAMETER (MMAX=10000)
CU USES covsrt,gaussj
INTEGER i,j,k,l,m,mfit
double precision sig2i,sum,wt,ym,afunc(MMAX),beta(MMAX)
mfit=0
ierr=1
do 11 j=1,ma
ia(j)=1
if(ia(j).ne.0) mfit=mfit+1
11 continue
if(mfit.eq.0) pause 'lfit: no parameters to be fitted'
do 13 j=1,mfit
do 12 k=1,mfit
covar(j,k)=0.0d0
12 continue
beta(j)=0.0d0
13 continue
do 17 i=1,ndat
call funcs(x(i),afunc,ma)
ym=y(i)
if(mfit.lt.ma) then
do 14 j=1,ma
if(ia(j).eq.0) ym=ym-a(j)*afunc(j)
14 continue
endif
sig2i=1.0d0/sig(i)**2
j=0
do 16 l=1,ma
if (ia(l).ne.0) then
j=j+1
wt=afunc(l)*sig2i
k=0
do 15 m=1,l
if (ia(m).ne.0) then
k=k+1
covar(j,k)=covar(j,k)+wt*afunc(m)
endif
15 continue
beta(j)=beta(j)+ym*wt
endif
16 continue
17 continue
do 19 j=2,mfit
do 18 k=1,j-1
covar(k,j)=covar(j,k)
18 continue
19 continue
call gaussj(covar,mfit,npc,beta,1,1,ierr)
if(ierr.eq.0)then
! singular matrix
return
endif
j=0
do 21 l=1,ma
if(ia(l).ne.0) then
j=j+1
a(l)=beta(j)
endif
21 continue
chisq=0.0d0
do 23 i=1,ndat
call funcs(x(i),afunc,ma)
sum=0.0d0
do 22 j=1,ma
sum=sum+a(j)*afunc(j)
22 continue
chisq=chisq+((y(i)-sum)/sig(i))**2
23 continue
call covsrt(covar,npc,ma,ia,mfit)
return
END
C (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.
SUBROUTINE covsrt(covar,npc,ma,ia,mfit)
implicit none
INTEGER ma,mfit,npc,ia(ma)
double precision covar(npc,npc)
INTEGER i,j,k
double precision swap
do 12 i=mfit+1,ma
do 11 j=1,i
covar(i,j)=0.0d0
covar(j,i)=0.0d0
11 continue
12 continue
k=mfit
do 15 j=ma,1,-1
if(ia(j).ne.0)then
do 13 i=1,ma
swap=covar(i,k)
covar(i,k)=covar(i,j)
covar(i,j)=swap
13 continue
do 14 i=1,ma
swap=covar(k,i)
covar(k,i)=covar(j,i)
covar(j,i)=swap
14 continue
k=k-1
endif
15 continue
return
END
C (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.
SUBROUTINE gaussj(a,n,np,b,m,mp,ierr)
implicit none
INTEGER m,mp,n,np,NMAX,ierr
double precision a(np,np),b(np,mp)
PARAMETER (NMAX=10000)
INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),
& ipiv(NMAX)
double precision big,dum,pivinv
ierr=1
do 11 j=1,n
ipiv(j)=0
11 continue
do 22 i=1,n
big=0.0d0
do 13 j=1,n
if(ipiv(j).ne.1)then
do 12 k=1,n
if (ipiv(k).eq.0) then
if (dabs(a(j,k)).ge.big)then
big=dabs(a(j,k))
irow=j
icol=k
endif
else if (ipiv(k).gt.1) then
! pause 'singular matrix in gaussj'
ierr=0
return
endif
12 continue
endif
13 continue
ipiv(icol)=ipiv(icol)+1
if (irow.ne.icol) then
do 14 l=1,n
dum=a(irow,l)
a(irow,l)=a(icol,l)
a(icol,l)=dum
14 continue
do 15 l=1,m
dum=b(irow,l)
b(irow,l)=b(icol,l)
b(icol,l)=dum
15 continue
endif
indxr(i)=irow
indxc(i)=icol
if (a(icol,icol).eq.0.0d0)then
! pause 'singular matrix in gaussj'
ierr=0
return
endif
pivinv=1.0d0/a(icol,icol)
a(icol,icol)=1.0d0
do 16 l=1,n
a(icol,l)=a(icol,l)*pivinv
16 continue
do 17 l=1,m
b(icol,l)=b(icol,l)*pivinv
17 continue
do 21 ll=1,n
if(ll.ne.icol)then
dum=a(ll,icol)
a(ll,icol)=0.0d0
do 18 l=1,n
a(ll,l)=a(ll,l)-a(icol,l)*dum
18 continue
do 19 l=1,m
b(ll,l)=b(ll,l)-b(icol,l)*dum
19 continue
endif
21 continue
22 continue
do 24 l=n,1,-1
if(indxr(l).ne.indxc(l))then
do 23 k=1,n
dum=a(k,indxr(l))
a(k,indxr(l))=a(k,indxc(l))
a(k,indxc(l))=dum
23 continue
endif
24 continue
return
END
C (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.
File diff suppressed because it is too large Load Diff
+19
View File
@@ -0,0 +1,19 @@
subroutine matrixsquare_up(m,n,A,B)
implicit none
! compute B=AA^T. B is symmetrical so only its upper triangle is computed.
integer m,n
double precision A(m,n),B(m,m)
integer i,j,k
do i=1,m
do j=i,m
B(i,j)=0.0d0
do k=1,n
B(i,j)=B(i,j)+A(i,k)*A(j,k)
enddo
enddo
enddo
return
end
+428
View File
@@ -0,0 +1,428 @@
program main
implicit none
double precision x(1000),y(1000),dx(1000),dy(1000),
* slope,fintcpt,rtmnsquare,xoutliers(100),youtliers(1000),
* x1(1000),y1(1000),k,b,sum1,sum2
integer nsamp,i,numoutliers,nsamp1
open(unit=1,file='testdata.txt')
i=1
10 read(1,*,end=100)x(i),y(i)
i=i+1
goto 10
100 nsamp=i-1
goto 200
nsamp=13
x(1)=-1.0d0
y(1)=3.0d0
x(2)=-3.0d0
y(2)=4.0d0
x(3)=-5.0d0
y(3)=5.0d0
x(4)=-7.0d0
y(4)=6.0d0
x(5)=-9.0d0
y(5)=7.0d0
x(6)=-11.0d0
y(6)=8.0d0
x(7)=-3.0d0
y(7)=1.0d0
x(8)=-5.0d0
y(8)=2.0d0
x(9)=-7.0d0
y(9)=3.0d0
x(10)=-9.0d0
y(10)=4.0d0
x(11)=-11.0d0
y(11)=5.0d0
x(12)=-4.0d0
y(12)=7.0d0
x(13)=-12.0d0
y(13)=1.0d0
200 slope=-2.0d0
fintcpt=0.0d0
call OrthSoilRespRegres(nsamp,x,y,slope,fintcpt)
write(*,*)slope,fintcpt
pause
slope=-2.0d0
fintcpt=0.0d0
call orthlinreg_outlier(nsamp,x,y,slope,
& fintcpt,dx,dy,rtmnsquare,xoutliers,youtliers,
& numoutliers)
write(*,*)slope/2.0d0,fintcpt,numoutliers
do i=1,numoutliers
write(*,*)xoutliers(i),youtliers(i)
enddo
end
subroutine orthlinreg_outlier(nsamp0,x0,y0,slope,
& fintcpt,dx,dy,rtmnsquare,xoutliers,youtliers,
& numoutliers)
implicit none
integer nsamp0,numoutliers
double precision x0(nsamp0),y0(nsamp0),slope,
& fintcpt,dx(nsamp0),dy(nsamp0),rtmnsquare,
& xoutliers(nsamp0),youtliers(nsamp0),xtest(nsamp0),
& ytest(nsamp0),slopetest,fintcpttest,dxtest(nsamp0),
& dytest(nsamp0),rtmnsquaretest,testmeasure(nsamp0),
& x(nsamp0),y(nsamp0)
integer iwhichside,nsamptest,isitoutlier,
& isoutlier_1side,i,j,nsamp
parameter (iwhichside=1)
numoutliers=0
nsamp=nsamp0
do i=1,nsamp
x(i)=x0(i)
y(i)=y0(i)
enddo
50 call orthlinreg(nsamp,x,y,slope,fintcpt,
& dx,dy,rtmnsquare)
write(*,*)slope,fintcpt,rtmnsquare
stop
nsamptest=nsamp-1
do i=1,nsamp
do j=1,nsamp
xtest(j)=x(j)
ytest(j)=y(j)
enddo
xtest(i)=x(nsamp)
ytest(i)=y(nsamp)
call orthlinreg(nsamptest,xtest,ytest,slopetest,
& fintcpttest,dxtest,dytest,rtmnsquaretest)
! write(*,*)i,slopetest,fintcpttest
! testmeasure(i)=(slopetest-slope)**2+
! & (fintcpttest-fintcpt)**2
testmeasure(i)=100.0d0*dabs(rtmnsquaretest-rtmnsquare)/
& rtmnsquare
! write(*,*)i,testmeasure(i)
enddo
isitoutlier=isoutlier_1side(nsamp,testmeasure,iwhichside)
if(isitoutlier.lt.1.or.isitoutlier.gt.nsamp)return
! outlier detected
numoutliers=numoutliers+1
xoutliers(numoutliers)=x(isitoutlier)
youtliers(numoutliers)=y(isitoutlier)
x(isitoutlier)=x(nsamp)
y(isitoutlier)=y(nsamp)
nsamp=nsamp-1
if(nsamp.le.2)then
write(*,*)'No enough good data left'
stop
endif
goto 50
return
end
! orthogonal linear regression
subroutine orthlinreg(nsamp,x,y,slope0,fintcpt0,
& dx,dy,rtmnsquare)
implicit none
integer nsamp
double precision x(nsamp),y(nsamp),dx1(nsamp),
& dy1(nsamp),slope(2),fintcpt(2),dx2(nsamp),
& dy2(nsamp),slope0,fintcpt0,dx(nsamp),dy(nsamp)
integer i,j
double precision w,u,v,xbar,ybar,root1,root2,
& a,b,c,rtmnsquare1,rtmnsquare2,rtmnsquare
xbar=0.0d0
ybar=0.0d0
w=0.0d0
u=0.0d0
v=0.0d0
do i=1,nsamp
xbar=xbar+x(i)
ybar=ybar+y(i)
w=w+x(i)*x(i)
u=u+y(i)*y(i)
v=v+x(i)*y(i)
enddo
xbar=xbar/dble(nsamp)
ybar=ybar/dble(nsamp)
w=w/dble(nsamp)
u=u/dble(nsamp)
v=v/dble(nsamp)
a=v-xbar*ybar
b=w-u-xbar*xbar+ybar*ybar
c=xbar*ybar-v
call quadraticroots(a,b,c,root1,root2)
slope(1)=root1
slope(2)=root2
fintcpt(1)=ybar-slope(1)*xbar
fintcpt(2)=ybar-slope(2)*xbar
rtmnsquare1=0.0d0
rtmnsquare2=0.0d0
do i=1,nsamp
dx1(i)=(y(i)-fintcpt(1)-x(i)*slope(1))*
& slope(1)/(1.0d0+slope(1)*slope(1))
dy1(i)=-(y(i)-fintcpt(1)-x(i)*slope(1))/
& (1.0d0+slope(1)*slope(1))
rtmnsquare1=rtmnsquare1+dx1(i)**2+dy1(i)**2
dx2(i)=(y(i)-fintcpt(2)-x(i)*slope(2))*
& slope(2)/(1.0d0+slope(2)*slope(2))
dy2(i)=-(y(i)-fintcpt(2)-x(i)*slope(2))/
& (1.0d0+slope(2)*slope(2))
rtmnsquare2=rtmnsquare2+dx2(i)**2+dy2(i)**2
enddo
rtmnsquare1=dsqrt(rtmnsquare1/dble(nsamp))
rtmnsquare2=dsqrt(rtmnsquare2/dble(nsamp))
if(rtmnsquare1.gt.rtmnsquare2)then
rtmnsquare=rtmnsquare2
slope0=slope(2)
fintcpt0=fintcpt(2)
do i=1,nsamp
dx(i)=dx2(i)
dy(i)=dy2(i)
enddo
else
rtmnsquare=rtmnsquare1
slope0=slope(1)
fintcpt0=fintcpt(1)
do i=1,nsamp
dx(i)=dx1(i)
dy(i)=dy1(i)
enddo
endif
return
end
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!
subroutine OrthSoilRespRegres(npoints,x0,y0,slope,fintcpt)
implicit none
c
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 PARAMETERS SPECIFYING MAXIMUM PROBLEM SIZES HANDLED BY THIS DRIVER
C MAXN MAXIMUM NUMBER OF OBSERVATIONS
C MAXM MAXIMUM NUMBER OF COLUMNS IN EXPLANATORY VARIABLE
C MAXNP MAXIMUM NUMBER OF FUNCTION PARAMETERS
C MAXNQ MAXIMUM NUMBER OF RESPONSES PER OBSERVATION
C PARAMETER DECLARATIONS AND SPECIFICATIONS
INTEGER LDIFX,LDSCLD,LDSTPD,LDWD,LDWE,LDX,LDY,LD2WD,LD2WE,
+ LIWORK,LWORK,MAXM,MAXN,MAXNP,MAXNQ
PARAMETER (MAXM=25,MAXN=50000,MAXNP=30,MAXNQ=1,
+ LDY=MAXN,LDX=MAXN,
+ LDWE=1,LD2WE=1,LDWD=1,LD2WD=1,
+ LDIFX=MAXN,LDSTPD=1,LDSCLD=1,
+ LWORK=18 + 11*MAXNP + MAXNP**2 + MAXM + MAXM**2 +
+ 4*MAXN*MAXNQ + 6*MAXN*MAXM + 2*MAXN*MAXNQ*MAXNP +
+ 2*MAXN*MAXNQ*MAXM + MAXNQ**2 +
+ 5*MAXNQ + MAXNQ*(MAXNP+MAXM) + LDWE*LD2WE*MAXNQ,
+ LIWORK=20+MAXNP+MAXNQ*(MAXNP+MAXM))
C VARIABLE DECLARATIONS
INTEGER I,INFO,IPRINT,J,JOB,L,LUNERR,LUNRPT,M,MAXIT,N,
+ NDIGIT,NP,NQ
INTEGER IFIXB(MAXNP),IFIXX(LDIFX,MAXM),IWORK(LIWORK)
DOUBLE PRECISION PARTOL,SSTOL,TAUFAC
DOUBLE PRECISION BETA(MAXNP),SCLB(MAXNP),SCLD(LDSCLD,MAXM),
+ STPB(MAXNP),STPD(LDSTPD,MAXM),
+ WD(LDWD,LD2WD,MAXM),WE(LDWE,LD2WE,MAXNQ),
+ WORK(LWORK),X(LDX,MAXM),Y(LDY,MAXNQ)
c
integer npoints,i1
double precision x0(npoints),y0(npoints),slope,fintcpt
EXTERNAL OrthRespFCN
c
C SPECIFY DEFAULT VALUES FOR DODRC ARGUMENTS
WE(1,1,1) = -1.0D0
WD(1,1,1) = -1.0D0
IFIXB(1) = -1
! IFIXX(1,1) = -1
! JOB = 00023
JOB=20
NDIGIT = -1
TAUFAC = -1.0D0
SSTOL = -1.0D0
PARTOL = -1.0D0
MAXIT = -1
! IPRINT = -1
! IPRINT=0
IPRINT=-1
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 = 9
LUNRPT = 9
c
N=npoints
M=1
NP=2
NQ=1
do I=1,N
do i1=1,M
X(I,i1)=x0(I)
enddo
Y(I,1)=y0(I)
enddo
BETA(1)=slope
BETA(2)=fintcpt
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(OrthRespFCN,
+ 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)
slope=BETA(1)
fintcpt=BETA(2)
return
END
c
SUBROUTINE OrthRespFCN(N,M,NP,NQ,
+ LDN,LDM,LDNP,
+ BETA,XPLUSD,
+ IFIXB,IFIXX,LDIFX,
+ IDEVAL,F,FJACB,FJACD,
+ ISTOP)
implicit none
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)
C CHECK FOR UNACCEPTABLE VALUES FOR THIS PROBLEM
c
!
IF (MOD(IDEVAL,10).GE.1) THEN
DO 110 L = 1,NQ
DO 100 I = 1,N
F(I,L)=BETA(2)+BETA(1)*XPLUSD(I,1)
100 CONTINUE
110 CONTINUE
END IF
C COMPUTE DERIVATIVES WITH RESPECT TO BETA
IF (MOD(IDEVAL/10,10).GE.1) THEN
DO 210 L = 1,NQ
DO 200 I = 1,N
FJACB(I,1,L)=XPLUSD(I,1)
FJACB(I,2,L)=1.0d0
200 CONTINUE
210 CONTINUE
ENDIF
C COMPUTE DERIVATIVES WITH RESPECT TO DELTA
IF (MOD(IDEVAL/100,10).GE.1) THEN
DO 310 L = 1,NQ
DO 300 I = 1,N
FJACD(I,1,L)=BETA(1)
300 CONTINUE
310 CONTINUE
END IF
RETURN
END
!
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
File diff suppressed because it is too large Load Diff
+111
View File
@@ -0,0 +1,111 @@
-0.21875 0.2
-0.205 0.44263776
-0.19555556 0.10590558
-0.18666667 0.36658142
-0.1765 0.06329114
-0.176 0.12291351
-0.17368421 0.22368421
-0.17142857 0.36669213
-0.16789474 0.09287926
-0.1675 0.12742718
-0.1573913 0.08733624
-0.15692308 0.11617673
-0.15 0.22163511
-0.149 -0.09330629
-0.14875 0.234375
-0.1475 0.23972603
-0.14714286 0.34413766
-0.14625 -0.0308642
-0.14125 0.27777778
-0.14 0.24025974
-0.13863636 0.08765339
-0.1375 0.29329609
-0.1355 0.27157895
-0.1337037 0.05646903
-0.13291667 0.01709402
-0.13090909 -0.08833272
-0.13074074 0.00829962
-0.13 0.31606027
-0.13 0.15448604
-0.13 0.10471204
-0.12888889 0.11695906
-0.12555556 0.14814815
-0.12545455 0.01196172
-0.12518519 0.05713427
-0.12444444 0.20643594
-0.12392857 0.08314967
-0.12344828 0.07705644
-0.121 0.10019268
-0.11652174 0.23768116
-0.11375 0.10877581
-0.1115625 0.08479021
-0.10757576 0.05117845
-0.10733333 0.2659176
-0.10647059 0.09332991
-0.105 0.0976423
-0.10441176 0.07311399
-0.10294118 0.04956306
-0.10285714 0.07193372
-0.10269231 0.25865701
-0.10269231 0.24972856
-0.09722222 0.12745539
-0.095 -0.00529661
-0.09466667 0.09460738
-0.09210526 0.14687882
-0.08705882 -0.05602241
-0.08541667 -0.03205128
-0.0821875 0.21426616
-0.08105263 0.00903546
-0.07888889 0.06196581
-0.07888889 -0.03878622
-0.075 0.17806268
-0.07315789 -0.02300236
-0.07166667 0.19302153
-0.07 0.33347404
-0.06736842 0.16951037
-0.06333333 0.20434227
-0.06219512 0.16561277
-0.06219512 0.17213638
-0.06071429 0.16789396
-0.05906977 0.15890265
-0.0575 0.16435011
-0.05678571 0.17027201
-0.05208333 0.1582618
-0.05194444 0.19674797
-0.05166667 0.01883239
-0.05148148 0.2589273
-0.0498 0.15882353
-0.04891892 0.1733227
-0.04833333 0.15214385
-0.04647059 0.17944798
-0.0462 0.14204426
-0.04588235 0.05317065
-0.0440625 0.23542945
-0.04358491 0.16895538
-0.04339623 0.13951771
-0.04285714 -0.08354219
-0.04243243 0.17635673
-0.04203704 0.14302166
-0.04132075 0.16332383
-0.0412963 0.15429157
-0.04109091 0.14965035
-0.04051282 0.10450334
-0.03738095 0.13354616
-0.03488372 0.17941003
-0.0347619 0.16658253
-0.0344186 0.19873627
-0.03352941 0.11351909
-0.03333333 0.00316106
-0.03275 0.16699411
-0.031 -0.03697479
-0.03090909 0.14141414
-0.03 0
-0.029375 0.12665198
-0.02823529 0.14436343
-0.02571429 0.26308866
-0.025 0.14719848
-0.02428571 0.1073493
-0.02384615 0.00624122
-0.02333333 0.10115891
-0.02136364 0.14080196
-0.02076923 0.19122664
@@ -0,0 +1,54 @@
!----------- Variables in likelihood common blocks -----------
!
!nlikepoints: the number of sampling points for determining the parameters
! in the likelihood function
!ndesignparam: the number of design parameters
!xsamplike(nlikepoints,ndesignparam): the coordinates of the sampling points
! of the design parameters
!ysamplike(nlikepoints): the actual values at the sampling points
!delvector: = R^(-1)(y-Fb)
!coeffbas: on exit, holds the coefficients of the basis functions
!work: on exit, holds the upper triangle of the matrix R^(-1)F[F^tR^(-1)F)^(-1)F^(t)R^(-1)-R^(-1)
!fbassamp: on exit, holds the transpose of [F^(t)R^(-1)F]^(-1)F^(t)R^(-1)
!FtRinvF_inv: holds the upper triangle of the symmetric matrix [F^(t)R^(-1)F]^(-1)
!plikemin: likelihood parameter lower bounds
!plikemax: likelihood parameter upper bounds
!isitlastcall: if it is the last call, calculates the matrices needed by the prediction
!maxnvars: the maximum number of cokriging variables
!nvars: the actual number of cokriging variables
!kbasis: the total number of basis functions for all variables
!kbasvar: the number of basis functions in each variable
!maxkbasvar: the maximum number of basis functions iin each variable
!indexvars: the integer vector that holds the number of observations for each
! cokriging variables.
!ithetap: =1, pist=2.0; =2, pdist to be optimized
!kpvt: an integer array used in different subroutines, defined in common blocks
! to save memory, does not transfer data through common blocks.
integer maxdesignparam,maxlikepoints,maxkbasvar,maxnvars
parameter(maxdesignparam=40,maxlikepoints=365*48,
& maxkbasvar=100,maxnvars=10)
integer nlikepoints,ndesignparam,kbasis,kbasvar,
& nvars,indexvars,ithetap,kpvt
common /likedim/nlikepoints,ndesignparam,kbasis,
& kbasvar(maxnvars),nvars,indexvars(maxnvars),
& ithetap,kpvt(maxlikepoints)
double precision xsamplike,ysamplike,coeffbas,delvector,
& work,fbassamp,FtRinvF_inv
common /likevar/xsamplike(maxlikepoints,maxdesignparam),
& ysamplike(maxlikepoints),delvector(maxlikepoints),
& coeffbas(maxkbasvar*maxnvars),
& work(maxlikepoints*(maxlikepoints+1)/2),
& fbassamp(maxlikepoints,maxkbasvar*maxnvars),
& FtRinvF_inv(maxkbasvar*maxnvars*(maxkbasvar*maxnvars+1)/2)
double precision plikemin,plikemax
common /likebounds/plikemin(maxnvars*(maxnvars+1)*
& (2*maxdesignparam+1)/2),
& plikemax(maxnvars*(maxnvars+1)*
& (2*maxdesignparam+1)/2)
logical isitlastcall
common /likelogic/isitlastcall
!----------------------------------------------------------
+977
View File
@@ -0,0 +1,977 @@
! Multivariate Universal Maximum Likelihood Estimation (MUMLE)
! Developed by
! Lianhong Gu
! Environmental Sciences Division
! Oak Ridge National Lab
! lianhong-gu@ornl.gov
!
! Version 1.0, July, 2006
subroutine coapproxyfunc(nkrigpoints,nxdim,xsampkrig,
& ysampkrig,nvarscp,indexvarscp,iupdatelikelihood,
& iwhattodo,ithetap0,ixnewpos,xnew,
& ypred,grad_xnew,rootmeanerror)
implicit none
include '../maxlikelihood/colikelihood.h'
!!!!!!!!!!!!!!!!!!!! Arguments !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-------------------Inputs--------------------------------------
!nkrigpoints: the number of available points for kriging.
!nxdim: the dimension in the design parameter vector x
!xsampkrig: the coordinates of the sampled points for kriging. Note that these
! points may or may not be the same as those points used for estimating
! the parameters in the likelihood function
!
! ****************************************************************************
! *Note different variables are stacked in the vector (xsampkrig,ysampkrig) *
! *in an orderly fashion. The number of different variables is given by *
! *nvarscp. The number of points in each variable is indicated by indexvarscp*
! ****************************************************************************
!
!ysampkrig: the value at xsampkrig
!nvarscp: the number of variables for cokriging
!indexvarscp: the number of samples in each variable for cokriging
!iupdatelikelihood: whether to update the likelihood function
! =0: the first entering, so initilize the parameters to be optimized
! and do likelihood function optimization
! =1: update the likelihood parameters using the saved values from
! the last optimization as the initial guesses
! =2: do not do likelihood optimization
!iwhattodo: =1, only make a predition at xnew
! =2, make a prediction and calculate the derivative at xnew
! =3, only calculate the derivative at xnew
!ithetap0: =ithetap: =1 set pdist to 2.0, =2: pdist to be optimized
!ixnewpos: which variable does xnew belong to.
!xnew: the coordinates of the point to be kriged
!
!---------------------- Outputs --------------------------------
!ypred: the predicted value at xnew
!rootmeanerror: the root mean square error of ypred
!grad_xnew: the derivative at xnew
integer nkrigpoints,nxdim,iupdatelikelihood,iwhattodo,
& nvarscp,indexvarscp(nvarscp),ithetap0,ixnewpos
double precision xsampkrig(nkrigpoints,nxdim),
& ysampkrig(nkrigpoints),xnew(nxdim),ypred,
& rootmeanerror,grad_xnew(nxdim)
!
!!!!!!!!!!!!!!!!!!! Local variables !!!!!!!!!!!!!!!!!!!!!!!!!
!
double precision ftol
parameter(ftol=1.0d-8)
external df1dimd_like,gradlikefunc,
& likelihoodfunc,f1dim_like,colikelihoodgrad
integer i,j,nplike,info,nii,index,irowinc,istart,
& icompactpostri,icompactposoff
integer nbd(maxnvars*(maxnvars+1)*
& (2*maxdesignparam+1)/2)
double precision plike(maxnvars*(maxnvars+1)*
& (2*maxdesignparam+1)/2),
& x(nxdim),flikemin,xsamp(nxdim),
& theta((maxnvars*(maxnvars+1))/2,maxdesignparam),
& pdist((maxnvars*(maxnvars+1))/2,maxdesignparam),
& theta1(nxdim),pdist1(nxdim),
& covars(maxnvars),correls(maxnvars*(maxnvars-1)/2),
& rx(nkrigpoints),fb(maxkbasvar),corr_like,
& xmean(maxdesignparam),
& gradient(maxkbasvar,nxdim)
double precision gradlike(maxnvars*(maxnvars+1)*
& (2*maxdesignparam+1)/2)
save nplike
save plike,theta,pdist,xmean,covars,correls
if(iupdatelikelihood.le.1)then
! update parameters and correlation functions in the likelihood function
!
! nbd is an INTEGER array of dimension n that must be set by the
! user to the type of bounds imposed on the variables:
! nbd(i)=0 if x(i) is unbounded,
! 1 if x(i) has only a lower bound,
! 2 if x(i) has both lower and upper bounds,
! 3 if x(i) has only an upper bound.
nlikepoints=nkrigpoints
nvars=nvarscp
ithetap=ithetap0
do i=1,nvars
indexvars(i)=indexvarscp(i)
enddo
ndesignparam=nxdim
do i=1,nlikepoints
ysamplike(i)=ysampkrig(i)
enddo
! Centralize the coordinates to stablize the linear system
do j=1,nxdim
xmean(j)=0.0d0
do i=1,nlikepoints
xmean(j)=xmean(j)+xsampkrig(i,j)
enddo
xmean(j)=xmean(j)/dble(nlikepoints)
do i=1,nlikepoints
xsamplike(i,j)=xsampkrig(i,j)-xmean(j)
enddo
enddo
call setinit_bounds(nvars,ndesignparam,ithetap,
& iupdatelikelihood,plikemin,plike,plikemax,nbd,nplike)
isitlastcall=.false.
! to use Lbfgsb_2_4, no need to call the cost function first
call Lbfgsb_2_4(nplike,plike,flikemin,plikemin,plikemax,
& nbd,colikelihoodgrad,ftol,info)
! to use frprmn, the cost function needs to be called first
! call likelihoodfunc(nplike,plike,flikemin)
! call frprmn(plike,nplike,ftol,flikemin,plikemin,
! & plikemax,gradlikefunc,f1dim_like,df1dimd_like)
! call dfpmin(plike,nplike,ftol,flikemin,likelihoodfunc,
! & plikemin,plikemax,gradlikefunc,f1dim_like,df1dimd_like)
write(*,*)flikemin,plike(1),plike(2),plike(3),info
call gradlikefunc(nplike,plike,flikemin,gradlike)
write(*,*)gradlike(1),gradlike(2),gradlike(3)
call nongradopt(nplike,likelihoodfunc,f1dim_like,plike,
& plikemin,plikemax,ftol,flikemin)
write(*,*)flikemin,plike(1),plike(2),plike(3)
call gradlikefunc(nplike,plike,flikemin,gradlike)
write(*,*)gradlike(1),gradlike(2),gradlike(3)
isitlastcall=.true.
call likelihoodfunc(nplike,plike,flikemin)
call colikeliparampart(nplike,plike,ithetap,
& ndesignparam,nvars,
& theta(1:(nvars*(nvars+1))/2,1:ndesignparam),
& pdist(1:(nvars*(nvars+1))/2,1:ndesignparam),
& covars,correls)
endif
!
do j=1,ndesignparam
x(j)=xnew(j)-xmean(j)
enddo
istart=0
do i=1,ixnewpos-1
istart=istart+kbasvar(i)
enddo
!
if(iwhattodo.eq.1.or.iwhattodo.eq.2)then
! make a prediction at x
call cofuncbasis(ixnewpos,index,ndesignparam,
& x,fb)
irowinc=0
do i=1,nvars
index=icompactpostri(ixnewpos,i)
do j=1,ndesignparam
theta1(j)=theta(index,j)
pdist1(j)=pdist(index,j)
enddo
if(ixnewpos.ne.i)then
index=icompactposoff(ixnewpos,i)
endif
do nii=1,indexvars(i)
irowinc=irowinc+1
do j=1,ndesignparam
xsamp(j)=xsamplike(irowinc,j)
enddo
if(ixnewpos.eq.i)then
rx(irowinc)=covars(i)*corr_like(x,xsamp,
& theta1,pdist1,ndesignparam)
else
rx(irowinc)=dsqrt(covars(i)*covars(ixnewpos))*
& correls(index)*corr_like(x,xsamp,
& theta1,pdist1,ndesignparam)
endif
enddo
enddo
ypred=0.0d0
do i=1,kbasvar(ixnewpos)
ypred=ypred+fb(i)*coeffbas(i+istart)
enddo
do i=1,nlikepoints
ypred=ypred+delvector(i)*rx(i)
enddo
rootmeanerror=covars(ixnewpos)
do i=1,kbasis
call symmatindex(kbasis,i,kpvt)
flikemin=0.0d0
do j=1,kbasis
if(j.gt.istart.and.j.le.(istart+kbasvar(ixnewpos)))then
flikemin=flikemin+FtRinvF_inv(kpvt(j))*fb(j-istart)
endif
enddo
if(i.gt.istart.and.i.le.(istart+kbasvar(ixnewpos)))then
rootmeanerror=rootmeanerror+fb(i-istart)*flikemin
endif
enddo
do i=1,kbasis
flikemin=0.0d0
do j=1,nlikepoints
flikemin=flikemin+fbassamp(j,i)*rx(j)
enddo
if(i.gt.istart.and.i.le.(istart+kbasvar(ixnewpos)))then
rootmeanerror=rootmeanerror-2.0d0*fb(i)*flikemin
endif
enddo
do i=1,nlikepoints
call symmatindex(nlikepoints,i,kpvt)
flikemin=0.0d0
do j=1,nlikepoints
flikemin=flikemin+work(kpvt(j))*rx(j)
enddo
rootmeanerror=rootmeanerror+rx(i)*flikemin
enddo
if(rootmeanerror.lt.0.0d0)then
! computational error
rootmeanerror=0.0d0
endif
rootmeanerror=dsqrt(rootmeanerror)
endif
!
if(iwhattodo.eq.2.or.iwhattodo.eq.3)then
! calculate the derivative of the approximated function with respect to x
call cogradbasx(ixnewpos,kbasvar(ixnewpos),ndesignparam,
& x,gradient(1:kbasvar(ixnewpos),1:ndesignparam))
do i=1,ndesignparam
grad_xnew(i)=0.0d0
do j=1,kbasvar(ixnewpos)
grad_xnew(i)=grad_xnew(i)+gradient(j,i)*coeffbas(j+istart)
enddo
enddo
irowinc=0
do i=1,nvars
if(i.eq.ixnewpos)then
flikemin=covars(ixnewpos)
else
index=icompactposoff(i,ixnewpos)
flikemin=dsqrt(covars(i)*covars(ixnewpos))*
& correls(index)
endif
index=icompactpostri(i,ixnewpos)
do j=1,ndesignparam
theta1(j)=theta(index,j)
pdist1(j)=pdist(index,j)
enddo
do j=1,indexvars(i)
irowinc=irowinc+1
do index=1,ndesignparam
xsamp(index)=xsamplike(irowinc,index)
enddo
call gradcorratx(x,xsamp,ndesignparam,theta1,
& pdist1,rx)
do index=1,ndesignparam
grad_xnew(index)=grad_xnew(index)+flikemin*
& rx(index)*delvector(irowinc)
enddo
enddo
enddo
endif
return
end
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
subroutine colikelihoodgrad(nplike,plike,dogradient,
& fvalue,gradlike,ierr)
implicit none
include '../maxlikelihood/colikelihood.h'
!
! calculate the transformed likelihood function value and derivatives
! with respect to input parameters if specified
!
!----------- Arguments ------------------------------------
integer nplike,ierr
logical dogradient
double precision plike(nplike),fvalue,gradlike(nplike)
!
!=> nplike: the number of likelihood parameters to be estimated
!=> plike: the likelihood parameters to be estimated
!=> dogradient: whether to do derivative calculations
!<= fvalue: the transformed likelihood function value
!<= gradlike: derivatives at plike
!<= ierr: error messages
! =0: every thing ok
! <0: parameter out of bounds. |ierr| denotes the out-of-bound parameter
! =11: the correlation matrix is not postitive definite (the determinant is negative)
!----------- End of list of arguments ---------------------
!------------- Local variables ----------------------------
integer i,j,n,t,inert(3),index,nii,
& irowinc,getwhichvar,ivar,jvar,icount,k,
& icompactposoff,icompactpostri,iloop,istart,jstart
double precision lnabsdet,signunit,delta,corr_like,
& covars(nvars),theta(nvars*(nvars+1)/2,ndesignparam),
& pdist(nvars*(nvars+1)/2,ndesignparam),
& correls(nvars*(nvars-1)/2),xi(ndesignparam),
& xj(ndesignparam),theta1(ndesignparam),
& pdist1(ndesignparam),derR(nlikepoints*(nlikepoints+1)/2),
& sumi,sumt,gradtheta(ndesignparam),gradp(ndesignparam)
do i=1,nplike
if(plike(i).lt.plikemin(i).or.
& plike(i).gt.plikemax(i))then
! penalize out-of-bounds parameters
fvalue=1.0d+100
ierr=-i
return
endif
enddo
call colikeliparampart(nplike,plike,ithetap,
& ndesignparam,nvars,
& theta(1:(nvars*(nvars+1))/2,1:ndesignparam),
& pdist(1:(nvars*(nvars+1))/2,1:ndesignparam),
& covars,correls)
t=0
do j=1,nlikepoints
jvar=getwhichvar(j,nvars,indexvars)
do n=1,ndesignparam
xi(n)=xsamplike(j,n)
enddo
do i=1,j
t=t+1
if(i.eq.j)then
work(t)=covars(jvar)
else
ivar=getwhichvar(i,nvars,indexvars)
k=icompactpostri(ivar,jvar)
do n=1,ndesignparam
xj(n)=xsamplike(i,n)
theta1(n)=theta(k,n)
pdist1(n)=pdist(k,n)
enddo
if(ivar.eq.jvar)then
work(t)=covars(jvar)*corr_like(xj,xi,
& theta1,pdist1,ndesignparam)
else
k=icompactposoff(ivar,jvar)
work(t)=dsqrt(covars(jvar)*covars(ivar))*
& correls(k)*corr_like(xj,xi,
& theta1,pdist1,ndesignparam)
endif
endif
enddo
enddo
call dspfa(work,nlikepoints,kpvt,i)
j=11
call dspdi(work,nlikepoints,kpvt,derR,inert,delvector,j)
! work now stores the upper triangle of the inverse covariance matrix
! determinant=derR(1)*10.0d0**derR(2)
signunit=dsign(1.0d0,derR(1))
lnabsdet=dlog(dabs(derR(1)))+derR(2)*dlog(10.0d0)
if(signunit.lt.0.0d0)then
! write(*,*)'Wrong! the covariance matrix has a negative'
! write(*,*)'determinant. Program stops in likelihoodfunc'
write(*,*)signunit,lnabsdet
ierr=11
return
endif
t=0
kbasis=0
do i=1,nvars
do j=1,indexvars(i)
t=t+1
do n=1,ndesignparam
xi(n)=xsamplike(t,n)
enddo
call cofuncbasis(i,kbasvar(i),ndesignparam,
& xi,delvector)
do n=1,kbasvar(i)
fbassamp(t,kbasis+n)=delvector(n)
enddo
enddo
if(i.gt.1)then
do j=t-indexvars(i)+1,t
do n=1,kbasis
fbassamp(j,n)=0.0d0
enddo
enddo
do j=1,t
do n=kbasis+1,kbasis+kbasvar(i)
fbassamp(j,n)=0.0d0
enddo
enddo
endif
kbasis=kbasis+kbasvar(i)
enddo
call solvebassystem(delta)
fvalue=delta+lnabsdet
!
! The actual likelihood function =exp(-(fvalue+nlikepoints*ln(2*pi))/2)
if(dogradient.eqv..false..or.dogradient.eqv..FALSE.)then
ierr=0
return
endif
if(isitlastcall.eqv..true..or.isitlastcall.eqv..TRUE.)then
ierr=0
return
endif
!
! section for derivatives starts
!
! delvector=R^(-1)(Y-FB) was computed in solvebassystem and will not be altered from now
!----Section for calculating derivatives with respect to theta and pdist-------
icount=0
do iloop=1,ithetap
do k=1,nvars*(nvars+1)/2
call ivarjvartri(k,nvars,indexvars,
& ivar,jvar,istart,jstart)
if(ivar.ne.jvar)then
nii=icompactposoff(ivar,jvar)
endif
do n=1,ndesignparam
theta1(n)=theta(k,n)
pdist1(n)=pdist(k,n)
enddo
ierr=0
do j=1,nlikepoints
do i=1,j
ierr=ierr+1
derR(ierr)=0.0d0
enddo
enddo
do t=1,ndesignparam
icount=icount+1
do j=jstart,jstart+indexvars(jvar)-1
do n=1,ndesignparam
xj(n)=xsamplike(j,n)
enddo
if(ivar.eq.jvar)then
kpvt(1)=j-1
else
kpvt(1)=istart+indexvars(ivar)-1
endif
do i=istart,kpvt(1)
do n=1,ndesignparam
xi(n)=xsamplike(i,n)
enddo
ierr=icompactpostri(i,j)
if(ivar.eq.jvar)then
if(iloop.eq.1)then
call gradcorr_liketheta(xi,xj,
& theta1,pdist1,ndesignparam,gradtheta)
derR(ierr)=covars(ivar)*gradtheta(t)
else
call gradcorr_likep(xi,xj,
& theta1,pdist1,ndesignparam,gradp)
derR(ierr)=covars(ivar)*gradp(t)
endif
else
if(iloop.eq.1)then
call gradcorr_liketheta(xi,xj,
& theta1,pdist1,ndesignparam,gradtheta)
derR(ierr)=dsqrt(covars(ivar)*covars(jvar))*
& correls(nii)*gradtheta(t)
else
call gradcorr_likep(xi,xj,
& theta1,pdist1,ndesignparam,gradp)
derR(ierr)=dsqrt(covars(ivar)*covars(jvar))*
& correls(nii)*gradp(t)
endif
endif
enddo
enddo
gradlike(icount)=0.0d0
do i=1,nlikepoints
call symmatindex(nlikepoints,i,kpvt)
sumi=0.0d0
do j=1,nlikepoints
sumi=sumi+derR(kpvt(j))*delvector(j)
enddo
gradlike(icount)=
& gradlike(icount)+delvector(i)*sumi
enddo
gradlike(icount)=-gradlike(icount)
!
!Now, the trace part
sumi=0.0d0
do i=1,nlikepoints
call symmatindex(nlikepoints,i,kpvt)
sumt=0.0d0
do j=1,nlikepoints
sumt=sumt+work(kpvt(j))*derR(kpvt(j))
enddo
sumi=sumi+sumt
enddo
! sumi is the trace
gradlike(icount)=gradlike(icount)+sumi
enddo
enddo
enddo
!-------------------------------------------------------------------
! ------- Section for calculating derivatives with respect to variance
do k=1,nvars
icount=icount+1
ierr=0
do j=1,nlikepoints
do i=1,j
ierr=ierr+1
derR(ierr)=0.0d0
enddo
enddo
! first, vertical elements
do i=1,k-1
ierr=icompactpostri(i,k)
do t=1,ndesignparam
theta1(t)=theta(ierr,t)
pdist1(t)=pdist(ierr,t)
enddo
ierr=icompactposoff(i,k)
call istartjstart(i,k,nvars,indexvars,istart,jstart)
do n=istart,istart+indexvars(i)-1
do t=1,ndesignparam
xi(t)=xsamplike(n,t)
enddo
do j=jstart,jstart+indexvars(k)-1
do t=1,ndesignparam
xj(t)=xsamplike(j,t)
enddo
ivar=icompactpostri(n,j)
derR(ivar)=0.5d0*dsqrt(covars(i)/covars(k))*
& correls(ierr)*corr_like(xi,xj,
& theta1,pdist1,ndesignparam)
enddo
enddo
enddo
ierr=icompactpostri(k,k)
do t=1,ndesignparam
theta1(t)=theta(ierr,t)
pdist1(t)=pdist(ierr,t)
enddo
call istartjstart(k,k,nvars,indexvars,istart,jstart)
do j=jstart,jstart+indexvars(k)-1
do t=1,ndesignparam
xi(t)=xsamplike(j,t)
enddo
do n=istart,j
ivar=icompactpostri(n,j)
if(n.eq.j)then
derR(ivar)=1.0d0
else
do t=1,ndesignparam
xj(t)=xsamplike(n,t)
enddo
derR(ivar)=corr_like(xi,xj,
& theta1,pdist1,ndesignparam)
endif
enddo
enddo
! then, horizontal elements
do i=k+1,nvars
ierr=icompactpostri(k,i)
do t=1,ndesignparam
theta1(t)=theta(ierr,t)
pdist1(t)=pdist(ierr,t)
enddo
ierr=icompactposoff(k,i)
call istartjstart(k,i,nvars,indexvars,istart,jstart)
do n=istart,istart+indexvars(i)-1
do t=1,ndesignparam
xi(t)=xsamplike(n,t)
enddo
do j=jstart,jstart+indexvars(k)-1
do t=1,ndesignparam
xj(t)=xsamplike(j,t)
enddo
ivar=icompactpostri(n,j)
derR(ivar)=0.5d0*dsqrt(covars(i)/covars(k))*
& correls(ierr)*corr_like(xi,xj,
& theta1,pdist1,ndesignparam)
enddo
enddo
enddo
gradlike(icount)=0.0d0
do i=1,nlikepoints
call symmatindex(nlikepoints,i,kpvt)
sumi=0.0d0
do j=1,nlikepoints
sumi=sumi+derR(kpvt(j))*delvector(j)
enddo
gradlike(icount)=
& gradlike(icount)+delvector(i)*sumi
enddo
gradlike(icount)=-gradlike(icount)
!
!Now, the trace part
sumi=0.0d0
do i=1,nlikepoints
call symmatindex(nlikepoints,i,kpvt)
sumt=0.0d0
do j=1,nlikepoints
sumt=sumt+work(kpvt(j))*derR(kpvt(j))
enddo
sumi=sumi+sumt
enddo
! sumi is the trace
gradlike(icount)=gradlike(icount)+sumi
enddo
!--------------------------------------------------------------------
! ------- Section for calculating derivatives with respect to correls
do k=1,nvars*(nvars-1)/2
icount=icount+1
ierr=0
do j=1,nlikepoints
do i=1,j
ierr=ierr+1
derR(ierr)=0.0d0
enddo
enddo
call ivarjvaroff(k,nvars,indexvars,
& ivar,jvar,istart,jstart)
ierr=icompactpostri(ivar,jvar)
do t=1,ndesignparam
theta1(t)=theta(ierr,t)
pdist1(t)=pdist(ierr,t)
enddo
do i=istart,istart+indexvars(ivar)-1
do t=1,ndesignparam
xi(t)=xsamplike(i,t)
enddo
do j=jstart,jstart+indexvars(jvar)-1
do t=1,ndesignparam
xj(t)=xsamplike(j,t)
enddo
ierr=icompactpostri(i,j)
derR(ierr)=dsqrt(covars(ivar)*covars(jvar))*
& corr_like(xi,xj,theta1,pdist1,ndesignparam)
enddo
enddo
gradlike(icount)=0.0d0
do i=1,nlikepoints
call symmatindex(nlikepoints,i,kpvt)
sumi=0.0d0
do j=1,nlikepoints
sumi=sumi+derR(kpvt(j))*delvector(j)
enddo
gradlike(icount)=
& gradlike(icount)+delvector(i)*sumi
enddo
gradlike(icount)=-gradlike(icount)
!Now, the trace part
sumi=0.0d0
do i=1,nlikepoints
call symmatindex(nlikepoints,i,kpvt)
sumt=0.0d0
do j=1,nlikepoints
sumt=sumt+work(kpvt(j))*derR(kpvt(j))
enddo
sumi=sumi+sumt
enddo
! sumi is the trace
gradlike(icount)=gradlike(icount)+sumi
enddo
!--------------------------------------------------------------
ierr=0
return
end
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
subroutine solvebassystem(delta)
implicit none
include '../maxlikelihood/colikelihood.h'
!
! Calculating the best linear estimates of the coefficients of the basis functions,
! delta (the residual square error), and matrices needed for prediction and estimation
! of prediction error
integer kcopy(kbasis),inert(3)
double precision FR(kbasis,nlikepoints),
& a(kbasis*(kbasis+1)/2),rhs(kbasis),sum,
& acopy(max0(kbasis*(kbasis+1)/2,nlikepoints)),
& amat(kbasis*(kbasis+1)/2),
& rhscopy(kbasis),delta
integer n,j,t,i,mark,irowinc,index,nii
do i=1,nlikepoints
call symmatindex(nlikepoints,i,kpvt)
do j=1,kbasis
FR(j,i)=0.0d0
do t=1,nlikepoints
FR(j,i)=FR(j,i)+work(kpvt(t))*fbassamp(t,j)
enddo
enddo
enddo
do i=1,kbasis
! right hand side of the linear system
rhs(i)=0.0d0
do j=1,nlikepoints
rhs(i)=rhs(i)+FR(i,j)*ysamplike(j)
enddo
! make a copy
rhscopy(i)=rhs(i)
enddo
! coefficient matrix, symmetric
t=0
do j=1,kbasis
do i=1,j
t=t+1
a(t)=0.0d0
do n=1,nlikepoints
a(t)=a(t)+FR(i,n)*fbassamp(n,j)
enddo
! make a copy
amat(t)=a(t)
enddo
enddo
! solve the basis function linear system
call dspfa(a,kbasis,kpvt,mark)
do i=1,kbasis*(kbasis+1)/2
acopy(i)=a(i)
enddo
do i=1,kbasis
kcopy(i)=kpvt(i)
enddo
if(isitlastcall.eqv..true..or.isitlastcall.eqv..TRUE.)then
! calculate the inverse only
do j=1,kbasis*(kbasis+1)/2
FtRinvF_inv(j)=a(j)
enddo
j=1
call dspdi(FtRinvF_inv,kbasis,kpvt,coeffbas,
& inert,delvector,j)
do i=1,kbasis
kpvt(i)=kcopy(i)
enddo
endif
call dspsl(a,kbasis,kpvt,rhs)
do i=1,kbasis
coeffbas(i)=rhs(i)
enddo
! one step improvement
do i=1,kbasis
call symmatindex(kbasis,i,kpvt)
sum=0.0d0
do t=1,kbasis
sum=sum+amat(kpvt(t))*coeffbas(t)
enddo
rhs(i)=(sum-rhscopy(i))*1.0d+9
enddo
call dspsl(acopy,kbasis,kcopy,rhs)
do i=1,kbasis
coeffbas(i)=coeffbas(i)-rhs(i)*1.0d-9
enddo
do i=1,nlikepoints
acopy(i)=0.0d0
do j=1,kbasis
acopy(i)=acopy(i)+coeffbas(j)*fbassamp(i,j)
enddo
acopy(i)=ysamplike(i)-acopy(i)
enddo
do t=1,nlikepoints
call symmatindex(nlikepoints,t,kpvt)
delvector(t)=0.0d0
do i=1,nlikepoints
delvector(t)=delvector(t)+work(kpvt(i))*acopy(i)
enddo
enddo
delta=0.0d0
do t=1,nlikepoints
delta=delta+delvector(t)*acopy(t)
enddo
if(isitlastcall.eqv..true..or.isitlastcall.eqv..TRUE.)then
! All parameters in the likelihood function have been optimized. Prepare
! matrices for MSE calculations
do i=1,kbasis
call symmatindex(kbasis,i,kpvt)
do j=1,nlikepoints
fbassamp(j,i)=0.0d0
do t=1,kbasis
fbassamp(j,i)=fbassamp(j,i)+FtRinvF_inv(kpvt(t))*
& FR(t,j)
enddo
enddo
enddo
t=0
do j=1,nlikepoints
do i=1,j
t=t+1
work(t)=-work(t)
do n=1,kbasis
work(t)=work(t)+fbassamp(i,n)*FR(n,j)
enddo
enddo
enddo
endif
return
end
!######################################################
DOUBLE PRECISION FUNCTION df1dimd_like(x)
INTEGER NMAX
double precision x
PARAMETER (NMAX=50)
CU USES gradlikefunc
INTEGER j,ncom
double precision df(NMAX),pcom(NMAX),xicom(NMAX),
& xt(NMAX)
,dummy
COMMON /f1com/ pcom,xicom,ncom
do 11 j=1,ncom
xt(j)=pcom(j)+x*xicom(j)
11 continue
call gradlikefunc(ncom,xt,dummy,df)
df1dimd_like=0.0d0
do 12 j=1,ncom
df1dimd_like=df1dimd_like+df(j)*xicom(j)
12 continue
return
END
C (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.
subroutine gradlikefunc(nplike,plike,fvalue,gradlike)
implicit none
integer nplike,ierr
double precision plike(nplike),fvalue,gradlike(nplike)
logical dogradient
dogradient=.true.
call colikelihoodgrad(nplike,plike,dogradient,
& fvalue,gradlike,ierr)
return
end
subroutine likelihoodfunc(nplike,plike,fvalue)
implicit none
integer nplike,ierr
double precision plike(nplike),fvalue,gradlike(nplike)
logical dogradient
dogradient=.false.
call colikelihoodgrad(nplike,plike,dogradient,
& fvalue,gradlike,ierr)
return
end
double precision function f1dim_like(x)
INTEGER NMAX
double precision x
PARAMETER (NMAX=50)
CU USES likelihood
INTEGER j,ncom
double precision pcom(NMAX),xicom(NMAX),xt(NMAX)
COMMON /f1com/ pcom,xicom,ncom
do 11 j=1,ncom
xt(j)=pcom(j)+x*xicom(j)
11 continue
call likelihoodfunc(ncom,xt,f1dim_like)
return
END
C (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.
!#################################################################
subroutine gradfunkmin(ndim,x,fatx,fjac)
implicit none
integer ndim,j
double precision x(ndim),fatx,fjac(ndim),h,eps,
& zero,temp,fxh
call likelihoodfunc(ndim,x,fatx)
eps = 1.0d-5
zero=0.0d0
do j = 1, ndim
temp = x(j)
h = eps*dabs(temp)
if (h .eq. zero) h = eps
x(j) = temp + h
call likelihoodfunc(ndim,x,fxh)
x(j) = temp
fjac(j) = (fxh - fatx)/h
enddo
return
end subroutine gradfunkmin
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
subroutine randpermut(npoints,ndim,x)
implicit none
!
! conduct random permutation
integer npoints,ndim,i,j,index
double precision x(npoints,ndim),xtemp(npoints),
& ran2
do i=1,ndim
do j=1,npoints
xtemp(j)=x(j,i)
enddo
do j=1,npoints
index=int(dble(npoints-j+1)*ran2()+1.0d0)
x(j,i)=xtemp(index)
xtemp(index)=xtemp(npoints-j+1)
enddo
enddo
return
end
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
c
c####################################################################
c random number generator
c
double precision function ran2()
implicit none
integer im1,im2,imm1,ia1,ia2,iq1,iq2,ir1,ir2,ntab,ndiv
double precision am,eps,rnmx
parameter(im1=2147483563,im2=2147483399,am=1.0d0/dble(im1),
&imm1=im1-1,ia1=40014,ia2=40692,iq1=53668,iq2=52774,ir1=
&12211,ir2=3791,ntab=32,ndiv=1+imm1/ntab,eps=1.2d-7,
&rnmx=1.0d0-eps)
integer idum2,j,k,iv(ntab),iy,idum
save iv,iy,idum2,idum
data idum2/123456789/,iv/ntab*0/,iy/0/,idum/-1/
if(idum.le.0) then
idum=max0(-idum,1)
idum2=idum
do 11 j=ntab+8,1,-1
k=idum/iq1
idum=ia1*(idum-k*iq1)-k*ir1
if(idum.lt.0) idum=idum+im1
if(j.le.ntab) iv(j)=idum
11 continue
iy=iv(1)
end if
k=idum/iq1
idum=ia1*(idum-k*iq1)-k*ir1
if(idum.lt.0)idum=idum+im1
k=idum2/iq2
idum2=ia2*(idum2-k*iq2)-k*ir2
if(idum2.lt.0) idum2=idum2+im2
j=1+iy/ndiv
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1)iy=iy+imm1
ran2=dmin1(am*dble(iy),rnmx)
return
end
@@ -0,0 +1,23 @@
program testnewton
implicit none
integer i
double precision x,f0,f1,func,recider
x=-1.5d0
do i=1,200
f0=func(x)
write(*,*)i,x,f0
pause
f1=func(x+f0)
recider=f0/(f1-f0)
x=x-recider*f0
enddo
end
double precision function func(x)
double precision x
func=x-(x*x+1.0d0)/2.0d0
!x^2-2*x+1=0
return
end
@@ -0,0 +1,66 @@
subroutine bookkeeping(nunknowns,xvar,fequ,
& icall,iflargest)
implicit none
include 'nslasystem.h'
integer nunknowns,icall,iflargest
double precision xvar(nunknowns),fequ(nunknowns)
integer iGuCall,i,j,k
parameter(iGuCall=49)
!--------------------------------------------------------------------------
iflargest=1
do j=2,nunknowns
if(dabs(fequ(j)).gt.dabs(fequ(iflargest)))
& iflargest=j
enddo
if(numeval.eq.maxeval)goto 100
if(numeval.eq.0.or.icall.eq.iGuCall)then
numeval=numeval+1
flargest(numeval)=dabs(fequ(iflargest))
do i=1,nunknowns
xevaluated(numeval,i)=xvar(i)
fevaluated(numeval,i)=fequ(i)
enddo
return
endif
100 do i=1,numeval
k=0
do j=1,nunknowns
if(dabs(xvar(j)-xevaluated(i,j)).gt.
& 1.0d-5*dabs(xvar(j)))k=1
enddo
if(k.eq.0)goto 500
enddo
if(numeval.lt.maxeval)then
numeval=numeval+1
flargest(numeval)=dabs(fequ(iflargest))
do i=1,nunknowns
xevaluated(numeval,i)=xvar(i)
fevaluated(numeval,i)=fequ(i)
enddo
return
endif
! replace a point
j=1
do i=2,numeval
if(flargest(j).lt.flargest(i))then
j=i
endif
enddo
if(dabs(fequ(iflargest)).lt.flargest(j))then
flargest(j)=dabs(fequ(iflargest))
do i=1,nunknowns
xevaluated(j,i)=xvar(i)
fevaluated(j,i)=fequ(i)
enddo
endif
return
! too close to the existing point i
500 if(dabs(fequ(iflargest)).lt.flargest(i))then
flargest(i)=dabs(fequ(iflargest))
do j=1,nunknowns
xevaluated(i,j)=xvar(j)
fevaluated(i,j)=fequ(j)
enddo
endif
return
end
+567
View File
@@ -0,0 +1,567 @@
C This file contains all the subroutines needed by the nonlinear solver broydn
C code modified based on on-line version in Numerical Recipes Website, Dec 28, 2004
SUBROUTINE broydn(x0min,x,x0max,STPMX,n,
& fveccopy,funcv,TOLF,ierr)
implicit none
INTEGER n,NP,MAXITS,ierr
double precision x(n),EPS,TOLF,TOLMIN,TOLX,STPMX
double precision x0min(n),x0max(n),fveccopy(n)
LOGICAL check
! PARAMETER (NP=1000,MAXITS=250,EPS=1.0d-7,TOLF=1.0d-4,
! & TOLMIN=1.d-6,TOLX=EPS)
PARAMETER(NP=1000,MAXITS=250,EPS=1.0d-7,TOLX=EPS)
CU USES fdjac,funcv,lnsrch,qrdcmp,qrupdt,rsolv
INTEGER i,its,j,k
double precision den,f,fold,stpmax,sum,temp,test,c(NP),
& d(NP),fvcold(NP),g(NP),p(NP),qt(NP,NP),r(NP,NP),
& s(NP),t(NP),w(NP),xold(NP),fvec(NP)
LOGICAL restrt,sing,skip
EXTERNAL funcv
TOLMIN=TOLF*0.01d0
call funcv(n,x,fvec,f)
test=0.0d0
do 11 i=1,n
if(dabs(fvec(i)).gt.test)test=dabs(fvec(i))
11 continue
if(test.lt..01d0*TOLF)then
ierr=0
! check=.false.
return
endif
sum=0.0d0
do 12 i=1,n
sum=sum+x(i)*x(i)
12 continue
stpmax=STPMX*dmax1(dsqrt(sum),dble(n))
restrt=.true.
do 42 its=1,MAXITS
if(restrt)then
do i=1,n
if(x(i).lt.x0min(i).or.x(i).gt.x0max(i))then
ierr=1
return
endif
enddo
call fdjac(n,x,fvec,NP,r,funcv)
call qrdcmp(r,n,NP,c,d,sing)
! if(sing) pause 'singular Jacobian in broydn'
if(sing)then
ierr=2
return
end if
do 14 i=1,n
do 13 j=1,n
qt(i,j)=0.0d0
13 continue
qt(i,i)=1.0d0
14 continue
do 18 k=1,n-1
if(c(k).ne.0.0d0)then
do 17 j=1,n
sum=0.0d0
do 15 i=k,n
sum=sum+r(i,k)*qt(i,j)
15 continue
sum=sum/c(k)
do 16 i=k,n
qt(i,j)=qt(i,j)-sum*r(i,k)
16 continue
17 continue
endif
18 continue
do 21 i=1,n
r(i,i)=d(i)
do 19 j=1,i-1
r(i,j)=0.0d0
19 continue
21 continue
else
do 22 i=1,n
s(i)=x(i)-xold(i)
22 continue
do 24 i=1,n
sum=0.0d0
do 23 j=i,n
sum=sum+r(i,j)*s(j)
23 continue
t(i)=sum
24 continue
skip=.true.
do 26 i=1,n
sum=0.0d0
do 25 j=1,n
sum=sum+qt(j,i)*t(j)
25 continue
w(i)=fvec(i)-fvcold(i)-sum
if(dabs(w(i)).ge.EPS*(dabs(fvec(i))+
& dabs(fvcold(i))))then
skip=.false.
else
w(i)=0.0d0
endif
26 continue
if(.not.skip)then
do 28 i=1,n
sum=0.0d0
do 27 j=1,n
sum=sum+qt(i,j)*w(j)
27 continue
t(i)=sum
28 continue
den=0.0d0
do 29 i=1,n
den=den+s(i)*s(i)
29 continue
do 31 i=1,n
s(i)=s(i)/den
31 continue
call qrupdt(r,qt,n,NP,t,s)
do 32 i=1,n
if(r(i,i).eq.0.0d0) then
write(*,*) 'r singular in broydn'
end if
d(i)=r(i,i)
32 continue
endif
endif
do 34 i=1,n
sum=0.0d0
do 33 j=1,n
sum=sum+qt(i,j)*fvec(j)
33 continue
p(i)=-sum
34 continue
do 36 i=n,1,-1
sum=0.0d0
do 35 j=1,i
sum=sum-r(j,i)*p(j)
35 continue
g(i)=sum
36 continue
do 37 i=1,n
xold(i)=x(i)
fvcold(i)=fvec(i)
37 continue
fold=f
call rsolv(r,n,NP,d,p)
! Gu modification starts
do 100 i=1,n
if(xold(i).lt.x0min(i).or.xold(i).gt.x0max(i))then
ierr=1
return
endif
100 continue
! Gu modification ends
call lnsrch(n,xold,fold,g,p,x,f,
& stpmax,check,funcv,fvec)
test=0.0d0
do 38 i=1,n
if(dabs(fvec(i)).gt.test)test=dabs(fvec(i))
fveccopy(i)=fvec(i)
38 continue
if(test.lt.TOLF)then
ierr=0
! check=.false.
return
endif
if(check)then
if(restrt)then
ierr=3
return
else
test=0.0d0
den=dmax1(f,.5d0*dble(n))
do 39 i=1,n
temp=dabs(g(i))*dmax1(dabs(x(i)),1.0d0)/den
if(temp.gt.test)test=temp
39 continue
if(test.lt.TOLMIN)then
ierr=4
return
else
restrt=.true.
endif
endif
else
restrt=.false.
test=0.0d0
do 41 i=1,n
temp=(dabs(x(i)-xold(i)))/dmax1(dabs(x(i)),1.0d0)
if(temp.gt.test)test=temp
41 continue
if(test.lt.TOLX)then
ierr=4
! check=.true.
return
endif
endif
42 continue
ierr=5
return
END
SUBROUTINE fdjac(n,x,fvec,np,df,funcv)
implicit none
INTEGER n,np
double precision df(np,np),fvec(n),x(n),EPS
PARAMETER (EPS=1.0d-4)
CU USES funcv
INTEGER i,j,k
double precision h,temp,f(n),fsqsum
external funcv
do 12 j=1,n
temp=x(j)
h=EPS*dabs(temp)
if(h.eq.0.0d0)h=EPS
x(j)=temp+h
h=x(j)-temp
call funcv(n,x,f,fsqsum)
x(j)=temp
do 11 i=1,n
df(i,j)=(f(i)-fvec(i))/h
11 continue
12 continue
return
END
c
SUBROUTINE lnsrch(n,xold,fold,g,p,x,f,
& stpmax,check,funcv,fvec)
implicit none
INTEGER n
LOGICAL check
DOUBLE PRECISION f,fold,stpmax,g(n),p(n),x(n),
*xold(n),ALF,TOLX,fvec(n)
PARAMETER (ALF=1.d-4,TOLX=1.d-7)
EXTERNAL funcv
CU USES funcv
INTEGER i
DOUBLE PRECISION a,alam,alam2,alamin,b,disc,
*f2,rhs1,rhs2,slope,sum,temp,test,tmplam
check=.false.
sum=0.0d0
do 11 i=1,n
sum=sum+p(i)*p(i)
11 continue
sum=dsqrt(sum)
if(sum.gt.stpmax)then
do 12 i=1,n
p(i)=p(i)*stpmax/sum
12 continue
endif
slope=0.0d0
do 13 i=1,n
slope=slope+g(i)*p(i)
13 continue
! if(slope.ge.0.0d0)pause 'roundoff problem in lnsrch'
test=0.0d0
do 14 i=1,n
temp=dabs(p(i))/dmax1(dabs(xold(i)),1.0d0)
if(temp.gt.test)test=temp
14 continue
alamin=TOLX/test
alam=1.0d0
1 continue
do 15 i=1,n
x(i)=xold(i)+alam*p(i)
15 continue
call funcv(n,x,fvec,f)
if(alam.lt.alamin)then
do 16 i=1,n
x(i)=xold(i)
16 continue
check=.true.
return
else if(f.le.fold+ALF*alam*slope)then
return
else
if(alam.eq.1.0d0)then
tmplam=-slope/(2.0d0*(f-fold-slope))
else
rhs1=f-fold-alam*slope
rhs2=f2-fold-alam2*slope
a=(rhs1/alam**2-rhs2/alam2**2)/(alam-alam2)
b=(-alam2*rhs1/alam**2+alam*rhs2/alam2**2)/
& (alam-alam2)
if(a.eq.0.0d0)then
tmplam=-slope/(2.0d0*b)
else
disc=b*b-3.0d0*a*slope
if(disc.lt.0.0d0) then
tmplam=0.5d0*alam
else if(b.le.0.0d0)then
tmplam=(-b+dsqrt(disc))/(3.0d0*a)
else
tmplam=-slope/(b+dsqrt(disc))
endif
endif
if(tmplam.gt..5d0*alam)tmplam=.5d0*alam
endif
endif
alam2=alam
f2=f
alam=dmax1(tmplam,.1d0*alam)
goto 1
END
c
SUBROUTINE qrdcmp(a,n,np,c,d,sing)
implicit none
INTEGER n,np
DOUBLE PRECISION a(np,np),c(n),d(n)
LOGICAL sing
INTEGER i,j,k
DOUBLE PRECISION scale,sigma,sum,tau
sing=.false.
do 17 k=1,n-1
scale=0.0d0
do 11 i=k,n
scale=dmax1(scale,dabs(a(i,k)))
11 continue
if(scale.eq.0.0d0)then
sing=.true.
c(k)=0.0d0
d(k)=0.0d0
else
do 12 i=k,n
a(i,k)=a(i,k)/scale
12 continue
sum=0.0d0
do 13 i=k,n
sum=sum+a(i,k)**2
13 continue
sigma=dsign(dsqrt(sum),a(k,k))
a(k,k)=a(k,k)+sigma
c(k)=sigma*a(k,k)
d(k)=-scale*sigma
do 16 j=k+1,n
sum=0.0d0
do 14 i=k,n
sum=sum+a(i,k)*a(i,j)
14 continue
tau=sum/c(k)
do 15 i=k,n
a(i,j)=a(i,j)-tau*a(i,k)
15 continue
16 continue
endif
17 continue
d(n)=a(n,n)
if(d(n).eq.0.0d0)sing=.true.
return
END
c
SUBROUTINE qrupdt(r,qt,n,np,u,v)
implicit none
INTEGER n,np
DOUBLE PRECISION r(np,np),qt(np,np),u(np),v(np)
CU USES rotate
INTEGER i,j,k
do 11 k=n,1,-1
if(u(k).ne.0.0d0)goto 1
11 continue
k=1
1 do 12 i=k-1,1,-1
call rotate(r,qt,n,np,i,u(i),-u(i+1))
if(u(i).eq.0.0d0)then
u(i)=dabs(u(i+1))
else if(dabs(u(i)).gt.dabs(u(i+1)))then
u(i)=dabs(u(i))*dsqrt(1.0d0+(u(i+1)/u(i))**2)
else
u(i)=dabs(u(i+1))*dsqrt(1.0d0+(u(i)/u(i+1))**2)
endif
12 continue
do 13 j=1,n
r(1,j)=r(1,j)+u(1)*v(j)
13 continue
do 14 i=1,k-1
call rotate(r,qt,n,np,i,r(i,i),-r(i+1,i))
14 continue
return
END
c
SUBROUTINE rsolv(a,n,np,d,b)
implicit none
INTEGER n,np
DOUBLE PRECISION a(np,np),b(n),d(n)
INTEGER i,j
DOUBLE PRECISION sum
b(n)=b(n)/d(n)
do 12 i=n-1,1,-1
sum=0.0d0
do 11 j=i+1,n
sum=sum+a(i,j)*b(j)
11 continue
b(i)=(b(i)-sum)/d(i)
12 continue
return
END
c
SUBROUTINE rotate(r,qt,n,np,i,a,b)
implicit none
INTEGER n,np,i
DOUBLE PRECISION a,b,r(np,np),qt(np,np)
INTEGER j
DOUBLE PRECISION c,fact,s,w,y
if(a.eq.0.0d0)then
c=0.0d0
s=dsign(1.0d0,b)
else if(dabs(a).gt.dabs(b))then
fact=b/a
c=dsign(1.0d0/dsqrt(1.0d0+fact**2),a)
s=fact*c
else
fact=a/b
s=dsign(1.0d0/dsqrt(1.0d0+fact**2),b)
c=fact*s
endif
do 11 j=i,n
y=r(i,j)
w=r(i+1,j)
r(i,j)=c*y-s*w
r(i+1,j)=s*y+c*w
11 continue
do 12 j=1,n
y=qt(i,j)
w=qt(i+1,j)
qt(i,j)=c*y-s*w
qt(i+1,j)=s*y+c*w
12 continue
return
END
subroutine xmprove(N,NP,a,b,x,mark)
implicit none
INTEGER i,j,idum,N,NP,indx(N),mark
double precision d,a(NP,NP),b(N),x(N),aa(NP,NP)
do 12 i=1,N
x(i)=b(i)
do 11 j=1,N
aa(i,j)=a(i,j)
11 continue
12 continue
call ludcmp(aa,N,NP,indx,d,mark)
if (mark .eq. 0) goto 20
call lubksb(aa,N,NP,indx,x)
call mprove(a,aa,N,NP,indx,b,x)
20 continue
return
END
SUBROUTINE mprove(a,alud,n,np,indx,b,x)
implicit none
INTEGER n,np,indx(n),NMAX
double precision a(np,np),alud(np,np),b(n),x(n)
PARAMETER (NMAX=500)
CU USES lubksb
INTEGER i,j
double precision r(NMAX)
DOUBLE PRECISION sdp
do 12 i=1,n
sdp=-b(i)
do 11 j=1,n
sdp=sdp+(a(i,j))*(x(j))
11 continue
r(i)=sdp
12 continue
call lubksb(alud,n,np,indx,r)
do 13 i=1,n
x(i)=x(i)-r(i)
13 continue
return
END
SUBROUTINE ludcmp(a,n,np,indx,d,mark)
implicit none
INTEGER n,np,indx(n),NMAX
double precision d,a(np,np),TINY
PARAMETER (NMAX=500,TINY=1.0d-20)
INTEGER i,imax,j,k,mark
double precision aamax,dum,sum,vv(NMAX)
mark=1
d=1.0d0
do 12 i=1,n
aamax=0.0d0
do 11 j=1,n
if (dabs(a(i,j)).gt.aamax) aamax=dabs(a(i,j))
11 continue
if (aamax.eq.0.0d0) then
! singular matrix
mark=0
return
end if
vv(i)=1.0d0/aamax
12 continue
do 19 j=1,n
do 14 i=1,j-1
sum=a(i,j)
do 13 k=1,i-1
sum=sum-a(i,k)*a(k,j)
13 continue
a(i,j)=sum
14 continue
aamax=0.0d0
do 16 i=j,n
sum=a(i,j)
do 15 k=1,j-1
sum=sum-a(i,k)*a(k,j)
15 continue
a(i,j)=sum
dum=vv(i)*dabs(sum)
if (dum.ge.aamax) then
imax=i
aamax=dum
endif
16 continue
if (j.ne.imax)then
do 17 k=1,n
dum=a(imax,k)
a(imax,k)=a(j,k)
a(j,k)=dum
17 continue
d=-d
vv(imax)=vv(j)
endif
indx(j)=imax
if(a(j,j).eq.0.0d0)a(j,j)=TINY
if(j.ne.n)then
dum=1.0d0/a(j,j)
do 18 i=j+1,n
a(i,j)=a(i,j)*dum
18 continue
endif
19 continue
return
END
SUBROUTINE lubksb(a,n,np,indx,b)
implicit none
INTEGER n,np,indx(n)
double precision a(np,np),b(n)
INTEGER i,ii,j,ll
double precision sum
ii=0
do 12 i=1,n
ll=indx(i)
sum=b(ll)
b(ll)=b(i)
if (ii.ne.0)then
do 11 j=ii,i-1
sum=sum-a(i,j)*b(j)
11 continue
else if (sum.ne.0.0d0) then
ii=i
endif
b(i)=sum
12 continue
do 14 i=n,1,-1
sum=b(i)
do 13 j=i+1,n
sum=sum-a(i,j)*b(j)
13 continue
b(i)=sum/a(i,i)
14 continue
return
END
@@ -0,0 +1,65 @@
subroutine cpbookkeeping(nunknowns,xvar,fequ,
& icall,iflargest)
implicit none
include 'cpnslasystem.h'
integer nunknowns,icall,iflargest
double precision xvar(nunknowns),fequ(nunknowns)
integer iGuCall,i,j,k
parameter(iGuCall=49)
!--------------------------------------------------------------------------
iflargest=1
do j=2,nunknowns
if(dabs(fequ(j)).gt.dabs(fequ(iflargest)))
& iflargest=j
enddo
if(numeval.eq.0.or.icall.eq.iGuCall)then
numeval=numeval+1
flargest(numeval)=dabs(fequ(iflargest))
do i=1,nunknowns
xevaluated(numeval,i)=xvar(i)
fevaluated(numeval,i)=fequ(i)
enddo
return
endif
do i=1,numeval
k=0
do j=1,nunknowns
if(dabs(xvar(j)-xevaluated(i,j)).gt.
& 1.0d-5*dabs(xvar(j)))k=1
enddo
if(k.eq.0)goto 500
enddo
if(numeval.lt.maxeval)then
numeval=numeval+1
flargest(numeval)=dabs(fequ(iflargest))
do i=1,nunknowns
xevaluated(numeval,i)=xvar(i)
fevaluated(numeval,i)=fequ(i)
enddo
return
endif
! replace a point
j=1
do i=2,numeval
if(flargest(j).lt.flargest(i))then
j=i
endif
enddo
if(dabs(fequ(iflargest)).lt.flargest(j))then
flargest(j)=dabs(fequ(iflargest))
do i=1,nunknowns
xevaluated(j,i)=xvar(i)
fevaluated(j,i)=fequ(i)
enddo
endif
return
! too close to the existing point i
500 if(dabs(fequ(iflargest)).lt.flargest(i))then
flargest(i)=dabs(fequ(iflargest))
do j=1,nunknowns
xevaluated(i,j)=xvar(j)
fevaluated(i,j)=fequ(j)
enddo
endif
return
end
+567
View File
@@ -0,0 +1,567 @@
C This file contains all the subroutines needed by the nonlinear solver broydn
C code modified based on on-line version in Numerical Recipes Website, Dec 28, 2004
SUBROUTINE cpbroydn(x0min,x,x0max,STPMX,n,
& fveccopy,funcv,TOLF,ierr)
implicit none
INTEGER n,NP,MAXITS,ierr
double precision x(n),EPS,TOLF,TOLMIN,TOLX,STPMX
double precision x0min(n),x0max(n),fveccopy(n)
LOGICAL check
! PARAMETER (NP=1000,MAXITS=250,EPS=1.0d-7,TOLF=1.0d-4,
! & TOLMIN=1.d-6,TOLX=EPS)
PARAMETER(NP=1000,MAXITS=250,EPS=1.0d-7,TOLX=EPS)
CU USES fdjac,funcv,lnsrch,qrdcmp,qrupdt,rsolv
INTEGER i,its,j,k
double precision den,f,fold,stpmax,sum,temp,test,c(NP),
& d(NP),fvcold(NP),g(NP),p(NP),qt(NP,NP),r(NP,NP),
& s(NP),t(NP),w(NP),xold(NP),fvec(NP)
LOGICAL restrt,sing,skip
EXTERNAL funcv
TOLMIN=TOLF*0.01d0
call funcv(n,x,fvec,f)
test=0.0d0
do 11 i=1,n
if(dabs(fvec(i)).gt.test)test=dabs(fvec(i))
11 continue
if(test.lt..01d0*TOLF)then
ierr=0
! check=.false.
return
endif
sum=0.0d0
do 12 i=1,n
sum=sum+x(i)*x(i)
12 continue
stpmax=STPMX*dmax1(dsqrt(sum),dble(n))
restrt=.true.
do 42 its=1,MAXITS
if(restrt)then
do i=1,n
if(x(i).lt.x0min(i).or.x(i).gt.x0max(i))then
ierr=1
return
endif
enddo
call cpfdjac(n,x,fvec,NP,r,funcv)
call cpqrdcmp(r,n,NP,c,d,sing)
! if(sing) pause 'singular Jacobian in broydn'
if(sing)then
ierr=2
return
end if
do 14 i=1,n
do 13 j=1,n
qt(i,j)=0.0d0
13 continue
qt(i,i)=1.0d0
14 continue
do 18 k=1,n-1
if(c(k).ne.0.0d0)then
do 17 j=1,n
sum=0.0d0
do 15 i=k,n
sum=sum+r(i,k)*qt(i,j)
15 continue
sum=sum/c(k)
do 16 i=k,n
qt(i,j)=qt(i,j)-sum*r(i,k)
16 continue
17 continue
endif
18 continue
do 21 i=1,n
r(i,i)=d(i)
do 19 j=1,i-1
r(i,j)=0.0d0
19 continue
21 continue
else
do 22 i=1,n
s(i)=x(i)-xold(i)
22 continue
do 24 i=1,n
sum=0.0d0
do 23 j=i,n
sum=sum+r(i,j)*s(j)
23 continue
t(i)=sum
24 continue
skip=.true.
do 26 i=1,n
sum=0.0d0
do 25 j=1,n
sum=sum+qt(j,i)*t(j)
25 continue
w(i)=fvec(i)-fvcold(i)-sum
if(dabs(w(i)).ge.EPS*(dabs(fvec(i))+
& dabs(fvcold(i))))then
skip=.false.
else
w(i)=0.0d0
endif
26 continue
if(.not.skip)then
do 28 i=1,n
sum=0.0d0
do 27 j=1,n
sum=sum+qt(i,j)*w(j)
27 continue
t(i)=sum
28 continue
den=0.0d0
do 29 i=1,n
den=den+s(i)*s(i)
29 continue
do 31 i=1,n
s(i)=s(i)/den
31 continue
call cpqrupdt(r,qt,n,NP,t,s)
do 32 i=1,n
if(r(i,i).eq.0.0d0) then
write(*,*) 'r singular in broydn'
end if
d(i)=r(i,i)
32 continue
endif
endif
do 34 i=1,n
sum=0.0d0
do 33 j=1,n
sum=sum+qt(i,j)*fvec(j)
33 continue
p(i)=-sum
34 continue
do 36 i=n,1,-1
sum=0.0d0
do 35 j=1,i
sum=sum-r(j,i)*p(j)
35 continue
g(i)=sum
36 continue
do 37 i=1,n
xold(i)=x(i)
fvcold(i)=fvec(i)
37 continue
fold=f
call cprsolv(r,n,NP,d,p)
! Gu modification starts
do 100 i=1,n
if(xold(i).lt.x0min(i).or.xold(i).gt.x0max(i))then
ierr=1
return
endif
100 continue
! Gu modification ends
call cplnsrch(n,xold,fold,g,p,x,f,
& stpmax,check,funcv,fvec)
test=0.0d0
do 38 i=1,n
if(dabs(fvec(i)).gt.test)test=dabs(fvec(i))
fveccopy(i)=fvec(i)
38 continue
if(test.lt.TOLF)then
ierr=0
! check=.false.
return
endif
if(check)then
if(restrt)then
ierr=3
return
else
test=0.0d0
den=dmax1(f,.5d0*dble(n))
do 39 i=1,n
temp=dabs(g(i))*dmax1(dabs(x(i)),1.0d0)/den
if(temp.gt.test)test=temp
39 continue
if(test.lt.TOLMIN)then
ierr=4
return
else
restrt=.true.
endif
endif
else
restrt=.false.
test=0.0d0
do 41 i=1,n
temp=(dabs(x(i)-xold(i)))/dmax1(dabs(x(i)),1.0d0)
if(temp.gt.test)test=temp
41 continue
if(test.lt.TOLX)then
ierr=4
! check=.true.
return
endif
endif
42 continue
ierr=5
return
END
SUBROUTINE cpfdjac(n,x,fvec,np,df,funcv)
implicit none
INTEGER n,np
double precision df(np,np),fvec(n),x(n),EPS
PARAMETER (EPS=1.0d-4)
CU USES funcv
INTEGER i,j,k
double precision h,temp,f(n),fsqsum
external funcv
do 12 j=1,n
temp=x(j)
h=EPS*dabs(temp)
if(h.eq.0.0d0)h=EPS
x(j)=temp+h
h=x(j)-temp
call funcv(n,x,f,fsqsum)
x(j)=temp
do 11 i=1,n
df(i,j)=(f(i)-fvec(i))/h
11 continue
12 continue
return
END
c
SUBROUTINE cplnsrch(n,xold,fold,g,p,x,f,
& stpmax,check,funcv,fvec)
implicit none
INTEGER n
LOGICAL check
DOUBLE PRECISION f,fold,stpmax,g(n),p(n),x(n),
*xold(n),ALF,TOLX,fvec(n)
PARAMETER (ALF=1.d-4,TOLX=1.d-7)
EXTERNAL funcv
CU USES funcv
INTEGER i
DOUBLE PRECISION a,alam,alam2,alamin,b,disc,
*f2,rhs1,rhs2,slope,sum,temp,test,tmplam
check=.false.
sum=0.0d0
do 11 i=1,n
sum=sum+p(i)*p(i)
11 continue
sum=dsqrt(sum)
if(sum.gt.stpmax)then
do 12 i=1,n
p(i)=p(i)*stpmax/sum
12 continue
endif
slope=0.0d0
do 13 i=1,n
slope=slope+g(i)*p(i)
13 continue
! if(slope.ge.0.0d0)pause 'roundoff problem in lnsrch'
test=0.0d0
do 14 i=1,n
temp=dabs(p(i))/dmax1(dabs(xold(i)),1.0d0)
if(temp.gt.test)test=temp
14 continue
alamin=TOLX/test
alam=1.0d0
1 continue
do 15 i=1,n
x(i)=xold(i)+alam*p(i)
15 continue
call funcv(n,x,fvec,f)
if(alam.lt.alamin)then
do 16 i=1,n
x(i)=xold(i)
16 continue
check=.true.
return
else if(f.le.fold+ALF*alam*slope)then
return
else
if(alam.eq.1.0d0)then
tmplam=-slope/(2.0d0*(f-fold-slope))
else
rhs1=f-fold-alam*slope
rhs2=f2-fold-alam2*slope
a=(rhs1/alam**2-rhs2/alam2**2)/(alam-alam2)
b=(-alam2*rhs1/alam**2+alam*rhs2/alam2**2)/
& (alam-alam2)
if(a.eq.0.0d0)then
tmplam=-slope/(2.0d0*b)
else
disc=b*b-3.0d0*a*slope
if(disc.lt.0.0d0) then
tmplam=0.5d0*alam
else if(b.le.0.0d0)then
tmplam=(-b+dsqrt(disc))/(3.0d0*a)
else
tmplam=-slope/(b+dsqrt(disc))
endif
endif
if(tmplam.gt..5d0*alam)tmplam=.5d0*alam
endif
endif
alam2=alam
f2=f
alam=dmax1(tmplam,.1d0*alam)
goto 1
END
c
SUBROUTINE cpqrdcmp(a,n,np,c,d,sing)
implicit none
INTEGER n,np
DOUBLE PRECISION a(np,np),c(n),d(n)
LOGICAL sing
INTEGER i,j,k
DOUBLE PRECISION scale,sigma,sum,tau
sing=.false.
do 17 k=1,n-1
scale=0.0d0
do 11 i=k,n
scale=dmax1(scale,dabs(a(i,k)))
11 continue
if(scale.eq.0.0d0)then
sing=.true.
c(k)=0.0d0
d(k)=0.0d0
else
do 12 i=k,n
a(i,k)=a(i,k)/scale
12 continue
sum=0.0d0
do 13 i=k,n
sum=sum+a(i,k)**2
13 continue
sigma=dsign(dsqrt(sum),a(k,k))
a(k,k)=a(k,k)+sigma
c(k)=sigma*a(k,k)
d(k)=-scale*sigma
do 16 j=k+1,n
sum=0.0d0
do 14 i=k,n
sum=sum+a(i,k)*a(i,j)
14 continue
tau=sum/c(k)
do 15 i=k,n
a(i,j)=a(i,j)-tau*a(i,k)
15 continue
16 continue
endif
17 continue
d(n)=a(n,n)
if(d(n).eq.0.0d0)sing=.true.
return
END
c
SUBROUTINE cpqrupdt(r,qt,n,np,u,v)
implicit none
INTEGER n,np
DOUBLE PRECISION r(np,np),qt(np,np),u(np),v(np)
CU USES rotate
INTEGER i,j,k
do 11 k=n,1,-1
if(u(k).ne.0.0d0)goto 1
11 continue
k=1
1 do 12 i=k-1,1,-1
call cprotate(r,qt,n,np,i,u(i),-u(i+1))
if(u(i).eq.0.0d0)then
u(i)=dabs(u(i+1))
else if(dabs(u(i)).gt.dabs(u(i+1)))then
u(i)=dabs(u(i))*dsqrt(1.0d0+(u(i+1)/u(i))**2)
else
u(i)=dabs(u(i+1))*dsqrt(1.0d0+(u(i)/u(i+1))**2)
endif
12 continue
do 13 j=1,n
r(1,j)=r(1,j)+u(1)*v(j)
13 continue
do 14 i=1,k-1
call cprotate(r,qt,n,np,i,r(i,i),-r(i+1,i))
14 continue
return
END
c
SUBROUTINE cprsolv(a,n,np,d,b)
implicit none
INTEGER n,np
DOUBLE PRECISION a(np,np),b(n),d(n)
INTEGER i,j
DOUBLE PRECISION sum
b(n)=b(n)/d(n)
do 12 i=n-1,1,-1
sum=0.0d0
do 11 j=i+1,n
sum=sum+a(i,j)*b(j)
11 continue
b(i)=(b(i)-sum)/d(i)
12 continue
return
END
c
SUBROUTINE cprotate(r,qt,n,np,i,a,b)
implicit none
INTEGER n,np,i
DOUBLE PRECISION a,b,r(np,np),qt(np,np)
INTEGER j
DOUBLE PRECISION c,fact,s,w,y
if(a.eq.0.0d0)then
c=0.0d0
s=dsign(1.0d0,b)
else if(dabs(a).gt.dabs(b))then
fact=b/a
c=dsign(1.0d0/dsqrt(1.0d0+fact**2),a)
s=fact*c
else
fact=a/b
s=dsign(1.0d0/dsqrt(1.0d0+fact**2),b)
c=fact*s
endif
do 11 j=i,n
y=r(i,j)
w=r(i+1,j)
r(i,j)=c*y-s*w
r(i+1,j)=s*y+c*w
11 continue
do 12 j=1,n
y=qt(i,j)
w=qt(i+1,j)
qt(i,j)=c*y-s*w
qt(i+1,j)=s*y+c*w
12 continue
return
END
subroutine cpxmprove(N,NP,a,b,x,mark)
implicit none
INTEGER i,j,idum,N,NP,indx(N),mark
double precision d,a(NP,NP),b(N),x(N),aa(NP,NP)
do 12 i=1,N
x(i)=b(i)
do 11 j=1,N
aa(i,j)=a(i,j)
11 continue
12 continue
call cpludcmp(aa,N,NP,indx,d,mark)
if (mark .eq. 0) goto 20
call cplubksb(aa,N,NP,indx,x)
call cpmprove(a,aa,N,NP,indx,b,x)
20 continue
return
END
SUBROUTINE cpmprove(a,alud,n,np,indx,b,x)
implicit none
INTEGER n,np,indx(n),NMAX
double precision a(np,np),alud(np,np),b(n),x(n)
PARAMETER (NMAX=500)
CU USES lubksb
INTEGER i,j
double precision r(NMAX)
DOUBLE PRECISION sdp
do 12 i=1,n
sdp=-b(i)
do 11 j=1,n
sdp=sdp+(a(i,j))*(x(j))
11 continue
r(i)=sdp
12 continue
call cplubksb(alud,n,np,indx,r)
do 13 i=1,n
x(i)=x(i)-r(i)
13 continue
return
END
SUBROUTINE cpludcmp(a,n,np,indx,d,mark)
implicit none
INTEGER n,np,indx(n),NMAX
double precision d,a(np,np),TINY
PARAMETER (NMAX=500,TINY=1.0d-20)
INTEGER i,imax,j,k,mark
double precision aamax,dum,sum,vv(NMAX)
mark=1
d=1.0d0
do 12 i=1,n
aamax=0.0d0
do 11 j=1,n
if (dabs(a(i,j)).gt.aamax) aamax=dabs(a(i,j))
11 continue
if (aamax.eq.0.0d0) then
! singular matrix
mark=0
return
end if
vv(i)=1.0d0/aamax
12 continue
do 19 j=1,n
do 14 i=1,j-1
sum=a(i,j)
do 13 k=1,i-1
sum=sum-a(i,k)*a(k,j)
13 continue
a(i,j)=sum
14 continue
aamax=0.0d0
do 16 i=j,n
sum=a(i,j)
do 15 k=1,j-1
sum=sum-a(i,k)*a(k,j)
15 continue
a(i,j)=sum
dum=vv(i)*dabs(sum)
if (dum.ge.aamax) then
imax=i
aamax=dum
endif
16 continue
if (j.ne.imax)then
do 17 k=1,n
dum=a(imax,k)
a(imax,k)=a(j,k)
a(j,k)=dum
17 continue
d=-d
vv(imax)=vv(j)
endif
indx(j)=imax
if(a(j,j).eq.0.0d0)a(j,j)=TINY
if(j.ne.n)then
dum=1.0d0/a(j,j)
do 18 i=j+1,n
a(i,j)=a(i,j)*dum
18 continue
endif
19 continue
return
END
SUBROUTINE cplubksb(a,n,np,indx,b)
implicit none
INTEGER n,np,indx(n)
double precision a(np,np),b(n)
INTEGER i,ii,j,ll
double precision sum
ii=0
do 12 i=1,n
ll=indx(i)
sum=b(ll)
b(ll)=b(i)
if (ii.ne.0)then
do 11 j=ii,i-1
sum=sum-a(i,j)*b(j)
11 continue
else if (sum.ne.0.0d0) then
ii=i
endif
b(i)=sum
12 continue
do 14 i=n,1,-1
sum=b(i)
do 13 j=i+1,n
sum=sum-a(i,j)*b(j)
13 continue
b(i)=sum/a(i,i)
14 continue
return
END
+315
View File
@@ -0,0 +1,315 @@
subroutine cpfixedpoint(funcnleq1,x0min,x0ori,xp,
& x0max,fequ,nunknowns,TOLF,stpmax,iwhichsolver)
implicit none
include 'cpnslasystem.h'
!-------- Inputs ---------------------------------------
! nunknowns: The number of unknowns to be solved
! x0ori(1:nunknowns): initial guess for the unknowns
! x0min(1:nunknowns): lower bound of the solution
! x0max(1:nunknowns): upper bound of the solution
! stpmax: the maximum length of the steps allowed to prevent search into
! undefined region.
! TOLF: Error tolerance
! funcnleq1: the subroutine name for the nonlinear system
integer nunknowns
double precision x0min(1:nunknowns),x0ori(1:nunknowns),
& x0max(1:nunknowns),TOLF,stpmax
! --------- Outputs -------------------------------------
! fequ(1:nunknowns): function values at the last step of iteration
! xp(1:nunknowns): final solutions or solutions not worse than x0ori
! iwhichsolver: =1,2,3,4 successful
! =-9999 failed, best solution returned
integer iwhichsolver
double precision fequ(1:nunknowns),xp(1:nunknowns)
! ---------Local variables --------------------------------
integer i,j,k,maxiter,notfound,ncount,ierr,
& ismallest,iGuCall
double precision swap,x1,x2,f1,f2,fsqsumold,
& fsqsumnew,xpold(nunknowns),fequold(nunknowns),
& gfuncsum(nunknowns),deltax(nunknowns),
& xpder(nunknowns),fjacob(nunknowns,nunknowns),
& fjacobcopy(nunknowns,nunknowns),fsqsum
logical check
parameter(maxiter=200,notfound=-9999,iGuCall=49)
integer iselect(300*maxiter)
logical resetran2
common /cpran2reset/resetran2
save /cpran2reset/
external funcnleq1
!-----------------------------------------------------------
resetran2=.true.
do i=1,nunknowns
xp(i)=x0ori(i)
enddo
iwhichsolver=notfound
numeval=0
!--------------------------------------------------------------
!Plain fixed-point method. Fixed-point method 1
do i=1,maxiter
call funcnleq1(nunknowns,xp,fequ,fsqsum)
call cpbookkeeping(nunknowns,xp,fequ,iGuCall,k)
if(dabs(fequ(k)).lt.TOLF)then
iwhichsolver=1
return
endif
do j=1,nunknowns
xp(j)=xp(j)-fequ(j)
if(xp(j).lt.x0min(j).or.xp(j).gt.x0max(j))then
call reinitialization(x0min(j),x0ori(j),
& x0max(j),xp(j),50000)
endif
enddo
enddo
!_____________________________________________________________________
!try fixed-point method 2
do j=1,nunknowns
call reinitialization(x0min(j),x0ori(j),
& x0max(j),xp(j),10000)
enddo
call funcnleq1(nunknowns,xp,fequ,fsqsum)
call cpbookkeeping(nunknowns,xp,fequ,iGuCall,k)
do i=1,nunknowns
xp(i)=xp(i)-fequ(i)
if(xp(i).lt.x0min(i).or.xp(i).gt.x0max(i))then
call reinitialization(x0min(i),x0ori(i),
& x0max(i),xp(i),2000)
endif
enddo
do i=1,maxiter
call funcnleq1(nunknowns,xp,fequ,fsqsum)
call cpbookkeeping(nunknowns,xp,fequ,iGuCall,k)
if(dabs(fequ(k)).lt.TOLF)then
iwhichsolver=2
return
endif
do j=1,nunknowns
ierr=0
x1=xevaluated(numeval-1,j)
f1=x1-fevaluated(numeval-1,j)
x2=xevaluated(numeval,j)
f2=x2-fevaluated(numeval,j)
if(dabs(f2-f1-x2+x1).gt.1.0d-20)then
ierr=1
xp(j)=(x1*(f2-f1)-f1*(x2-x1))/
& (f2-f1-x2+x1)
if(xp(j).le.x0min(j).or.xp(j)
& .ge.x0max(j))then
ierr=0
endif
endif
if(ierr.le.0.and.numeval.ge.3)then
! haven't found a usable new point yet, first try the opposite sign point
ncount=0
do k=1,numeval-2
if((fevaluated(k,j)*fevaluated(numeval,j))
& .lt.0.0d0)then
ncount=ncount+1
iselect(ncount)=k
endif
enddo
if(ncount.gt.0)then
! there are points at different sides of the zero.
ismallest=1
do k=2,ncount
if(dabs(xevaluated(iselect(k),j)-x2).lt.
& dabs(xevaluated(iselect(ismallest),j)-x2))then
ismallest=k
endif
enddo
ierr=1
x1=xevaluated(iselect(ismallest),j)
f1=x1-fevaluated(iselect(ismallest),j)
xp(j)=(x1*(f2-f1)-f1*(x2-x1))/
& (f2-f1-x2+x1)
else
! all at the same sides of the zero.
do k=1,numeval-2
x1=xevaluated(k,j)
f1=x1-fevaluated(k,j)
if(dabs(f2-f1-x2+x1).gt.1.0d-10)then
xp(j)=(x1*(f2-f1)-f1*(x2-x1))/
& (f2-f1-x2+x1)
if(xp(j).gt.x0min(j).and.xp(j).lt.x0max(j))then
ierr=1
endif
endif
if(ierr.eq.1)goto 10
enddo
10 continue
endif
endif
if(ierr.eq.0)then
call reinitialization(x0min(j),
& xevaluated(numeval,j),x0max(j),xp(j),1000)
endif
enddo
ierr=0
do k=1,nunknowns
if(xp(k).ne.xevaluated(numeval,k))ierr=1
enddo
if(ierr.eq.0)then
do k=1,nunknowns
call reinitialization(x0min(k),
& xevaluated(numeval,k),x0max(k),xp(k),25000)
enddo
endif
enddo
!__________________________________________________________________
!Try fixed-point method 3
do i=1,nunknowns
xp(i)=x0ori(i)+1.0d-6
if(xp(i).lt.x0min(i).or.xp(i).gt.x0max(i))then
call reinitialization(x0min(i),x0ori(i),
& x0max(i),xp(i),250910)
endif
enddo
do j=1,maxiter
call funcnleq1(nunknowns,xp,fequ,fsqsum)
call cpbookkeeping(nunknowns,xp,fequ,iGuCall,k)
if(dabs(fequ(k)).lt.TOLF)then
iwhichsolver=3
return
endif
do i=1,nunknowns
xp(i)=xp(i)-fequ(i)
if(xp(i).lt.x0min(i).or.xp(i).gt.x0max(i))then
call reinitialization(x0min(i),x0ori(i),
& x0max(i),xp(i),25500)
endif
enddo
call funcnleq1(nunknowns,xp,fequ,fsqsum)
call cpbookkeeping(nunknowns,xp,fequ,iGuCall,k)
if(dabs(fequ(k)).lt.TOLF)then
iwhichsolver=3
return
endif
do i=1,nunknowns
if(fevaluated(numeval,i).eq.
& fevaluated(numeval-1,i))then
x1=(xevaluated(numeval,i)+
& xevaluated(numeval-1,i))/2.0d0
call reinitialization(x0min(i),x1,
& x0max(i),xp(i),35678)
else
xp(i)=(xevaluated(numeval,i)*fevaluated(numeval-1,i)
& -xevaluated(numeval-1,i)*fevaluated(numeval,i))/
& (fevaluated(numeval-1,i)-fevaluated(numeval,i))
if(xp(i).lt.x0min(i).or.xp(i).gt.x0max(i))then
call reinitialization(x0min(i),x0ori(i),
& x0max(i),xp(i),45678)
endif
endif
enddo
enddo
!------------------------------------------------------------
!Try fixed-point method 4
!11 call funcnleq1(nunknowns,xp,fequ,fsqsumold)
! call cpbookkeeping(nunknowns,xp,fequ,iGuCall,i)
fsqsumold=0.0d0
do i=1,nunknowns
xpold(i)=xevaluated(1,i)
fequold(i)=fevaluated(1,i)
fsqsumold=fsqsumold+fequold(i)*fequold(i)
enddo
fsqsumold=0.5d0*fsqsumold
do k=1,maxiter
do j=1,nunknowns
do i=1,nunknowns
xpder(i)=xpold(i)
enddo
if(dabs(fequold(j)).lt.1.0d-10)then
xpder(j)=xpold(j)+1.0d-5
else
xpder(j)=xpold(j)-fequold(j)
endif
if(xpder(j).lt.x0min(j).or.xpder(j).
& gt.x0max(j))then
call reinitialization(x0min(j),xpold(j),
& x0max(j),xpder(j),89000)
endif
call funcnleq1(nunknowns,xpder,fequ,fsqsumnew)
call cpbookkeeping(nunknowns,xpder,fequ,iGuCall,i)
if(dabs(fequ(i)).lt.TOLF)then
iwhichsolver=4
return
endif
do i=1,nunknowns
fjacob(i,j)=(fequ(i)-fequold(i))/
& (xpder(j)-xpold(j))
fjacobcopy(i,j)=fjacob(i,j)
enddo
gfuncsum(j)=(fsqsumnew-fsqsumold)/
& (xpder(j)-xpold(j))
enddo
call cpxmprove(nunknowns,nunknowns,
& fjacob,fequold,deltax,ierr)
!if ierr = 0, matrix is singular. ierr = 1, everything is ok.
if(ierr.eq.0)then
call adsor(fjacobcopy,nunknowns,nunknowns,
& fequold,deltax,ierr)
if(ierr.ne.1)ierr=0
endif
if(ierr.ne.0)then
do i=1,nunknowns
deltax(i)=-deltax(i)
enddo
call cplnsrch(nunknowns,xpold,fsqsumold,
& gfuncsum,deltax,xp,fsqsumnew,stpmax,
& check,funcnleq1,fequ)
if(check.eq..true..or.check.eq..TRUE.)then
do i=1,nunknowns
call reinitialization(x0min(i),xpold(i),
& x0max(i),xp(i),6678)
enddo
endif
do i=1,nunknowns
if(xp(i).lt.x0min(i).or.xp(i).gt.x0max(i))then
call reinitialization(x0min(i),x0ori(i),
& x0max(i),xp(i),678)
endif
enddo
else
do i=1,nunknowns
call reinitialization(x0min(i),x0ori(i),
& x0max(i),xp(i),75678)
enddo
endif
call funcnleq1(nunknowns,xp,fequ,fsqsumold)
call cpbookkeeping(nunknowns,xp,fequ,iGuCall,i)
if(dabs(fequ(i)).lt.TOLF)then
iwhichsolver=4
return
endif
do i=1,nunknowns
xpold(i)=xp(i)
fequold(i)=fequ(i)
enddo
enddo
!_____________________________________________________________
!If all four methods failed, choose the best xp
do i=1,numeval
do k=i+1,numeval
if(flargest(k).lt.flargest(i))then
swap=flargest(k)
flargest(k)=flargest(i)
flargest(i)=swap
do ncount=1,nunknowns
swap=xevaluated(k,ncount)
xevaluated(k,ncount)=xevaluated(i,ncount)
xevaluated(i,ncount)=swap
swap=fevaluated(k,ncount)
fevaluated(k,ncount)=fevaluated(i,ncount)
fevaluated(i,ncount)=swap
enddo
endif
enddo
enddo
! Best solution found so far
do i=1,nunknowns
xp(i)=xevaluated(1,i)
fequ(i)=fevaluated(1,i)
enddo
return
end subroutine cpfixedpoint
@@ -0,0 +1,116 @@
subroutine cpnonsyssolver(funcnleq1,fmin_funcnleq1,
& f1dim_funcnleq1,x0min,x0ori,xp,x0max,fp,
& nunknowns,iwhichsolver)
implicit none
integer nunknowns,iwhichsolver
double precision x0min(nunknowns),x0ori(nunknowns),
& xp(nunknowns),x0max(nunknowns),fp(nunknowns)
!-------- Specified values ---------------------------------------
!funcnleq1: the subroutine that calculates the functional values of the
! the nonlinear system in the following form:
! funcnleq1(nunknowns,xp,fp,fsqsum)
!fmin_funcnleq1: the subroutine that calls funcnleq1 and returns fsqsum (half
! of the sum of the squared functional values of the nonlinear system)
! fmin_funcnleq1(nunknowns,xp,fsqsum)
!f1dim_funcnleq1: a function subroutine that returns fsqsum
! f1dim_funcnleq1(xp)
! nunknowns: The number of unknowns to be solved
! x0ori(1:nunknowns): initial guess for the unknowns
! x0min(1:nunknowns): lower bound of the solution
! x0max(1:nunknowns): upper bound of the solution
! --------- Calculated values -------------------------------------
! fp(1:nunknowns): function values at the last step of iteration
! xp(1:nunknowns): final solutions
! iwhichsolver:
! =1 solved by plain fixed point method 1
! =2 solved by fixed point method 2
! =3 solved by fixed point method 3
! =4 solved by fixed point method 4
! =6 solved by broydn
! =7 Solved by multiobjective minimization.
! =-9999 Best approximation returned. Solution may not be accurate.
! --------- Local variables ---------------------------------------
double precision x0(nunknowns),TOLF,stpmax,scldstpmax,
& sum,tb,tp,xb(nunknowns),fb(nunknowns),fsqsum,
& f1dim_funcnleq1
integer i,irepeat,maxrepeats,IERR,notfound
intrinsic dble
parameter(maxrepeats=100,notfound=-9999,TOLF=1.0d-10)
external funcnleq1,fmin_funcnleq1,f1dim_funcnleq1
!-------------------------------------------------------------------
stpmax=0.0d0
sum=0.0d0
do i=1, nunknowns
x0(i)=x0ori(i)
sum=sum+x0ori(i)*x0ori(i)
stpmax=stpmax+
& (x0min(i)-x0max(i))*(x0min(i)-x0max(i))
enddo
stpmax=dsqrt(stpmax)/4.0d0
scldstpmax=stpmax/dmax1(dsqrt(sum),dble(nunknowns))
! In Numerical Recipes, scldstpmax (STPMX) is 100
scldstpmax=dmax1(100.0d0,scldstpmax)
iwhichsolver=notfound
do irepeat=1,maxrepeats
call cpfixedpoint(funcnleq1,x0min,x0,xp,
& x0max,fp,nunknowns,TOLF,stpmax,iwhichsolver)
if(iwhichsolver.ne.notfound)return
tp=dabs(fp(1))
xb(1)=xp(1)
do i=2,nunknowns
if(dabs(fp(i)).gt.tp)tp=dabs(fp(i))
xb(i)=xp(i)
enddo
call cpbroydn(x0min,xb,x0max,scldstpmax,nunknowns,
& fb,funcnleq1,TOLF,IERR)
call funcnleq1(nunknowns,xb,fb,fsqsum)
tb=dabs(fb(1))
do i=2,nunknowns
if(dabs(fb(i)).gt.tb)tb=dabs(fb(i))
enddo
do i=1,nunknowns
if(xb(i).lt.x0min(i).or.xb(i).gt.x0max(i))then
tb=1.0d+100
endif
enddo
if(tb.lt.tp)then
do i=1,nunknowns
xp(i)=xb(i)
fp(i)=fb(i)
enddo
if(tb.lt.TOLF)then
iwhichsolver=6
return
endif
endif
fsqsum=0.0d0
do i=1,nunknowns
fsqsum=fsqsum+fp(i)*fp(i)
enddo
tp=fsqsum
call cpnongradopt(nunknowns,fmin_funcnleq1,
& f1dim_funcnleq1,xp,x0min,x0max,TOLF,fsqsum)
if(dabs(tp-fsqsum).gt.TOLF)then
call cpRepeatCompassSearch(nunknowns,xp,fsqsum,
& x0min,x0max,fmin_funcnleq1,f1dim_funcnleq1,
& TOLF)
endif
call funcnleq1(nunknowns,xp,fp,fsqsum)
tp=dabs(fp(1))
do i=2,nunknowns
if(dabs(fp(i)).gt.tp)tp=dabs(fp(i))
enddo
if(tp.lt.TOLF)then
iwhichsolver=7
return
endif
IERR=0
do i=1,nunknowns
if(dabs(xp(i)-x0(i)).gt.TOLF)IERR=1
enddo
if(IERR.eq.0)return
do i=1,nunknowns
x0(i)=xp(i)
enddo
enddo
end subroutine cpnonsyssolver
@@ -0,0 +1,18 @@
!------------------ Common Blocks -------------------------
integer numeval,maxndim,maxeval
parameter(maxndim=1000,maxeval=15000)
double precision xevaluated,fevaluated,flargest
common /cpFuncvRegresInteg/numeval
common /cpFuncvRegresDble/xevaluated(1:maxeval,1:maxndim),
& fevaluated(1:maxeval,1:maxndim),
& flargest(1:maxeval)
save /cpFuncvRegresInteg/,/cpFuncvRegresDble/
! numeval: the number of times that the system is evaluated so far
! iflargest: the index of the largest function for the latest evaluation
! xevaluated: the positions where the system is evaluated
! fevaluated: the function values at xevaluated
! flargest: the largest absolute function value
! maxndim: the maximum allowable dimensions of the system
! maxeval: the maximum allowable number of function evaluations
!--------------------------------------------------------------
+370
View File
@@ -0,0 +1,370 @@
subroutine fixedpoint(funcnleq1,x0min,x0ori,xp,
& x0max,fequ,nunknowns,TOLF,stpmax,iwhichsolver)
implicit none
include 'nslasystem.h'
!-------- Inputs ---------------------------------------
! nunknowns: The number of unknowns to be solved
! x0ori(1:nunknowns): initial guess for the unknowns
! x0min(1:nunknowns): lower bound of the solution
! x0max(1:nunknowns): upper bound of the solution
! stpmax: the maximum length of the steps allowed to prevent search into
! undefined region.
! TOLF: Error tolerance
! funcnleq1: the subroutine name for the nonlinear system
integer nunknowns
double precision x0min(1:nunknowns),x0ori(1:nunknowns),
& x0max(1:nunknowns),TOLF,stpmax
! --------- Outputs -------------------------------------
! fequ(1:nunknowns): function values at the last step of iteration
! xp(1:nunknowns): final solutions or solutions not worse than x0ori
! iwhichsolver: =0,1,2,3,4 successful
! =-9999 failed, best solution returned
integer iwhichsolver
double precision fequ(1:nunknowns),xp(1:nunknowns)
! ---------Local variables --------------------------------
integer i,j,k,n,maxiter,notfound,ncount,ierr,
& ismallest,iGuCall
double precision swap,x1,x2,f1,f2,fsqsumold,
& fsqsumnew,xpold(nunknowns),fequold(nunknowns),
& gfuncsum(nunknowns),deltax(nunknowns),
& xpder(nunknowns),fjacob(nunknowns,nunknowns),
& fjacobcopy(nunknowns,nunknowns),fsqsum,term
logical check
parameter(maxiter=200,notfound=-9999,iGuCall=49)
integer iselect(300*maxiter)
logical resetran2
common /ran2reset/resetran2
save /ran2reset/
external funcnleq1
!-----------------------------------------------------------
resetran2=.true.
do i=1,nunknowns
xp(i)=x0ori(i)
enddo
iwhichsolver=notfound
numeval=0
!--------------------------------------------------------------
!Plain fixed-point method. Fixed-point method 1
do i=1,maxiter
call funcnleq1(nunknowns,xp,fequ,fsqsum)
call bookkeeping(nunknowns,xp,fequ,iGuCall,k)
if(dabs(fequ(k)).lt.TOLF)then
iwhichsolver=1
return
endif
do j=1,nunknowns
xp(j)=xp(j)-fequ(j)
if(xp(j).lt.x0min(j).or.xp(j).gt.x0max(j))then
call reinitialization(x0min(j),x0ori(j),
& x0max(j),xp(j),50000)
endif
enddo
enddo
!_____________________________________________________________________
!try approximation to the newton method, iwhichsolver=0. this would work
!if the equations are independent
1 do i=1,nunknowns
xp(i)=x0ori(i)
enddo
do i=1,maxiter
call funcnleq1(nunknowns,xp,fequ,fsqsum)
n=0
do j=1,nunknowns
if(dabs(fequ(j)).gt.0.0d0)then
else
n=1
call reinitialization(x0min(j),x0ori(j),
& x0max(j),xp(j),50000)
endif
enddo
if(n.ne.0)goto 2
call bookkeeping(nunknowns,xp,fequ,iGuCall,k)
if(dabs(fequ(k)).lt.TOLF)then
iwhichsolver=0
return
endif
do j=1,nunknowns
xpold(j)=xp(j)
fequold(j)=fequ(j)
xp(j)=xp(j)+fequ(j)
enddo
call funcnleq1(nunknowns,xp,fequ,fsqsum)
do j=1,nunknowns
if(dabs(fequ(j)).gt.0.0d0)then
else
n=1
call reinitialization(x0min(j),x0ori(j),
& x0max(j),xp(j),50000)
endif
enddo
if(n.ne.0)goto 2
call bookkeeping(nunknowns,xp,fequ,iGuCall,k)
do j=1,nunknowns
if(fequ(j).ne.fequold(j))then
xpder(j)=fequold(j)/(fequ(j)-fequold(j))
else
xpder(j)=0.0d0
endif
xp(j)=xpold(j)-xpder(j)*fequold(j)
if(xp(j).lt.x0min(j).or.xp(j).gt.x0max(j))then
call reinitialization(x0min(j),x0ori(j),
& x0max(j),xp(j),50000)
endif
enddo
2 continue
enddo
!_____________________________________________________________________
!try fixed-point method 2
do j=1,nunknowns
call reinitialization(x0min(j),x0ori(j),
& x0max(j),xp(j),10000)
enddo
call funcnleq1(nunknowns,xp,fequ,fsqsum)
call bookkeeping(nunknowns,xp,fequ,iGuCall,k)
do i=1,nunknowns
xp(i)=xp(i)-fequ(i)
if(xp(i).lt.x0min(i).or.xp(i).gt.x0max(i))then
call reinitialization(x0min(i),x0ori(i),
& x0max(i),xp(i),2000)
endif
enddo
do i=1,maxiter
call funcnleq1(nunknowns,xp,fequ,fsqsum)
call bookkeeping(nunknowns,xp,fequ,iGuCall,k)
if(dabs(fequ(k)).lt.TOLF)then
iwhichsolver=2
return
endif
do j=1,nunknowns
ierr=0
x1=xevaluated(numeval-1,j)
f1=x1-fevaluated(numeval-1,j)
x2=xevaluated(numeval,j)
f2=x2-fevaluated(numeval,j)
if(dabs(f2-f1-x2+x1).gt.1.0d-20)then
ierr=1
xp(j)=(x1*(f2-f1)-f1*(x2-x1))/
& (f2-f1-x2+x1)
if(xp(j).le.x0min(j).or.xp(j)
& .ge.x0max(j))then
ierr=0
endif
endif
if(ierr.le.0.and.numeval.ge.3)then
! haven't found a usable new point yet, first try the opposite sign point
ncount=0
do k=1,numeval-2
if((fevaluated(k,j)*fevaluated(numeval,j))
& .lt.0.0d0)then
ncount=ncount+1
iselect(ncount)=k
endif
enddo
if(ncount.gt.0)then
! there are points at different sides of the zero.
ismallest=1
do k=2,ncount
if(dabs(xevaluated(iselect(k),j)-x2).lt.
& dabs(xevaluated(iselect(ismallest),j)-x2))then
ismallest=k
endif
enddo
ierr=1
x1=xevaluated(iselect(ismallest),j)
f1=x1-fevaluated(iselect(ismallest),j)
xp(j)=(x1*(f2-f1)-f1*(x2-x1))/
& (f2-f1-x2+x1)
else
! all at the same sides of the zero.
do k=1,numeval-2
x1=xevaluated(k,j)
f1=x1-fevaluated(k,j)
if(dabs(f2-f1-x2+x1).gt.1.0d-10)then
xp(j)=(x1*(f2-f1)-f1*(x2-x1))/
& (f2-f1-x2+x1)
if(xp(j).gt.x0min(j).and.xp(j).lt.x0max(j))then
ierr=1
endif
endif
if(ierr.eq.1)goto 10
enddo
10 continue
endif
endif
if(ierr.eq.0)then
call reinitialization(x0min(j),
& xevaluated(numeval,j),x0max(j),xp(j),1000)
endif
enddo
ierr=0
do k=1,nunknowns
if(xp(k).ne.xevaluated(numeval,k))ierr=1
enddo
if(ierr.eq.0)then
do k=1,nunknowns
call reinitialization(x0min(k),
& xevaluated(numeval,k),x0max(k),xp(k),25000)
enddo
endif
enddo
!__________________________________________________________________
!Try fixed-point method 3
do i=1,nunknowns
xp(i)=x0ori(i)+1.0d-6
if(xp(i).lt.x0min(i).or.xp(i).gt.x0max(i))then
call reinitialization(x0min(i),x0ori(i),
& x0max(i),xp(i),250910)
endif
enddo
do j=1,maxiter
call funcnleq1(nunknowns,xp,fequ,fsqsum)
call bookkeeping(nunknowns,xp,fequ,iGuCall,k)
if(dabs(fequ(k)).lt.TOLF)then
iwhichsolver=3
return
endif
do i=1,nunknowns
xp(i)=xp(i)-fequ(i)
if(xp(i).lt.x0min(i).or.xp(i).gt.x0max(i))then
call reinitialization(x0min(i),x0ori(i),
& x0max(i),xp(i),25500)
endif
enddo
call funcnleq1(nunknowns,xp,fequ,fsqsum)
call bookkeeping(nunknowns,xp,fequ,iGuCall,k)
if(dabs(fequ(k)).lt.TOLF)then
iwhichsolver=3
return
endif
do i=1,nunknowns
if(fevaluated(numeval,i).eq.
& fevaluated(numeval-1,i))then
x1=(xevaluated(numeval,i)+
& xevaluated(numeval-1,i))/2.0d0
call reinitialization(x0min(i),x1,
& x0max(i),xp(i),35678)
else
xp(i)=(xevaluated(numeval,i)*fevaluated(numeval-1,i)
& -xevaluated(numeval-1,i)*fevaluated(numeval,i))/
& (fevaluated(numeval-1,i)-fevaluated(numeval,i))
if(xp(i).lt.x0min(i).or.xp(i).gt.x0max(i))then
call reinitialization(x0min(i),x0ori(i),
& x0max(i),xp(i),45678)
endif
endif
enddo
enddo
!------------------------------------------------------------
!Try fixed-point method 4
!11 call funcnleq1(nunknowns,xp,fequ,fsqsumold)
! call bookkeeping(nunknowns,xp,fequ,iGuCall,i)
fsqsumold=0.0d0
do i=1,nunknowns
xpold(i)=xevaluated(numeval,i)
fequold(i)=fevaluated(numeval,i)
fsqsumold=fsqsumold+fequold(i)*fequold(i)
enddo
term=fsqsumold
do k=1,maxiter/5
do j=1,nunknowns
do i=1,nunknowns
xpder(i)=xpold(i)
enddo
if(dabs(fequold(j)).lt.1.0d-10)then
xpder(j)=xpold(j)+1.0d-5
else
xpder(j)=xpold(j)-fequold(j)
endif
if(xpder(j).lt.x0min(j).or.xpder(j).
& gt.x0max(j))then
call reinitialization(x0min(j),xpold(j),
& x0max(j),xpder(j),89000)
endif
call funcnleq1(nunknowns,xpder,fequ,fsqsumnew)
call bookkeeping(nunknowns,xpder,fequ,iGuCall,i)
if(dabs(fequ(i)).lt.TOLF)then
iwhichsolver=4
return
endif
do i=1,nunknowns
fjacob(i,j)=(fequ(i)-fequold(i))/
& (xpder(j)-xpold(j))
fjacobcopy(i,j)=fjacob(i,j)
enddo
gfuncsum(j)=(fsqsumnew-fsqsumold)/
& (xpder(j)-xpold(j))
enddo
call xmprove(nunknowns,nunknowns,
& fjacob,fequold,deltax,ierr)
!if ierr = 0, matrix is singular. ierr = 1, everything is ok.
if(ierr.eq.0)then
call adsor(fjacobcopy,nunknowns,nunknowns,
& fequold,deltax,ierr)
if(ierr.ne.1)ierr=0
endif
if(ierr.ne.0)then
do i=1,nunknowns
deltax(i)=-deltax(i)
enddo
call lnsrch(nunknowns,xpold,fsqsumold,
& gfuncsum,deltax,xp,fsqsumnew,stpmax,
& check,funcnleq1,fequ)
if(check.eq..true..or.check.eq..TRUE.)then
do i=1,nunknowns
call reinitialization(x0min(i),xpold(i),
& x0max(i),xp(i),6678)
enddo
endif
do i=1,nunknowns
if(xp(i).lt.x0min(i).or.xp(i).gt.x0max(i))then
call reinitialization(x0min(i),x0ori(i),
& x0max(i),xp(i),678)
endif
enddo
else
do i=1,nunknowns
call reinitialization(x0min(i),x0ori(i),
& x0max(i),xp(i),75678)
enddo
endif
call funcnleq1(nunknowns,xp,fequ,fsqsumold)
call bookkeeping(nunknowns,xp,fequ,iGuCall,i)
if(dabs(fequ(i)).lt.TOLF)then
iwhichsolver=4
return
endif
do i=1,nunknowns
xpold(i)=xp(i)
fequold(i)=fequ(i)
enddo
if(fsqsumold.ge.term)goto 30
term=fsqsumold
enddo
!_____________________________________________________________
!If all four methods failed, choose the best xp
30 do i=1,numeval
do k=i+1,numeval
if(flargest(k).lt.flargest(i))then
swap=flargest(k)
flargest(k)=flargest(i)
flargest(i)=swap
do ncount=1,nunknowns
swap=xevaluated(k,ncount)
xevaluated(k,ncount)=xevaluated(i,ncount)
xevaluated(i,ncount)=swap
swap=fevaluated(k,ncount)
fevaluated(k,ncount)=fevaluated(i,ncount)
fevaluated(i,ncount)=swap
enddo
endif
enddo
enddo
! Best solution found so far
do i=1,nunknowns
xp(i)=xevaluated(1,i)
fequ(i)=fevaluated(1,i)
enddo
return
end subroutine fixedpoint
+115
View File
@@ -0,0 +1,115 @@
subroutine nonsyssolver(funcnleq1,fmin_funcnleq1,
& f1dim_funcnleq1,x0min,x0ori,xp,x0max,fp,
& nunknowns,iwhichsolver)
implicit none
integer nunknowns,iwhichsolver
double precision x0min(nunknowns),x0ori(nunknowns),
& xp(nunknowns),x0max(nunknowns),fp(nunknowns)
!-------- Specified values ---------------------------------------
!funcnleq1: the subroutine that calculates the functional values of the
! the nonlinear system in the following form:
! funcnleq1(nunknowns,xp,fp,fsqsum)
!fmin_funcnleq1: the subroutine that calls funcnleq1 and returns fsqsum
! fmin_funcnleq1(nunknowns,xp,fsqsum)
!f1dim_funcnleq1: a function subroutine that returns fsqsum
! f1dim_funcnleq1(xp)
! nunknowns: The number of unknowns to be solved
! x0ori(1:nunknowns): initial guess for the unknowns
! x0min(1:nunknowns): lower bound of the solution
! x0max(1:nunknowns): upper bound of the solution
! --------- Calculated values -------------------------------------
! fp(1:nunknowns): function values at the last step of iteration
! xp(1:nunknowns): final solutions
! iwhichsolver:
! =1 solved by plain fixed point method 1
! =2 solved by fixed point method 2
! =3 solved by fixed point method 3
! =4 solved by fixed point method 4
! =6 solved by broydn
! =7 Solved by multiobjective minimization.
! =-9999 Best approximation returned. Solution may not be accurate.
! --------- Local variables ---------------------------------------
double precision x0(nunknowns),TOLF,stpmax,scldstpmax,
& sum,tb,tp,xb(nunknowns),fb(nunknowns),fsqsum,
& f1dim_funcnleq1
integer i,irepeat,maxrepeats,IERR,notfound
intrinsic dble
parameter(maxrepeats=100,notfound=-9999,TOLF=1.0d-7)
external funcnleq1,fmin_funcnleq1,f1dim_funcnleq1
!-------------------------------------------------------------------
stpmax=0.0d0
sum=0.0d0
do i=1, nunknowns
x0(i)=x0ori(i)
sum=sum+x0ori(i)*x0ori(i)
stpmax=stpmax+
& (x0min(i)-x0max(i))*(x0min(i)-x0max(i))
enddo
stpmax=dsqrt(stpmax)/4.0d0
scldstpmax=stpmax/dmax1(dsqrt(sum),dble(nunknowns))
! In Numerical Recipes, scldstpmax (STPMX) is 100
scldstpmax=dmax1(100.0d0,scldstpmax)
iwhichsolver=notfound
do irepeat=1,maxrepeats
call fixedpoint(funcnleq1,x0min,x0,xp,
& x0max,fp,nunknowns,TOLF,stpmax,iwhichsolver)
if(iwhichsolver.ne.notfound)return
tp=dabs(fp(1))
xb(1)=xp(1)
do i=2,nunknowns
if(dabs(fp(i)).gt.tp)tp=dabs(fp(i))
xb(i)=xp(i)
enddo
call broydn(x0min,xb,x0max,scldstpmax,nunknowns,
& fb,funcnleq1,TOLF,IERR)
call funcnleq1(nunknowns,xb,fb,fsqsum)
tb=dabs(fb(1))
do i=2,nunknowns
if(dabs(fb(i)).gt.tb)tb=dabs(fb(i))
enddo
do i=1,nunknowns
if(xb(i).lt.x0min(i).or.xb(i).gt.x0max(i))then
tb=1.0d+100
endif
enddo
if(tb.lt.tp)then
do i=1,nunknowns
xp(i)=xb(i)
fp(i)=fb(i)
enddo
if(tb.lt.TOLF)then
iwhichsolver=6
return
endif
endif
fsqsum=0.0d0
do i=1,nunknowns
fsqsum=fsqsum+fp(i)*fp(i)
enddo
tp=fsqsum
call nongradopt(nunknowns,fmin_funcnleq1,
& f1dim_funcnleq1,xp,x0min,x0max,TOLF,fsqsum)
! if(dabs(tp-fsqsum).gt.TOLF)then
! call RepeatCompassSearch(nunknowns,xp,fsqsum,
! & x0min,x0max,fmin_funcnleq1,f1dim_funcnleq1,
! & TOLF)
! endif
call funcnleq1(nunknowns,xp,fp,fsqsum)
tp=dabs(fp(1))
do i=2,nunknowns
if(dabs(fp(i)).gt.tp)tp=dabs(fp(i))
enddo
if(tp.lt.TOLF)then
iwhichsolver=7
return
endif
IERR=0
do i=1,nunknowns
if(dabs(xp(i)-x0(i)).gt.TOLF)IERR=1
enddo
if(IERR.eq.0)return
do i=1,nunknowns
x0(i)=xp(i)
enddo
enddo
end subroutine nonsyssolver
+18
View File
@@ -0,0 +1,18 @@
!------------------ Common Blocks -------------------------
integer numeval,maxndim,maxeval
parameter(maxndim=1000,maxeval=15000)
double precision xevaluated,fevaluated,flargest
common /FuncvRegresInteg/numeval
common /FuncvRegresDble/xevaluated(1:maxeval,1:maxndim),
& fevaluated(1:maxeval,1:maxndim),
& flargest(1:maxeval)
save /FuncvRegresInteg/,/FuncvRegresDble/
! numeval: the number of times that the system is evaluated so far
! iflargest: the index of the largest function for the latest evaluation
! xevaluated: the positions where the system is evaluated
! fevaluated: the function values at xevaluated
! flargest: the largest absolute function value
! maxndim: the maximum allowable dimensions of the system
! maxeval: the maximum allowable number of function evaluations
!--------------------------------------------------------------
@@ -0,0 +1,85 @@
program test
implicit none
integer nunknowns,iwhichsolver,i,j
double precision x0min(11),x0ori(11),xp(11),
& x0max(11),fequ(11),f1dim_funcsys
external funcsys,fsqsum_funcsys,f1dim_funcsys
nunknowns=11
do i=1,nunknowns
x0min(i)=-0.00001d0
x0ori(i)=1.0d0
x0max(i)=100.0d0
enddo
nunknowns=2
x0ori(1)=0.0d0
x0ori(2)=3.0d0
call nonsyssolver(funcsys,fsqsum_funcsys,
& f1dim_funcsys,x0min,x0ori,xp,x0max,
& fequ,nunknowns,iwhichsolver)
do i=1,nunknowns
write(*,*)fequ(i),xp(i),iwhichsolver
enddo
end
subroutine funcsys(nunknowns,x,f,fsqsum)
implicit none
integer nunknowns,i
double precision x(nunknowns),f(nunknowns),
& fsqsum
double precision R,p,K5,K6,K7,K8,K9,K10
parameter(R=10.0d0,p=40.0d0,
& K5=1.0d0,K6=1.0d0,
& K7=1.0d0,K8=0.1d0,
& K9=1.0d0,K10=0.1d0)
f(1)=x(1)-(1.0d0+0.5d0*dsin(x(1)))
f(2)=x(2)-(3.0d0+2.0d0*dsin(x(2)))
! Combustion of propane problem
! f(1)=x(1)-(3.0d0-x(4))
! f(2)=x(2)-(R-2.0d0*x(1)-x(4)-x(7)-
! & x(8)-x(9)-2.0d0*x(10))
! f(3)=x(3)-(2.0d0*R-0.5d0*x(9))
! f(4)=x(4)-x(1)*x(5)/(K5*x(2))
! f(5)=x(5)-(4-x(2)-0.5d0*x(6)-0.5d0*x(7))
! f(6)=x(6)-K6*dsqrt(x(2)*x(4)*x(11)/(p*x(1)))
! f(7)=x(7)-K7*dsqrt(x(1)*x(2)*x(11)/(p*x(4)))
! f(8)=x(8)-K8*x(1)*x(11)/(p*x(4))
! f(9)=x(9)-K9*(x(1)/x(4))*dsqrt(x(3)*x(11)/p)
! f(10)=x(10)-K10*x(1)*x(1)*x(11)/(p*x(4)*x(4))
! f(11)=x(11)-(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)
! & +x(7)+x(8)+x(9)+x(10))
fsqsum=0.0d0
do i=1,nunknowns
fsqsum=fsqsum+f(i)*f(i)
enddo
fsqsum=0.5d0*fsqsum
return
end
subroutine fsqsum_funcsys(nunknowns,xp,fsqsum)
implicit none
integer nunknowns
double precision xp(nunknowns),fsqsum,
& fequ(nunknowns)
call funcsys(nunknowns,xp,fequ,fsqsum)
return
end
double precision function f1dim_funcsys(x)
INTEGER NMAX
double precision x
PARAMETER (NMAX=1000)
CU USES funcsys
INTEGER j,ncom
double precision pcom(NMAX),xicom(NMAX),
& xt(NMAX),fequ(NMAX)
COMMON /f1com/ pcom,xicom,ncom
save /f1com/
do 11 j=1,ncom
xt(j)=pcom(j)+x*xicom(j)
11 continue
call funcsys(ncom,xt,fequ,f1dim_funcsys)
return
END
Binary file not shown.
@@ -0,0 +1,13 @@
SUBROUTINE addint(uf,uc,res,nf)
INTEGER nf
DOUBLE PRECISION res(nf,nf),uc(nf/2+1,nf/2+1),uf(nf,nf)
CU USES interp
INTEGER i,j
call interp(res,uc,nf)
do 12 j=1,nf
do 11 i=1,nf
uf(i,j)=uf(i,j)+res(i,j)
11 continue
12 continue
return
END
@@ -0,0 +1,32 @@
SUBROUTINE airy(x,ai,bi,aip,bip)
REAL ai,aip,bi,bip,x
CU USES bessik,bessjy
REAL absx,ri,rip,rj,rjp,rk,rkp,rootx,ry,ryp,z,PI,THIRD,TWOTHR,
*ONOVRT
PARAMETER (PI=3.1415927,THIRD=1./3.,TWOTHR=2.*THIRD,
*ONOVRT=.57735027)
absx=abs(x)
rootx=sqrt(absx)
z=TWOTHR*absx*rootx
if(x.gt.0.)then
call bessik(z,THIRD,ri,rk,rip,rkp)
ai=rootx*ONOVRT*rk/PI
bi=rootx*(rk/PI+2.*ONOVRT*ri)
call bessik(z,TWOTHR,ri,rk,rip,rkp)
aip=-x*ONOVRT*rk/PI
bip=x*(rk/PI+2.*ONOVRT*ri)
else if(x.lt.0.)then
call bessjy(z,THIRD,rj,ry,rjp,ryp)
ai=.5*rootx*(rj-ONOVRT*ry)
bi=-.5*rootx*(ry+ONOVRT*rj)
call bessjy(z,TWOTHR,rj,ry,rjp,ryp)
aip=.5*absx*(ONOVRT*ry+rj)
bip=.5*absx*(ONOVRT*rj-ry)
else
ai=.35502805
bi=ai/ONOVRT
aip=-.25881940
bip=-aip/ONOVRT
endif
return
END
@@ -0,0 +1,85 @@
SUBROUTINE amebsa(p,y,mp,np,ndim,pb,yb,ftol,funk,iter,temptr)
INTEGER iter,mp,ndim,np,NMAX
REAL ftol,temptr,yb,p(mp,np),pb(np),y(mp),funk
PARAMETER (NMAX=200)
EXTERNAL funk
CU USES amotsa,funk,ran1
INTEGER i,idum,ihi,ilo,inhi,j,m,n
REAL rtol,sum,swap,tt,yhi,ylo,ynhi,ysave,yt,ytry,psum(NMAX),
*amotsa,ran1
COMMON /ambsa/ tt,idum
tt=-temptr
1 do 12 n=1,ndim
sum=0.
do 11 m=1,ndim+1
sum=sum+p(m,n)
11 continue
psum(n)=sum
12 continue
2 ilo=1
inhi=1
ihi=2
ylo=y(1)+tt*log(ran1(idum))
ynhi=ylo
yhi=y(2)+tt*log(ran1(idum))
if (ylo.gt.yhi) then
ihi=1
inhi=2
ilo=2
ynhi=yhi
yhi=ylo
ylo=ynhi
endif
do 13 i=3,ndim+1
yt=y(i)+tt*log(ran1(idum))
if(yt.le.ylo) then
ilo=i
ylo=yt
endif
if(yt.gt.yhi) then
inhi=ihi
ynhi=yhi
ihi=i
yhi=yt
else if(yt.gt.ynhi) then
inhi=i
ynhi=yt
endif
13 continue
rtol=2.*abs(yhi-ylo)/(abs(yhi)+abs(ylo))
if (rtol.lt.ftol.or.iter.lt.0) then
swap=y(1)
y(1)=y(ilo)
y(ilo)=swap
do 14 n=1,ndim
swap=p(1,n)
p(1,n)=p(ilo,n)
p(ilo,n)=swap
14 continue
return
endif
iter=iter-2
ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,-1.0)
if (ytry.le.ylo) then
ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,2.0)
else if (ytry.ge.ynhi) then
ysave=yhi
ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,0.5)
if (ytry.ge.ysave) then
do 16 i=1,ndim+1
if(i.ne.ilo)then
do 15 j=1,ndim
psum(j)=0.5*(p(i,j)+p(ilo,j))
p(i,j)=psum(j)
15 continue
y(i)=funk(psum)
endif
16 continue
iter=iter-ndim
goto 1
endif
else
iter=iter+1
endif
goto 2
END
@@ -0,0 +1,71 @@
SUBROUTINE amoeba(p,y,mp,np,ndim,ftol,funk,iter)
INTEGER iter,mp,ndim,np,NMAX,ITMAX
REAL ftol,p(mp,np),y(mp),funk
PARAMETER (NMAX=20,ITMAX=5000)
EXTERNAL funk
CU USES amotry,funk
INTEGER i,ihi,ilo,inhi,j,m,n
REAL rtol,sum,swap,ysave,ytry,psum(NMAX),amotry
iter=0
1 do 12 n=1,ndim
sum=0.
do 11 m=1,ndim+1
sum=sum+p(m,n)
11 continue
psum(n)=sum
12 continue
2 ilo=1
if (y(1).gt.y(2)) then
ihi=1
inhi=2
else
ihi=2
inhi=1
endif
do 13 i=1,ndim+1
if(y(i).le.y(ilo)) ilo=i
if(y(i).gt.y(ihi)) then
inhi=ihi
ihi=i
else if(y(i).gt.y(inhi)) then
if(i.ne.ihi) inhi=i
endif
13 continue
rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo)))
if (rtol.lt.ftol) then
swap=y(1)
y(1)=y(ilo)
y(ilo)=swap
do 14 n=1,ndim
swap=p(1,n)
p(1,n)=p(ilo,n)
p(ilo,n)=swap
14 continue
return
endif
if (iter.ge.ITMAX) pause 'ITMAX exceeded in amoeba'
iter=iter+2
ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,-1.0)
if (ytry.le.y(ilo)) then
ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,2.0)
else if (ytry.ge.y(inhi)) then
ysave=y(ihi)
ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,0.5)
if (ytry.ge.ysave) then
do 16 i=1,ndim+1
if(i.ne.ilo)then
do 15 j=1,ndim
psum(j)=0.5*(p(i,j)+p(ilo,j))
p(i,j)=psum(j)
15 continue
y(i)=funk(psum)
endif
16 continue
iter=iter+ndim
goto 1
endif
else
iter=iter-1
endif
goto 2
END
@@ -0,0 +1,24 @@
FUNCTION amotry(p,y,psum,mp,np,ndim,funk,ihi,fac)
INTEGER ihi,mp,ndim,np,NMAX
REAL amotry,fac,p(mp,np),psum(np),y(mp),funk
PARAMETER (NMAX=20)
EXTERNAL funk
CU USES funk
INTEGER j
REAL fac1,fac2,ytry,ptry(NMAX)
fac1=(1.-fac)/ndim
fac2=fac1-fac
do 11 j=1,ndim
ptry(j)=psum(j)*fac1-p(ihi,j)*fac2
11 continue
ytry=funk(ptry)
if (ytry.lt.y(ihi)) then
y(ihi)=ytry
do 12 j=1,ndim
psum(j)=psum(j)-p(ihi,j)+ptry(j)
p(ihi,j)=ptry(j)
12 continue
endif
amotry=ytry
return
END
@@ -0,0 +1,33 @@
FUNCTION amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,fac)
INTEGER ihi,mp,ndim,np,NMAX
REAL amotsa,fac,yb,yhi,p(mp,np),pb(np),psum(np),y(mp),funk
PARAMETER (NMAX=200)
EXTERNAL funk
CU USES funk,ran1
INTEGER idum,j
REAL fac1,fac2,tt,yflu,ytry,ptry(NMAX),ran1
COMMON /ambsa/ tt,idum
fac1=(1.-fac)/ndim
fac2=fac1-fac
do 11 j=1,ndim
ptry(j)=psum(j)*fac1-p(ihi,j)*fac2
11 continue
ytry=funk(ptry)
if (ytry.le.yb) then
do 12 j=1,ndim
pb(j)=ptry(j)
12 continue
yb=ytry
endif
yflu=ytry-tt*log(ran1(idum))
if (yflu.lt.yhi) then
y(ihi)=ytry
yhi=yflu
do 13 j=1,ndim
psum(j)=psum(j)-p(ihi,j)+ptry(j)
p(ihi,j)=ptry(j)
13 continue
endif
amotsa=yflu
return
END
@@ -0,0 +1,62 @@
SUBROUTINE anneal(x,y,iorder,ncity)
INTEGER ncity,iorder(ncity)
REAL x(ncity),y(ncity)
CU USES irbit1,metrop,ran3,revcst,revers,trncst,trnspt
INTEGER i,i1,i2,idec,idum,iseed,j,k,nlimit,nn,nover,nsucc,n(6),
*irbit1
REAL de,path,t,tfactr,ran3,alen,x1,x2,y1,y2
LOGICAL ans
alen(x1,x2,y1,y2)=sqrt((x2-x1)**2+(y2-y1)**2)
nover=100*ncity
nlimit=10*ncity
tfactr=0.9
path=0.0
t=0.5
do 11 i=1,ncity-1
i1=iorder(i)
i2=iorder(i+1)
path=path+alen(x(i1),x(i2),y(i1),y(i2))
11 continue
i1=iorder(ncity)
i2=iorder(1)
path=path+alen(x(i1),x(i2),y(i1),y(i2))
idum=-1
iseed=111
do 13 j=1,100
nsucc=0
do 12 k=1,nover
1 n(1)=1+int(ncity*ran3(idum))
n(2)=1+int((ncity-1)*ran3(idum))
if (n(2).ge.n(1)) n(2)=n(2)+1
nn=1+mod((n(1)-n(2)+ncity-1),ncity)
if (nn.lt.3) goto 1
idec=irbit1(iseed)
if (idec.eq.0) then
n(3)=n(2)+int(abs(nn-2)*ran3(idum))+1
n(3)=1+mod(n(3)-1,ncity)
call trncst(x,y,iorder,ncity,n,de)
call metrop(de,t,ans)
if (ans) then
nsucc=nsucc+1
path=path+de
call trnspt(iorder,ncity,n)
endif
else
call revcst(x,y,iorder,ncity,n,de)
call metrop(de,t,ans)
if (ans) then
nsucc=nsucc+1
path=path+de
call revers(iorder,ncity,n)
endif
endif
if (nsucc.ge.nlimit) goto 2
12 continue
2 write(*,*)
write(*,*) 'T =',t,' Path Length =',path
write(*,*) 'Successful Moves: ',nsucc
t=t*tfactr
if (nsucc.eq.0) return
13 continue
return
END
@@ -0,0 +1,14 @@
DOUBLE PRECISION FUNCTION anorm2(a,n)
INTEGER n
DOUBLE PRECISION a(n,n)
INTEGER i,j
DOUBLE PRECISION sum
sum=0.d0
do 12 j=1,n
do 11 i=1,n
sum=sum+a(i,j)**2
11 continue
12 continue
anorm2=sqrt(sum)/n
return
END
@@ -0,0 +1,20 @@
SUBROUTINE arcmak(nfreq,nchh,nradd)
INTEGER nchh,nradd,nfreq(nchh),MC,NWK,MAXINT
PARAMETER (MC=512,NWK=20,MAXINT=2147483647)
INTEGER j,jdif,minint,nc,nch,nrad,ncum,ncumfq(MC+2),ilob(NWK),
*iupb(NWK)
COMMON /arccom/ ncumfq,iupb,ilob,nch,nrad,minint,jdif,nc,ncum
SAVE /arccom/
if(nchh.gt.MC)pause 'MC too small in arcmak'
if(nradd.gt.256)pause 'nradd may not exceed 256 in arcmak'
minint=MAXINT/nradd
nch=nchh
nrad=nradd
ncumfq(1)=0
do 11 j=2,nch+1
ncumfq(j)=ncumfq(j-1)+max(nfreq(j-1),1)
11 continue
ncumfq(nch+2)=ncumfq(nch+1)+1
ncum=ncumfq(nch+2)
return
END
@@ -0,0 +1,75 @@
SUBROUTINE arcode(ich,code,lcode,lcd,isign)
INTEGER ich,isign,lcd,lcode,MC,NWK
CHARACTER*1 code(lcode)
PARAMETER (MC=512,NWK=20)
CU USES arcsum
INTEGER ihi,j,ja,jdif,jh,jl,k,m,minint,nc,nch,nrad,ilob(NWK),
*iupb(NWK),ncumfq(MC+2),ncum,JTRY
COMMON /arccom/ ncumfq,iupb,ilob,nch,nrad,minint,jdif,nc,ncum
SAVE /arccom/
JTRY(j,k,m)=int((dble(k)*dble(j))/dble(m))
if (isign.eq.0) then
jdif=nrad-1
do 11 j=NWK,1,-1
iupb(j)=nrad-1
ilob(j)=0
nc=j
if(jdif.gt.minint)return
jdif=(jdif+1)*nrad-1
11 continue
pause 'NWK too small in arcode'
else
if (isign.gt.0) then
if(ich.gt.nch.or.ich.lt.0)pause 'bad ich in arcode'
else
ja=ichar(code(lcd))-ilob(nc)
do 12 j=nc+1,NWK
ja=ja*nrad+(ichar(code(j+lcd-nc))-ilob(j))
12 continue
ich=0
ihi=nch+1
1 if(ihi-ich.gt.1) then
m=(ich+ihi)/2
if (ja.ge.JTRY(jdif,ncumfq(m+1),ncum)) then
ich=m
else
ihi=m
endif
goto 1
endif
if(ich.eq.nch)return
endif
jh=JTRY(jdif,ncumfq(ich+2),ncum)
jl=JTRY(jdif,ncumfq(ich+1),ncum)
jdif=jh-jl
call arcsum(ilob,iupb,jh,NWK,nrad,nc)
call arcsum(ilob,ilob,jl,NWK,nrad,nc)
do 13 j=nc,NWK
if(ich.ne.nch.and.iupb(j).ne.ilob(j))goto 2
if(lcd.gt.lcode)pause 'lcode too small in arcode'
if(isign.gt.0) code(lcd)=char(ilob(j))
lcd=lcd+1
13 continue
return
2 nc=j
j=0
3 if (jdif.lt.minint) then
j=j+1
jdif=jdif*nrad
goto 3
endif
if (nc-j.lt.1) pause 'NWK too small in arcode'
if(j.ne.0)then
do 14 k=nc,NWK
iupb(k-j)=iupb(k)
ilob(k-j)=ilob(k)
14 continue
endif
nc=nc-j
do 15 k=NWK-j+1,NWK
iupb(k)=0
ilob(k)=0
15 continue
endif
return
END
@@ -0,0 +1,18 @@
SUBROUTINE arcsum(iin,iout,ja,nwk,nrad,nc)
INTEGER ja,nc,nrad,nwk,iin(*),iout(*)
INTEGER j,jtmp,karry
karry=0
do 11 j=nwk,nc+1,-1
jtmp=ja
ja=ja/nrad
iout(j)=iin(j)+(jtmp-ja*nrad)+karry
if (iout(j).ge.nrad) then
iout(j)=iout(j)-nrad
karry=1
else
karry=0
endif
11 continue
iout(nc)=iin(nc)+ja+karry
return
END
@@ -0,0 +1,10 @@
SUBROUTINE asolve(n,b,x,itrnsp)
INTEGER n,itrnsp,ija,NMAX,i
DOUBLE PRECISION x(n),b(n),sa
PARAMETER (NMAX=1000)
COMMON /mat/ sa(NMAX),ija(NMAX)
do 11 i=1,n
x(i)=b(i)/sa(i)
11 continue
return
END
@@ -0,0 +1,13 @@
SUBROUTINE atimes(n,x,r,itrnsp)
INTEGER n,itrnsp,ija,NMAX
DOUBLE PRECISION x(n),r(n),sa
PARAMETER (NMAX=1000)
COMMON /mat/ sa(NMAX),ija(NMAX)
CU USES dsprsax,dsprstx
if (itrnsp.eq.0) then
call dsprsax(sa,ija,x,r,n)
else
call dsprstx(sa,ija,x,r,n)
endif
return
END
@@ -0,0 +1,20 @@
SUBROUTINE avevar(data,n,ave,var)
INTEGER n
REAL ave,var,data(n)
INTEGER j
REAL s,ep
ave=0.0
do 11 j=1,n
ave=ave+data(j)
11 continue
ave=ave/n
var=0.0
ep=0.0
do 12 j=1,n
s=data(j)-ave
ep=ep+s
var=var+s*s
12 continue
var=(var-ep**2/n)/(n-1)
return
END
@@ -0,0 +1,44 @@
PROGRAM badluk
INTEGER ic,icon,idwk,ifrac,im,iybeg,iyend,iyyy,jd,jday,n,julday
REAL TIMZON,frac
PARAMETER (TIMZON=-5./24.)
DATA iybeg,iyend /1900,2000/
CU USES flmoon,julday
write (*,'(1x,a,i5,a,i5)') 'Full moons on Friday the 13th from',
*iybeg,' to',iyend
do 12 iyyy=iybeg,iyend
do 11 im=1,12
jday=julday(im,13,iyyy)
idwk=mod(jday+1,7)
if(idwk.eq.5) then
n=12.37*(iyyy-1900+(im-0.5)/12.)
icon=0
1 call flmoon(n,2,jd,frac)
ifrac=nint(24.*(frac+TIMZON))
if(ifrac.lt.0)then
jd=jd-1
ifrac=ifrac+24
endif
if(ifrac.gt.12)then
jd=jd+1
ifrac=ifrac-12
else
ifrac=ifrac+12
endif
if(jd.eq.jday)then
write (*,'(/1x,i2,a,i2,a,i4)') im,'/',13,'/',iyyy
write (*,'(1x,a,i2,a)') 'Full moon ',ifrac,
*' hrs after midnight (EST).'
goto 2
else
ic=isign(1,jday-jd)
if(ic.eq.-icon) goto 2
icon=ic
n=n+ic
endif
goto 1
2 continue
endif
11 continue
12 continue
END
@@ -0,0 +1,47 @@
SUBROUTINE balanc(a,n,np)
INTEGER n,np
REAL a(np,np),RADIX,SQRDX
PARAMETER (RADIX=2.,SQRDX=RADIX**2)
INTEGER i,j,last
REAL c,f,g,r,s
1 continue
last=1
do 14 i=1,n
c=0.
r=0.
do 11 j=1,n
if(j.ne.i)then
c=c+abs(a(j,i))
r=r+abs(a(i,j))
endif
11 continue
if(c.ne.0..and.r.ne.0.)then
g=r/RADIX
f=1.
s=c+r
2 if(c.lt.g)then
f=f*RADIX
c=c*SQRDX
goto 2
endif
g=r*RADIX
3 if(c.gt.g)then
f=f/RADIX
c=c/SQRDX
goto 3
endif
if((c+r)/f.lt.0.95*s)then
last=0
g=1./f
do 12 j=1,n
a(i,j)=a(i,j)*g
12 continue
do 13 j=1,n
a(j,i)=a(j,i)*f
13 continue
endif
endif
14 continue
if(last.eq.0)goto 1
return
END
@@ -0,0 +1,31 @@
SUBROUTINE banbks(a,n,m1,m2,np,mp,al,mpl,indx,b)
INTEGER m1,m2,mp,mpl,n,np,indx(n)
REAL a(np,mp),al(np,mpl),b(n)
INTEGER i,k,l,mm
REAL dum
mm=m1+m2+1
if(mm.gt.mp.or.m1.gt.mpl.or.n.gt.np) pause 'bad args in banbks'
l=m1
do 12 k=1,n
i=indx(k)
if(i.ne.k)then
dum=b(k)
b(k)=b(i)
b(i)=dum
endif
if(l.lt.n)l=l+1
do 11 i=k+1,l
b(i)=b(i)-al(k,i-k)*b(k)
11 continue
12 continue
l=1
do 14 i=n,1,-1
dum=b(i)
do 13 k=2,l
dum=dum-a(i,k)*b(k+i-1)
13 continue
b(i)=dum/a(i,1)
if(l.lt.mm) l=l+1
14 continue
return
END
@@ -0,0 +1,51 @@
SUBROUTINE bandec(a,n,m1,m2,np,mp,al,mpl,indx,d)
INTEGER m1,m2,mp,mpl,n,np,indx(n)
REAL d,a(np,mp),al(np,mpl),TINY
PARAMETER (TINY=1.e-20)
INTEGER i,j,k,l,mm
REAL dum
mm=m1+m2+1
if(mm.gt.mp.or.m1.gt.mpl.or.n.gt.np) pause 'bad args in bandec'
l=m1
do 13 i=1,m1
do 11 j=m1+2-i,mm
a(i,j-l)=a(i,j)
11 continue
l=l-1
do 12 j=mm-l,mm
a(i,j)=0.
12 continue
13 continue
d=1.
l=m1
do 18 k=1,n
dum=a(k,1)
i=k
if(l.lt.n)l=l+1
do 14 j=k+1,l
if(abs(a(j,1)).gt.abs(dum))then
dum=a(j,1)
i=j
endif
14 continue
indx(k)=i
if(dum.eq.0.) a(k,1)=TINY
if(i.ne.k)then
d=-d
do 15 j=1,mm
dum=a(k,j)
a(k,j)=a(i,j)
a(i,j)=dum
15 continue
endif
do 17 i=k+1,l
dum=a(i,1)/a(k,1)
al(k,i-k)=dum
do 16 j=2,mm
a(i,j-1)=a(i,j)-dum*a(k,j)
16 continue
a(i,mm)=0.
17 continue
18 continue
return
END
@@ -0,0 +1,13 @@
SUBROUTINE banmul(a,n,m1,m2,np,mp,x,b)
INTEGER m1,m2,mp,n,np
REAL a(np,mp),b(n),x(n)
INTEGER i,j,k
do 12 i=1,n
b(i)=0.
k=i-m1-1
do 11 j=max(1,1-k),min(m1+m2+1,n-k)
b(i)=b(i)+a(i,j)*x(j+k)
11 continue
12 continue
return
END
@@ -0,0 +1,35 @@
SUBROUTINE bcucof(y,y1,y2,y12,d1,d2,c)
REAL d1,d2,c(4,4),y(4),y1(4),y12(4),y2(4)
INTEGER i,j,k,l
REAL d1d2,xx,cl(16),wt(16,16),x(16)
SAVE wt
DATA wt/1,0,-3,2,4*0,-3,0,9,-6,2,0,-6,4,8*0,3,0,-9,6,-2,0,6,-4,10*
*0,9,-6,2*0,-6,4,2*0,3,-2,6*0,-9,6,2*0,6,-4,4*0,1,0,-3,2,-2,0,6,-4,
*1,0,-3,2,8*0,-1,0,3,-2,1,0,-3,2,10*0,-3,2,2*0,3,-2,6*0,3,-2,2*0,
*-6,4,2*0,3,-2,0,1,-2,1,5*0,-3,6,-3,0,2,-4,2,9*0,3,-6,3,0,-2,4,-2,
*10*0,-3,3,2*0,2,-2,2*0,-1,1,6*0,3,-3,2*0,-2,2,5*0,1,-2,1,0,-2,4,
*-2,0,1,-2,1,9*0,-1,2,-1,0,1,-2,1,10*0,1,-1,2*0,-1,1,6*0,-1,1,2*0,
*2,-2,2*0,-1,1/
d1d2=d1*d2
do 11 i=1,4
x(i)=y(i)
x(i+4)=y1(i)*d1
x(i+8)=y2(i)*d2
x(i+12)=y12(i)*d1d2
11 continue
do 13 i=1,16
xx=0.
do 12 k=1,16
xx=xx+wt(i,k)*x(k)
12 continue
cl(i)=xx
13 continue
l=0
do 15 i=1,4
do 14 j=1,4
l=l+1
c(i,j)=cl(l)
14 continue
15 continue
return
END
@@ -0,0 +1,23 @@
SUBROUTINE bcuint(y,y1,y2,y12,x1l,x1u,x2l,x2u,x1,x2,ansy,ansy1,
*ansy2)
REAL ansy,ansy1,ansy2,x1,x1l,x1u,x2,x2l,x2u,y(4),y1(4),y12(4),
*y2(4)
CU USES bcucof
INTEGER i
REAL t,u,c(4,4)
call bcucof(y,y1,y2,y12,x1u-x1l,x2u-x2l,c)
if(x1u.eq.x1l.or.x2u.eq.x2l)pause 'bad input in bcuint'
t=(x1-x1l)/(x1u-x1l)
u=(x2-x2l)/(x2u-x2l)
ansy=0.
ansy2=0.
ansy1=0.
do 11 i=4,1,-1
ansy=t*ansy+((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
ansy2=t*ansy2+(3.*c(i,4)*u+2.*c(i,3))*u+c(i,2)
ansy1=u*ansy1+(3.*c(4,i)*t+2.*c(3,i))*t+c(2,i)
11 continue
ansy1=ansy1/(x1u-x1l)
ansy2=ansy2/(x2u-x2l)
return
END
@@ -0,0 +1,19 @@
SUBROUTINE beschb(x,gam1,gam2,gampl,gammi)
INTEGER NUSE1,NUSE2
DOUBLE PRECISION gam1,gam2,gammi,gampl,x
PARAMETER (NUSE1=5,NUSE2=5)
CU USES chebev
REAL xx,c1(7),c2(8),chebev
SAVE c1,c2
DATA c1/-1.142022680371168d0,6.5165112670737d-3,3.087090173086d-4,
*-3.4706269649d-6,6.9437664d-9,3.67795d-11,-1.356d-13/
DATA c2/1.843740587300905d0,-7.68528408447867d-2,
*1.2719271366546d-3,-4.9717367042d-6,-3.31261198d-8,2.423096d-10,
*-1.702d-13,-1.49d-15/
xx=8.d0*x*x-1.d0
gam1=chebev(-1.,1.,c1,NUSE1,xx)
gam2=chebev(-1.,1.,c2,NUSE2,xx)
gampl=gam2-x*gam1
gammi=gam2+x*gam1
return
END
@@ -0,0 +1,32 @@
FUNCTION bessi(n,x)
INTEGER n,IACC
REAL bessi,x,BIGNO,BIGNI
PARAMETER (IACC=40,BIGNO=1.0e10,BIGNI=1.0e-10)
CU USES bessi0
INTEGER j,m
REAL bi,bim,bip,tox,bessi0
if (n.lt.2) pause 'bad argument n in bessi'
if (x.eq.0.) then
bessi=0.
else
tox=2.0/abs(x)
bip=0.0
bi=1.0
bessi=0.
m=2*((n+int(sqrt(float(IACC*n)))))
do 11 j=m,1,-1
bim=bip+float(j)*tox*bi
bip=bi
bi=bim
if (abs(bi).gt.BIGNO) then
bessi=bessi*BIGNI
bi=bi*BIGNI
bip=bip*BIGNI
endif
if (j.eq.n) bessi=bip
11 continue
bessi=bessi*bessi0(x)/bi
if (x.lt.0..and.mod(n,2).eq.1) bessi=-bessi
endif
return
END
@@ -0,0 +1,21 @@
FUNCTION bessi0(x)
REAL bessi0,x
REAL ax
DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0,
*1.2067492d0,0.2659732d0,0.360768d-1,0.45813d-2/
DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1,
*0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1,0.2635537d-1,
*-0.1647633d-1,0.392377d-2/
if (abs(x).lt.3.75) then
y=(x/3.75)**2
bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))
else
ax=abs(x)
y=3.75/ax
bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
*(q7+y*(q8+y*q9))))))))
endif
return
END
@@ -0,0 +1,22 @@
FUNCTION bessi1(x)
REAL bessi1,x
REAL ax
DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,
*0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/
DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,-0.3988024d-1,
*-0.362018d-2,0.163801d-2,-0.1031555d-1,0.2282967d-1,-0.2895312d-1,
*0.1787654d-1,-0.420059d-2/
if (abs(x).lt.3.75) then
y=(x/3.75)**2
bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
else
ax=abs(x)
y=3.75/ax
bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
*(q7+y*(q8+y*q9))))))))
if(x.lt.0.)bessi1=-bessi1
endif
return
END
@@ -0,0 +1,129 @@
SUBROUTINE bessik(x,xnu,ri,rk,rip,rkp)
INTEGER MAXIT
REAL ri,rip,rk,rkp,x,xnu,XMIN
DOUBLE PRECISION EPS,FPMIN,PI
PARAMETER (EPS=1.e-10,FPMIN=1.e-30,MAXIT=10000,XMIN=2.,
*PI=3.141592653589793d0)
CU USES beschb
INTEGER i,l,nl
DOUBLE PRECISION a,a1,b,c,d,del,del1,delh,dels,e,f,fact,fact2,ff,
*gam1,gam2,gammi,gampl,h,p,pimu,q,q1,q2,qnew,ril,ril1,rimu,rip1,
*ripl,ritemp,rk1,rkmu,rkmup,rktemp,s,sum,sum1,x2,xi,xi2,xmu,xmu2
if(x.le.0..or.xnu.lt.0.) pause 'bad arguments in bessik'
nl=int(xnu+.5d0)
xmu=xnu-nl
xmu2=xmu*xmu
xi=1.d0/x
xi2=2.d0*xi
h=xnu*xi
if(h.lt.FPMIN)h=FPMIN
b=xi2*xnu
d=0.d0
c=h
do 11 i=1,MAXIT
b=b+xi2
d=1.d0/(b+d)
c=b+1.d0/c
del=c*d
h=del*h
if(abs(del-1.d0).lt.EPS)goto 1
11 continue
pause 'x too large in bessik; try asymptotic expansion'
1 continue
ril=FPMIN
ripl=h*ril
ril1=ril
rip1=ripl
fact=xnu*xi
do 12 l=nl,1,-1
ritemp=fact*ril+ripl
fact=fact-xi
ripl=fact*ritemp+ril
ril=ritemp
12 continue
f=ripl/ril
if(x.lt.XMIN) then
x2=.5d0*x
pimu=PI*xmu
if(abs(pimu).lt.EPS)then
fact=1.d0
else
fact=pimu/sin(pimu)
endif
d=-log(x2)
e=xmu*d
if(abs(e).lt.EPS)then
fact2=1.d0
else
fact2=sinh(e)/e
endif
call beschb(xmu,gam1,gam2,gampl,gammi)
ff=fact*(gam1*cosh(e)+gam2*fact2*d)
sum=ff
e=exp(e)
p=0.5d0*e/gampl
q=0.5d0/(e*gammi)
c=1.d0
d=x2*x2
sum1=p
do 13 i=1,MAXIT
ff=(i*ff+p+q)/(i*i-xmu2)
c=c*d/i
p=p/(i-xmu)
q=q/(i+xmu)
del=c*ff
sum=sum+del
del1=c*(p-i*ff)
sum1=sum1+del1
if(abs(del).lt.abs(sum)*EPS)goto 2
13 continue
pause 'bessk series failed to converge'
2 continue
rkmu=sum
rk1=sum1*xi2
else
b=2.d0*(1.d0+x)
d=1.d0/b
delh=d
h=delh
q1=0.d0
q2=1.d0
a1=.25d0-xmu2
c=a1
q=c
a=-a1
s=1.d0+q*delh
do 14 i=2,MAXIT
a=a-2*(i-1)
c=-a*c/i
qnew=(q1-b*q2)/a
q1=q2
q2=qnew
q=q+c*qnew
b=b+2.d0
d=1.d0/(b+a*d)
delh=(b*d-1.d0)*delh
h=h+delh
dels=q*delh
s=s+dels
if(abs(dels/s).lt.EPS)goto 3
14 continue
pause 'bessik: failure to converge in cf2'
3 continue
h=a1*h
rkmu=sqrt(PI/(2.d0*x))*exp(-x)/s
rk1=rkmu*(xmu+x+.5d0-h)*xi
endif
rkmup=xmu*xi*rkmu-rk1
rimu=xi/(f*rkmu-rkmup)
ri=(rimu*ril1)/ril
rip=(rimu*rip1)/ril
do 15 i=1,nl
rktemp=(xmu+i)*xi2*rk1+rkmu
rkmu=rk1
rk1=rktemp
15 continue
rk=rkmu
rkp=xnu*xi*rkmu-rk1
return
END
@@ -0,0 +1,49 @@
FUNCTION bessj(n,x)
INTEGER n,IACC
REAL bessj,x,BIGNO,BIGNI
PARAMETER (IACC=40,BIGNO=1.e10,BIGNI=1.e-10)
CU USES bessj0,bessj1
INTEGER j,jsum,m
REAL ax,bj,bjm,bjp,sum,tox,bessj0,bessj1
if(n.lt.2)pause 'bad argument n in bessj'
ax=abs(x)
if(ax.eq.0.)then
bessj=0.
else if(ax.gt.float(n))then
tox=2./ax
bjm=bessj0(ax)
bj=bessj1(ax)
do 11 j=1,n-1
bjp=j*tox*bj-bjm
bjm=bj
bj=bjp
11 continue
bessj=bj
else
tox=2./ax
m=2*((n+int(sqrt(float(IACC*n))))/2)
bessj=0.
jsum=0
sum=0.
bjp=0.
bj=1.
do 12 j=m,1,-1
bjm=j*tox*bj-bjp
bjp=bj
bj=bjm
if(abs(bj).gt.BIGNO)then
bj=bj*BIGNI
bjp=bjp*BIGNI
bessj=bessj*BIGNI
sum=sum*BIGNI
endif
if(jsum.ne.0)sum=sum+bj
jsum=1-jsum
if(j.eq.n)bessj=bjp
12 continue
sum=2.*sum-bj
bessj=bessj/sum
endif
if(x.lt.0..and.mod(n,2).eq.1)bessj=-bessj
return
END
@@ -0,0 +1,28 @@
FUNCTION bessj0(x)
REAL bessj0,x
REAL ax,xx,z
DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
*s1,s2,s3,s4,s5,s6,y
SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
*s5,s6
DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
*-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1,
*.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
DATA r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,
*651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/,s1,s2,
*s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0,
*59272.64853d0,267.8532712d0,1.d0/
if(abs(x).lt.8.)then
y=x**2
bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y*
*(s4+y*(s5+y*s6)))))
else
ax=abs(x)
z=8./ax
y=z**2
xx=ax-.785398164
bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
*p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
endif
return
END
@@ -0,0 +1,28 @@
FUNCTION bessj1(x)
REAL bessj1,x
REAL ax,xx,z
DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
*s1,s2,s3,s4,s5,s6,y
SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
*s5,s6
DATA r1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0,
*242396853.1d0,-2972611.439d0,15704.48260d0,-30.16036606d0/,s1,s2,
*s3,s4,s5,s6/144725228442.d0,2300535178.d0,18583304.74d0,
*99447.43394d0,376.9991397d0,1.d0/
DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4,
*.2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/.04687499995d0,
*-.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/
if(abs(x).lt.8.)then
y=x**2
bessj1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+
*y*(s4+y*(s5+y*s6)))))
else
ax=abs(x)
z=8./ax
y=z**2
xx=ax-2.356194491
bessj1=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
*p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))*sign(1.,x)
endif
return
END
@@ -0,0 +1,162 @@
SUBROUTINE bessjy(x,xnu,rj,ry,rjp,ryp)
INTEGER MAXIT
REAL rj,rjp,ry,ryp,x,xnu,XMIN
DOUBLE PRECISION EPS,FPMIN,PI
PARAMETER (EPS=1.e-10,FPMIN=1.e-30,MAXIT=10000,XMIN=2.,
*PI=3.141592653589793d0)
CU USES beschb
INTEGER i,isign,l,nl
DOUBLE PRECISION a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,
*f,fact,fact2,fact3,ff,gam,gam1,gam2,gammi,gampl,h,p,pimu,pimu2,q,
*r,rjl,rjl1,rjmu,rjp1,rjpl,rjtemp,ry1,rymu,rymup,rytemp,sum,sum1,
*temp,w,x2,xi,xi2,xmu,xmu2
if(x.le.0..or.xnu.lt.0.) pause 'bad arguments in bessjy'
if(x.lt.XMIN)then
nl=int(xnu+.5d0)
else
nl=max(0,int(xnu-x+1.5d0))
endif
xmu=xnu-nl
xmu2=xmu*xmu
xi=1.d0/x
xi2=2.d0*xi
w=xi2/PI
isign=1
h=xnu*xi
if(h.lt.FPMIN)h=FPMIN
b=xi2*xnu
d=0.d0
c=h
do 11 i=1,MAXIT
b=b+xi2
d=b-d
if(abs(d).lt.FPMIN)d=FPMIN
c=b-1.d0/c
if(abs(c).lt.FPMIN)c=FPMIN
d=1.d0/d
del=c*d
h=del*h
if(d.lt.0.d0)isign=-isign
if(abs(del-1.d0).lt.EPS)goto 1
11 continue
pause 'x too large in bessjy; try asymptotic expansion'
1 continue
rjl=isign*FPMIN
rjpl=h*rjl
rjl1=rjl
rjp1=rjpl
fact=xnu*xi
do 12 l=nl,1,-1
rjtemp=fact*rjl+rjpl
fact=fact-xi
rjpl=fact*rjtemp-rjl
rjl=rjtemp
12 continue
if(rjl.eq.0.d0)rjl=EPS
f=rjpl/rjl
if(x.lt.XMIN) then
x2=.5d0*x
pimu=PI*xmu
if(abs(pimu).lt.EPS)then
fact=1.d0
else
fact=pimu/sin(pimu)
endif
d=-log(x2)
e=xmu*d
if(abs(e).lt.EPS)then
fact2=1.d0
else
fact2=sinh(e)/e
endif
call beschb(xmu,gam1,gam2,gampl,gammi)
ff=2.d0/PI*fact*(gam1*cosh(e)+gam2*fact2*d)
e=exp(e)
p=e/(gampl*PI)
q=1.d0/(e*PI*gammi)
pimu2=0.5d0*pimu
if(abs(pimu2).lt.EPS)then
fact3=1.d0
else
fact3=sin(pimu2)/pimu2
endif
r=PI*pimu2*fact3*fact3
c=1.d0
d=-x2*x2
sum=ff+r*q
sum1=p
do 13 i=1,MAXIT
ff=(i*ff+p+q)/(i*i-xmu2)
c=c*d/i
p=p/(i-xmu)
q=q/(i+xmu)
del=c*(ff+r*q)
sum=sum+del
del1=c*p-i*del
sum1=sum1+del1
if(abs(del).lt.(1.d0+abs(sum))*EPS)goto 2
13 continue
pause 'bessy series failed to converge'
2 continue
rymu=-sum
ry1=-sum1*xi2
rymup=xmu*xi*rymu-ry1
rjmu=w/(rymup-f*rymu)
else
a=.25d0-xmu2
p=-.5d0*xi
q=1.d0
br=2.d0*x
bi=2.d0
fact=a*xi/(p*p+q*q)
cr=br+q*fact
ci=bi+p*fact
den=br*br+bi*bi
dr=br/den
di=-bi/den
dlr=cr*dr-ci*di
dli=cr*di+ci*dr
temp=p*dlr-q*dli
q=p*dli+q*dlr
p=temp
do 14 i=2,MAXIT
a=a+2*(i-1)
bi=bi+2.d0
dr=a*dr+br
di=a*di+bi
if(abs(dr)+abs(di).lt.FPMIN)dr=FPMIN
fact=a/(cr*cr+ci*ci)
cr=br+cr*fact
ci=bi-ci*fact
if(abs(cr)+abs(ci).lt.FPMIN)cr=FPMIN
den=dr*dr+di*di
dr=dr/den
di=-di/den
dlr=cr*dr-ci*di
dli=cr*di+ci*dr
temp=p*dlr-q*dli
q=p*dli+q*dlr
p=temp
if(abs(dlr-1.d0)+abs(dli).lt.EPS)goto 3
14 continue
pause 'cf2 failed in bessjy'
3 continue
gam=(p-f)/q
rjmu=sqrt(w/((p-f)*gam+q))
rjmu=sign(rjmu,rjl)
rymu=rjmu*gam
rymup=rymu*(p+q/gam)
ry1=xmu*xi*rymu-rymup
endif
fact=rjmu/rjl
rj=rjl1*fact
rjp=rjp1*fact
do 15 i=1,nl
rytemp=(xmu+i)*xi2*ry1-rymu
rymu=ry1
ry1=rytemp
15 continue
ry=rymu
ryp=xnu*xi*rymu-ry1
return
END
@@ -0,0 +1,18 @@
FUNCTION bessk(n,x)
INTEGER n
REAL bessk,x
CU USES bessk0,bessk1
INTEGER j
REAL bk,bkm,bkp,tox,bessk0,bessk1
if (n.lt.2) pause 'bad argument n in bessk'
tox=2.0/x
bkm=bessk0(x)
bk=bessk1(x)
do 11 j=1,n-1
bkp=bkm+j*tox*bk
bkm=bk
bk=bkp
11 continue
bessk=bk
return
END
@@ -0,0 +1,21 @@
FUNCTION bessk0(x)
REAL bessk0,x
CU USES bessi0
REAL bessi0
DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y
SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7
DATA p1,p2,p3,p4,p5,p6,p7/-0.57721566d0,0.42278420d0,0.23069756d0,
*0.3488590d-1,0.262698d-2,0.10750d-3,0.74d-5/
DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,-0.7832358d-1,0.2189568d-1,
*-0.1062446d-1,0.587872d-2,-0.251540d-2,0.53208d-3/
if (x.le.2.0) then
y=x*x/4.0
bessk0=(-log(x/2.0)*bessi0(x))+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*
*(p6+y*p7))))))
else
y=(2.0/x)
bessk0=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
*q7))))))
endif
return
END
@@ -0,0 +1,21 @@
FUNCTION bessk1(x)
REAL bessk1,x
CU USES bessi1
REAL bessi1
DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y
SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7
DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,0.15443144d0,-0.67278579d0,
*-0.18156897d0,-0.1919402d-1,-0.110404d-2,-0.4686d-4/
DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,0.23498619d0,-0.3655620d-1,
*0.1504268d-1,-0.780353d-2,0.325614d-2,-0.68245d-3/
if (x.le.2.0) then
y=x*x/4.0
bessk1=(log(x/2.0)*bessi1(x))+(1.0/x)*(p1+y*(p2+y*(p3+y*(p4+y*
*(p5+y*(p6+y*p7))))))
else
y=2.0/x
bessk1=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
*q7))))))
endif
return
END
@@ -0,0 +1,18 @@
FUNCTION bessy(n,x)
INTEGER n
REAL bessy,x
CU USES bessy0,bessy1
INTEGER j
REAL by,bym,byp,tox,bessy0,bessy1
if(n.lt.2)pause 'bad argument n in bessy'
tox=2./x
by=bessy1(x)
bym=bessy0(x)
do 11 j=1,n-1
byp=j*tox*by-bym
bym=by
by=byp
11 continue
bessy=by
return
END
@@ -0,0 +1,28 @@
FUNCTION bessy0(x)
REAL bessy0,x
CU USES bessj0
REAL xx,z,bessj0
DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
*s1,s2,s3,s4,s5,s6,y
SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
*s5,s6
DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
*-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1,
*.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
DATA r1,r2,r3,r4,r5,r6/-2957821389.d0,7062834065.d0,
*-512359803.6d0,10879881.29d0,-86327.92757d0,228.4622733d0/,s1,s2,
*s3,s4,s5,s6/40076544269.d0,745249964.8d0,7189466.438d0,
*47447.26470d0,226.1030244d0,1.d0/
if(x.lt.8.)then
y=x**2
bessy0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y*
*(s4+y*(s5+y*s6)))))+.636619772*bessj0(x)*log(x)
else
z=8./x
y=z**2
xx=x-.785398164
bessy0=sqrt(.636619772/x)*(sin(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
*p5))))+z*cos(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
endif
return
END
@@ -0,0 +1,28 @@
FUNCTION bessy1(x)
REAL bessy1,x
CU USES bessj1
REAL xx,z,bessj1
DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
*s1,s2,s3,s4,s5,s6,s7,y
SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
*s5,s6,s7
DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4,
*.2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/.04687499995d0,
*-.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/
DATA r1,r2,r3,r4,r5,r6/-.4900604943d13,.1275274390d13,
*-.5153438139d11,.7349264551d9,-.4237922726d7,.8511937935d4/,s1,s2,
*s3,s4,s5,s6,s7/.2499580570d14,.4244419664d12,.3733650367d10,
*.2245904002d8,.1020426050d6,.3549632885d3,1.d0/
if(x.lt.8.)then
y=x**2
bessy1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+
*y*(s4+y*(s5+y*(s6+y*s7))))))+.636619772*(bessj1(x)*log(x)-1./x)
else
z=8./x
y=z**2
xx=x-2.356194491
bessy1=sqrt(.636619772/x)*(sin(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
*p5))))+z*cos(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
endif
return
END
@@ -0,0 +1,7 @@
FUNCTION beta(z,w)
REAL beta,w,z
CU USES gammln
REAL gammln
beta=exp(gammln(z)+gammln(w)-gammln(z+w))
return
END
@@ -0,0 +1,37 @@
FUNCTION betacf(a,b,x)
INTEGER MAXIT
REAL betacf,a,b,x,EPS,FPMIN
PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30)
INTEGER m,m2
REAL aa,c,d,del,h,qab,qam,qap
qab=a+b
qap=a+1.
qam=a-1.
c=1.
d=1.-qab*x/qap
if(abs(d).lt.FPMIN)d=FPMIN
d=1./d
h=d
do 11 m=1,MAXIT
m2=2*m
aa=m*(b-m)*x/((qam+m2)*(a+m2))
d=1.+aa*d
if(abs(d).lt.FPMIN)d=FPMIN
c=1.+aa/c
if(abs(c).lt.FPMIN)c=FPMIN
d=1./d
h=h*d*c
aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2))
d=1.+aa*d
if(abs(d).lt.FPMIN)d=FPMIN
c=1.+aa/c
if(abs(c).lt.FPMIN)c=FPMIN
d=1./d
del=d*c
h=h*del
if(abs(del-1.).lt.EPS)goto 1
11 continue
pause 'a or b too big, or MAXIT too small in betacf'
1 betacf=h
return
END
@@ -0,0 +1,18 @@
FUNCTION betai(a,b,x)
REAL betai,a,b,x
CU USES betacf,gammln
REAL bt,betacf,gammln
if(x.lt.0..or.x.gt.1.)pause 'bad argument x in betai'
if(x.eq.0..or.x.eq.1.)then
bt=0.
else
bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.-x))
endif
if(x.lt.(a+1.)/(a+b+2.))then
betai=bt*betacf(a,b,x)/a
return
else
betai=1.-bt*betacf(b,a,1.-x)/b
return
endif
END
@@ -0,0 +1,8 @@
FUNCTION bico(n,k)
INTEGER k,n
REAL bico
CU USES factln
REAL factln
bico=nint(exp(factln(n)-factln(k)-factln(n-k)))
return
END
@@ -0,0 +1,28 @@
SUBROUTINE bksub(ne,nb,jf,k1,k2,c,nci,ncj,nck)
INTEGER jf,k1,k2,nb,nci,ncj,nck,ne
REAL c(nci,ncj,nck)
INTEGER i,im,j,k,kp,nbf
REAL xx
nbf=ne-nb
im=1
do 13 k=k2,k1,-1
if (k.eq.k1) im=nbf+1
kp=k+1
do 12 j=1,nbf
xx=c(j,jf,kp)
do 11 i=im,ne
c(i,jf,k)=c(i,jf,k)-c(i,j,k)*xx
11 continue
12 continue
13 continue
do 16 k=k1,k2
kp=k+1
do 14 i=1,nb
c(i,1,k)=c(i+nbf,jf,k)
14 continue
do 15 i=1,nbf
c(i+nb,1,k)=c(i,jf,kp)
15 continue
16 continue
return
END
@@ -0,0 +1,54 @@
FUNCTION bnldev(pp,n,idum)
INTEGER idum,n
REAL bnldev,pp,PI
CU USES gammln,ran1
PARAMETER (PI=3.141592654)
INTEGER j,nold
REAL am,em,en,g,oldg,p,pc,pclog,plog,pold,sq,t,y,gammln,ran1
SAVE nold,pold,pc,plog,pclog,en,oldg
DATA nold /-1/, pold /-1./
if(pp.le.0.5)then
p=pp
else
p=1.-pp
endif
am=n*p
if (n.lt.25)then
bnldev=0.
do 11 j=1,n
if(ran1(idum).lt.p)bnldev=bnldev+1.
11 continue
else if (am.lt.1.) then
g=exp(-am)
t=1.
do 12 j=0,n
t=t*ran1(idum)
if (t.lt.g) goto 1
12 continue
j=n
1 bnldev=j
else
if (n.ne.nold) then
en=n
oldg=gammln(en+1.)
nold=n
endif
if (p.ne.pold) then
pc=1.-p
plog=log(p)
pclog=log(pc)
pold=p
endif
sq=sqrt(2.*am*pc)
2 y=tan(PI*ran1(idum))
em=sq*y+am
if (em.lt.0..or.em.ge.en+1.) goto 2
em=int(em)
t=1.2*sq*(1.+y**2)*exp(oldg-gammln(em+1.)-gammln(en-em+1.)+em*
*plog+(en-em)*pclog)
if (ran1(idum).gt.t) goto 2
bnldev=em
endif
if (p.ne.pp) bnldev=n-bnldev
return
END
@@ -0,0 +1,83 @@
FUNCTION brent(ax,bx,cx,f,tol,xmin)
INTEGER ITMAX
REAL brent,ax,bx,cx,tol,xmin,f,CGOLD,ZEPS
EXTERNAL f
PARAMETER (ITMAX=100,CGOLD=.3819660,ZEPS=1.0e-10)
INTEGER iter
REAL a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
a=min(ax,cx)
b=max(ax,cx)
v=bx
w=v
x=v
e=0.
fx=f(x)
fv=fx
fw=fx
do 11 iter=1,ITMAX
xm=0.5*(a+b)
tol1=tol*abs(x)+ZEPS
tol2=2.*tol1
if(abs(x-xm).le.(tol2-.5*(b-a))) goto 3
if(abs(e).gt.tol1) then
r=(x-w)*(fx-fv)
q=(x-v)*(fx-fw)
p=(x-v)*q-(x-w)*r
q=2.*(q-r)
if(q.gt.0.) p=-p
q=abs(q)
etemp=e
e=d
if(abs(p).ge.abs(.5*q*etemp).or.p.le.q*(a-x).or.p.ge.q*(b-x))
*goto 1
d=p/q
u=x+d
if(u-a.lt.tol2 .or. b-u.lt.tol2) d=sign(tol1,xm-x)
goto 2
endif
1 if(x.ge.xm) then
e=a-x
else
e=b-x
endif
d=CGOLD*e
2 if(abs(d).ge.tol1) then
u=x+d
else
u=x+sign(tol1,d)
endif
fu=f(u)
if(fu.le.fx) then
if(u.ge.x) then
a=x
else
b=x
endif
v=w
fv=fw
w=x
fw=fx
x=u
fx=fu
else
if(u.lt.x) then
a=u
else
b=u
endif
if(fu.le.fw .or. w.eq.x) then
v=w
fv=fw
w=u
fw=fu
else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then
v=u
fv=fu
endif
endif
11 continue
pause 'brent exceed maximum iterations'
3 xmin=x
brent=fx
return
END
@@ -0,0 +1,170 @@
SUBROUTINE broydn(x,n,check)
INTEGER n,nn,NP,MAXITS
REAL x(n),fvec,EPS,TOLF,TOLMIN,TOLX,STPMX
LOGICAL check
PARAMETER (NP=40,MAXITS=200,EPS=1.e-7,TOLF=1.e-4,TOLMIN=1.e-6,
*TOLX=EPS,STPMX=100.)
COMMON /newtv/ fvec(NP),nn
CU USES fdjac,fmin,lnsrch,qrdcmp,qrupdt,rsolv
INTEGER i,its,j,k
REAL den,f,fold,stpmax,sum,temp,test,c(NP),d(NP),fvcold(NP),g(NP),
*p(NP),qt(NP,NP),r(NP,NP),s(NP),t(NP),w(NP),xold(NP),fmin
LOGICAL restrt,sing,skip
EXTERNAL fmin
nn=n
f=fmin(x)
test=0.
do 11 i=1,n
if(abs(fvec(i)).gt.test)test=abs(fvec(i))
11 continue
if(test.lt..01*TOLF)then
check=.false.
return
endif
sum=0.
do 12 i=1,n
sum=sum+x(i)**2
12 continue
stpmax=STPMX*max(sqrt(sum),float(n))
restrt=.true.
do 44 its=1,MAXITS
if(restrt)then
call fdjac(n,x,fvec,NP,r)
call qrdcmp(r,n,NP,c,d,sing)
if(sing) pause 'singular Jacobian in broydn'
do 14 i=1,n
do 13 j=1,n
qt(i,j)=0.
13 continue
qt(i,i)=1.
14 continue
do 18 k=1,n-1
if(c(k).ne.0.)then
do 17 j=1,n
sum=0.
do 15 i=k,n
sum=sum+r(i,k)*qt(i,j)
15 continue
sum=sum/c(k)
do 16 i=k,n
qt(i,j)=qt(i,j)-sum*r(i,k)
16 continue
17 continue
endif
18 continue
do 21 i=1,n
r(i,i)=d(i)
do 19 j=1,i-1
r(i,j)=0.
19 continue
21 continue
else
do 22 i=1,n
s(i)=x(i)-xold(i)
22 continue
do 24 i=1,n
sum=0.
do 23 j=i,n
sum=sum+r(i,j)*s(j)
23 continue
t(i)=sum
24 continue
skip=.true.
do 26 i=1,n
sum=0.
do 25 j=1,n
sum=sum+qt(j,i)*t(j)
25 continue
w(i)=fvec(i)-fvcold(i)-sum
if(abs(w(i)).ge.EPS*(abs(fvec(i))+abs(fvcold(i))))then
skip=.false.
else
w(i)=0.
endif
26 continue
if(.not.skip)then
do 28 i=1,n
sum=0.
do 27 j=1,n
sum=sum+qt(i,j)*w(j)
27 continue
t(i)=sum
28 continue
den=0.
do 29 i=1,n
den=den+s(i)**2
29 continue
do 31 i=1,n
s(i)=s(i)/den
31 continue
call qrupdt(r,qt,n,NP,t,s)
do 32 i=1,n
if(r(i,i).eq.0.) pause 'r singular in broydn'
d(i)=r(i,i)
32 continue
endif
endif
do 34 i=1,n
sum=0.
do 33 j=1,n
sum=sum+qt(i,j)*fvec(j)
33 continue
g(i)=sum
34 continue
do 36 i=n,1,-1
sum=0.
do 35 j=1,i
sum=sum+r(j,i)*g(j)
35 continue
g(i)=sum
36 continue
do 37 i=1,n
xold(i)=x(i)
fvcold(i)=fvec(i)
37 continue
fold=f
do 39 i=1,n
sum=0.
do 38 j=1,n
sum=sum+qt(i,j)*fvec(j)
38 continue
p(i)=-sum
39 continue
call rsolv(r,n,NP,d,p)
call lnsrch(n,xold,fold,g,p,x,f,stpmax,check,fmin)
test=0.
do 41 i=1,n
if(abs(fvec(i)).gt.test)test=abs(fvec(i))
41 continue
if(test.lt.TOLF)then
check=.false.
return
endif
if(check)then
if(restrt)then
return
else
test=0.
den=max(f,.5*n)
do 42 i=1,n
temp=abs(g(i))*max(abs(x(i)),1.)/den
if(temp.gt.test)test=temp
42 continue
if(test.lt.TOLMIN)then
return
else
restrt=.true.
endif
endif
else
restrt=.false.
test=0.
do 43 i=1,n
temp=(abs(x(i)-xold(i)))/max(abs(x(i)),1.)
if(temp.gt.test)test=temp
43 continue
if(test.lt.TOLX)return
endif
44 continue
pause 'MAXITS exceeded in broydn'
END
@@ -0,0 +1,109 @@
SUBROUTINE bsstep(y,dydx,nv,x,htry,eps,yscal,hdid,hnext,derivs)
INTEGER nv,NMAX,KMAXX,IMAX
REAL eps,hdid,hnext,htry,x,dydx(nv),y(nv),yscal(nv),SAFE1,SAFE2,
*REDMAX,REDMIN,TINY,SCALMX
PARAMETER (NMAX=50,KMAXX=8,IMAX=KMAXX+1,SAFE1=.25,SAFE2=.7,
*REDMAX=1.e-5,REDMIN=.7,TINY=1.e-30,SCALMX=.1)
CU USES derivs,mmid,pzextr
INTEGER i,iq,k,kk,km,kmax,kopt,nseq(IMAX)
REAL eps1,epsold,errmax,fact,h,red,scale,work,wrkmin,xest,xnew,
*a(IMAX),alf(KMAXX,KMAXX),err(KMAXX),yerr(NMAX),ysav(NMAX),
*yseq(NMAX)
LOGICAL first,reduct
SAVE a,alf,epsold,first,kmax,kopt,nseq,xnew
EXTERNAL derivs
DATA first/.true./,epsold/-1./
DATA nseq /2,4,6,8,10,12,14,16,18/
if(eps.ne.epsold)then
hnext=-1.e29
xnew=-1.e29
eps1=SAFE1*eps
a(1)=nseq(1)+1
do 11 k=1,KMAXX
a(k+1)=a(k)+nseq(k+1)
11 continue
do 13 iq=2,KMAXX
do 12 k=1,iq-1
alf(k,iq)=eps1**((a(k+1)-a(iq+1))/((a(iq+1)-a(1)+1.)*(2*k+
*1)))
12 continue
13 continue
epsold=eps
do 14 kopt=2,KMAXX-1
if(a(kopt+1).gt.a(kopt)*alf(kopt-1,kopt))goto 1
14 continue
1 kmax=kopt
endif
h=htry
do 15 i=1,nv
ysav(i)=y(i)
15 continue
if(h.ne.hnext.or.x.ne.xnew)then
first=.true.
kopt=kmax
endif
reduct=.false.
2 do 17 k=1,kmax
xnew=x+h
if(xnew.eq.x)pause 'step size underflow in bsstep'
call mmid(ysav,dydx,nv,x,h,nseq(k),yseq,derivs)
xest=(h/nseq(k))**2
call pzextr(k,xest,yseq,y,yerr,nv)
if(k.ne.1)then
errmax=TINY
do 16 i=1,nv
errmax=max(errmax,abs(yerr(i)/yscal(i)))
16 continue
errmax=errmax/eps
km=k-1
err(km)=(errmax/SAFE1)**(1./(2*km+1))
endif
if(k.ne.1.and.(k.ge.kopt-1.or.first))then
if(errmax.lt.1.)goto 4
if(k.eq.kmax.or.k.eq.kopt+1)then
red=SAFE2/err(km)
goto 3
else if(k.eq.kopt)then
if(alf(kopt-1,kopt).lt.err(km))then
red=1./err(km)
goto 3
endif
else if(kopt.eq.kmax)then
if(alf(km,kmax-1).lt.err(km))then
red=alf(km,kmax-1)*SAFE2/err(km)
goto 3
endif
else if(alf(km,kopt).lt.err(km))then
red=alf(km,kopt-1)/err(km)
goto 3
endif
endif
17 continue
3 red=min(red,REDMIN)
red=max(red,REDMAX)
h=h*red
reduct=.true.
goto 2
4 x=xnew
hdid=h
first=.false.
wrkmin=1.e35
do 18 kk=1,km
fact=max(err(kk),SCALMX)
work=fact*a(kk+1)
if(work.lt.wrkmin)then
scale=fact
wrkmin=work
kopt=kk+1
endif
18 continue
hnext=h/scale
if(kopt.ge.k.and.kopt.ne.kmax.and..not.reduct)then
fact=max(scale/alf(kopt-1,kopt),SCALMX)
if(a(kopt+1)*fact.le.wrkmin)then
hnext=h/fact
kopt=kopt+1
endif
endif
return
END
@@ -0,0 +1,22 @@
SUBROUTINE caldat(julian,mm,id,iyyy)
INTEGER id,iyyy,julian,mm,IGREG
PARAMETER (IGREG=2299161)
INTEGER ja,jalpha,jb,jc,jd,je
if(julian.ge.IGREG)then
jalpha=int(((julian-1867216)-0.25)/36524.25)
ja=julian+1+jalpha-int(0.25*jalpha)
else
ja=julian
endif
jb=ja+1524
jc=int(6680.+((jb-2439870)-122.1)/365.25)
jd=365*jc+int(0.25*jc)
je=int((jb-jd)/30.6001)
id=jb-jd-int(30.6001*je)
mm=je-1
if(mm.gt.12)mm=mm-12
iyyy=jc-4715
if(mm.gt.2)iyyy=iyyy-1
if(iyyy.le.0)iyyy=iyyy-1
return
END
@@ -0,0 +1,16 @@
SUBROUTINE chder(a,b,c,cder,n)
INTEGER n
REAL a,b,c(n),cder(n)
INTEGER j
REAL con
cder(n)=0.
cder(n-1)=2*(n-1)*c(n)
do 11 j=n-2,1,-1
cder(j)=cder(j+2)+2*j*c(j+1)
11 continue
con=2./(b-a)
do 12 j=1,n
cder(j)=cder(j)*con
12 continue
return
END
@@ -0,0 +1,18 @@
FUNCTION chebev(a,b,c,m,x)
INTEGER m
REAL chebev,a,b,x,c(m)
INTEGER j
REAL d,dd,sv,y,y2
if ((x-a)*(x-b).gt.0.) pause 'x not in range in chebev'
d=0.
dd=0.
y=(2.*x-a-b)/(b-a)
y2=2.*y
do 11 j=m,2,-1
sv=d
d=y2*d-dd+c(j)
dd=sv
11 continue
chebev=y*d-dd+0.5*c(1)
return
END
@@ -0,0 +1,24 @@
SUBROUTINE chebft(a,b,c,n,func)
INTEGER n,NMAX
REAL a,b,c(n),func,PI
EXTERNAL func
PARAMETER (NMAX=50, PI=3.141592653589793d0)
INTEGER j,k
REAL bma,bpa,fac,y,f(NMAX)
DOUBLE PRECISION sum
bma=0.5*(b-a)
bpa=0.5*(b+a)
do 11 k=1,n
y=cos(PI*(k-0.5)/n)
f(k)=func(y*bma+bpa)
11 continue
fac=2./n
do 13 j=1,n
sum=0.d0
do 12 k=1,n
sum=sum+f(k)*cos((PI*(j-1))*((k-0.5d0)/n))
12 continue
c(j)=fac*sum
13 continue
return
END
@@ -0,0 +1,27 @@
SUBROUTINE chebpc(c,d,n)
INTEGER n,NMAX
REAL c(n),d(n)
PARAMETER (NMAX=50)
INTEGER j,k
REAL sv,dd(NMAX)
do 11 j=1,n
d(j)=0.
dd(j)=0.
11 continue
d(1)=c(n)
do 13 j=n-1,2,-1
do 12 k=n-j+1,2,-1
sv=d(k)
d(k)=2.*d(k-1)-dd(k)
dd(k)=sv
12 continue
sv=d(1)
d(1)=-dd(1)+c(j)
dd(1)=sv
13 continue
do 14 j=n,2,-1
d(j)=d(j-1)-dd(j)
14 continue
d(1)=-dd(1)+0.5*c(1)
return
END
@@ -0,0 +1,18 @@
SUBROUTINE chint(a,b,c,cint,n)
INTEGER n
REAL a,b,c(n),cint(n)
INTEGER j
REAL con,fac,sum
con=0.25*(b-a)
sum=0.
fac=1.
do 11 j=2,n-1
cint(j)=con*(c(j-1)-c(j+1))/(j-1)
sum=sum+fac*cint(j)
fac=-fac
11 continue
cint(n)=con*c(n-1)/(n-1)
sum=sum+fac*cint(n)
cint(1)=2.*sum
return
END
@@ -0,0 +1,32 @@
FUNCTION chixy(bang)
REAL chixy,bang,BIG
INTEGER NMAX
PARAMETER (NMAX=1000,BIG=1.E30)
INTEGER nn,j
REAL xx(NMAX),yy(NMAX),sx(NMAX),sy(NMAX),ww(NMAX),aa,offs,avex,
*avey,sumw,b
COMMON /fitxyc/ xx,yy,sx,sy,ww,aa,offs,nn
b=tan(bang)
avex=0.
avey=0.
sumw=0.
do 11 j=1,nn
ww(j)=(b*sx(j))**2+sy(j)**2
if(ww(j).lt.1./BIG) then
ww(j)=BIG
else
ww(j)=1./ww(j)
endif
sumw=sumw+ww(j)
avex=avex+ww(j)*xx(j)
avey=avey+ww(j)*yy(j)
11 continue
avex=avex/sumw
avey=avey/sumw
aa=avey-b*avex
chixy=-offs
do 12 j=1,nn
chixy=chixy+ww(j)*(yy(j)-aa-b*xx(j))**2
12 continue
return
END
@@ -0,0 +1,21 @@
SUBROUTINE choldc(a,n,np,p)
INTEGER n,np
REAL a(np,np),p(n)
INTEGER i,j,k
REAL sum
do 13 i=1,n
do 12 j=i,n
sum=a(i,j)
do 11 k=i-1,1,-1
sum=sum-a(i,k)*a(j,k)
11 continue
if(i.eq.j)then
if(sum.le.0.)pause 'choldc failed'
p(i)=sqrt(sum)
else
a(j,i)=sum/p(i)
endif
12 continue
13 continue
return
END
@@ -0,0 +1,21 @@
SUBROUTINE cholsl(a,n,np,p,b,x)
INTEGER n,np
REAL a(np,np),b(n),p(n),x(n)
INTEGER i,k
REAL sum
do 12 i=1,n
sum=b(i)
do 11 k=i-1,1,-1
sum=sum-a(i,k)*x(k)
11 continue
x(i)=sum/p(i)
12 continue
do 14 i=n,1,-1
sum=x(i)
do 13 k=i+1,n
sum=sum-a(k,i)*x(k)
13 continue
x(i)=sum/p(i)
14 continue
return
END
@@ -0,0 +1,15 @@
SUBROUTINE chsone(bins,ebins,nbins,knstrn,df,chsq,prob)
INTEGER knstrn,nbins
REAL chsq,df,prob,bins(nbins),ebins(nbins)
CU USES gammq
INTEGER j
REAL gammq
df=nbins-knstrn
chsq=0.
do 11 j=1,nbins
if(ebins(j).le.0.)pause 'bad expected number in chsone'
chsq=chsq+(bins(j)-ebins(j))**2/ebins(j)
11 continue
prob=gammq(0.5*df,0.5*chsq)
return
END
@@ -0,0 +1,18 @@
SUBROUTINE chstwo(bins1,bins2,nbins,knstrn,df,chsq,prob)
INTEGER knstrn,nbins
REAL chsq,df,prob,bins1(nbins),bins2(nbins)
CU USES gammq
INTEGER j
REAL gammq
df=nbins-knstrn
chsq=0.
do 11 j=1,nbins
if(bins1(j).eq.0..and.bins2(j).eq.0.)then
df=df-1.
else
chsq=chsq+(bins1(j)-bins2(j))**2/(bins1(j)+bins2(j))
endif
11 continue
prob=gammq(0.5*df,0.5*chsq)
return
END
@@ -0,0 +1,70 @@
SUBROUTINE cisi(x,ci,si)
INTEGER MAXIT
REAL ci,si,x,EPS,EULER,PIBY2,FPMIN,TMIN
PARAMETER (EPS=6.e-8,EULER=.57721566,MAXIT=100,PIBY2=1.5707963,
*FPMIN=1.e-30,TMIN=2.)
INTEGER i,k
REAL a,err,fact,sign,sum,sumc,sums,t,term,absc
COMPLEX h,b,c,d,del
LOGICAL odd
absc(h)=abs(real(h))+abs(aimag(h))
t=abs(x)
if(t.eq.0.)then
si=0.
ci=-1./FPMIN
return
endif
if(t.gt.TMIN)then
b=cmplx(1.,t)
c=1./FPMIN
d=1./b
h=d
do 11 i=2,MAXIT
a=-(i-1)**2
b=b+2.
d=1./(a*d+b)
c=b+a/c
del=c*d
h=h*del
if(absc(del-1.).lt.EPS)goto 1
11 continue
pause 'cf failed in cisi'
1 continue
h=cmplx(cos(t),-sin(t))*h
ci=-real(h)
si=PIBY2+aimag(h)
else
if(t.lt.sqrt(FPMIN))then
sumc=0.
sums=t
else
sum=0.
sums=0.
sumc=0.
sign=1.
fact=1.
odd=.true.
do 12 k=1,MAXIT
fact=fact*t/k
term=fact/k
sum=sum+sign*term
err=term/abs(sum)
if(odd)then
sign=-sign
sums=sum
sum=sumc
else
sumc=sum
sum=sums
endif
if(err.lt.EPS)goto 2
odd=.not.odd
12 continue
pause 'maxits exceeded in cisi'
endif
2 si=sums
ci=sumc+log(t)+EULER
endif
if(x.lt.0.)si=-si
return
END
@@ -0,0 +1,38 @@
SUBROUTINE cntab1(nn,ni,nj,chisq,df,prob,cramrv,ccc)
INTEGER ni,nj,nn(ni,nj),MAXI,MAXJ
REAL ccc,chisq,cramrv,df,prob,TINY
PARAMETER (MAXI=100,MAXJ=100,TINY=1.e-30)
CU USES gammq
INTEGER i,j,nni,nnj
REAL expctd,sum,sumi(MAXI),sumj(MAXJ),gammq
sum=0
nni=ni
nnj=nj
do 12 i=1,ni
sumi(i)=0.
do 11 j=1,nj
sumi(i)=sumi(i)+nn(i,j)
sum=sum+nn(i,j)
11 continue
if(sumi(i).eq.0.)nni=nni-1
12 continue
do 14 j=1,nj
sumj(j)=0.
do 13 i=1,ni
sumj(j)=sumj(j)+nn(i,j)
13 continue
if(sumj(j).eq.0.)nnj=nnj-1
14 continue
df=nni*nnj-nni-nnj+1
chisq=0.
do 16 i=1,ni
do 15 j=1,nj
expctd=sumj(j)*sumi(i)/sum
chisq=chisq+(nn(i,j)-expctd)**2/(expctd+TINY)
15 continue
16 continue
prob=gammq(0.5*df,0.5*chisq)
cramrv=sqrt(chisq/(sum*min(nni-1,nnj-1)))
ccc=sqrt(chisq/(chisq+sum))
return
END
@@ -0,0 +1,50 @@
SUBROUTINE cntab2(nn,ni,nj,h,hx,hy,hygx,hxgy,uygx,uxgy,uxy)
INTEGER ni,nj,nn(ni,nj),MAXI,MAXJ
REAL h,hx,hxgy,hy,hygx,uxgy,uxy,uygx,TINY
PARAMETER (MAXI=100,MAXJ=100,TINY=1.e-30)
INTEGER i,j
REAL p,sum,sumi(MAXI),sumj(MAXJ)
sum=0
do 12 i=1,ni
sumi(i)=0.0
do 11 j=1,nj
sumi(i)=sumi(i)+nn(i,j)
sum=sum+nn(i,j)
11 continue
12 continue
do 14 j=1,nj
sumj(j)=0.
do 13 i=1,ni
sumj(j)=sumj(j)+nn(i,j)
13 continue
14 continue
hx=0.
do 15 i=1,ni
if(sumi(i).ne.0.)then
p=sumi(i)/sum
hx=hx-p*log(p)
endif
15 continue
hy=0.
do 16 j=1,nj
if(sumj(j).ne.0.)then
p=sumj(j)/sum
hy=hy-p*log(p)
endif
16 continue
h=0.
do 18 i=1,ni
do 17 j=1,nj
if(nn(i,j).ne.0)then
p=nn(i,j)/sum
h=h-p*log(p)
endif
17 continue
18 continue
hygx=h-hx
hxgy=h-hy
uygx=(hy-hygx)/(hy+TINY)
uxgy=(hx-hxgy)/(hx+TINY)
uxy=2.*(hx+hy-h)/(hx+hy+TINY)
return
END
@@ -0,0 +1,31 @@
SUBROUTINE convlv(data,n,respns,m,isign,ans)
INTEGER isign,m,n,NMAX
REAL data(n),respns(n)
COMPLEX ans(n)
PARAMETER (NMAX=4096)
CU USES realft,twofft
INTEGER i,no2
COMPLEX fft(NMAX)
do 11 i=1,(m-1)/2
respns(n+1-i)=respns(m+1-i)
11 continue
do 12 i=(m+3)/2,n-(m-1)/2
respns(i)=0.0
12 continue
call twofft(data,respns,fft,ans,n)
no2=n/2
do 13 i=1,no2+1
if (isign.eq.1) then
ans(i)=fft(i)*ans(i)/no2
else if (isign.eq.-1) then
if (abs(ans(i)).eq.0.0) pause
*'deconvolving at response zero in convlv'
ans(i)=fft(i)/ans(i)/no2
else
pause 'no meaning for isign in convlv'
endif
13 continue
ans(1)=cmplx(real(ans(1)),real(ans(no2+1)))
call realft(ans,n,-1)
return
END
@@ -0,0 +1,11 @@
SUBROUTINE copy(aout,ain,n)
INTEGER n
DOUBLE PRECISION ain(n,n),aout(n,n)
INTEGER i,j
do 12 i=1,n
do 11 j=1,n
aout(j,i)=ain(j,i)
11 continue
12 continue
return
END
@@ -0,0 +1,17 @@
SUBROUTINE correl(data1,data2,n,ans)
INTEGER n,NMAX
REAL data1(n),data2(n)
COMPLEX ans(n)
PARAMETER (NMAX=4096)
CU USES realft,twofft
INTEGER i,no2
COMPLEX fft(NMAX)
call twofft(data1,data2,fft,ans,n)
no2=n/2
do 11 i=1,no2+1
ans(i)=fft(i)*conjg(ans(i))/float(no2)
11 continue
ans(1)=cmplx(real(ans(1)),real(ans(no2+1)))
call realft(ans,n,-1)
return
END
@@ -0,0 +1,33 @@
SUBROUTINE cosft1(y,n)
INTEGER n
REAL y(n+1)
CU USES realft
INTEGER j
REAL sum,y1,y2
DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
theta=3.141592653589793d0/n
wr=1.0d0
wi=0.0d0
wpr=-2.0d0*sin(0.5d0*theta)**2
wpi=sin(theta)
sum=0.5*(y(1)-y(n+1))
y(1)=0.5*(y(1)+y(n+1))
do 11 j=1,n/2-1
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
y1=0.5*(y(j+1)+y(n-j+1))
y2=(y(j+1)-y(n-j+1))
y(j+1)=y1-wi*y2
y(n-j+1)=y1+wi*y2
sum=sum+wr*y2
11 continue
call realft(y,n,+1)
y(n+1)=y(2)
y(2)=sum
do 12 j=4,n,2
sum=sum+y(j)
y(j)=sum
12 continue
return
END
@@ -0,0 +1,69 @@
SUBROUTINE cosft2(y,n,isign)
INTEGER isign,n
REAL y(n)
CU USES realft
INTEGER i
REAL sum,sum1,y1,y2,ytemp
DOUBLE PRECISION theta,wi,wi1,wpi,wpr,wr,wr1,wtemp,PI
PARAMETER (PI=3.141592653589793d0)
theta=0.5d0*PI/n
wr=1.0d0
wi=0.0d0
wr1=cos(theta)
wi1=sin(theta)
wpr=-2.0d0*wi1**2
wpi=sin(2.d0*theta)
if(isign.eq.1)then
do 11 i=1,n/2
y1=0.5*(y(i)+y(n-i+1))
y2=wi1*(y(i)-y(n-i+1))
y(i)=y1+y2
y(n-i+1)=y1-y2
wtemp=wr1
wr1=wr1*wpr-wi1*wpi+wr1
wi1=wi1*wpr+wtemp*wpi+wi1
11 continue
call realft(y,n,1)
do 12 i=3,n,2
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
y1=y(i)*wr-y(i+1)*wi
y2=y(i+1)*wr+y(i)*wi
y(i)=y1
y(i+1)=y2
12 continue
sum=0.5*y(2)
do 13 i=n,2,-2
sum1=sum
sum=sum+y(i)
y(i)=sum1
13 continue
else if(isign.eq.-1)then
ytemp=y(n)
do 14 i=n,4,-2
y(i)=y(i-2)-y(i)
14 continue
y(2)=2.0*ytemp
do 15 i=3,n,2
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
y1=y(i)*wr+y(i+1)*wi
y2=y(i+1)*wr-y(i)*wi
y(i)=y1
y(i+1)=y2
15 continue
call realft(y,n,-1)
do 16 i=1,n/2
y1=y(i)+y(n-i+1)
y2=(0.5/wi1)*(y(i)-y(n-i+1))
y(i)=0.5*(y1+y2)
y(n-i+1)=0.5*(y1-y2)
wtemp=wr1
wr1=wr1*wpr-wi1*wpi+wr1
wi1=wi1*wpr+wtemp*wpi+wi1
16 continue
endif
return
END
@@ -0,0 +1,29 @@
SUBROUTINE covsrt(covar,npc,ma,ia,mfit)
INTEGER ma,mfit,npc,ia(ma)
REAL covar(npc,npc)
INTEGER i,j,k
REAL swap
do 12 i=mfit+1,ma
do 11 j=1,i
covar(i,j)=0.
covar(j,i)=0.
11 continue
12 continue
k=mfit
do 15 j=ma,1,-1
if(ia(j).ne.0)then
do 13 i=1,ma
swap=covar(i,k)
covar(i,k)=covar(i,j)
covar(i,j)=swap
13 continue
do 14 i=1,ma
swap=covar(k,i)
covar(k,i)=covar(j,i)
covar(j,i)=swap
14 continue
k=k-1
endif
15 continue
return
END
@@ -0,0 +1,29 @@
SUBROUTINE crank(n,w,s)
INTEGER n
REAL s,w(n)
INTEGER j,ji,jt
REAL rank,t
s=0.
j=1
1 if(j.lt.n)then
if(w(j+1).ne.w(j))then
w(j)=j
j=j+1
else
do 11 jt=j+1,n
if(w(jt).ne.w(j))goto 2
11 continue
jt=n+1
2 rank=0.5*(j+jt-1)
do 12 ji=j,jt-1
w(ji)=rank
12 continue
t=jt-j
s=s+t**3-t
j=jt
endif
goto 1
endif
if(j.eq.n)w(n)=n
return
END
@@ -0,0 +1,28 @@
SUBROUTINE cyclic(a,b,c,alpha,beta,r,x,n)
INTEGER n,NMAX
REAL alpha,beta,a(n),b(n),c(n),r(n),x(n)
PARAMETER (NMAX=500)
CU USES tridag
INTEGER i
REAL fact,gamma,bb(NMAX),u(NMAX),z(NMAX)
if(n.le.2)pause 'n too small in cyclic'
if(n.gt.NMAX)pause 'NMAX too small in cyclic'
gamma=-b(1)
bb(1)=b(1)-gamma
bb(n)=b(n)-alpha*beta/gamma
do 11 i=2,n-1
bb(i)=b(i)
11 continue
call tridag(a,bb,c,r,x,n)
u(1)=gamma
u(n)=alpha
do 12 i=2,n-1
u(i)=0.
12 continue
call tridag(a,bb,c,u,z,n)
fact=(x(1)+beta*x(n)/gamma)/(1.+z(1)+beta*z(n)/gamma)
do 13 i=1,n
x(i)=x(i)-fact*z(i)
13 continue
return
END
@@ -0,0 +1,35 @@
SUBROUTINE daub4(a,n,isign)
INTEGER n,isign,NMAX
REAL a(n),C3,C2,C1,C0
PARAMETER (C0=0.4829629131445341,C1=0.8365163037378079,
*C2=0.2241438680420134,C3=-0.1294095225512604,NMAX=1024)
REAL wksp(NMAX)
INTEGER nh,nh1,i,j
if(n.lt.4)return
if(n.gt.NMAX) pause 'wksp too small in daub4'
nh=n/2
nh1=nh+1
if (isign.ge.0) then
i=1
do 11 j=1,n-3,2
wksp(i)=C0*a(j)+C1*a(j+1)+C2*a(j+2)+C3*a(j+3)
wksp(i+nh)=C3*a(j)-C2*a(j+1)+C1*a(j+2)-C0*a(j+3)
i=i+1
11 continue
wksp(i)=C0*a(n-1)+C1*a(n)+C2*a(1)+C3*a(2)
wksp(i+nh)=C3*a(n-1)-C2*a(n)+C1*a(1)-C0*a(2)
else
wksp(1)=C2*a(nh)+C1*a(n)+C0*a(1)+C3*a(nh1)
wksp(2)=C3*a(nh)-C0*a(n)+C1*a(1)-C2*a(nh1)
j=3
do 12 i=1,nh-1
wksp(j)=C2*a(i)+C1*a(i+nh)+C0*a(i+1)+C3*a(i+nh1)
wksp(j+1)=C3*a(i)-C0*a(i+nh)+C1*a(i+1)-C2*a(i+nh1)
j=j+2
12 continue
endif
do 13 i=1,n
a(i)=wksp(i)
13 continue
return
END
@@ -0,0 +1,36 @@
FUNCTION dawson(x)
INTEGER NMAX
REAL dawson,x,H,A1,A2,A3
PARAMETER (NMAX=6,H=0.4,A1=2./3.,A2=0.4,A3=2./7.)
INTEGER i,init,n0
REAL d1,d2,e1,e2,sum,x2,xp,xx,c(NMAX)
SAVE init,c
DATA init/0/
if(init.eq.0)then
init=1
do 11 i=1,NMAX
c(i)=exp(-((2.*float(i)-1.)*H)**2)
11 continue
endif
if(abs(x).lt.0.2)then
x2=x**2
dawson=x*(1.-A1*x2*(1.-A2*x2*(1.-A3*x2)))
else
xx=abs(x)
n0=2*nint(0.5*xx/H)
xp=xx-float(n0)*H
e1=exp(2.*xp*H)
e2=e1**2
d1=float(n0+1)
d2=d1-2.
sum=0.
do 12 i=1,NMAX
sum=sum+c(i)*(e1/d1+1./(d2*e1))
d1=d1+2.
d2=d2-2.
e1=e2*e1
12 continue
dawson=0.5641895835*sign(exp(-xp**2),x)*sum
endif
return
END
@@ -0,0 +1,110 @@
FUNCTION dbrent(ax,bx,cx,f,df,tol,xmin)
INTEGER ITMAX
REAL dbrent,ax,bx,cx,tol,xmin,df,f,ZEPS
EXTERNAL df,f
PARAMETER (ITMAX=100,ZEPS=1.0e-10)
INTEGER iter
REAL a,b,d,d1,d2,du,dv,dw,dx,e,fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,
*v,w,x,xm
LOGICAL ok1,ok2
a=min(ax,cx)
b=max(ax,cx)
v=bx
w=v
x=v
e=0.
fx=f(x)
fv=fx
fw=fx
dx=df(x)
dv=dx
dw=dx
do 11 iter=1,ITMAX
xm=0.5*(a+b)
tol1=tol*abs(x)+ZEPS
tol2=2.*tol1
if(abs(x-xm).le.(tol2-.5*(b-a))) goto 3
if(abs(e).gt.tol1) then
d1=2.*(b-a)
d2=d1
if(dw.ne.dx) d1=(w-x)*dx/(dx-dw)
if(dv.ne.dx) d2=(v-x)*dx/(dx-dv)
u1=x+d1
u2=x+d2
ok1=((a-u1)*(u1-b).gt.0.).and.(dx*d1.le.0.)
ok2=((a-u2)*(u2-b).gt.0.).and.(dx*d2.le.0.)
olde=e
e=d
if(.not.(ok1.or.ok2))then
goto 1
else if (ok1.and.ok2)then
if(abs(d1).lt.abs(d2))then
d=d1
else
d=d2
endif
else if (ok1)then
d=d1
else
d=d2
endif
if(abs(d).gt.abs(0.5*olde))goto 1
u=x+d
if(u-a.lt.tol2 .or. b-u.lt.tol2) d=sign(tol1,xm-x)
goto 2
endif
1 if(dx.ge.0.) then
e=a-x
else
e=b-x
endif
d=0.5*e
2 if(abs(d).ge.tol1) then
u=x+d
fu=f(u)
else
u=x+sign(tol1,d)
fu=f(u)
if(fu.gt.fx)goto 3
endif
du=df(u)
if(fu.le.fx) then
if(u.ge.x) then
a=x
else
b=x
endif
v=w
fv=fw
dv=dw
w=x
fw=fx
dw=dx
x=u
fx=fu
dx=du
else
if(u.lt.x) then
a=u
else
b=u
endif
if(fu.le.fw .or. w.eq.x) then
v=w
fv=fw
dv=dw
w=u
fw=fu
dw=du
else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then
v=u
fv=fu
dv=du
endif
endif
11 continue
pause 'dbrent exceeded maximum iterations'
3 xmin=x
dbrent=fx
return
END
@@ -0,0 +1,23 @@
SUBROUTINE ddpoly(c,nc,x,pd,nd)
INTEGER nc,nd
REAL x,c(nc),pd(nd)
INTEGER i,j,nnd
REAL const
pd(1)=c(nc)
do 11 j=2,nd
pd(j)=0.
11 continue
do 13 i=nc-1,1,-1
nnd=min(nd,nc+1-i)
do 12 j=nnd,2,-1
pd(j)=pd(j)*x+pd(j-1)
12 continue
pd(1)=pd(1)*x+c(i)
13 continue
const=2.
do 14 i=3,nd
pd(i)=const*pd(i)
const=const*i
14 continue
return
END
@@ -0,0 +1,27 @@
LOGICAL FUNCTION decchk(string,n,ch)
INTEGER n
CHARACTER string*(*),ch*1
INTEGER ij(10,10),ip(10,8),i,j,k,m
SAVE ij,ip
DATA ip/0,1,2,3,4,5,6,7,8,9,1,5,7,6,2,8,3,0,9,4,5,8,0,3,7,9,6,1,4,
*2,8,9,1,6,0,4,3,5,2,7,9,4,5,3,1,2,6,8,7,0,4,2,8,6,5,7,3,9,0,1,2,7,
*9,3,8,0,6,4,1,5,7,0,4,6,9,1,3,2,5,8/,ij/0,1,2,3,4,5,6,7,8,9,1,2,3,
*4,0,9,5,6,7,8,2,3,4,0,1,8,9,5,6,7,3,4,0,1,2,7,8,9,5,6,4,0,1,2,3,6,
*7,8,9,5,5,6,7,8,9,0,1,2,3,4,6,7,8,9,5,4,0,1,2,3,7,8,9,5,6,3,4,0,1,
*2,8,9,5,6,7,2,3,4,0,1,9,5,6,7,8,1,2,3,4,0/
k=0
m=0
do 11 j=1,n
i=ichar(string(j:j))
if (i.ge.48.and.i.le.57)then
k=ij(k+1,ip(mod(i+2,10)+1,mod(m,8)+1)+1)
m=m+1
endif
11 continue
decchk=(k.eq.0)
do 12 i=0,9
if (ij(k+1,ip(i+1,mod(m,8)+1)+1).eq.0) goto 1
12 continue
1 ch=char(i+48)
return
end
@@ -0,0 +1,18 @@
FUNCTION df1dim(x)
INTEGER NMAX
REAL df1dim,x
PARAMETER (NMAX=50)
CU USES dfunc
INTEGER j,ncom
REAL df(NMAX),pcom(NMAX),xicom(NMAX),xt(NMAX)
COMMON /f1com/ pcom,xicom,ncom
do 11 j=1,ncom
xt(j)=pcom(j)+x*xicom(j)
11 continue
call dfunc(xt,df)
df1dim=0.
do 12 j=1,ncom
df1dim=df1dim+df(j)*xicom(j)
12 continue
return
END
@@ -0,0 +1,52 @@
SUBROUTINE dfour1(data,nn,isign)
INTEGER isign,nn
DOUBLE PRECISION data(2*nn)
INTEGER i,istep,j,m,mmax,n
DOUBLE PRECISION tempi,tempr
DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
n=2*nn
j=1
do 11 i=1,n,2
if(j.gt.i)then
tempr=data(j)
tempi=data(j+1)
data(j)=data(i)
data(j+1)=data(i+1)
data(i)=tempr
data(i+1)=tempi
endif
m=n/2
1 if ((m.ge.2).and.(j.gt.m)) then
j=j-m
m=m/2
goto 1
endif
j=j+m
11 continue
mmax=2
2 if (n.gt.mmax) then
istep=2*mmax
theta=6.28318530717959d0/(isign*mmax)
wpr=-2.d0*sin(0.5d0*theta)**2
wpi=sin(theta)
wr=1.d0
wi=0.d0
do 13 m=1,mmax,2
do 12 i=m,n,istep
j=i+mmax
tempr=wr*data(j)-wi*data(j+1)
tempi=wr*data(j+1)+wi*data(j)
data(j)=data(i)-tempr
data(j+1)=data(i+1)-tempi
data(i)=data(i)+tempr
data(i+1)=data(i+1)+tempi
12 continue
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
13 continue
mmax=istep
goto 2
endif
return
END
@@ -0,0 +1,89 @@
SUBROUTINE dfpmin(p,n,gtol,iter,fret,func,dfunc)
INTEGER iter,n,NMAX,ITMAX
REAL fret,gtol,p(n),func,EPS,STPMX,TOLX
PARAMETER (NMAX=50,ITMAX=200,STPMX=100.,EPS=3.e-8,TOLX=4.*EPS)
EXTERNAL dfunc,func
CU USES dfunc,func,lnsrch
INTEGER i,its,j
LOGICAL check
REAL den,fac,fad,fae,fp,stpmax,sum,sumdg,sumxi,temp,test,dg(NMAX),
*g(NMAX),hdg(NMAX),hessin(NMAX,NMAX),pnew(NMAX),xi(NMAX)
fp=func(p)
call dfunc(p,g)
sum=0.
do 12 i=1,n
do 11 j=1,n
hessin(i,j)=0.
11 continue
hessin(i,i)=1.
xi(i)=-g(i)
sum=sum+p(i)**2
12 continue
stpmax=STPMX*max(sqrt(sum),float(n))
do 27 its=1,ITMAX
iter=its
call lnsrch(n,p,fp,g,xi,pnew,fret,stpmax,check,func)
fp=fret
do 13 i=1,n
xi(i)=pnew(i)-p(i)
p(i)=pnew(i)
13 continue
test=0.
do 14 i=1,n
temp=abs(xi(i))/max(abs(p(i)),1.)
if(temp.gt.test)test=temp
14 continue
if(test.lt.TOLX)return
do 15 i=1,n
dg(i)=g(i)
15 continue
call dfunc(p,g)
test=0.
den=max(fret,1.)
do 16 i=1,n
temp=abs(g(i))*max(abs(p(i)),1.)/den
if(temp.gt.test)test=temp
16 continue
if(test.lt.gtol)return
do 17 i=1,n
dg(i)=g(i)-dg(i)
17 continue
do 19 i=1,n
hdg(i)=0.
do 18 j=1,n
hdg(i)=hdg(i)+hessin(i,j)*dg(j)
18 continue
19 continue
fac=0.
fae=0.
sumdg=0.
sumxi=0.
do 21 i=1,n
fac=fac+dg(i)*xi(i)
fae=fae+dg(i)*hdg(i)
sumdg=sumdg+dg(i)**2
sumxi=sumxi+xi(i)**2
21 continue
if(fac**2.gt.EPS*sumdg*sumxi)then
fac=1./fac
fad=1./fae
do 22 i=1,n
dg(i)=fac*xi(i)-fad*hdg(i)
22 continue
do 24 i=1,n
do 23 j=1,n
hessin(i,j)=hessin(i,j)+fac*xi(i)*xi(j)-fad*hdg(i)*hdg(j)+
*fae*dg(i)*dg(j)
23 continue
24 continue
endif
do 26 i=1,n
xi(i)=0.
do 25 j=1,n
xi(i)=xi(i)-hessin(i,j)*g(j)
25 continue
26 continue
27 continue
pause 'too many iterations in dfpmin'
return
END

Some files were not shown because too many files have changed in this diff Show More