[cig-commits] r8598 - seismo/2D/SPECFEM2D/trunk
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:59:46 PST 2007
Author: walter
Date: 2007-12-07 15:59:45 -0800 (Fri, 07 Dec 2007)
New Revision: 8598
Added:
seismo/2D/SPECFEM2D/trunk/netlib_specfun_erf.f90
Modified:
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
replaced erf() from Numerical Recipes with erf() from Netlib.
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2007-11-20 22:09:36 UTC (rev 8597)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2007-12-07 23:59:45 UTC (rev 8598)
@@ -33,7 +33,7 @@
#LIB = /opt/metis-4.0/gcc64/lib/libmetis.a /opt/scotch-4.0/gcc64/lib/libscotch.a /opt/scotch-4.0/gcc64/lib/libscotcherr.a
LIB =
-OBJS_MESHFEM2D = $O/part_unstruct.o $O/meshfem2D.o $O/read_value_parameters.o
+OBJS_MESHFEM2D = $O/part_unstruct.o $O/meshfem2D.o $O/read_value_parameters.o $O/numerical_recipes.o
OBJS_SPECFEM2D = $O/checkgrid.o $O/datim.o $O/enforce_acoustic_free_surface.o\
$O/compute_forces_acoustic.o $O/compute_forces_elastic.o\
@@ -41,7 +41,7 @@
$O/plotpost.o $O/locate_receivers.o $O/locate_source_force.o $O/compute_gradient_attenuation.o\
$O/specfem2D.o $O/write_seismograms.o $O/define_external_model.o $O/createnum_fast.o $O/createnum_slow.o\
$O/define_shape_functions.o $O/attenuation_model.o $O/create_color_image.o $O/compute_vector_field.o $O/compute_pressure.o\
- $O/recompute_jacobian.o $O/compute_arrays_source.o $O/locate_source_moment_tensor.o $O/numerical_recipes.o\
+ $O/recompute_jacobian.o $O/compute_arrays_source.o $O/locate_source_moment_tensor.o $O/netlib_specfun_erf.o\
$O/construct_acoustic_surface.o $O/assemble_MPI.o $O/compute_energy.o\
$O/attenuation_compute_param.o
@@ -157,6 +157,9 @@
$O/numerical_recipes.o: numerical_recipes.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/numerical_recipes.o numerical_recipes.f90
+$O/netlib_specfun_erf.o: netlib_specfun_erf.f90
+ ${F90} $(FLAGS_CHECK) -c -o $O/netlib_specfun_erf.o netlib_specfun_erf.f90
+
$O/define_external_model.o: define_external_model.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/define_external_model.o define_external_model.f90
@@ -173,4 +176,4 @@
${F90} $(FLAGS_CHECK) -c -o $O/assemble_MPI.o assemble_MPI.F90
$O/attenuation_compute_param.o: attenuation_compute_param.c
- ${CC} -c -o $O/attenuation_compute_param.o attenuation_compute_param.c
\ No newline at end of file
+ ${CC} -c -o $O/attenuation_compute_param.o attenuation_compute_param.c
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-11-20 22:09:36 UTC (rev 8597)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-12-07 23:59:45 UTC (rev 8598)
@@ -1324,81 +1324,3 @@
end function value_spline
-! --------------------------------------
-
-! compute spline coefficients (Numerical Recipes)
-! modified to use dynamic allocation
-
- subroutine spline(x,y,n,yp1,ypn,y2)
-
- implicit none
-
- integer n
- double precision, dimension(n) :: x,y,y2
- double precision, dimension(:), allocatable :: u
- double precision yp1,ypn
-
- integer i,k
- double precision sig,p,qn,un
-
- allocate(u(n))
-
- y2(1)=-0.5d0
- u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
-
- do i=2,n-1
- sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
- p=sig*y2(i-1)+2.d0
- y2(i)=(sig-1.d0)/p
- u(i)=(6.d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
- enddo
-
- qn=0.5d0
- un=(3.d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
- y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)
-
- do k=n-1,1,-1
- y2(k)=y2(k)*y2(k+1)+u(k)
- enddo
-
- deallocate(u)
-
- end subroutine spline
-
-! --------------
-
-! evaluate spline (adapted from Numerical Recipes)
-
- subroutine splint(xa,ya,y2a,n,x,y)
-
- implicit none
-
- integer n
- double precision, dimension(n) :: XA,YA,Y2A
- double precision x,y
-
- integer k,klo,khi
- double precision h,a,b
-
- KLO = 1
- KHI = N
-
- do while (KHI-KLO > 1)
- K=(KHI+KLO)/2
- if(XA(K) > X) then
- KHI=K
- else
- KLO=K
- endif
- enddo
-
- H = XA(KHI) - XA(KLO)
- IF (H == 0.d0) stop 'bad input in spline evaluation'
-
- A = (XA(KHI)-X) / H
- B = (X-XA(KLO)) / H
-
- Y = A*YA(KLO) + B*YA(KHI) + ((A**3-A)*Y2A(KLO) + (B**3-B)*Y2A(KHI))*(H**2)/6.d0
-
- end subroutine splint
-
Added: seismo/2D/SPECFEM2D/trunk/netlib_specfun_erf.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/netlib_specfun_erf.f90 2007-11-20 22:09:36 UTC (rev 8597)
+++ seismo/2D/SPECFEM2D/trunk/netlib_specfun_erf.f90 2007-12-07 23:59:45 UTC (rev 8598)
@@ -0,0 +1,257 @@
+
+ subroutine calerf(ARG,RESULT,JINT)
+
+!------------------------------------------------------------------
+!
+! This routine can be freely obtained from Netlib
+! at http://www.netlib.org/specfun/erf
+!
+! Most Netlib software packages have no restrictions on their use
+! but Netlib recommends that you check with the authors to be sure.
+! See http://www.netlib.org/misc/faq.html#2.3 for details.
+!
+!------------------------------------------------------------------
+!
+! This packet evaluates erf(x) for a real argument x.
+! It contains one FUNCTION type subprogram: ERF,
+! and one SUBROUTINE type subprogram, CALERF. The calling
+! statements for the primary entries are:
+!
+! Y = ERF(X)
+!
+! The routine CALERF is intended for internal packet use only,
+! all computations within the packet being concentrated in this
+! routine. The function subprograms invoke CALERF with the
+! statement
+!
+! call CALERF(ARG,RESULT,JINT)
+!
+! where the parameter usage is as follows
+!
+! Function Parameters for CALERF
+! call ARG Result JINT
+!
+! ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0
+!
+! The main computation evaluates near-minimax approximations
+! from "Rational Chebyshev approximations for the error function"
+! by William J. Cody, Math. Comp., 1969, PP. 631-638. This
+! transportable program uses rational functions that theoretically
+! approximate erf(x) and erfc(x) to at least 18 significant
+! decimal digits. The accuracy achieved depends on the arithmetic
+! system, the compiler, the intrinsic functions, and proper
+! selection of the machine-dependent constants.
+!
+!*******************************************************************
+!*******************************************************************
+!
+! Explanation of machine-dependent constants
+!
+! XMIN = the smallest positive floating-point number.
+! XINF = the largest positive finite floating-point number.
+! XNEG = the largest negative argument acceptable to ERFCX;
+! the negative of the solution to the equation
+! 2*exp(x*x) = XINF.
+! XSMALL = argument below which erf(x) may be represented by
+! 2*x/sqrt(pi) and above which x*x will not underflow.
+! A conservative value is the largest machine number X
+! such that 1.0 + X = 1.0 to machine precision.
+! XBIG = largest argument acceptable to ERFC; solution to
+! the equation: W(x) * (1-0.5/x**2) = XMIN, where
+! W(x) = exp(-x*x)/[x*sqrt(pi)].
+! XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to
+! machine precision. A conservative value is
+! 1/[2*sqrt(XSMALL)]
+! XMAX = largest acceptable argument to ERFCX; the minimum
+! of XINF and 1/[sqrt(pi)*XMIN].
+!
+! Approximate IEEE double precision values are defined below.
+!
+!*******************************************************************
+!*******************************************************************
+!
+! Error returns
+!
+! The program returns ERFC = 0 for ARG >= XBIG;
+!
+! Author: William J. Cody
+! Mathematics and Computer Science Division
+! Argonne National Laboratory
+! Argonne, IL 60439, USA
+!
+! Latest modification: March 19, 1990
+!
+! Converted to Fortran90 and slightly modified by
+! Dimitri Komatitsch, University of Pau, France, November 2007.
+!
+!------------------------------------------------------------------
+
+ implicit none
+
+ integer I,JINT
+ double precision A,ARG,B,C,D,DEL,FOUR,HALF,P,ONE,Q,RESULT,SIXTEEN,SQRPI, &
+ TWO,THRESHOLD,X,XBIG,XDEN,XHUGE,XINF,XMAX,XNEG,XNUM,XSMALL, &
+ Y,YSQ,ZERO
+ dimension A(5),B(4),C(9),D(8),P(6),Q(5)
+
+!------------------------------------------------------------------
+! Mathematical constants
+!------------------------------------------------------------------
+ data FOUR,ONE,HALF,TWO,ZERO/4.0D0,1.0D0,0.5D0,2.0D0,0.0D0/, &
+ SQRPI/5.6418958354775628695D-1/,THRESHOLD/0.46875D0/, &
+ SIXTEEN/16.0D0/
+
+!------------------------------------------------------------------
+! Machine-dependent constants
+!------------------------------------------------------------------
+ data XINF,XNEG,XSMALL/1.79D308,-26.628D0,1.11D-16/, &
+ XBIG,XHUGE,XMAX/26.543D0,6.71D7,2.53D307/
+
+!------------------------------------------------------------------
+! Coefficients for approximation to erf in first interval
+!------------------------------------------------------------------
+ data A/3.16112374387056560D00,1.13864154151050156D02, &
+ 3.77485237685302021D02,3.20937758913846947D03, &
+ 1.85777706184603153D-1/
+ data B/2.36012909523441209D01,2.44024637934444173D02, &
+ 1.28261652607737228D03,2.84423683343917062D03/
+
+!------------------------------------------------------------------
+! Coefficients for approximation to erfc in second interval
+!------------------------------------------------------------------
+ data C/5.64188496988670089D-1,8.88314979438837594D0, &
+ 6.61191906371416295D01,2.98635138197400131D02, &
+ 8.81952221241769090D02,1.71204761263407058D03, &
+ 2.05107837782607147D03,1.23033935479799725D03, &
+ 2.15311535474403846D-8/
+ data D/1.57449261107098347D01,1.17693950891312499D02, &
+ 5.37181101862009858D02,1.62138957456669019D03, &
+ 3.29079923573345963D03,4.36261909014324716D03, &
+ 3.43936767414372164D03,1.23033935480374942D03/
+
+!------------------------------------------------------------------
+! Coefficients for approximation to erfc in third interval
+!------------------------------------------------------------------
+ data P/3.05326634961232344D-1,3.60344899949804439D-1, &
+ 1.25781726111229246D-1,1.60837851487422766D-2, &
+ 6.58749161529837803D-4,1.63153871373020978D-2/
+ data Q/2.56852019228982242D00,1.87295284992346047D00, &
+ 5.27905102951428412D-1,6.05183413124413191D-2, &
+ 2.33520497626869185D-3/
+
+ X = ARG
+ Y = ABS(X)
+ if (Y <= THRESHOLD) then
+
+!------------------------------------------------------------------
+! Evaluate erf for |X| <= 0.46875
+!------------------------------------------------------------------
+ YSQ = ZERO
+ if (Y > XSMALL) YSQ = Y * Y
+ XNUM = A(5)*YSQ
+ XDEN = YSQ
+
+ do I = 1, 3
+ XNUM = (XNUM + A(I)) * YSQ
+ XDEN = (XDEN + B(I)) * YSQ
+ enddo
+
+ RESULT = X * (XNUM + A(4)) / (XDEN + B(4))
+ if (JINT /= 0) RESULT = ONE - RESULT
+ if (JINT == 2) RESULT = EXP(YSQ) * RESULT
+ goto 800
+
+!------------------------------------------------------------------
+! Evaluate erfc for 0.46875 <= |X| <= 4.0
+!------------------------------------------------------------------
+ else if (Y <= FOUR) then
+ XNUM = C(9)*Y
+ XDEN = Y
+
+ do I = 1, 7
+ XNUM = (XNUM + C(I)) * Y
+ XDEN = (XDEN + D(I)) * Y
+ enddo
+
+ RESULT = (XNUM + C(8)) / (XDEN + D(8))
+ if (JINT /= 2) then
+ YSQ = AINT(Y*SIXTEEN)/SIXTEEN
+ DEL = (Y-YSQ)*(Y+YSQ)
+ RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
+ endif
+
+!------------------------------------------------------------------
+! Evaluate erfc for |X| > 4.0
+!------------------------------------------------------------------
+ else
+ RESULT = ZERO
+ if (Y >= XBIG) then
+ if (JINT /= 2 .OR. Y >= XMAX) goto 300
+ if (Y >= XHUGE) then
+ RESULT = SQRPI / Y
+ goto 300
+ endif
+ endif
+ YSQ = ONE / (Y * Y)
+ XNUM = P(6)*YSQ
+ XDEN = YSQ
+
+ do I = 1, 4
+ XNUM = (XNUM + P(I)) * YSQ
+ XDEN = (XDEN + Q(I)) * YSQ
+ enddo
+
+ RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5))
+ RESULT = (SQRPI - RESULT) / Y
+ if (JINT /= 2) then
+ YSQ = AINT(Y*SIXTEEN)/SIXTEEN
+ DEL = (Y-YSQ)*(Y+YSQ)
+ RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
+ endif
+ endif
+
+!------------------------------------------------------------------
+! Fix up for negative argument, erf, etc.
+!------------------------------------------------------------------
+ 300 if (JINT == 0) then
+ RESULT = (HALF - RESULT) + HALF
+ if (X < ZERO) RESULT = -RESULT
+ else if (JINT == 1) then
+ if (X < ZERO) RESULT = TWO - RESULT
+ else
+ if (X < ZERO) then
+ if (X < XNEG) then
+ RESULT = XINF
+ else
+ YSQ = AINT(X*SIXTEEN)/SIXTEEN
+ DEL = (X-YSQ)*(X+YSQ)
+ Y = EXP(YSQ*YSQ) * EXP(DEL)
+ RESULT = (Y+Y) - RESULT
+ endif
+ endif
+ endif
+
+ 800 return
+
+ end subroutine calerf
+
+!--------------------------------------------------------------------
+
+ double precision function netlib_specfun_erf(X)
+
+! This subprogram computes approximate values for erf(x).
+! (see comments heading CALERF).
+!
+! Author/date: William J. Cody, January 8, 1985
+
+ implicit none
+
+ integer JINT
+ double precision X, RESULT
+
+ JINT = 0
+ call calerf(X,RESULT,JINT)
+ netlib_specfun_erf = RESULT
+
+ end function netlib_specfun_erf
+
Modified: seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90 2007-11-20 22:09:36 UTC (rev 8597)
+++ seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90 2007-12-07 23:59:45 UTC (rev 8598)
@@ -6,165 +6,81 @@
!
!=====================================================================
-! double precision routines
+! --------------------------------------
- double precision function erf(x)
+! compute spline coefficients (Numerical Recipes)
+! modified to use dynamic allocation
- implicit none
+ subroutine spline(x,y,n,yp1,ypn,y2)
- double precision x
-
-! this routine uses routine gammp
- double precision gammp
-
- if(x<0.)then
- erf=-gammp(0.5d0,x**2)
- else
- erf=gammp(0.5d0,x**2)
- endif
-
- end function erf
-
-! ---------------------------------
-
- double precision function gammp(a,x)
-
implicit none
- double precision a,x
+ integer n
+ double precision, dimension(n) :: x,y,y2
+ double precision, dimension(:), allocatable :: u
+ double precision yp1,ypn
-! this routine uses routines gcf and gser
- double precision gammcf,gamser,gln
+ integer i,k
+ double precision sig,p,qn,un
- if(x<0.d0 .or. a <= 0.d0) call exit_MPI('bad arguments in gammp')
+ allocate(u(n))
- if(x<a+1.d0)then
- call gser(gamser,a,x,gln)
- gammp=gamser
- else
- call gcf(gammcf,a,x,gln)
- gammp=1.d0-gammcf
- endif
+ y2(1)=-0.5d0
+ u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
- end function gammp
+ do i=2,n-1
+ sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
+ p=sig*y2(i-1)+2.d0
+ y2(i)=(sig-1.d0)/p
+ u(i)=(6.d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
+ enddo
-! ---------------------------------
+ qn=0.5d0
+ un=(3.d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
+ y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)
- subroutine gcf(gammcf,a,x,gln)
-
- implicit none
-
- double precision a,gammcf,gln,x
-
- double precision, parameter :: EPS=3.d-7,FPMIN=1.d-30
- integer, parameter :: ITMAX=100
-
-! this routine uses routine gammln
-
- integer i
- double precision an,b,c,d,del,h
-
- double precision, external :: gammln
-
- gln=gammln(a)
- b=x+1.d0-a
- c=1.d0/FPMIN
- d=1.d0/b
- h=d
- do i=1,ITMAX
- an=-i*(i-a)
- b=b+2.d0
- d=an*d+b
- if(dabs(d)<FPMIN)d=FPMIN
- c=b+an/c
- if(dabs(c)<FPMIN)c=FPMIN
- d=1.d0/d
- del=d*c
- h=h*del
- if(dabs(del-1.d0)<EPS) then
- gammcf=exp(-x+a*log(x)-gln)*h
- return
- endif
+ do k=n-1,1,-1
+ y2(k)=y2(k)*y2(k+1)+u(k)
enddo
- call exit_MPI('a too large, ITMAX too small in gcf')
+ deallocate(u)
- end subroutine gcf
+ end subroutine spline
-! ---------------------------------
+! --------------
- subroutine gser(gamser,a,x,gln)
+! evaluate spline (adapted from Numerical Recipes)
+ subroutine splint(xa,ya,y2a,n,x,y)
+
implicit none
- double precision a,gamser,gln,x
-
- integer, parameter :: ITMAX=100
- double precision, parameter :: EPS=3.d-7
-
-! this routine uses routine gammln
-
integer n
- double precision ap,del,sumval
+ double precision, dimension(n) :: XA,YA,Y2A
+ double precision x,y
- double precision, external :: gammln
+ integer k,klo,khi
+ double precision h,a,b
- gln=gammln(a)
+ KLO = 1
+ KHI = N
- if(x <= 0.d0)then
- if(x<0.d0) call exit_MPI('x < 0 in gser')
- gamser=0.d0
- return
- endif
-
- ap=a
- sumval=1.d0/a
- del=sumval
-
- do n=1,ITMAX
- ap=ap+1.d0
- del=del*x/ap
- sumval=sumval+del
- if(dabs(del)<dabs(sumval)*EPS) then
- gamser=sumval*exp(-x+a*log(x)-gln)
- return
+ do while (KHI-KLO > 1)
+ K=(KHI+KLO)/2
+ if(XA(K) > X) then
+ KHI=K
+ else
+ KLO=K
endif
enddo
- call exit_MPI('a too large, ITMAX too small in gser')
+ H = XA(KHI) - XA(KLO)
+ IF (H == 0.d0) stop 'bad input in spline evaluation'
- end subroutine gser
+ A = (XA(KHI)-X) / H
+ B = (X-XA(KLO)) / H
-! ---------------------------------
+ Y = A*YA(KLO) + B*YA(KHI) + ((A**3-A)*Y2A(KLO) + (B**3-B)*Y2A(KHI))*(H**2)/6.d0
- double precision function gammln(xx)
+ end subroutine splint
- implicit none
-
- double precision xx
-
- integer j
- double precision ser,stp,tmp,x,y,cof(6)
-
- cof(1) = 76.18009172947146d0
- cof(2) = -86.50532032941677d0
- cof(3) = 24.01409824083091d0
- cof(4) = -1.231739572450155d0
- cof(5) = 0.1208650973866179d-2
- cof(6) = -0.5395239384953d-5
-
- stp = 2.5066282746310005d0
-
- x=xx
- y=x
- tmp=x+5.5d0
- tmp=(x+0.5d0)*log(tmp)-tmp
- ser=1.000000000190015d0
- do j=1,6
- y=y+1.d0
- ser=ser+cof(j)/y
- enddo
- gammln=tmp+log(stp*ser/x)
-
- end function gammln
-
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-11-20 22:09:36 UTC (rev 8597)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-07 23:59:45 UTC (rev 8598)
@@ -1471,7 +1471,7 @@
else if(time_function_type == 5) then
hdur = 1.d0 / f0
hdur_gauss = hdur * 5.d0 / 3.d0
- source_time_function(it) = factor * 0.5d0*(1.0d0 + erf(SOURCE_DECAY_MIMIC_TRIANGLE*(time-t0)/hdur_gauss))
+ source_time_function(it) = factor * 0.5d0*(1.0d0 + netlib_specfun_erf(SOURCE_DECAY_MIMIC_TRIANGLE*(time-t0)/hdur_gauss))
else
call exit_MPI('unknown source time function')
More information about the cig-commits
mailing list