C ALGORITHM 824, COLLECTED ALGORITHMS FROM ACM.
C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C VOL. 29,NO. 3, September, 2003, P. 287--296.
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# Doc/
# Doc/framework_abstract
# Doc/index.html
# Doc/simplex_abstract
# Fortran90/
# Fortran90/Drivers/
# Fortran90/Drivers/Makefile
# Fortran90/Drivers/details.f90
# Fortran90/Drivers/details.out
# Fortran90/Drivers/ex_cutet.f90
# Fortran90/Drivers/ex_cutet.out
# Fortran90/Drivers/ex_decuhr2d.f90
# Fortran90/Drivers/ex_decuhr2d.out
# Fortran90/Drivers/ex_decuhr3d.f90
# Fortran90/Drivers/ex_decuhr3d.out
# Fortran90/Drivers/ex_qag.f90
# Fortran90/Drivers/ex_qag.out
# Fortran90/Drivers/ex_qags.f90
# Fortran90/Drivers/ex_qags.out
# Fortran90/Drivers/ex_triex.f90
# Fortran90/Drivers/ex_triex.out
# Fortran90/Drivers/simplexpapertest.f90
# Fortran90/Drivers/simplexpapertest.out
# Fortran90/Src/
# Fortran90/Src/Makefile
# Fortran90/Src/buckley.f90
# Fortran90/Src/check.f90
# Fortran90/Src/cui.f90
# Fortran90/Src/divide.f90
# Fortran90/Src/ds_routines.f90
# Fortran90/Src/error_handling.f90
# Fortran90/Src/global_all.f90
# Fortran90/Src/internal_types.f90
# Fortran90/Src/region_processor.f90
# Fortran90/Src/rule_1.f90
# Fortran90/Src/rule_c2.f90
# Fortran90/Src/rule_c3.f90
# Fortran90/Src/rule_cn.f90
# Fortran90/Src/rule_general.f90
# Fortran90/Src/rule_t2.f90
# Fortran90/Src/rule_t3.f90
# Fortran90/Src/rule_tn.f90
# Fortran90/Src/volume.f90
# This archive created: Fri Oct 24 17:54:24 2003
export PATH; PATH=/bin:$PATH
if test ! -d 'Doc'
then
mkdir 'Doc'
fi
cd 'Doc'
if test -f 'framework_abstract'
then
echo shar: will not over-write existing file "'framework_abstract'"
else
cat << "SHAR_EOF" > 'framework_abstract'
Algorithm 8xx:
CUBPACK: a package for automatic cubature; framework description
Ronald Cools and Ann Haegemans
Abstract:
CUBPACK aims to offer a collection of re-usable code for automatic
n-dimensional (n >= 1) numerical integration of functions over a
collection of regions, i.e., quadrature and cubature. The current
version allows this region to consist of a union of n-simplices and
n-parellellepids. The framework of CUBPACK is described as well as
its user interface. The functionality of several well known routines
is embedded. New features include integration algorithms using the
epsilon-algorithm for extrapolation for regions other than triangles
and the implementation of a new type of subdivision for 3-cubes.
SHAR_EOF
fi # end of overwriting check
if test -f 'index.html'
then
echo shar: will not over-write existing file "'index.html'"
else
cat << "SHAR_EOF" > 'index.html'
CUBPACK documentation page
CUBPACK documentation page
Structure of the files
You will see 3 directories:
- CUBPACK: This contains the source code
- Examples: Some example programs
- Doc: The directory where you found this file
Installation (Unix/Linux)
- Place yourself in the subdirectory CUBPACK and open the
file Makefile
in an editor.
- Assign the name of your Fortran (at least 90) compiler
to the variable FC.
- Assign your favorite compiler options
to the variable CFLAGS.
- Assign the suffix of modules
to the variable MSUFF.
(This is not very important.)
- Execute the command
'make'.
You should not get any error messages...
- Remove redundant files by executing
'make clean'.
Using CUBPACK - try some of the examples first
- Place yourself in the subdirectory Examples and open
the file Makefile in an editor.
- Make the same changes you made to the previous Makefile.
Pay special attention to the variable CFLAGS.
It is here also used to tell the compiler where to look
for CUBPACK modules.
The file you received uses
'-I../CUBPACK'
for this. This works for NAGWare f95.
For SUN's f95 use
'-M../CUBPACK'.
Consult the manual page of your compiler otherwise.
- You will see some example programs mentioned in this Makefile.
Select one and assign its name (no suffixes!) to the
variable MAIN.
- Execute the command
'make'.
You should not get any error messages...
- You will now have an executable with a name equal to what
you selected. Enjoy.
- You can clean up all created files using
'make clean'.
First you might want to know more about the version you installed.
The package contains a subroutine which you can call by
CALL CUBPACK_INFO()
that prints such details. A
sample program and its
output
(obtained on
the platform that produced all sample output) are included.
Here follows a brief description of the examples included.
- 1-dimensional examples
Observe that this uses the reduced calling sequence for scalar
functions and single regions.
- ex_qag:
Simulates QAG from QUADPACK and uses the same test examples.
Look here for the output.
- ex_qags
Simulates QAGS from QUADPACK and uses the same test examples.
Look here for the output.
- 2-dimensional examples
- ex_triex:
Uses the test examples presented in the ACM TOMS paper
describing TRIEX. It uses the same strategy as TRIEX but a
different integration rule and error estimator.
Look here for the output.
- ex_decuhr2d:
Uses a test example presented in the Numerical Algorithms paper
describing DECUHR for the square. It demonstrates the restart
feature with the default integration routine.
Look here for the output.
- 3-dimensional example
- ex_cutet:
Uses one of the test programs included in the distribtion of DCUTET.
It forces the routine to subdivide a tetrahedron in 8 congruent
subregions.
Look here for the output.
- ex_decuhr3d:
Uses the two test examples presented in the Numerical Algorithms paper
describing DECUHR for the 3-dimensional cube.
It demonstrates the use of different integration strategies
(global adaptive with subdivision in 2 (JOB=12) and with a dynamic
subdivision in 2,4 or 8 (JOB=1). This exploits the RESTART feature.
It also uses the extrapolation strategy previously only available
for finite intervals (QAGS) and triangles (TRIEX).
In addition it demonstrates the use of the default (KEY=0) and other
available integration rules. (Clearly showing that the default is
the best available.)
Look here for the output.
- 5-dimensional examples
- simplexpapertest:
This is the test program given in the appendix of the paper
describing a new subdivision strategy for simplices, as part
of CUBPACK.
Look here for the output.
Using CUBPACK - do it yourself
Only one subroutine name is available to the users:
CUBATR.
The most general interface is presented first.
Another interface is provided for those who need to compute an integral
of a scalar function over a single basic region.
A third way to call CUBATR is provided to clear all memory saved by
the previous call.
Index
Full calling sequence
CALL CUBATR &
(DIMENS,NumFun,Integrand,NumRgn,Vertices,RgType,Value,AbsErr, &
! and optional parameters
IFAIL,Neval,EpsAbs,EpsRel,Restart,MinPts,MaxPts,Key,Job,Tune)
!-----------------------------------------------------------------------
! Input parameters
! ----------------
!
! DIMENS Integer.
! The dimension of the region of integration.
!
! NumFun Integer.
! Number of components of the integrand.
!
! Integrand
! Externally declared function for computing all components
! of the integrand at the given evaluation point.
! It must have input parameter X:
! X(1) The x-coordinate of the evaluation point.
! X(2) The y-coordinate of the evaluation point.
! ...
! X(DIMENS) The z-coordinate of the evaluation point.
! and NumFun, the number of components of the integrand.
! It must be compatible with the following interface:
! INTERFACE
! FUNCTION Integrand(NUMFUN,X) RESULT(Value)
! USE Precision_Model
! INTEGER, INTENT(IN) :: NUMFUN
! REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
! REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
! END FUNCTION Integrand
! END INTERFACE
!
! NumRgn Integer.
! The number of given regions.
!
! Vertices
! Real array of dimension (DIMENS,DIMENS+1,NumRgn).
! Vertices(1:DIMENS,K,L) are the x, y, ... coordinates
! respectively of vertex K of region L, where
! K = 1,...,DIMENS+1 and L = 1,...,NumRgn.
!
! RgType Integer array of dimension (NumRgn).
! RgType(L) describes the type of region L.
!
! Value Real array of dimension NumFun.
! Approximations to all components of the integral if
! the procedure is restarted.
!
! AbsErr Real array of dimension NumFun.
! Estimates of absolute errors if the procedure is restarted.
!
! IFAIL Optional integer argument.
! This follows the NAG convention:
! IFAIL = 1 : soft silent error
! Control returned to calling program.
! IFAIL = -1: soft noisy error
! Error message is printed.
! Control returned to calling program.
! IFAIL = 0 : hard noisy error
! Error message is printed and program is stopped.
! Default IFAIL = -1.
!
! EpsAbs Optional real argument.
! Requested absolute error.
! Default EpsAbs = 0.
!
! EpsRel Optional real argument.
! Requested relative error.
! Default EpsRel = sqrt(machine precision).
!
! Restart Optional boolean argument.
! If Restart = FALSE, this is the first attempt to compute
! the integral.
! If Restart = TRUE, then we restart a previous attempt.
! In this case the only parameters for CUBATR that may
! be changed (with respect to the previous call of CUBATR)
! are MinPts, MaxPts, EpsAbs, EpsRel, Key and Restart.
! Default Restart = FALSE.
!
! MinPts Optional integer argument.
! The minimum allowed number of integrand evaluations.
! Default MinPts = 0.
!
! MaxPts Optional integer argument.
! The maximum allowed number of integrand evaluations.
! Default MaxPts = enough to do 500 subdivisions.
!
! Key Optional integer argument.
! Can be used by Rule_General to choose between several
! local integration rules.
! Default Key = 2 if Dimension=1 and extrapolation is used
! (This corresponds to QAGS)
! Default Key = 0 otherwise
!
! Job Optional integer argument.
! If |Job| = 0, then nothing will be done except freeing all
! allocated memory.
! This is usefull after a call of CUBATR if no
! Restart will be done later and memory usage
! might become an issue later.
! Equivalently, one can call CUBATR()
! without any arguments.
! = 1, the global adaptive algorithm is called
! = 2, extrapolation using the epsilon algorithm is used.
! = 11, a region will be divided in 2**DIMENS subregions
! and the global adaptive algorithm is called.
! In combination with Key=0, this resembles DUCTRI and DCUTET.
! = 12, a region will be divided in 2 subregions
! and the global adaptive algorithm is called.
! In combination with Key=3 or 4, this resembles DCUHRE.
! If Job < 0, then an overview of the Region Collection is dumped.
! This will create the files tmp_integerstore and tmp_realstore.
! Default Job = 1.
!
! Tune Optional real argument.
! Can be used by Global_Adapt or the local error estimators
! to influence the reliability. 0 <= Tune <= 1.
! Tune = 1 is the most reliable available.
! Default Tune = 1.
! Note that this is an experimental and controversial parameter.
! In this version, only Tune = 1 is supported for all regions.
!
! Output parameters
! -----------------
!
! Value Real array of dimension NumFun.
! Approximations to all components of the integral
!
! AbsErr Real array of dimension NumFun.
! Estimates of absolute errors.
!
! NEval Optional Integer.
! Number of integrand evaluations used by CUBATR for this call.
!
! IFAIL Optional Integer.
! IFAIL = 0 for normal exit.
!
! AbsErr(K) <= EpsAbs or
! AbsErr(K) <= ABS(Value(K))*EpsRel with MaxPts or less
! function evaluations for all values of K,
! 1 <= K <= NumFun .
!
! IFAIL = 1 if MaxPts was too small to obtain the required
! accuracy. In this case Global_Adapt returns values of
! Value with estimated absolute errors AbsErr.
!
! IFAIL > 1 in more serious case of trouble.
!-----------------------------------------------------------------------
Calling sequence for scalar functions and single regions
CALL CUBATR &
(DIMENS,Integrand,SVertices,SRgType,SValue,SAbsErr, &
! and optional parameters &
IFAIL,Neval,EpsAbs,EpsRel,Restart,MaxPts,Key,Job)
!-----------------------------------------------------------------------
! Input parameters
! ----------------
!
! DIMENS Integer.
! The dimension of the region of integration.
!
! Integrand
! Externally declared function for computing all components
! of the integrand at the given evaluation point.
! It must have input parameter X:
! X(1) The x-coordinate of the evaluation point.
! X(2) The y-coordinate of the evaluation point.
! ...
! X(DIMENS) The z-coordinate of the evaluation point.
! and NumFun, the number of components of the integrand.
! It must be compatible with the following interface:
! INTERFACE
! FUNCTION Integrand(NUMFUN,X) RESULT(Value)
! USE Precision_Model
! INTEGER, INTENT(IN) :: NUMFUN
! REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
! REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
! END FUNCTION Integrand
! END INTERFACE
!
! SVertices
! Real array of dimension (DIMENS,DIMENS+1).
! Vertices(1:DIMENS,K) are the x, y, ... coordinates
! respectively of vertex K of the region, where
! K = 1,...,DIMENS+1.
!
! SRgType Integer.
! RgType describes the type of region L.
!
! SValue Real.
! Approximation to the integral if the procedure is restarted.
!
! SAbsErr Real.
! Estimate of the absolute error if the procedure is restarted.
!
! IFAIL Optional integer argument.
! This follows the NAG convention:
! IFAIL = 1 : soft silent error
! Control returned to calling program.
! IFAIL = -1: soft noisy error
! Error message is printed.
! Control returned to calling program.
! IFAIL = 0 : hard noisy error
! Error message is printed and program is stopped.
! Default IFAIL = -1.
!
! EpsAbs Optional real argument.
! Requested absolute error.
! Default EpsAbs = 0.
!
! EpsRel Optional real argument.
! Requested relative error.
! Default EpsRel = sqrt(machine precision).
!
! Restart Optional boolean argument.
! If Restart = FALSE, this is the first attempt to compute
! the integral.
! If Restart = TRUE, then we restart a previous attempt.
! In this case the only parameters for CUBATR that may
! be changed (with respect to the previous call of CUBATR)
! are MinPts, MaxPts, EpsAbs, EpsRel, Key and Restart.
! Default Restart = FALSE.
!
! MaxPts Optional integer argument.
! The maximum allowed number of integrand evaluations.
! Default MaxPts = enough to do 500 subdivisions.
!
! Key Optional integer argument.
! Can be used by Rule_General to choose between several
! local integration rules.
! Default Key = 2 if Dimension=1 and extrapolation is used
! (This corresponds to QAGS)
! Default Key = 0 otherwise
!
! Job Optional integer argument.
! If |Job| = 0, then nothing will be done except freeing all
! allocated memory.
! This is usefull after a call of CUBATR if no
! Restart will be done later and memory usage
! might become an issue later.
! Equivalently, one can call CUBATR()
! without any arguments.
! = 1, the global adaptive algorithm is called
! = 2, extrapolation using the epsilon algorithm is used.
! = 11, a region will be divided in 2**DIMENS subregions
! and the global adaptive algorithm is called.
! In combination with Key=0, this resembles DUCTRI and DCUTET.
! = 12, a region will be divided in 2 subregions
! and the global adaptive algorithm is called.
! In combination with Key=3 or 4, this resembles DCUHRE.
! If Job < 0, then an overview of the Region Collection is dumped.
! This will create the files tmp_integerstore and tmp_realstore.
! Default Job = 1.
!
! Output parameters
! -----------------
!
! SValue Real.
! Approximation to the integral
!
! AbsErr Real.
! Estimate of the absolute error.
!
! NEval Optional Integer.
! Number of integrand evaluations used by CUBATR for this call.
!
! IFAIL Optional Integer.
! IFAIL = 0 for normal exit.
!
! AbsErr(K) <= EpsAbs or
! AbsErr(K) <= ABS(Value(K))*EpsRel with MaxPts or less
! function evaluations for all values of K,
! 1 <= K <= NumFun .
!
! IFAIL = 1 if MaxPts was too small to obtain the required
! accuracy. In this case Global_Adapt returns values of
! Value with estimated absolute errors AbsErr.
!
! IFAIL > 1 in more serious case of trouble.
!-----------------------------------------------------------------------
Clearing saved memory
Because CUBPACK has a restart feature, information is saved by
the module. To clean this up explicitely, use an additional
CALL CUBATR()
This is equivalent, but much more convenient, than using the general
calling sequence with JOB=0.
This file is maintained by
Ronald Cools.
This page was last modified on Friday 19 July 2002, 10:38:48 CEST.
CUBPACK's home is
here.
SHAR_EOF
fi # end of overwriting check
if test -f 'simplex_abstract'
then
echo shar: will not over-write existing file "'simplex_abstract'"
else
cat << "SHAR_EOF" > 'simplex_abstract'
An Adaptive Numerical Cubature Algorithm for Simplices
Alan Genz and Ronald Cools
Abstract:
A globally adaptive algorithm for numerical cubature of a vector of
functions over a collection of n-dimensional simplices is described.
The algorithm is based on a subdivision strategy that chooses for
subdivision at each stage the subregion (of the input simplices) with
the largest estimated error. This subregion is divided into two, three
or four equal volume subregions by cutting selected edges. These edges
are selected using information about the smoothness of the integrands
in the edge directions. The algorithm allows a choice from several
embedded cubature rule sequences for approximate integration and error
estimation. A Fortran 95 implementation as a part of CUBPACK is also
discussed. Testing of the algorithm is described.
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Fortran90'
then
mkdir 'Fortran90'
fi
cd 'Fortran90'
if test ! -d 'Drivers'
then
mkdir 'Drivers'
fi
cd 'Drivers'
if test -f 'Makefile'
then
echo shar: will not over-write existing file "'Makefile'"
else
cat << "SHAR_EOF" > 'Makefile'
# Suffix for Fortran 90 programs
.SUFFIXES:
.SUFFIXES : .f90 $(SUFFIXES)
#-----------------------------------------------------------------------
# Name of Fortran 90 compiler
FC = f95
# Select your compiler flags
# Example for SUN's f95
#CFLAGS = -M../Src -fast
# Example for NAGWare f95
CFLAGS = -ieee=full -nan -C=all -I../Src
#CFLAGS = -I../Src
# Suffix for Modules (created by compiler)
MSUFF = mod
#-----------------------------------------------------------------------
# Select main program (do not give a suffix!)
# With this distribution the following test programs are distributed:
EXAMPLES= details ex_qag ex_qags ex_triex ex_decuhr2d ex_cutet ex_decuhr3d simplexpapertest
# details : prints info on current distribution
# ex_qag : 1-dimensional integration, tests of QAG from QUADPACK
# ex_qags : 1-dimensional integration, tests of QAGS from QUADPACK
# ex_triex : 2-dimensional integration (triangle), tests of TRIEX
# ex_decuhr2d: 2-dimensional integration (square), test of DECUHR
# ex_cutet : 3-dimensional integration (tetrahedron), test of DCUTET
# ex_decuhr3d: 3-dimensional integration (cube), tests of DECUHR
# simplexpapertest : 5-dimension integration (simplex)
MAIN = simplexpapertest
#-----------------------------------------------------------------------
# The target (executable) will have the same name as the main
# program, without suffix.
TARGET = $(MAIN)
#-----------------------------------------------------------------------
.f90.o:
$(FC) -c $(CFLAGS) $<
#-----------------------------------------------------------------------
all: $(MAIN).o
$(FC) -o $(TARGET) $(MAIN).o $(libcui)
#-----------------------------------------------------------------------
# Place library in current directory (or change the -L../Src)
libcui = -L../Src -lcubpack
clean:
/bin/rm -f *.o *.$(MSUFF) $(EXAMPLES)
SHAR_EOF
fi # end of overwriting check
if test -f 'details.f90'
then
echo shar: will not over-write existing file "'details.f90'"
else
cat << "SHAR_EOF" > 'details.f90'
! This routine only prints details on the current distribtion.
!
PROGRAM Details
USE Precision_Model
USE CUI ! Cubpack User Interface
CALL CUBPACK_INFO()
STOP
END PROGRAM Details
SHAR_EOF
fi # end of overwriting check
if test -f 'details.out'
then
echo shar: will not over-write existing file "'details.out'"
else
cat << "SHAR_EOF" > 'details.out'
---------------------------------------------------------------
CUBPACK information
-------------------
The model for real numbers in the current installed version,
obtained with the declaration REAL(KIND=stnd), has the
following characteristics:
base = 2
digits in this base = 53
This implies:
machine epsilon = 2.2204460492503131E-16
largest real number = 1.7976931348623157E+308
smallest normalized number = 2.2250738585072014E-308
The lowest relative error that may be obtained with this
version is about 0.11E-13
Asking for lower error will push the routine to use the
maximal number of function evaluations it is allowed.
This version accepts a collection of hyper-rectangles
(and parallelepipeds) and simplices as integration regions.
Extrapolation using the epsilon-algorithm is available
for dimensions 1, 2 and 3.
The following values of KEY give different integration rules:
- finite interval: KEY = 1, 2, 3, 4, 5.
KEY < 1 defaults to 1; KEY > 5 defaults to 5.
- n-cube: KEY = 3, 4 uses rule of degree 2*KEY+1
otherwise, uses for a square a rule of degree 13
3-cube a rule of degree 11
a rule of degree 7
- n-simplex: KEY = 1, 2, 3, 4 uses rule of degree 2*KEY+1
otherwise, uses for a triangle a rule of degree 13
tetrahedron a rule of degree 8
a rule of degree 7
KEY = 0 corresponds to our preferred choice.
---------------------------------------------------------------
SHAR_EOF
fi # end of overwriting check
if test -f 'ex_cutet.f90'
then
echo shar: will not over-write existing file "'ex_cutet.f90'"
else
cat << "SHAR_EOF" > 'ex_cutet.f90'
! This file contains the first test program of DCUTET,
! a routine for integration over tetrahedrons.
! See Berntsen, Cools & Espelid, ACM TOMS Vol 19, 1993.
!
MODULE Integrand
IMPLICIT NONE
PUBLIC :: F
CONTAINS
FUNCTION F(NUMFUN,X) RESULT(Value)
USE Precision_Model
INTEGER, INTENT(IN) :: NUMFUN
REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
Value(1) = EXP(X(1)*X(1)+X(2)*X(2)+X(3)*X(3))
RETURN
END FUNCTION F
END MODULE Integrand
PROGRAM Example_CUTET
USE Precision_Model
USE CUI ! Cubpack User Interface
USE Integrand
IMPLICIT NONE
INTEGER, PARAMETER :: n=3, & ! the dimension
m=1, & ! the number of simple regions
l=1, & ! the length of the integrand vector
Tetrahedron = 1, &
Parallelogram = 2
INTEGER, DIMENSION(1:m) :: RgType
INTEGER :: NEval
REAL(kind=stnd), DIMENSION(1:n,0:n,1:m) :: Vertices
REAL(kind=stnd), DIMENSION(1:l) :: IntegralValue, AbsErr
REAL(kind=stnd) :: epsrel
RgType(1) = Tetrahedron
Vertices(1:n,0,1) = (/0 , 0, 0 /)
Vertices(1:n,1,1) = (/1 , 0, 0 /)
Vertices(1:n,2,1) = (/0 , 1, 0 /)
Vertices(1:n,3,1) = (/0 , 0, 1 /)
WRITE(unit=*,fmt=*) "Simulation of DCUTET:"
epsrel = 1.0e-5_stnd
WRITE(unit=*,fmt=*) " CUBATR will now be called with epsrel = ",epsrel
CALL CUBATR(n,l,F,m,Vertices,RgType,IntegralValue,AbsErr,Epsrel=epsrel,NEval=NEval,JOB=11)
WRITE(unit=*,fmt=*) " Integral = ",IntegralValue
WRITE(unit=*,fmt=*) " with estimated error < ",AbsErr
WRITE(unit=*,fmt=*) " The number of integrand evaluations used = ",NEval
epsrel = 1.0e-8_stnd
WRITE(unit=*,fmt=*) " CUBATR will now be called with epsrel = ",epsrel
CALL CUBATR(n,l,F,m,Vertices,RgType,IntegralValue,AbsErr,Epsrel=epsrel,NEval=NEval,Restart=.true.,JOB=11)
WRITE(unit=*,fmt=*) " Integral = ",IntegralValue
WRITE(unit=*,fmt=*) " with estimated error < ",AbsErr
WRITE(unit=*,fmt=*) " The number of additional integrand evaluations used = ",NEval
STOP
END PROGRAM Example_CUTET
SHAR_EOF
fi # end of overwriting check
if test -f 'ex_cutet.out'
then
echo shar: will not over-write existing file "'ex_cutet.out'"
else
cat << "SHAR_EOF" > 'ex_cutet.out'
Simulation of DCUTET:
CUBATR will now be called with epsrel = 1.0000000000000001E-05
Integral = 0.2277999969898576
with estimated error < 1.9752829160297815E-06
The number of integrand evaluations used = 43
CUBATR will now be called with epsrel = 1.0000000000000000E-08
Integral = 0.2277998982512967
with estimated error < 2.1808179040984215E-09
The number of additional integrand evaluations used = 688
SHAR_EOF
fi # end of overwriting check
if test -f 'ex_decuhr2d.f90'
then
echo shar: will not over-write existing file "'ex_decuhr2d.f90'"
else
cat << "SHAR_EOF" > 'ex_decuhr2d.f90'
! This file contains a test example for the square
! that appears in the DECUHR paper:
! Espelid & Genz, Numerical Algorithms Vol 8, 1994.
! The restart feature is used with the default integration routine.
MODULE Integrand
IMPLICIT NONE
PUBLIC :: F
PRIVATE
CONTAINS
FUNCTION F(NUMFUN,X) RESULT(Value)
USE Precision_Model
INTEGER, INTENT(IN) :: NUMFUN
REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
REAL(kind=stnd), DIMENSION(1:NUMFUN) :: Value
! The integrand is scaled so that the exact result is 1.
Value(1) = exp(2*x(1)+x(2)*(1-x(1))) * (1-x(1)) / sqrt(x(1))
Value(1) = Value(1)/3.22289153891637_stnd
RETURN
END FUNCTION F
END MODULE Integrand
PROGRAM Example2D
USE Precision_Model
USE CUI ! Cubpack User Interface
USE Integrand
INTEGER, PARAMETER :: n=2, & ! the dimension
m=1, & ! the number of simple regions
l=1, & ! the length of the integrand vector
Cube=2
INTEGER, DIMENSION(1:m) :: RgType
INTEGER :: NEval, j
REAL(kind=stnd), DIMENSION(1:n,0:n,1:m) :: Vertices
REAL(kind=stnd), DIMENSION(1:l) :: IntegralValue, AbsErr
REAL(kind=stnd) :: EpsRel
LOGICAL :: Restart
RgType(1) = Cube
Vertices(1:n,0,1) = (/0 , 0 /)
Vertices(1:n,1,1) = (/1 , 0 /)
Vertices(1:n,2,1) = (/0 , 1 /)
epsrel = 0.01_stnd
Restart= .false.
do j = 1,4
Print *,"Request relative accuracy = ",epsrel
CALL CUBATR(n,l,F,m,Vertices(:,:,1:m),RgType,IntegralValue,AbsErr, &
NEval=NEval,EpsRel=epsrel,Restart=Restart)
Print *,"-> Integral approximation = ",IntegralValue
Print *,"-> with estimated error < ",AbsErr
Print *,"-> The number of integrand evaluations used = ",NEval
Print *,""
epsrel=epsrel*0.01_stnd
Restart = .true.
end do
STOP
END PROGRAM Example2D
SHAR_EOF
fi # end of overwriting check
if test -f 'ex_decuhr2d.out'
then
echo shar: will not over-write existing file "'ex_decuhr2d.out'"
else
cat << "SHAR_EOF" > 'ex_decuhr2d.out'
Request relative accuracy = 1.0000000000000000E-02
-> Integral approximation = 0.9997541078696911
-> with estimated error < 9.7141400501658589E-03
-> The number of integrand evaluations used = 1073
Request relative accuracy = 1.0000000000000000E-04
-> Integral approximation = 0.9999980789642501
-> with estimated error < 7.5892903270461572E-05
-> The number of integrand evaluations used = 1036
Request relative accuracy = 1.0000000000000002E-06
-> Integral approximation = 0.9999999787718864
-> with estimated error < 8.3967209823977500E-07
-> The number of integrand evaluations used = 962
Request relative accuracy = 1.0000000000000002E-08
-> Integral approximation = 0.9999999998306773
-> with estimated error < 7.7298122812463085E-09
-> The number of integrand evaluations used = 1036
SHAR_EOF
fi # end of overwriting check
if test -f 'ex_decuhr3d.f90'
then
echo shar: will not over-write existing file "'ex_decuhr3d.f90'"
else
cat << "SHAR_EOF" > 'ex_decuhr3d.f90'
! This file contains the test examples for the 3-dimensional cube
! that appear in the DECUHR paper:
! Espelid & Genz, Numerical Algorithms Vol 8, 1994.
!
MODULE Integrand
IMPLICIT NONE
PUBLIC :: F, INIT
PRIVATE
INTEGER, PRIVATE :: Fnr
CONTAINS
FUNCTION F(NUMFUN,X) RESULT(Value)
USE Precision_Model
INTEGER, INTENT(IN) :: NUMFUN
REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
REAL(kind=stnd), DIMENSION(1:NUMFUN) :: Value
REAL(kind=stnd) :: r
SELECT CASE (Fnr)
CASE(1)
Value(1) = exp(x(1)+x(1)*x(2)+x(3)/3) /(sqrt(x(1)+x(2))*2.7878925361_stnd)
CASE(2)
r = sqrt(x(1)*x(1)+x(2)*x(2)+x(3)*x(3))
Value(1) = -log(r)*exp(x(1)*x(2)+x(3))/(sqrt(r)*0.1176364548_stnd)
END SELECT
RETURN
END FUNCTION F
SUBROUTINE INIT (invoer)
INTEGER, INTENT(IN) :: invoer
Fnr = invoer
RETURN
END SUBROUTINE INIT
END MODULE Integrand
PROGRAM Example3D
USE Precision_Model
USE CUI ! Cubpack User Interface
USE Integrand
INTEGER, PARAMETER :: n=3, & ! the dimension
m=1, & ! the number of simple regions
l=1, & ! the length of the integrand vector
Simplex = 1, &
Cube=2
INTEGER, DIMENSION(1:m) :: RgType
INTEGER :: NEval, i,j,k, testfunction, Job, Key
INTEGER, DIMENSION(1:3) :: jobtypes = (/ 1,12,2 /)
INTEGER, DIMENSION(1:3) :: keytypes = (/ 0,3,4 /)
REAL(kind=stnd), DIMENSION(1:n,0:n,1:m) :: Vertices
REAL(kind=stnd), DIMENSION(1:l) :: IntegralValue, AbsErr
REAL(kind=stnd) :: EpsRel
LOGICAL :: RESTART
RgType(1) = Cube
Vertices(1:n,0,1) = (/0 , 0 , 0 /)
Vertices(1:n,1,1) = (/1 , 0 , 0 /)
Vertices(1:n,2,1) = (/0 , 1 , 0 /)
Vertices(1:n,3,1) = (/0 , 0 , 1 /)
do testfunction = 1,2
Print *,"Testfunction ",testfunction
Print *,"==============="
CALL INIT(testfunction)
do i = 1,size(jobtypes)
Job = jobtypes(i)
do k = 1,size(keytypes)
Key = keytypes(k)
Print *,"JOB = ",Job," and Key=",Key
Print *," Err Req Integral Error est Cost"
epsrel = 0.1_stnd
RESTART= .false.
do j = 1,4
CALL CUBATR(n,l,F,m,Vertices(:,:,1:m),RgType,IntegralValue,&
AbsErr,NEval=NEval,EpsRel=epsrel, &
JOB=Job,KEY=Key,Restart=RESTART,MaxPTS=100000)
Print '(E10.3,F22.15,E10.3,I8)',EpsRel, IntegralValue, AbsErr, NEval
epsrel=epsrel*0.01_stnd
if (Job /= 2) then
RESTART = .true.
end if
end do
CALL CUBATR()
Print *,""
end do
Print *,"----------------------------------------------------------"
end do
end do
STOP
END PROGRAM Example3D
SHAR_EOF
fi # end of overwriting check
if test -f 'ex_decuhr3d.out'
then
echo shar: will not over-write existing file "'ex_decuhr3d.out'"
else
cat << "SHAR_EOF" > 'ex_decuhr3d.out'
Testfunction 1
===============
JOB = 1 and Key= 0
Err Req Integral Error est Cost
0.100E+00 1.001163725242669 0.108E-01 89
0.100E-02 1.000051697523749 0.484E-03 1068
0.100E-04 1.000000810182230 0.768E-05 1424
0.100E-06 1.000000005542309 0.979E-07 2492
JOB = 1 and Key= 3
Err Req Integral Error est Cost
0.100E+00 1.000985711929658 0.652E-01 195
0.100E-02 1.000002133902706 0.993E-03 1482
0.100E-04 1.000000022023047 0.922E-05 2028
0.100E-06 0.999999999727435 0.973E-07 5148
JOB = 1 and Key= 4
Err Req Integral Error est Cost
0.100E+00 1.001940586463901 0.196E-01 77
0.100E-02 1.000123252414042 0.656E-03 924
0.100E-04 1.000000279816102 0.892E-05 3542
0.100E-06 1.000000012150892 0.863E-07 5236
----------------------------------------------------------
JOB = 12 and Key= 0
Err Req Integral Error est Cost
0.100E+00 1.001163725242669 0.108E-01 89
0.100E-02 1.000068703503122 0.745E-03 890
0.100E-04 1.000000816891973 0.791E-05 1602
0.100E-06 1.000000002954057 0.941E-07 2848
JOB = 12 and Key= 3
Err Req Integral Error est Cost
0.100E+00 1.000985711929658 0.652E-01 195
0.100E-02 1.000002133902706 0.993E-03 1482
0.100E-04 1.000000022023047 0.922E-05 2028
0.100E-06 0.999999999727435 0.973E-07 5148
JOB = 12 and Key= 4
Err Req Integral Error est Cost
0.100E+00 1.001940586463901 0.196E-01 77
0.100E-02 1.000123252414042 0.656E-03 924
0.100E-04 1.000000279816102 0.892E-05 3542
0.100E-06 1.000000012150892 0.863E-07 5236
----------------------------------------------------------
JOB = 2 and Key= 0
Err Req Integral Error est Cost
0.100E+00 1.001163725242669 0.108E-01 89
0.100E-02 0.999998761741008 0.415E-03 2225
0.100E-04 0.999999837660737 0.117E-05 5073
0.100E-06 0.999999999775124 0.125E-07 34265
JOB = 2 and Key= 3
Err Req Integral Error est Cost
0.100E+00 0.999993132571583 0.101E-02 975
0.100E-02 0.999999250124814 0.137E-04 4719
0.100E-04 0.999999135849675 0.721E-05 6279
0.100E-06 0.999999997856459 0.213E-07 86463
JOB = 2 and Key= 4
Err Req Integral Error est Cost
0.100E+00 1.001940586463901 0.196E-01 77
0.100E-02 0.999999525248721 0.689E-03 3157
0.100E-04 0.999999823055542 0.324E-06 10549
0.100E-06 1.000000002448551 0.807E-07 56133
----------------------------------------------------------
Testfunction 2
===============
JOB = 1 and Key= 0
Err Req Integral Error est Cost
0.100E+00 1.001219527452425 0.500E-01 801
0.100E-02 1.000009844244836 0.392E-03 2136
0.100E-04 1.000000513199047 0.323E-05 2136
0.100E-06 1.000000344478394 0.926E-07 2848
JOB = 1 and Key= 3
Err Req Integral Error est Cost
0.100E+00 0.999939124769939 0.988E-01 5577
0.100E-02 0.999999121249005 0.986E-03 30420
0.100E-04 0.999999961128279 0.998E-05 81042
-> Allowed number of function evaluations reached.
0.100E-06 0.999999991528949 0.569E-06 99996
JOB = 1 and Key= 4
Err Req Integral Error est Cost
0.100E+00 1.013854234672262 0.921E-01 693
0.100E-02 1.000003579641635 0.989E-03 19250
-> Allowed number of function evaluations reached.
0.100E-04 1.000000025953699 0.316E-04 99946
-> Allowed number of function evaluations reached.
0.100E-06 1.000000002122953 0.546E-05 99946
----------------------------------------------------------
JOB = 12 and Key= 0
Err Req Integral Error est Cost
0.100E+00 0.995926356362354 0.756E-01 445
0.100E-02 0.999987084932052 0.879E-03 1602
0.100E-04 0.999999961038715 0.610E-05 3382
0.100E-06 1.000000206555651 0.982E-07 5696
JOB = 12 and Key= 3
Err Req Integral Error est Cost
0.100E+00 0.999939124769939 0.988E-01 5577
0.100E-02 0.999999121249005 0.986E-03 30420
0.100E-04 0.999999961128279 0.998E-05 81042
-> Allowed number of function evaluations reached.
0.100E-06 0.999999991528949 0.569E-06 99996
JOB = 12 and Key= 4
Err Req Integral Error est Cost
0.100E+00 1.013854234672262 0.921E-01 693
0.100E-02 1.000003579641635 0.989E-03 19250
-> Allowed number of function evaluations reached.
0.100E-04 1.000000025953699 0.316E-04 99946
-> Allowed number of function evaluations reached.
0.100E-06 1.000000002122953 0.546E-05 99946
----------------------------------------------------------
JOB = 2 and Key= 0
Err Req Integral Error est Cost
0.100E+00 1.000021050889802 0.120E-02 1513
0.100E-02 1.000001617541351 0.198E-04 2225
0.100E-04 1.000000369271170 0.167E-05 2937
0.100E-06 0.999999996279877 0.545E-07 11481
JOB = 2 and Key= 3
Err Req Integral Error est Cost
0.100E+00 1.000143707726010 0.633E-02 2847
0.100E-02 1.000010592078562 0.134E-03 15015
-> Allowed number of function evaluations reached.
0.100E-04 0.999999640754199 0.118E-04 99879
-> Allowed number of function evaluations reached.
0.100E-06 1.000010575989407 0.133E-03 99879
JOB = 2 and Key= 4
Err Req Integral Error est Cost
0.100E+00 1.006357416426450 0.206E-01 4389
0.100E-02 0.999984668237667 0.125E-03 24101
-> Allowed number of function evaluations reached.
0.100E-04 0.999984729148196 0.321E-04 99869
-> Allowed number of function evaluations reached.
0.100E-06 0.999984729091079 0.317E-04 99869
----------------------------------------------------------
SHAR_EOF
fi # end of overwriting check
if test -f 'ex_qag.f90'
then
echo shar: will not over-write existing file "'ex_qag.f90'"
else
cat << "SHAR_EOF" > 'ex_qag.f90'
MODULE Integrand
PRIVATE
PUBLIC :: F
CONTAINS
FUNCTION F(NUMFUN,X) RESULT(Value)
USE Precision_Model
INTEGER, INTENT(IN) :: NUMFUN
REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
Value(1) = (abs(x(1)-1.0_stnd/3.0_stnd))**(0.8_stnd)
RETURN
END FUNCTION F
END MODULE Integrand
PROGRAM Example_QAG
!
! This is equivalent to the test program for QAG from QUADPACK
! The output is identical to the orginal Fortran 77 programs.
! 25 May 1999
USE Precision_Model
USE CUI ! Cubpack User Interface
USE Integrand
IMPLICIT NONE
INTEGER, PARAMETER :: n=1, & ! the dimension
Finite_interval = 1
INTEGER :: RgType, NEval, j, key
REAL(kind=stnd), DIMENSION(1:n,0:n) :: Vertices
REAL(kind=stnd) :: IntegralValue, AbsErr, epsrel
epsrel = 0.1_stnd
do j=1,10
epsrel = epsrel*0.1_stnd
do key = 1,6
RgType = Finite_interval
Vertices(1,:) = (/ 0 , 1 /)
CALL CUBATR(n,F,Vertices,RgType,IntegralValue,AbsErr, &
MAXPTS=50000,EpsRel=epsrel,Key=key,NEval=NEval,JOB=1)
write(unit=*,fmt="( "" Results for key = "",I3 )") KEY
write(unit=*,fmt="( "" Results for epsrel = "",es9.2 )") EPSREL
write(unit=*,fmt="( "" INTEGRAL APPROXIMATION = "",es15.8 )") IntegralValue
write(unit=*,fmt="( "" ESTIMATE OF ABSOLUTE ERROR = "",es9.2 )") ABSERR
write(unit=*,fmt="( "" NUMBER OF FUNCTION EVALATIONS = "",i5 )") NEVAL
end do
end do
STOP
END PROGRAM Example_QAG
SHAR_EOF
fi # end of overwriting check
if test -f 'ex_qag.out'
then
echo shar: will not over-write existing file "'ex_qag.out'"
else
cat << "SHAR_EOF" > 'ex_qag.out'
Results for key = 1
Results for epsrel = 1.00E-02
INTEGRAL APPROXIMATION = 3.44676993E-01
ESTIMATE OF ABSOLUTE ERROR = 1.08E-03
NUMBER OF FUNCTION EVALATIONS = 135
Results for key = 2
Results for epsrel = 1.00E-02
INTEGRAL APPROXIMATION = 3.44670156E-01
ESTIMATE OF ABSOLUTE ERROR = 1.10E-03
NUMBER OF FUNCTION EVALATIONS = 189
Results for key = 3
Results for epsrel = 1.00E-02
INTEGRAL APPROXIMATION = 3.44675548E-01
ESTIMATE OF ABSOLUTE ERROR = 3.05E-03
NUMBER OF FUNCTION EVALATIONS = 217
Results for key = 4
Results for epsrel = 1.00E-02
INTEGRAL APPROXIMATION = 3.44688417E-01
ESTIMATE OF ABSOLUTE ERROR = 1.56E-03
NUMBER OF FUNCTION EVALATIONS = 205
Results for key = 5
Results for epsrel = 1.00E-02
INTEGRAL APPROXIMATION = 3.44845059E-01
ESTIMATE OF ABSOLUTE ERROR = 6.05E-04
NUMBER OF FUNCTION EVALATIONS = 51
Results for key = 6
Results for epsrel = 1.00E-02
INTEGRAL APPROXIMATION = 3.44700813E-01
ESTIMATE OF ABSOLUTE ERROR = 2.75E-03
NUMBER OF FUNCTION EVALATIONS = 183
Results for key = 1
Results for epsrel = 1.00E-03
INTEGRAL APPROXIMATION = 3.44670425E-01
ESTIMATE OF ABSOLUTE ERROR = 3.10E-04
NUMBER OF FUNCTION EVALATIONS = 165
Results for key = 2
Results for epsrel = 1.00E-03
INTEGRAL APPROXIMATION = 3.44668462E-01
ESTIMATE OF ABSOLUTE ERROR = 3.16E-04
NUMBER OF FUNCTION EVALATIONS = 231
Results for key = 3
Results for epsrel = 1.00E-03
INTEGRAL APPROXIMATION = 3.44668420E-01
ESTIMATE OF ABSOLUTE ERROR = 2.51E-04
NUMBER OF FUNCTION EVALATIONS = 341
Results for key = 4
Results for epsrel = 1.00E-03
INTEGRAL APPROXIMATION = 3.44669481E-01
ESTIMATE OF ABSOLUTE ERROR = 1.29E-04
NUMBER OF FUNCTION EVALATIONS = 369
Results for key = 5
Results for epsrel = 1.00E-03
INTEGRAL APPROXIMATION = 3.44718690E-01
ESTIMATE OF ABSOLUTE ERROR = 1.74E-04
NUMBER OF FUNCTION EVALATIONS = 153
Results for key = 6
Results for epsrel = 1.00E-03
INTEGRAL APPROXIMATION = 3.44670504E-01
ESTIMATE OF ABSOLUTE ERROR = 2.27E-04
NUMBER OF FUNCTION EVALATIONS = 427
Results for key = 1
Results for epsrel = 1.00E-04
INTEGRAL APPROXIMATION = 3.44667997E-01
ESTIMATE OF ABSOLUTE ERROR = 2.56E-05
NUMBER OF FUNCTION EVALATIONS = 225
Results for key = 2
Results for epsrel = 1.00E-04
INTEGRAL APPROXIMATION = 3.44667836E-01
ESTIMATE OF ABSOLUTE ERROR = 2.60E-05
NUMBER OF FUNCTION EVALATIONS = 315
Results for key = 3
Results for epsrel = 1.00E-04
INTEGRAL APPROXIMATION = 3.44667832E-01
ESTIMATE OF ABSOLUTE ERROR = 2.07E-05
NUMBER OF FUNCTION EVALATIONS = 465
Results for key = 4
Results for epsrel = 1.00E-04
INTEGRAL APPROXIMATION = 3.44667920E-01
ESTIMATE OF ABSOLUTE ERROR = 1.06E-05
NUMBER OF FUNCTION EVALATIONS = 533
Results for key = 5
Results for epsrel = 1.00E-04
INTEGRAL APPROXIMATION = 3.44671978E-01
ESTIMATE OF ABSOLUTE ERROR = 1.43E-05
NUMBER OF FUNCTION EVALATIONS = 357
Results for key = 6
Results for epsrel = 1.00E-04
INTEGRAL APPROXIMATION = 3.44668004E-01
ESTIMATE OF ABSOLUTE ERROR = 1.87E-05
NUMBER OF FUNCTION EVALATIONS = 671
Results for key = 1
Results for epsrel = 1.00E-05
INTEGRAL APPROXIMATION = 3.44667797E-01
ESTIMATE OF ABSOLUTE ERROR = 2.11E-06
NUMBER OF FUNCTION EVALATIONS = 285
Results for key = 2
Results for epsrel = 1.00E-05
INTEGRAL APPROXIMATION = 3.44667784E-01
ESTIMATE OF ABSOLUTE ERROR = 2.15E-06
NUMBER OF FUNCTION EVALATIONS = 399
Results for key = 3
Results for epsrel = 1.00E-05
INTEGRAL APPROXIMATION = 3.44667784E-01
ESTIMATE OF ABSOLUTE ERROR = 1.71E-06
NUMBER OF FUNCTION EVALATIONS = 589
Results for key = 4
Results for epsrel = 1.00E-05
INTEGRAL APPROXIMATION = 3.44667820E-01
ESTIMATE OF ABSOLUTE ERROR = 3.06E-06
NUMBER OF FUNCTION EVALATIONS = 615
Results for key = 5
Results for epsrel = 1.00E-05
INTEGRAL APPROXIMATION = 3.44668126E-01
ESTIMATE OF ABSOLUTE ERROR = 1.18E-06
NUMBER OF FUNCTION EVALATIONS = 561
Results for key = 6
Results for epsrel = 1.00E-05
INTEGRAL APPROXIMATION = 3.44667798E-01
ESTIMATE OF ABSOLUTE ERROR = 1.54E-06
NUMBER OF FUNCTION EVALATIONS = 915
Results for key = 1
Results for epsrel = 1.00E-06
INTEGRAL APPROXIMATION = 3.44667781E-01
ESTIMATE OF ABSOLUTE ERROR = 1.74E-07
NUMBER OF FUNCTION EVALATIONS = 345
Results for key = 2
Results for epsrel = 1.00E-06
INTEGRAL APPROXIMATION = 3.44667780E-01
ESTIMATE OF ABSOLUTE ERROR = 1.77E-07
NUMBER OF FUNCTION EVALATIONS = 483
Results for key = 3
Results for epsrel = 1.00E-06
INTEGRAL APPROXIMATION = 3.44667780E-01
ESTIMATE OF ABSOLUTE ERROR = 1.41E-07
NUMBER OF FUNCTION EVALATIONS = 713
Results for key = 4
Results for epsrel = 1.00E-06
INTEGRAL APPROXIMATION = 3.44667783E-01
ESTIMATE OF ABSOLUTE ERROR = 2.52E-07
NUMBER OF FUNCTION EVALATIONS = 779
Results for key = 5
Results for epsrel = 1.00E-06
INTEGRAL APPROXIMATION = 3.44667879E-01
ESTIMATE OF ABSOLUTE ERROR = 3.39E-07
NUMBER OF FUNCTION EVALATIONS = 663
Results for key = 6
Results for epsrel = 1.00E-06
INTEGRAL APPROXIMATION = 3.44667781E-01
ESTIMATE OF ABSOLUTE ERROR = 1.27E-07
NUMBER OF FUNCTION EVALATIONS = 1159
Results for key = 1
Results for epsrel = 1.00E-07
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.44E-08
NUMBER OF FUNCTION EVALATIONS = 405
Results for key = 2
Results for epsrel = 1.00E-07
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.46E-08
NUMBER OF FUNCTION EVALATIONS = 567
Results for key = 3
Results for epsrel = 1.00E-07
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.16E-08
NUMBER OF FUNCTION EVALATIONS = 837
Results for key = 4
Results for epsrel = 1.00E-07
INTEGRAL APPROXIMATION = 3.44667780E-01
ESTIMATE OF ABSOLUTE ERROR = 2.08E-08
NUMBER OF FUNCTION EVALATIONS = 943
Results for key = 5
Results for epsrel = 1.00E-07
INTEGRAL APPROXIMATION = 3.44667787E-01
ESTIMATE OF ABSOLUTE ERROR = 2.80E-08
NUMBER OF FUNCTION EVALATIONS = 867
Results for key = 6
Results for epsrel = 1.00E-07
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.05E-08
NUMBER OF FUNCTION EVALATIONS = 1403
Results for key = 1
Results for epsrel = 1.00E-08
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.20E-09
NUMBER OF FUNCTION EVALATIONS = 465
Results for key = 2
Results for epsrel = 1.00E-08
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.20E-09
NUMBER OF FUNCTION EVALATIONS = 651
Results for key = 3
Results for epsrel = 1.00E-08
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 3.34E-09
NUMBER OF FUNCTION EVALATIONS = 899
Results for key = 4
Results for epsrel = 1.00E-08
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.71E-09
NUMBER OF FUNCTION EVALATIONS = 1107
Results for key = 5
Results for epsrel = 1.00E-08
INTEGRAL APPROXIMATION = 3.44667780E-01
ESTIMATE OF ABSOLUTE ERROR = 2.31E-09
NUMBER OF FUNCTION EVALATIONS = 1071
Results for key = 6
Results for epsrel = 1.00E-08
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 3.02E-09
NUMBER OF FUNCTION EVALATIONS = 1525
Results for key = 1
Results for epsrel = 1.00E-09
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.13E-10
NUMBER OF FUNCTION EVALATIONS = 525
Results for key = 2
Results for epsrel = 1.00E-09
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 9.93E-11
NUMBER OF FUNCTION EVALATIONS = 735
Results for key = 3
Results for epsrel = 1.00E-09
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 2.75E-10
NUMBER OF FUNCTION EVALATIONS = 1023
Results for key = 4
Results for epsrel = 1.00E-09
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.41E-10
NUMBER OF FUNCTION EVALATIONS = 1271
Results for key = 5
Results for epsrel = 1.00E-09
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.90E-10
NUMBER OF FUNCTION EVALATIONS = 1275
Results for key = 6
Results for epsrel = 1.00E-09
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 2.49E-10
NUMBER OF FUNCTION EVALATIONS = 1769
Results for key = 1
Results for epsrel = 1.00E-10
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 2.34E-11
NUMBER OF FUNCTION EVALATIONS = 585
Results for key = 2
Results for epsrel = 1.00E-10
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 2.85E-11
NUMBER OF FUNCTION EVALATIONS = 777
Results for key = 3
Results for epsrel = 1.00E-10
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 2.27E-11
NUMBER OF FUNCTION EVALATIONS = 1147
Results for key = 4
Results for epsrel = 1.00E-10
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.17E-11
NUMBER OF FUNCTION EVALATIONS = 1435
Results for key = 5
Results for epsrel = 1.00E-10
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.57E-11
NUMBER OF FUNCTION EVALATIONS = 1479
Results for key = 6
Results for epsrel = 1.00E-10
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 2.05E-11
NUMBER OF FUNCTION EVALATIONS = 2013
Results for key = 1
Results for epsrel = 1.00E-11
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.93E-12
NUMBER OF FUNCTION EVALATIONS = 705
Results for key = 2
Results for epsrel = 1.00E-11
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 2.35E-12
NUMBER OF FUNCTION EVALATIONS = 861
Results for key = 3
Results for epsrel = 1.00E-11
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.88E-12
NUMBER OF FUNCTION EVALATIONS = 1271
Results for key = 4
Results for epsrel = 1.00E-11
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 3.35E-12
NUMBER OF FUNCTION EVALATIONS = 1517
Results for key = 5
Results for epsrel = 1.00E-11
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.30E-12
NUMBER OF FUNCTION EVALATIONS = 1683
Results for key = 6
Results for epsrel = 1.00E-11
INTEGRAL APPROXIMATION = 3.44667779E-01
ESTIMATE OF ABSOLUTE ERROR = 1.70E-12
NUMBER OF FUNCTION EVALATIONS = 2257
SHAR_EOF
fi # end of overwriting check
if test -f 'ex_qags.f90'
then
echo shar: will not over-write existing file "'ex_qags.f90'"
else
cat << "SHAR_EOF" > 'ex_qags.f90'
MODULE Integrand
PRIVATE
PUBLIC :: F
CONTAINS
FUNCTION F(NUMFUN,X) RESULT(Value)
USE Precision_Model
INTEGER, INTENT(IN) :: NUMFUN
REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
Value(1) = log(x(1))*(x(1)**0.2_stnd)
RETURN
END FUNCTION F
END MODULE Integrand
PROGRAM Example_QAGS
!
! This is equivalent to the test program for QAGS from QUADPACK
! The output is identical to the orginal Fortran 77 programs
! (at least on Linux using NAGWare f95 and f77).
! Minor differences for small numbers possible on e.g. Sun.
USE Precision_Model
USE CUI ! Cubpack User Interface
USE Integrand
INTEGER, PARAMETER :: n=1, & ! the dimension
Finite_interval = 1
INTEGER :: RgType, NEval, j, key
REAL(kind=stnd):: IntegralValue, AbsErr, epsrel
REAL(kind=stnd), DIMENSION(1:n,0:n) :: Vertices
epsrel = 0.1_stnd
do j=1,10
epsrel = epsrel*0.1_stnd
RgType = Finite_interval
Vertices(1,:) = (/ 0 , 1 /)
CALL CUBATR(n,F,Vertices,RgType,IntegralValue,AbsErr, &
EpsRel=epsrel,NEval=NEval,JOB=2)
! EpsRel=epsrel,Key=2,NEval=NEval,JOB=2)
! ^^/^^ ^^^/^
! QAGS uses dqk21 only ---------/ /
! (Key=2 is therefore made default.) /
! uses extrapolation ------------------------/
!
! write(unit=*,fmt="( "" Results for key = "",I3 )") KEY
write(unit=*,fmt="( "" Resuls for epsrel = "",es9.2 )") EPSREL
write(unit=*,fmt="( "" INTEGRAL APPROXIMATION = "",es15.8 )") IntegralValue
write(unit=*,fmt="( "" ESTIMATE OF ABSOLUTE ERROR = "",es9.2 )") ABSERR
write(unit=*,fmt="( "" NUMBER OF FUNCTION EVALATIONS = "",i5 )") NEVAL
write(unit=*,fmt="( "" ERROR CODE = 0"" )") ! IER
end do
STOP
END PROGRAM Example_QAGS
SHAR_EOF
fi # end of overwriting check
if test -f 'ex_qags.out'
then
echo shar: will not over-write existing file "'ex_qags.out'"
else
cat << "SHAR_EOF" > 'ex_qags.out'
Resuls for epsrel = 1.00E-02
INTEGRAL APPROXIMATION = -6.94451546E-01
ESTIMATE OF ABSOLUTE ERROR = 4.30E-03
NUMBER OF FUNCTION EVALATIONS = 189
ERROR CODE = 0
Resuls for epsrel = 1.00E-03
INTEGRAL APPROXIMATION = -6.94444444E-01
ESTIMATE OF ABSOLUTE ERROR = 4.20E-06
NUMBER OF FUNCTION EVALATIONS = 231
ERROR CODE = 0
Resuls for epsrel = 1.00E-04
INTEGRAL APPROXIMATION = -6.94444444E-01
ESTIMATE OF ABSOLUTE ERROR = 4.20E-06
NUMBER OF FUNCTION EVALATIONS = 231
ERROR CODE = 0
Resuls for epsrel = 1.00E-05
INTEGRAL APPROXIMATION = -6.94444444E-01
ESTIMATE OF ABSOLUTE ERROR = 4.20E-06
NUMBER OF FUNCTION EVALATIONS = 231
ERROR CODE = 0
Resuls for epsrel = 1.00E-06
INTEGRAL APPROXIMATION = -6.94444444E-01
ESTIMATE OF ABSOLUTE ERROR = 8.84E-15
NUMBER OF FUNCTION EVALATIONS = 315
ERROR CODE = 0
Resuls for epsrel = 1.00E-07
INTEGRAL APPROXIMATION = -6.94444444E-01
ESTIMATE OF ABSOLUTE ERROR = 8.84E-15
NUMBER OF FUNCTION EVALATIONS = 315
ERROR CODE = 0
Resuls for epsrel = 1.00E-08
INTEGRAL APPROXIMATION = -6.94444444E-01
ESTIMATE OF ABSOLUTE ERROR = 8.84E-15
NUMBER OF FUNCTION EVALATIONS = 315
ERROR CODE = 0
Resuls for epsrel = 1.00E-09
INTEGRAL APPROXIMATION = -6.94444444E-01
ESTIMATE OF ABSOLUTE ERROR = 8.84E-15
NUMBER OF FUNCTION EVALATIONS = 315
ERROR CODE = 0
Resuls for epsrel = 1.00E-10
INTEGRAL APPROXIMATION = -6.94444444E-01
ESTIMATE OF ABSOLUTE ERROR = 8.84E-15
NUMBER OF FUNCTION EVALATIONS = 315
ERROR CODE = 0
Resuls for epsrel = 1.00E-11
INTEGRAL APPROXIMATION = -6.94444444E-01
ESTIMATE OF ABSOLUTE ERROR = 8.84E-15
NUMBER OF FUNCTION EVALATIONS = 315
ERROR CODE = 0
SHAR_EOF
fi # end of overwriting check
if test -f 'ex_triex.f90'
then
echo shar: will not over-write existing file "'ex_triex.f90'"
else
cat << "SHAR_EOF" > 'ex_triex.f90'
!
! This file is equivalent to the tests mentioned in the TRIEX paper.
! The results are NOT fully identical because CUBPACK uses another
! integration rule and error estimator for a triangle.
!
MODULE Integrand
IMPLICIT NONE
PUBLIC :: F, INIT
PRIVATE
INTEGER, PRIVATE :: Fnr
CONTAINS
FUNCTION F(NUMFUN,X) RESULT(Value)
USE Precision_Model
INTEGER, INTENT(IN) :: NUMFUN
REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
INTEGER, DIMENSION(1:10), PARAMETER :: &
P = (/ 1,0,1,1,0,0,1,1,1,1 /), &
Q = (/ 0,0,0,0,0,0,1,1,1,0 /), &
R = (/ 0,1,0,0,1,1,0,1,0,1 /)
REAL(kind=stnd), DIMENSION(1:10), PARAMETER :: &
A=(/0.0_stnd,0.0_stnd,0.5_stnd,1.0/3.0_stnd,0.0_stnd,0.0_stnd,1.0_stnd,0.0_stnd,0.5_stnd,0.0_stnd /), &
B = 0, &
C = (/0.0_stnd,1.0_stnd,0.0_stnd,0.0_stnd,0.5_stnd,1.0/3.0_stnd,0.0_stnd,0.0_stnd,0.0_stnd,1.0/3.0_stnd/)
SELECT CASE (Fnr)
CASE (1:10)
Value(1) = p(Fnr)*sqrt(abs(x(1)-a(Fnr)))+q(Fnr)*sqrt(abs(x(2)-b(Fnr))) &
+ r(Fnr)*sqrt(abs(x(1)+x(2)-c(Fnr)))
CASE(11)
Value(1) = log(x(1)+x(2))
CASE (12)
Value(1) = 1/sqrt(x(1)*x(1)+x(2)*x(2))
CASE (13)
Value(1) = log(sqrt(x(1)*x(1)+x(2)*x(2)))/sqrt(x(1)*x(1)+x(2)*x(2))
CASE (14)
Value(1) = sin(x(1)) * cos(5*x(2))
CASE (15)
! I guess there is a typo in the triex paper and t must be 1.
Value(1) = sin(11*x(1)) * cos(x(2))
CASE (16)
Value(1) = x(1)**(-0.2_stnd)*9.0_stnd/6.25_stnd
CASE (17)
Value(1) = 1.0_stnd/sqrt(x(1)) + 1.0_stnd/sqrt(x(2)) &
+ 1.0_stnd/sqrt(x(1) + x(2))
Value(1) = Value(1)*3.0_stnd/10.0_stnd
CASE (18)
Value(1) = 3.0_stnd/(sqrt(1 - x(1) - x(2))*4.0_stnd)
CASE (19)
Value(1) = (x(1)*x(2))**(-0.2_stnd)/0.9481026454955768_stnd
CASE (20)
Value(1) = -2.0_stnd * log(x(1)*x(2))/ 3.0_stnd
CASE (21)
Value(1) = (1.0_stnd/sqrt(abs(x(1)-0.25_stnd)) + &
1.0_stnd/sqrt(abs(x(2) - 0.5_stnd)))/3.11357229949_stnd
END SELECT
RETURN
END FUNCTION F
SUBROUTINE INIT (invoer)
INTEGER, INTENT(IN) :: invoer
Fnr = invoer
RETURN
END SUBROUTINE INIT
END MODULE Integrand
! ---- Now follows the main program ---
PROGRAM Example_Triex
USE Precision_Model
USE CUI ! Cubpack User Interface
USE Integrand
IMPLICIT NONE
INTEGER, PARAMETER :: n=2, & ! the dimension
m=1, & ! the number of simple regions
l=1, & ! the length of the integrand vector
Triangle = 1, & ! increase readability
Parallelogram = 2
INTEGER, DIMENSION(1:n) :: RgType
INTEGER :: NEval,i
REAL(kind=stnd), DIMENSION(1:n,0:n,1:m) :: Vertices
REAL(kind=stnd), DIMENSION(1:l) :: IntegralValue, AbsErr
REAL(kind=stnd) :: epsrel
RgType(1) = Triangle
Vertices(1:n,0,1) = (/0 , 0 /)
Vertices(1:n,1,1) = (/1 , 0 /)
Vertices(1:n,2,1) = (/0 , 1 /)
epsrel = 1.0e-6_stnd
do i = 1,21
WRITE(unit=*,fmt=*)
WRITE(unit=*,fmt=*) "Testfunction ",i
CALL INIT(i)
WRITE(unit=*,fmt=*) "CUBATR will now be called with epsrel = ",real(epsrel)
CALL CUBATR(n,l,F,m,Vertices(:,:,1:m),RgType,IntegralValue, &
AbsErr,Epsrel=epsrel,NEval=NEval,JOB=2,MaxPts=100000)
WRITE(unit=*,fmt=*) "Integral = ",IntegralValue
WRITE(unit=*,fmt=*) "with estimated error < ",real(AbsErr)
WRITE(unit=*,fmt=*) "The number of integrand evaluations used = ",NEval
CALL CUBATR()
end do
STOP
END PROGRAM Example_Triex
SHAR_EOF
fi # end of overwriting check
if test -f 'ex_triex.out'
then
echo shar: will not over-write existing file "'ex_triex.out'"
else
cat << "SHAR_EOF" > 'ex_triex.out'
Testfunction 1
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.2666666176679433
with estimated error < 2.3370785E-07
The number of integrand evaluations used = 1665
Testfunction 2
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.2666666176679433
with estimated error < 2.3370785E-07
The number of integrand evaluations used = 1665
Testfunction 3
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.2357022998775372
with estimated error < 5.8716449E-08
The number of integrand evaluations used = 4181
Testfunction 4
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.2079633504260309
with estimated error < 4.4174242E-10
The number of integrand evaluations used = 12469
Testfunction 5
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.2357022998775372
with estimated error < 5.8716449E-08
The number of integrand evaluations used = 4181
Testfunction 6
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.2832240788968500
with estimated error < 4.4854476E-10
The number of integrand evaluations used = 12617
Testfunction 7
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.6666666666666838
with estimated error < 8.2636426E-14
The number of integrand evaluations used = 1665
Testfunction 8
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.9333332848090127
with estimated error < 2.2865441E-07
The number of integrand evaluations used = 2553
Testfunction 9
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.5023689976522531
with estimated error < 1.5647734E-07
The number of integrand evaluations used = 6105
Testfunction 10
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.5498908051048528
with estimated error < 6.5724684E-08
The number of integrand evaluations used = 32005
Testfunction 11
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = -0.2499999999997612
with estimated error < 1.9270631E-12
The number of integrand evaluations used = 481
Testfunction 12
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 1.2464504802856151
with estimated error < 8.6946741E-08
The number of integrand evaluations used = 1369
Testfunction 13
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = -1.5280234547382086
with estimated error < 2.9076742E-07
The number of integrand evaluations used = 2405
Testfunction 14
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 4.3052326655855164E-02
with estimated error < 4.1217030E-15
The number of integrand evaluations used = 777
Testfunction 15
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 8.5468091995326595E-02
with estimated error < 1.5776196E-10
The number of integrand evaluations used = 777
Testfunction 16
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.9999999999997983
with estimated error < 1.6637984E-12
The number of integrand evaluations used = 8473
Testfunction 17
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.9999999999994391
with estimated error < 4.7042812E-12
The number of integrand evaluations used = 15577
Testfunction 18
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.9999999999988229
with estimated error < 8.9261949E-12
The number of integrand evaluations used = 8473
Testfunction 19
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.9999997392694022
with estimated error < 4.2480073E-07
The number of integrand evaluations used = 71077
Testfunction 20
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.9999999999997948
with estimated error < 1.3740056E-12
The number of integrand evaluations used = 15577
Testfunction 21
CUBATR will now be called with epsrel = 1.0000000E-06
Integral = 0.9999999999992429
with estimated error < 5.0617345E-12
The number of integrand evaluations used = 18093
SHAR_EOF
fi # end of overwriting check
if test -f 'simplexpapertest.f90'
then
echo shar: will not over-write existing file "'simplexpapertest.f90'"
else
cat << "SHAR_EOF" > 'simplexpapertest.f90'
! This file contains the full example of the paper
! A. Genz & R. Cools
! An adaptive numerical cubature algorithm for simplices
!
MODULE INTEGRAND
USE PRECISION_MODEL
IMPLICIT NONE
PRIVATE
PUBLIC :: F
CONTAINS
FUNCTION F( L, X ) RESULT(FUN)
INTEGER, INTENT(IN) :: L
REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
REAL(kind=stnd), DIMENSION(0:L-1) :: FUN
!
REAL(kind=stnd) :: S
INTEGER :: I
INTEGER, PARAMETER :: N = 5
S = 0
DO I = 1, N
S = S + ( I*X(I) )**2
END DO
FUN(0) = EXP(-S)
FUN(1:N) = X(1:N)*FUN(0)
END FUNCTION F
!
END MODULE INTEGRAND
!
PROGRAM Simplex_Example
USE Precision_Model
USE CUI ! Cubpack User Interface
USE INTEGRAND
INTEGER, PARAMETER :: n = 5, & ! the dimension
nf = n + 1 ! number of integrand functions
INTEGER :: NEval, i, Inform
REAL(kind=stnd), DIMENSION(n,0:n,1) :: Simplex
REAL(kind=stnd), DIMENSION(n,0:n,2) :: TwoSimplices
REAL(kind=stnd), DIMENSION(0:n) :: Value, AbsErr
Inform = -1 ! Soft noisy errors
Simplex = 0 ! Build the unit simplex
DO i = 1, n
Simplex(i,i,1) = 1
END DO
TwoSimplices(:,:,1:1) = Simplex ! Split it into 2 parts.
TwoSimplices(1,0,1) = 0.5_stnd ! Therefore one has to change one
TwoSimplices(:,:,2:2) = Simplex ! coordinate of the unit simplex.
TwoSimplices(1,1,2) = 0.5_stnd
! The two subregions | are simplices | |
! v v v
CALL CUBATR( n, nf, F, 2, TwoSimplices, (/1,1/), Value, AbsErr, &
Inform, NEval, EpsRel=5.0e-4_stnd, MaxPts=200000 )
Print "(""Expected values are "" /5F12.8)", Value(1:n)/Value(0)
Print "(""with estimated errors < ""/5F12.8)", AbsErr(1:n)/Value(0)
Print *,"The number of integrand evaluations used was ", NEval
Inform = -1 ! Soft noisy errors
! The single region | is a simplex |
! v v
CALL CUBATR( n, nf, F, 1, Simplex, (/1/), Value, AbsErr, &
Inform, NEval, EpsRel=5.0e-4_stnd, MaxPts=200000 )
Print "(""Expected values are "" /5F12.8)", Value(1:n)/Value(0)
Print "(""with estimated errors < ""/5F12.8)", AbsErr(1:n)/Value(0)
Print *,"The number of integrand evaluations used was ", NEval
END PROGRAM Simplex_Example
SHAR_EOF
fi # end of overwriting check
if test -f 'simplexpapertest.out'
then
echo shar: will not over-write existing file "'simplexpapertest.out'"
else
cat << "SHAR_EOF" > 'simplexpapertest.out'
Expected values are
0.22419120 0.17813751 0.14013325 0.11372875 0.09523524
with estimated errors <
0.00009069 0.00006152 0.00006977 0.00005160 0.00004529
The number of integrand evaluations used was 178627
Expected values are
0.22419087 0.17813878 0.14013454 0.11372938 0.09523551
with estimated errors <
0.00006262 0.00006472 0.00006016 0.00005668 0.00004474
The number of integrand evaluations used was 76735
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Src'
then
mkdir 'Src'
fi
cd 'Src'
if test -f 'Makefile'
then
echo shar: will not over-write existing file "'Makefile'"
else
cat << "SHAR_EOF" > 'Makefile'
# Suffix for Fortran 90 programs
.SUFFIXES:
.SUFFIXES : .f90 $(SUFFIXES)
#-----------------------------------------------------------------------
# The following are platform dependent
# Name of Fortran 90 (or 95) compiler
FC = f95
# Select your compiler flags
CFLAGS = -ieee=full -nan -C=all -gline
#CFLAGS = -nan -gline
#CFLAGS = -fast
# Suffix for Modules (created by the Fortran compiler)
MSUFF = mod
AR= ar
#-----------------------------------------------------------------------
.f90.o:
$(FC) -c $(CFLAGS) $<
#-----------------------------------------------------------------------
OBJS1= internal_types.o ds_routines.o divide.o \
rule_tn.o rule_t3.o rule_t2.o rule_c2.o rule_c3.o rule_cn.o \
rule_general.o rule_1.o region_processor.o volume.o \
check.o global_all.o error_handling.o
OBJS2=buckley.o cui.o
MOD1= internal_types.$(MSUFF) ds_routines.$(MSUFF) divide.$(MSUFF) \
rule_t3.$(MSUFF) rule_t2.$(MSUFF) rule_c2.$(MSUFF) rule_c3.$(MSUFF) \
rule_general.$(MSUFF) rule_1.$(MSUFF) region_processor.$(MSUFF) volume.$(MSUFF) \
check.$(MSUFF) global_all.$(MSUFF) error_handling.$(MSUFF)
SRC=cui.f90 check.f90 global_all.f90 error_handling.f90 \
buckley.f90 internal_types.f90 ds_routines.f90 divide.f90 \
rule_tn.f90 rule_t3.f90 rule_t2.f90 rule_c2.f90 rule_c3.f90 rule_cn.f90 \
rule_general.f90 rule_1.f90 region_processor.f90 volume.f90
#-----------------------------------------------------------------------
all: libcubpack.a
#-----------------------------------------------------------------------
# dependencies
rule_1.o: buckley.o
rule_c2.o: buckley.o
rule_t2.o: buckley.o
rule_c3.o: buckley.o
rule_t3.o: buckley.o
rule_tn.o: buckley.o
rule_cn.o: buckley.o
volume.o: buckley.o
divide.o: buckley.o internal_types.o
rule_general.o: rule_t3.o rule_t2.o rule_c2.o rule_c3.o rule_1.o \
internal_types.o
region_processor.o : divide.o rule_general.o
global_all.o: region_processor.o volume.o ds_routines.o internal_types.o
cui.o: global_all.o check.o error_handling.o ds_routines.o internal_types.o
internal_types.o: buckley.o
#-----------------------------------------------------------------------
libcubpack.a: $(OBJS1) $(OBJS2)
$(AR) rvu libcubpack.a $(OBJS1) $(OBJS2)
veryclean:
/bin/rm -f $(OBJS1) $(OBJS2) *.$(MSUFF) libcubpack.a
clean:
/bin/rm -f $(OBJS1)
SHAR_EOF
fi # end of overwriting check
if test -f 'buckley.f90'
then
echo shar: will not over-write existing file "'buckley.f90'"
else
cat << "SHAR_EOF" > 'buckley.f90'
! This file is F-compatible, except for upper/lower case conventions.
!--------------------------------------------------------------------
MODULE Extended ! precision specification for real computations
! This requests the processor to use a real implementation 'stnd'
! which provides at least 20 decimal digits of precision and an
! exponent range of at least 10 ^ +- 80. It is expected that this
! precision may not be available on all machines.
! In July 2002, we found this available on
! SUN (Solaris) with f95
! IBM (AIX) with xlf90
! DEC alpha with f90
IMPLICIT NONE
Integer, PUBLIC, PARAMETER :: stnd = Selected_Real_Kind ( 20, 80 )
!-------------------------
! A few computations are preferably done in higher precision 'extd'. The
! numbers chosen here should be such that the underlying hardware will
! select a higher precision for kind 'extd' than for kind 'stnd', if
! this is feasible. If a higher precision is not readily available,
! the same values may be used as are given above for 'stnd'. It is
! anticipated that on many machines, such an even higher precision may
! not be available.
Integer, PUBLIC, PARAMETER :: extd = Selected_Real_Kind ( 30, 80 )
!-------------------------
end Module Extended
MODULE Integers ! precision specification for integer computations
! This is provided for those machines where using short integers
! has some advantage. The range here is from -999 to +999. Note that
! 999 = 10^3 -1. No harm will be done if short integers are made
! the same as long integers.
IMPLICIT NONE
Integer, PUBLIC, PARAMETER :: short = Selected_Int_Kind ( 3 )
!-------------------------
! The range here is at least from -9 999 999 to +9 999 999. Note
! that 10^7 - 1 = 9,999,999. This may limit the largest possible value
! of the dimension n of problems which can be solved. If n is to
! be larger, the 7 should be replaced with a number k so that
! n is considerably less than 10^k.
Integer, PUBLIC, PARAMETER :: long = Selected_Int_Kind ( 7 )
!-------------------------
end Module Integers
MODULE Low ! precision specification for real computations
! This requests the processor to use a real implementation 'stnd'
! which provides at least 6 decimal digits of precision and an
! exponent range of at least 10 ^ +- 35. This would be suitable for
! low accuracy computations. It is expected that this precision will
! be available on all machines.
IMPLICIT NONE
Integer, PUBLIC, PARAMETER :: stnd = Selected_Real_Kind ( 6, 35 )
!-------------------------
! A few computations are preferably done in higher precision 'extd'. The
! numbers chosen here should be such that the underlying hardware will
! select a higher precision for kind 'extd' than for kind 'stnd', if
! this is feasible. If a higher precision is not readily available,
! the same values may be used as are given above for 'stnd'. It is
! anticipated that on most machines this higher precision will also
! be available.
Integer, PUBLIC, PARAMETER :: extd = Selected_Real_Kind ( 12, 35 )
!-------------------------
end Module Low
MODULE Normal ! precision specification for real computations
! This requests the processor to use a real implementation 'stnd'
! which provides at least 12 decimal digits of precision and an
! exponent range of at least 10 ^ +- 50. It is expected that this
! precision will be available on all machines.
IMPLICIT NONE
Integer, PUBLIC, PARAMETER :: stnd = Selected_Real_Kind ( 12, 50 )
!-------------------------
! A few computations are preferably done in higher precision 'extd'. The
! numbers chosen here should be such that the underlying hardware will
! select a higher precision for kind 'extd' than for kind 'stnd', if
! this is feasible. If a higher precision is not readily available,
! the same values may be used as are given above for 'stnd'. It is
! anticipated that on many machines this higher precision may
! not be available.
!Integer, PUBLIC, PARAMETER :: extd = Selected_Real_Kind ( 20, 50 ) ! preferred
Integer, PUBLIC, PARAMETER :: extd = Selected_Real_Kind ( 12, 50 ) ! NAG f95
!-------------------------
end Module Normal
MODULE Precision_Model
! This provides a convenient way of selecting the precision
! required for a computation. By simply ensuring that a leading '!'
! appears on all but exactly one of the following USE statements,
! and then recompiling all routines, the precision of an entire
! computation can be altered.
! USE Low
USE Normal
! USE Extended
USE Integers
! This is the original F90 code
! PRIVATE
! PUBLIC :: stnd, extd, short, long
! This is the F code
PUBLIC
end Module Precision_Model
SHAR_EOF
fi # end of overwriting check
if test -f 'check.f90'
then
echo shar: will not over-write existing file "'check.f90'"
else
cat << "SHAR_EOF" > 'check.f90'
! This file is F-compatible, except for upper/lower case conventions.
!--------------------------------------------------------------------
Module Check_Input
USE Precision_Model, ONLY: stnd
USE internal_types
Implicit NONE
PRIVATE
PUBLIC :: CHECK
PRIVATE :: CHECK_RESTART, CHECK_NORESTART
INTERFACE CHECK
MODULE PROCEDURE CHECK_RESTART, CHECK_NORESTART
END INTERFACE
CONTAINS
SUBROUTINE CHECK_RESTART(DIMENS,NUMFUN,BOTTRH,BOTTIH,&
ISTORE,INFORM,JOB,IFAIL,EPSABS,EPSREL, &
MINPTS,MAXPTS,TUNE)
!***BEGIN PROLOGUE CHECK_RESTART
!***REVISION DATE 950503 (YYMMDD)
!***REVISION DATE 980407 (YYMMDD)
!***REVISION DATE 990531 (YYMMDD)
!***REVISION DATE 990624 (YYMMDD)
!***REVISION DATE 010719 (YYMMDD)
!***AUTHOR
! Ronald Cools, Dept. of Computer Science,
! Katholieke Universiteit Leuven, Celestijnenlaan 200A,
! B-3001 Heverlee, Belgium
! Email: Ronald.Cools@cs.kuleuven.ac.be
!***PURPOSE CHECK_RESTART checks the validity of the
! input parameters to CUBATR on restart.
!***DESCRIPTION
!
!
! ON ENTRY
!
! DIMENS Integer.
! Number of variables.
! NUMFUN Integer.
! Number of components of the integral.
! NUMRGN Integer.
! NUMRGN should contain the initial number of regions.
! MINPTS Integer.
! Minimum number of FUNSUB calls.
! MAXPTS Integer.
! Maximum number of FUNSUB calls.
! TUNE Real.
! Requested reliability.
! JOB Integer.
! Describes what algorithm CUBATR must use.
! EPSABS Real.
! Requested absolute accuracy.
! EPSREL Real.
! Requested relative accuracy.
! BOTTRH Integer.
! Defines the length of the real array containing
! the region collection.
! BOTTIH Integer.
! Defines the length of the integer array containing
! the region collection.
! RGTYPE Integer array of dimension (NUMRGN).
! RGTYPE(L) describes the type of region L.
!
! ON RETURN
!
! IFAIL Integer.
! If IFAIL has an illegal value on entry, it is reset to 0.
!
! INFORM Integer.
! INFORM = 0 for normal exit.
! INFORM = 8 if DIMENS <= 1 ( RESTART = FALSE )
! DIMENS used inconsistently ( RESTAR = TRUE ).
! INFORM = 16 if NUMFUN < 1. ( RESTART = FALSE )
! NUMFUN used inconsistently ( RESTART = TRUE ).
! INFORM = 32 if NUMRGN < 1. ( RESTART = FALSE )
! BOTTIH used inconsistently ( RESTART = TRUE ).
! INFORM = 64 if there is no support for this region (RESTART = FALSE )
! BOTTRH used inconsistently ( RESTART = TRUE ).
! INFORM = 128 if ISTORE is strangely small (RESTART = TRUE)
! INFORM = 256 if MAXPTS <1 or MINPTS > MAXPTS.
! INFORM = 512 if TUNE > 1 OR TUNE < 0.
! INFORM = 1024 if EPSABS < 0 or EPSREL < 0
! INFORM = 2048 if IFAIL not in {-1,0,1}
! INFORM = 4096 if ABS(JOB) not in {0,1,2,11,12}
! If more than one input parameter is wrong,
! then INFORM is set to the sum of the values mentioned above.
! If errors occured during a restart, INFORM = INFORM + 8192.
!
!***END PROLOGUE CHECK_RESTART
!
! Global variables.
!
INTEGER, INTENT(IN) :: DIMENS,NUMFUN,BOTTRH,BOTTIH,JOB
INTEGER, INTENT(OUT) :: INFORM
INTEGER, DIMENSION(:), INTENT(IN) :: ISTORE
INTEGER, OPTIONAL, INTENT(IN OUT) :: IFAIL
INTEGER, OPTIONAL, INTENT(IN) :: MINPTS,MAXPTS
REAL(kind=stnd), OPTIONAL, INTENT(IN) :: EPSABS,EPSREL,TUNE
!
!***FIRST EXECUTABLE STATEMENT CHECK_RESTART
!
INFORM = 0
IF ( SIZE(ISTORE) <= 15) THEN
INFORM = 128
ELSE
! We assume that istore is present and that BOTTIH and
! BORTRH have the right values.
!
! Check valid DIMENS.
!
IF (DIMENS /= ISTORE(1)) THEN
INFORM = 8
END IF
!
! Check positive NUMFUN.
!
IF (NUMFUN /= ISTORE(5)) THEN
INFORM = INFORM + 16
END IF
!
! Check workspace.
!
IF ((BOTTIH /= ISTORE(7)) .OR. (SIZE(ISTORE) /= BOTTIH)) THEN
INFORM = INFORM + 32
END IF
IF (BOTTRH /= ISTORE(8)) THEN
INFORM = INFORM + 64
END IF
END IF
IF (INFORM /= 0) THEN
INFORM = INFORM + 8192
END IF
!
! Check valid limits on allowed number of function evaluations.
!
IF (PRESENT(MAXPTS)) THEN
IF (PRESENT(MINPTS)) THEN
IF (MINPTS > MAXPTS) THEN
INFORM = INFORM + 256
END IF
ELSE
IF (MAXPTS < 1) THEN
INFORM = INFORM + 256
END IF
END IF
END IF
!
! Check valid requested reliablitiy.
!
IF (PRESENT(TUNE)) THEN
IF ((TUNE < 0) .OR. (TUNE > 1)) THEN
INFORM = INFORM + 512
END IF
END IF
!
! Check valid accuracy requests.
!
IF (PRESENT(EPSABS)) THEN
IF (EPSABS < 0) THEN
INFORM = INFORM + 1024
END IF
ELSE IF (PRESENT(EPSREL)) THEN
IF (EPSREL < 0) THEN
INFORM = INFORM + 1024
END IF
END IF
!
! Check valid IFAIL
!
IF (PRESENT(IFAIL)) THEN
IF ((IFAIL < -1) .OR. (IFAIL > 1)) THEN
INFORM = INFORM + 2048
IFAIL = 0
END IF
END IF
!
! Check valid JOB. JOB = 0 cannot occur at this stage
!
IF ((ABS(JOB) /= 1) .AND. (ABS(JOB) /= 2) .AND. (ABS(JOB) /= 11) .AND. (ABS(JOB) /= 12)) THEN
INFORM = INFORM + 4096
END IF
!
RETURN
END SUBROUTINE CHECK_RESTART
SUBROUTINE CHECK_NORESTART(DIMENS,NUMFUN,NUMRGN,RGTYPE, &
INFORM,JOB,IFAIL,EPSABS,EPSREL, &
MINPTS,MAXPTS,TUNE)
!***BEGIN PROLOGUE CHECK_NORESTART
!***REVISION DATE 950503 (YYMMDD)
!***REVISION DATE 980407 (YYMMDD)
!***REVISION DATE 990531 (YYMMDD)
!***AUTHOR
! Ronald Cools, Dept. of Computer Science,
! Katholieke Universiteit Leuven, Celestijnenlaan 200A,
! B-3001 Heverlee, Belgium
! Email: Ronald.Cools@cs.kuleuven.ac.be
!***PURPOSE CHECK_NORESTART checks the validity of the
! input parameters to CUBATR.
!***DESCRIPTION
!
!
! ON ENTRY
!
! DIMENS Integer.
! Number of variables.
! NUMFUN Integer.
! Number of components of the integral.
! NUMRGN Integer.
! NUMRGN should contain the initial number of regions.
! MINPTS Integer.
! Minimum number of FUNSUB calls.
! MAXPTS Integer.
! Maximum number of FUNSUB calls.
! TUNE Real.
! Requested reliability.
! JOB Integer.
! Describes what algorithm CUBATR must use.
! EPSABS Real.
! Requested absolute accuracy.
! EPSREL Real.
! Requested relative accuracy.
! BOTTRH Integer.
! Defines the length of the real array containing
! the region collection.
! BOTTIH Integer.
! Defines the length of the integer array containing
! the region collection.
! RGTYPE Integer array of dimension (NUMRGN).
! RGTYPE(L) describes the type of region L.
!
! ON RETURN
!
! IFAIL Integer.
! If IFAIL has an illegal value on entry, it is reset to 0.
!
! INFORM Integer.
! INFORM = 0 for normal exit.
! INFORM = 8 if DIMENS <= 1 ( RESTART = FALSE )
! DIMENS used inconsistently ( RESTAR = TRUE ).
! INFORM = 16 if NUMFUN < 1. ( RESTART = FALSE )
! NUMFUN used inconsistently ( RESTART = TRUE ).
! INFORM = 32 if NUMRGN < 1. ( RESTART = FALSE )
! BOTTIH used inconsistently ( RESTART = TRUE ).
! INFORM = 64 if there is no support for this region (RESTART = FALSE )
! BOTTRH used inconsistently ( RESTART = TRUE ).
! INFORM = 128 if JOB={2,11} and DIMENS > 3
! INFORM = 256 if MAXPTS <1 or MINPTS > MAXPTS.
! INFORM = 512 if TUNE > 1 OR TUNE < 0.
! INFORM = 1024 if EPSABS < 0 or EPSREL < 0
! INFORM = 2048 if IFAIL not in {-1,0,1}
! INFORM = 4096 if ABS(JOB) not in {0,1,2,11,12}
! If more than one input parameter is wrong,
! then INFORM is set to the sum of the values mentioned above.
! If errors occured during a restart, INFORM = INFORM + 8192.
!
!***END PROLOGUE CHECK_NORESTART
!
! Global variables.
!
INTEGER, INTENT(IN) :: DIMENS,NUMFUN,NUMRGN,JOB
INTEGER, INTENT(OUT) :: INFORM
INTEGER, DIMENSION(:), INTENT(IN) :: RGTYPE
INTEGER, OPTIONAL, INTENT(IN OUT) :: IFAIL
INTEGER, OPTIONAL, INTENT(IN) :: MINPTS,MAXPTS
REAL(kind=stnd), OPTIONAL, INTENT(IN) :: EPSABS,EPSREL,TUNE
!
! Local variables
!
INTEGER :: I
!
!***FIRST EXECUTABLE STATEMENT CHECK_NORESTART
!
INFORM = 0
!
! Check valid DIMENS.
!
IF (DIMENS <= 0) THEN
INFORM = 8
END IF
!
! Check positive NUMFUN.
!
IF (NUMFUN < 1) THEN
INFORM = INFORM + 16
END IF
!
! Check valid NUMRGN.
!
IF (NUMRGN <= 0) THEN
INFORM = INFORM + 32
END IF
!
! Check valid region type
!
DO I = 1,NUMRGN
IF ((RGTYPE(I) /= Simplex) .and. (RGTYPE(I) /= Hyperrectangle)) THEN
INFORM = INFORM + 64
EXIT
END IF
END DO
!
! Check if DIMENS and JOB or in agreement for special values of JOB
!
IF ( ((JOB == 2) .OR. (JOB == 11)) .AND. (DIMENS > 3)) THEN
INFORM = INFORM + 128
END IF
!
! Check valid limits on allowed number of function evaluations.
!
IF (PRESENT(MAXPTS)) THEN
IF (PRESENT(MINPTS)) THEN
IF (MINPTS > MAXPTS) THEN
INFORM = INFORM + 256
END IF
ELSE
IF (MAXPTS < 1) THEN
INFORM = INFORM + 256
END IF
END IF
END IF
!
! Check valid requested reliablitiy.
!
IF (PRESENT(TUNE)) THEN
IF ((TUNE < 0) .OR. (TUNE > 1)) THEN
INFORM = INFORM + 512
END IF
END IF
!
! Check valid accuracy requests.
!
IF (PRESENT(EPSABS)) THEN
IF (EPSABS < 0) THEN
INFORM = INFORM + 1024
END IF
ELSE IF (PRESENT(EPSREL)) THEN
IF (EPSREL < 0) THEN
INFORM = INFORM + 1024
END IF
END IF
!
! Check valid IFAIL
!
IF (PRESENT(IFAIL)) THEN
IF ((IFAIL < -1) .OR. (IFAIL > 1)) THEN
INFORM = INFORM + 2048
IFAIL = 0
END IF
END IF
!
! Check valid JOB. JOB = 0 cannot occur at this stage
!
IF ((ABS(JOB) /= 1) .AND. (ABS(JOB) /= 2) .AND. (ABS(JOB) /= 11) .AND. (ABS(JOB) /= 12)) THEN
INFORM = INFORM + 4096
END IF
!
RETURN
END SUBROUTINE CHECK_NORESTART
END MODULE Check_Input
SHAR_EOF
fi # end of overwriting check
if test -f 'cui.f90'
then
echo shar: will not over-write existing file "'cui.f90'"
else
cat << "SHAR_EOF" > 'cui.f90'
!------------------------!
! Cubpack User Interface !
!------------------------!
Module CUI
USE Precision_Model
USE internal_types
Implicit NONE
PRIVATE
PUBLIC :: CUBATR, CUBPACK_INFO
!-----------------------------------------------------------------------
!***BEGIN PROLOGUE CUBATR
!***DATE WRITTEN 901114 (YYMMDD)
!***REVISION DATE 970620 (YYMMDD)
!***REVISION DATE 980406 (YYMMDD) (MDIV removed)
!***REVISION DATE 000809 (YYMMDD)
!***REVISION DATE 010719 (YYMMDD)
!***REVISION DATE 020715 (YYMMDD) (CUBPACK_INFO added)
!***AUTHOR
! Ronald Cools, Dept. of Computer Science,
! Katholieke Universiteit Leuven, Celestijnenlaan 200A,
! B-3001 Heverlee, Belgium
! Email: Ronald.Cools@cs.kuleuven.ac.be
!
!***PURPOSE Computation of integrals over a collection of regions.
!
!***DESCRIPTION
! CUBATR is the driver routine for CUBPACK and the only
! routine that a user has to deal with (at the moment).
!
!-----------------------------------------------------------------------
PRIVATE :: CUBATR_X, CUBATR_1, CUBATR_CLEAR
!
! Module variables
!
INTEGER, SAVE, PRIVATE :: PreJob=0
INTEGER, PRIVATE :: BOTTIH,BOTTRH
INTEGER, DIMENSION(:), PRIVATE, ALLOCATABLE :: IWork
REAL(kind=stnd), DIMENSION(:), PRIVATE, ALLOCATABLE :: RWork
TYPE(EPSALG_MEM), PRIVATE :: M
INTERFACE CUBATR
MODULE PROCEDURE CUBATR_X, CUBATR_1, CUBATR_CLEAR
END INTERFACE
CONTAINS
SUBROUTINE CUBATR_X &
(DIMENS,NumFun,Integrand,NumRgn,Vertices,RgType,Value,AbsErr, &
! and optional parameters
IFAIL,Neval,EpsAbs,EpsRel,Restart,MinPts,MaxPts,Key,Job,Tune)
!-----------------------------------------------------------------------
! Input parameters
! ----------------
!
! DIMENS Integer.
! The dimension of the region of integration.
!
! NumFun Integer.
! Number of components of the integrand.
!
! Integrand
! Externally declared function for computing all components
! of the integrand at the given evaluation point.
! It must have input parameter X:
! X(1) The x-coordinate of the evaluation point.
! X(2) The y-coordinate of the evaluation point.
! ...
! X(DIMENS) The z-coordinate of the evaluation point.
! and NumFun, the number of components of the integrand.
! It must be compatible with the following interface:
! INTERFACE
! FUNCTION Integrand(NUMFUN,X) RESULT(Value)
! USE Precision_Model
! INTEGER, INTENT(IN) :: NUMFUN
! REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
! REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
! END FUNCTION Integrand
! END INTERFACE
!
! NumRgn Integer.
! The number of given regions.
!
! Vertices
! Real array of dimension (DIMENS,DIMENS+1,NumRgn).
! Vertices(1:DIMENS,K,L) are the x, y, ... coordinates
! respectively of vertex K of region L, where
! K = 1,...,DIMENS+1 and L = 1,...,NumRgn.
!
! RgType Integer array of dimension (NumRgn).
! RgType(L) describes the type of region L.
!
! Value Real array of dimension NumFun.
! Approximations to all components of the integral if
! the procedure is restarted.
!
! AbsErr Real array of dimension NumFun.
! Estimates of absolute errors if the procedure is restarted.
!
! IFAIL Optional integer argument.
! This follows the NAG convention:
! IFAIL = 1 : soft silent error
! Control returned to calling program.
! IFAIL = -1: soft noisy error
! Error message is printed.
! Control returned to calling program.
! IFAIL = 0 : hard noisy error
! Error message is printed and program is stopped.
! Default IFAIL = -1.
!
! EpsAbs Optional real argument.
! Requested absolute error.
! Default EpsAbs = 0.
!
! EpsRel Optional real argument.
! Requested relative error.
! Default EpsRel = sqrt(machine precision).
!
! Restart Optional boolean argument.
! If Restart = FALSE, this is the first attempt to compute
! the integral.
! If Restart = TRUE, then we restart a previous attempt.
! In this case the only parameters for CUBATR that may
! be changed (with respect to the previous call of CUBATR)
! are MinPts, MaxPts, EpsAbs, EpsRel, Key and Restart.
! Default Restart = FALSE.
!
! MinPts Optional integer argument.
! The minimum allowed number of integrand evaluations.
! Default MinPts = 0.
!
! MaxPts Optional integer argument.
! The maximum allowed number of integrand evaluations.
! Default MaxPts = enough to do 500 subdivisions.
!
! Key Optional integer argument.
! Can be used by Rule_General to choose between several
! local integration rules.
! Default Key = 2 if Dimension=1 and extrapolation is used
! (This corresponds to QAGS)
! Default Key = 0 otherwise
!
! Job Optional integer argument.
! If |Job| = 0, then nothing will be done except freeing all
! allocated memory.
! This is usefull after a call of CUBATR if no
! Restart will be done later and memory usage
! might become an issue later.
! Equivalently, one can call CUBATR()
! without any arguments.
! = 1, the global adaptive algorithm is called
! = 2, extrapolation using the epsilon algorithm is used.
! = 11, a region will be divided in 2**DIMENS subregions
! and the global adaptive algorithm is called.
! In combination with Key=0, this resembles DUCTRI and DCUTET.
! = 12, a region will be divided in 2 subregions
! and the global adaptive algorithm is called.
! In combination with Key=3 or 4, this resembles DCUHRE.
! If Job < 0, then an overview of the Region Collection is dumped.
! This will create the files tmp_integerstore and tmp_realstore.
! Default Job = 1.
!
! Tune Optional real argument.
! Can be used by Global_Adapt or the local error estimators
! to influence the reliability. 0 <= Tune <= 1.
! Tune = 1 is the most reliable available.
! Default Tune = 1.
! Note that this is an experimental and controversial parameter.
! In this version, only Tune = 1 is supported for all regions.
!
! Output parameters
! -----------------
!
! Value Real array of dimension NumFun.
! Approximations to all components of the integral
!
! AbsErr Real array of dimension NumFun.
! Estimates of absolute errors.
!
! NEval Optional Integer.
! Number of integrand evaluations used by CUBATR for this call.
!
! IFAIL Optional Integer.
! IFAIL = 0 for normal exit.
!
! AbsErr(K) <= EpsAbs or
! AbsErr(K) <= ABS(Value(K))*EpsRel with MaxPts or less
! function evaluations for all values of K,
! 1 <= K <= NumFun .
!
! IFAIL = 1 if MaxPts was too small to obtain the required
! accuracy. In this case Global_Adapt returns values of
! Value with estimated absolute errors AbsErr.
!
! IFAIL > 1 in more serious case of trouble.
!-----------------------------------------------------------------------
! MODULES USED
USE Check_Input
USE Error_Handling
USE DS_ROUTINES, ONLY: DSCOPY, DSINIT, DSUSED, DSSTAT, DSSUM, DSPINT
USE CubatureRule_General, ONLY: Rule_Cost
USE Global_Adaptive_Algorithm
!***END PROLOGUE CUBATR
!-----------------------------------------------------------------------
!
! Global variables
!
INTERFACE
FUNCTION Integrand(NUMFUN,X) RESULT(Value)
USE Precision_Model
INTEGER, INTENT(IN) :: NUMFUN
REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
END FUNCTION Integrand
END INTERFACE
INTEGER, INTENT(IN) :: DIMENS,NumFun,NumRgn
INTEGER, DIMENSION(:), INTENT(IN) :: RgType
LOGICAL, INTENT(IN), OPTIONAL :: Restart
INTEGER, INTENT(OUT), OPTIONAL :: NEval
INTEGER, INTENT(IN), OPTIONAL :: Job,Key,MaxPts,MinPts
INTEGER, INTENT(IN OUT), OPTIONAL :: IFAIL
REAL(kind=stnd), INTENT(IN), OPTIONAL :: Tune,EpsAbs,EpsRel
REAL(kind=stnd), INTENT(IN), DIMENSION(:,:,:) :: Vertices
REAL(kind=stnd), INTENT(IN OUT), DIMENSION(:) :: AbsErr, Value
!
! Named constants
!
INTEGER, PARAMETER :: NRINFO=1, NIINFO=5
!
! Local variables
!
INTEGER :: BLOCK,i,Inform,Leval,LJob,LMaxPts,LMinPts,MinCost, &
NRVERT,NrSub,NRVACA,MAXRGN,RULCLS,STATUS,Tmp
LOGICAL :: EpsAlg,LRestart
REAL(kind=stnd) :: LEpsAbs, LEpsRel
REAL(kind=stnd), DIMENSION(:), ALLOCATABLE :: TmpRWork
INTEGER, DIMENSION(:), ALLOCATABLE :: TmpIWork
TYPE(INTEGRATOR_INFO) :: CINFO
TYPE(USER_INFO) :: UINFO
!-----------------------------------------------------------------------
!
! Check array sizes
! Array size mismatch results in hard error.
Inform = 0
IF (size(rgtype) < numrgn) THEN
write(unit=*,fmt=*) "Error: size(rgtype) < numrgn"
Inform = Inform + 1
END IF
IF (size(abserr) < numfun) THEN
write(unit=*,fmt=*) "Error: size(abserr) < numfun"
Inform = Inform + 1
END IF
IF (size(Value) < numfun) THEN
write(unit=*,fmt=*) "Error: size(Value) < numfun"
Inform = Inform + 1
END IF
IF ((size(vertices,1) /= dimens) .or. (size(vertices,2) /= dimens+1) &
.or. (size(vertices,3) < numrgn)) THEN
Inform = Inform + 1
write(unit=*,fmt=*)"Error: size(vertices) /= (/dimens,dimens+1,numrgn/)"
END IF
IF (Inform /= 0) THEN
WRITE(unit=*,fmt=*) "Array size mismatch results in hard error."
STOP ! "Array size mismatch results in hard error."
END IF
!-----------------------------------------------------------------------
IF (PRESENT(NEval)) THEN
NEval = 0
END IF
!-----------------------------------------------------------------------
IF (PRESENT(Job)) THEN
LJob = Job
IF (Job == 0) THEN
CALL CUBATR_CLEAR()
RETURN
END IF
ELSE
LJob = 1
END IF
!-----------------------------------------------------------------------
!
! Set optional arguments
!
IF ( PRESENT(Restart)) THEN
LRestart = Restart
ELSE
LRestart = .FALSE.
END IF
IF ( PRESENT(Key)) THEN
CINFO%Key = Key
ELSE
IF ((ABS(LJob) == 2) .AND. (DIMENS == 1)) THEN
CINFO%Key = 2 ! simulate QAGS
ELSE
CINFO%Key = 0
END IF
END IF
!-----------------------------------------------------------------------
!
! Check input parameters
!
IF ( .NOT. LRestart) THEN
CALL CHECK(DIMENS,NumFun,NumRgn,RgType,Inform, &
LJob,IFAIL,EpsAbs,EpsRel,MinPts,MaxPts,Tune)
ELSE IF ( ALLOCATED(IWork) ) THEN
CALL CHECK(DIMENS,NumFun,BOTTRH,BOTTIH,IWork, &
Inform,LJob,IFAIL,EpsAbs,EpsRel,MinPts,MaxPts,Tune)
ELSE
Inform = 4096 ! There is nothing to restart from
END IF
IF (Inform /= 0) THEN
CALL Handle_Error(Inform,IFAIL)
RETURN
END IF
!-----------------------------------------------------------------------
RULCLS = Rule_Cost( DIMENS, RgType(1), CINFO%Key )
MinCost = RULCLS
DO i = 2,NumRgn
Tmp = Rule_Cost( DIMENS, RgType(i), CINFO%Key)
RULCLS = max(RULCLS,Tmp)
MinCost = MinCost + Tmp
END DO
!-----------------------------------------------------------------------
!
! Set optional arguments
!
IF ( PRESENT(MinPts)) THEN
LMinPts = MinPts
ELSE
LMinPts = 0
END IF
IF ( PRESENT(MaxPts)) THEN
LMaxPts = MaxPts
ELSE
LMaxPts = 500*RULCLS
END IF
IF ( PRESENT(Tune)) THEN
CINFO%Tune = Tune
ELSE
CINFO%Tune = 1
END IF
IF ( PRESENT(EpsAbs)) THEN
LEpsAbs = EpsAbs
ELSE
LEpsAbs = 0
END IF
IF ( PRESENT(EpsRel)) THEN
LEpsRel = EpsRel
ELSE
LEpsRel = SQRT(EPSILON(LEpsRel))
END IF
!-----------------------------------------------------------------------
!
! Set other parameters of the Global Adaptive algorithm
!
!
! NrSub is an upper limit for the number of subregions after subdivision.
! This influence memory managment, so don't exagerate here.
IF (DIMENS <= 3) THEN
NrSub = 2**DIMENS
ELSE
NrSub = 4
END IF
!
EpsAlg = ( ABS(LJob) == 2 )
CINFO%UNIFORM_SUBDIV = EpsAlg
CINFO%NrSub = NrSub
IF (( ABS(LJob) == 11) .AND. (DIMENS <= 3)) THEN
! simulate dcutri and dcutet ; NrSub = 2**DIMENS ; EpsAlg = .FALSE.
CINFO%UNIFORM_SUBDIV = .TRUE.
END IF
IF ( ABS(LJob) == 12 ) THEN
! simulate dcuhre; NrSub = 2 ; EpsAlg = .FALSE.
CINFO%NrSub = 2
END IF
NRVERT = DIMENS + 1 ! Only cubes and simplices are implemented here.
UINFO%NumFun = NumFun
UINFO%NumRgn = NumRgn
UINFO%MinPts = LMinPts
UINFO%MaxPts = LMaxPts
UINFO%EpsAbs = LEpsAbs
UINFO%EpsRel = LEpsRel
!-----------------------------------------------------------------------
IF (LRestart) THEN
! This requires allocating larger arrays and copying
! the region collection.
ALLOCATE(TmpRWork(SIZE(RWork)),STAT=status)
IF (status /= 0) THEN
WRITE(unit=*,fmt=*) "Problem allocating real workspace."
STOP ! "Problem allocating real workspace."
END IF
ALLOCATE(TmpIWork(SIZE(IWork)),STAT=status)
IF (status /= 0) THEN
WRITE(unit=*,fmt=*) "Problem allocating integer workspace."
STOP ! "Problem allocating integer workspace."
END IF
MAXRGN = DSUSED(IWork)
CALL DSCOPY(IWork,RWork,TmpIWork,TmpRWork)
ELSE
MAXRGN = NumRgn
END IF
!-----------------------------------------------------------------------
! NRVACA is the number of regions the global adaptive algorithm
! removes from the data structure for further processing.
! In some routines for shared memory parallel machines
! this is the variable MDIV
NRVACA = 1
! MAXRGN depends on the number of function evalutions
MAXRGN = MAXRGN + 1 + (NrSub-1)*(LMaxPts - RULCLS*NumRgn)/(RULCLS*NrSub)
!
! Compute length of workspace needed.
!
BOTTIH = MAXRGN*(1+NIINFO) + 15 + NRVACA
BLOCK = NRINFO+NRVERT*DIMENS+2*NumFun
IF (NumFun > 1) THEN
BLOCK = BLOCK + 1
END IF
BOTTRH = MAXRGN*BLOCK
!
! Allocate space for the region collection
!
IF (ALLOCATED(RWork)) THEN
DEALLOCATE(RWork)
END IF
ALLOCATE(RWork(BOTTRH),STAT=status)
IF (status /= 0) THEN
WRITE(unit=*,fmt=*) "Problem allocating real workspace."
STOP ! "Problem allocating real workspace."
END IF
IF (ALLOCATED(IWork)) THEN
DEALLOCATE(IWork)
END IF
ALLOCATE(IWork(BOTTIH),STAT=status)
IF (status /= 0) THEN
WRITE(unit=*,fmt=*) "Problem allocating integer workspace."
STOP ! "Problem allocating integer workspace."
END IF
!
! Initialise region collection
!
IF ( LRestart) THEN
CALL DSCOPY(TmpIWork,TmpRWork,IWork,RWork)
DEALLOCATE(TmpIWork,TmpRWork)
ELSE
IF ( MinCost > LMaxPts ) THEN
Inform = 128 ! Dit nummer werd al gebruikt !
ELSE
CALL DSINIT(DIMENS,NRVERT,NIINFO,NRINFO,NumFun,NRVACA, &
BOTTIH,BOTTRH,IWork,Inform)
END IF
IF (Inform /= 0) THEN
CALL Handle_Error(Inform,IFAIL)
RETURN
END IF
END IF
!-----------------------------------------------------------------------
!
! Call integration routine
!
If (EpsAlg) THEN
IF ( LRestart .AND. (PreJob /= ABS(LJob))) THEN
Inform = 3
ELSE
! Observe that only relevant array sections are passed !
CALL Global_Adapt_Extrap(DIMENS,CINFO,UINFO,NRVERT,NIINFO, &
NRINFO, Vertices(1:DIMENS,1:NRVERT,1:NUMRGN), &
RgType(1:NUMRGN),Integrand,LRestart, &
Value(1:NUMFUN),AbsErr(1:NUMFUN),LEval,Inform, &
RWork,IWork,M)
END IF
ELSE
IF ( LRestart .AND. (PreJob /= LJob)) THEN
IF ( ASSOCIATED(M%RESLA)) THEN
DEALLOCATE(M%RESLA,M%ERLARG,M%RESULT1,M%ABSERR1,M%RCOPY)
END IF
CALL DSPINT(IWork,RWork)
CALL DSSUM(Value,Abserr,IWork,RWork,Inform)
END IF
! Observe that only relevant array sections are passed !
CALL Global_Adapt(DIMENS,CINFO,UINFO,NRVERT,NIINFO,NRINFO, &
Vertices(1:DIMENS,1:NRVERT,1:NUMRGN), &
RgType(1:NUMRGN),Integrand,LRestart, &
Value(1:NUMFUN),AbsErr(1:NUMFUN),LEval,Inform, &
RWork,IWork)
END IF
IF (PRESENT(NEval)) THEN
NEval = LEval
END IF
!-----------------------------------------------------------------------
IF (LJob < 0) THEN
WRITE(unit=*,fmt=*) "Debug mode: dumping region collection overview."
CALL DSSTAT(IWork(:),RWork(:)) ! For debugging.
END IF
!-----------------------------------------------------------------------
! IF ((Inform >= 8) .or. (Inform == 3)) THEN
! Something went wrong but the data structure remains untouched
! and so this call can be ignored.
IF ((Inform < 8) .AND. (Inform /= 3)) THEN
PreJob = ABS(LJob)
END IF
CALL Handle_Error(Inform,IFAIL)
RETURN
END SUBROUTINE CUBATR_X
SUBROUTINE CUBATR_1 &
(DIMENS,Integrand,SVertices,SRgType,SValue,SAbsErr, &
! and optional parameters &
IFAIL,Neval,EpsAbs,EpsRel,Restart,MaxPts,Key,Job)
!-----------------------------------------------------------------------
! Input parameters
! ----------------
!
! DIMENS Integer.
! The dimension of the region of integration.
!
! Integrand
! Externally declared function for computing all components
! of the integrand at the given evaluation point.
! It must have input parameter X:
! X(1) The x-coordinate of the evaluation point.
! X(2) The y-coordinate of the evaluation point.
! ...
! X(DIMENS) The z-coordinate of the evaluation point.
! and NumFun, the number of components of the integrand.
! It must be compatible with the following interface:
! INTERFACE
! FUNCTION Integrand(NUMFUN,X) RESULT(Value)
! USE Precision_Model
! INTEGER, INTENT(IN) :: NUMFUN
! REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
! REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
! END FUNCTION Integrand
! END INTERFACE
!
! SVertices
! Real array of dimension (DIMENS,DIMENS+1).
! Vertices(1:DIMENS,K) are the x, y, ... coordinates
! respectively of vertex K of the region, where
! K = 1,...,DIMENS+1.
!
! SRgType Integer.
! RgType describes the type of region L.
!
! SValue Real.
! Approximation to the integral if the procedure is restarted.
!
! SAbsErr Real.
! Estimate of the absolute error if the procedure is restarted.
!
! IFAIL Optional integer argument.
! This follows the NAG convention:
! IFAIL = 1 : soft silent error
! Control returned to calling program.
! IFAIL = -1: soft noisy error
! Error message is printed.
! Control returned to calling program.
! IFAIL = 0 : hard noisy error
! Error message is printed and program is stopped.
! Default IFAIL = -1.
!
! EpsAbs Optional real argument.
! Requested absolute error.
! Default EpsAbs = 0.
!
! EpsRel Optional real argument.
! Requested relative error.
! Default EpsRel = sqrt(machine precision).
!
! Restart Optional boolean argument.
! If Restart = FALSE, this is the first attempt to compute
! the integral.
! If Restart = TRUE, then we restart a previous attempt.
! In this case the only parameters for CUBATR that may
! be changed (with respect to the previous call of CUBATR)
! are MinPts, MaxPts, EpsAbs, EpsRel, Key and Restart.
! Default Restart = FALSE.
!
! MaxPts Optional integer argument.
! The maximum allowed number of integrand evaluations.
! Default MaxPts = enough to do 500 subdivisions.
!
! Key Optional integer argument.
! Can be used by Rule_General to choose between several
! local integration rules.
! Default Key = 2 if Dimension=1 and extrapolation is used
! (This corresponds to QAGS)
! Default Key = 0 otherwise
!
! Job Optional integer argument.
! If |Job| = 0, then nothing will be done except freeing all
! allocated memory.
! This is usefull after a call of CUBATR if no
! Restart will be done later and memory usage
! might become an issue later.
! Equivalently, one can call CUBATR()
! without any arguments.
! = 1, the global adaptive algorithm is called
! = 2, extrapolation using the epsilon algorithm is used.
! = 11, a region will be divided in 2**DIMENS subregions
! and the global adaptive algorithm is called.
! In combination with Key=0, this resembles DUCTRI and DCUTET.
! = 12, a region will be divided in 2 subregions
! and the global adaptive algorithm is called.
! In combination with Key=3 or 4, this resembles DCUHRE.
! If Job < 0, then an overview of the Region Collection is dumped.
! This will create the files tmp_integerstore and tmp_realstore.
! Default Job = 1.
!
! Output parameters
! -----------------
!
! SValue Real.
! Approximation to the integral
!
! AbsErr Real.
! Estimate of the absolute error.
!
! NEval Optional Integer.
! Number of integrand evaluations used by CUBATR for this call.
!
! IFAIL Optional Integer.
! IFAIL = 0 for normal exit.
!
! AbsErr(K) <= EpsAbs or
! AbsErr(K) <= ABS(Value(K))*EpsRel with MaxPts or less
! function evaluations for all values of K,
! 1 <= K <= NumFun .
!
! IFAIL = 1 if MaxPts was too small to obtain the required
! accuracy. In this case Global_Adapt returns values of
! Value with estimated absolute errors AbsErr.
!
! IFAIL > 1 in more serious case of trouble.
!-----------------------------------------------------------------------
!
! Global variables
!
INTERFACE
FUNCTION Integrand(NUMFUN,X) RESULT(Value)
USE Precision_Model
INTEGER, INTENT(IN) :: NUMFUN
REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
END FUNCTION Integrand
END INTERFACE
LOGICAL, OPTIONAL, INTENT(IN) :: Restart
INTEGER, INTENT(IN) :: DIMENS,SRgType
INTEGER, INTENT(OUT), OPTIONAL :: NEval
INTEGER, INTENT(IN), OPTIONAL :: Key,MaxPts,Job
INTEGER, INTENT(IN OUT), OPTIONAL :: IFAIL
REAL(kind=stnd), INTENT(IN), OPTIONAL :: EpsAbs,EpsRel
REAL(kind=stnd), DIMENSION(:,:), INTENT(IN) :: SVertices
REAL(kind=stnd), INTENT(IN OUT) :: SValue,SAbsErr
!
! Local variables
!
INTEGER, DIMENSION(1) :: RgType
REAL(kind=stnd), DIMENSION(1) :: Value, AbsErr
REAL(kind=stnd), DIMENSION(DIMENS,DIMENS+1,1) :: Vertices
!-------------------
RgType(1) = SRgType
Vertices(:,:,1) = SVertices
IF (PRESENT(Restart)) THEN
IF ( Restart ) THEN
Value(1) = SValue
AbsErr(1) = SAbsErr
END IF
END IF
CALL CUBATR &
(DIMENS,1,Integrand,1,Vertices,RgType,Value,AbsErr, &
! and optional parameters &
ifail=IFAIL,neval=Neval,epsabs=EpsAbs,epsrel=EpsRel, &
restart=Restart,maxpts=MaxPts,key=key,job=Job)
SValue = Value(1)
SAbsErr = AbsErr(1)
RETURN
END SUBROUTINE CUBATR_1
SUBROUTINE CUBATR_CLEAR()
IF ( ALLOCATED(Iwork) ) THEN
DEALLOCATE(RWork,IWork)
END IF
IF ( ASSOCIATED(M%RESLA)) THEN
DEALLOCATE(M%RESLA,M%ERLARG,M%RESULT1,M%ABSERR1,M%RCOPY)
END IF
PreJob = 0
RETURN
END SUBROUTINE CUBATR_CLEAR
SUBROUTINE CUBPACK_INFO()
REAL(kind=stnd) :: x=1.0e-30 ! lowest accuracy of cubature formula constants
print *," ---------------------------------------------------------------"
print *," CUBPACK information"
print *," -------------------"
print *," The model for real numbers in the current installed version,"
print *," obtained with the declaration REAL(KIND=stnd), has the"
print *," following characteristics:"
print *," base = ",radix(x)
print *," digits in this base = ", digits(x)
! print *," highest exponent = ", maxexponent(x)
! print *," lowest exponent = "minexponent(x), "(normalized numbers)
print *," This implies:"
print *," machine epsilon = ",epsilon(x)
print *," largest real number = ", huge(x)
print *," smallest normalized number = ", tiny(x)
print *
print *," The lowest relative error that may be obtained with this"
print "("" version is about "",G8.2)",max(50*epsilon(x),x)
print *," Asking for lower error will push the routine to use the"
print *," maximal number of function evaluations it is allowed."
print *
print *
print *," This version accepts a collection of hyper-rectangles"
print *," (and parallelepipeds) and simplices as integration regions."
print *," Extrapolation using the epsilon-algorithm is available"
print *," for dimensions 1, 2 and 3."
print *," The following values of KEY give different integration rules:"
print *," - finite interval: KEY = 1, 2, 3, 4, 5."
print *," KEY < 1 defaults to 1; KEY > 5 defaults to 5."
print *," - n-cube: KEY = 3, 4 uses rule of degree 2*KEY+1"
print *," otherwise, uses for a square a rule of degree 13"
print *," 3-cube a rule of degree 11"
print *," a rule of degree 7"
print *," - n-simplex: KEY = 1, 2, 3, 4 uses rule of degree 2*KEY+1"
print *," otherwise, uses for a triangle a rule of degree 13"
print *," tetrahedron a rule of degree 8"
print *," a rule of degree 7"
print *," KEY = 0 corresponds to our preferred choice."
print *," ---------------------------------------------------------------"
RETURN
END SUBROUTINE CUBPACK_INFO
END Module CUI
SHAR_EOF
fi # end of overwriting check
if test -f 'divide.f90'
then
echo shar: will not over-write existing file "'divide.f90'"
else
cat << "SHAR_EOF" > 'divide.f90'
! This file is F-compatible, except for upper/lower case conventions.
!--------------------------------------------------------------------
MODULE Subdivisions
USE Precision_Model, ONLY: stnd
USE internal_types
Implicit NONE
PRIVATE
PUBLIC :: DIVIDE
PRIVATE :: DIVSMP
CONTAINS
SUBROUTINE DIVIDE(DIMENS,NRVERT,MAXSUB,UNIFORM_SUBDIV,NUMFUN,VEROLD, &
INFOLD,RINFOL,Integrand,OUTSUB,NUM,VERNEW,INFNEW,RINFNE,IFAIL)
!***BEGIN PROLOGUE DIVIDE
!***DATE WRITTEN 900615 (YYMMDD)
!***REVISION DATE 910506 (YYMMDD)
!***REVISION DATE 970401 (YYMMDD) (subdivision modifications)
!***REVISION DATE 980331 (YYMMDD) (1D added)
!***REVISION DATE 980406 (YYMMDD) (2div for nCube added)
!***REVISION DATE 980408 (YYMMDD) (F conversion)
!***REVISION DATE 990525 (YYMMDD) (re-organising)
!***REVISION DATE 990602 (YYMMDD) (2/4-division for T2 activated)
!***REVISION DATE 990624 (YYMMDD) (2/4/8-division for C3 activated)
!***REVISION DATE 010919 (YYMMDD) (infold(4) code changed elsewhere)
!***REVISION DATE 020716 (YYMMDD) (replace MOD intrinsic by MODULO)
!***AUTHOR
! Ronald Cools, Dept. of Computer Science,
! Katholieke Universiteit Leuven, Celestijnenlaan 200A,
! B-3001 Heverlee, Belgium
! Email: Ronald.Cools@cs.kuleuven.ac.be
!
! Alan Genz
! Department of Mathematics
! Washington State University
! Pullman, WA 99164-3113, USA
! Email: AlanGenz@wsu.edu
!
!***PURPOSE To divide a given region in OUTSUB subregions
! with equal volume.
!
!***DESCRIPTION
! Input parameters
! ----------------
! DIMENS Integer, dimension of the regions
! NRVERT Integer, number of vertices that describe a region
! MAXSUB Integer.
! The given region is divided in at most MAXSUB subregions
! UNIFORM_SUBDIV Logical.
! If true, this routine does a 2**dim subdivision.
! If false, this routine may decide to divide into less than
! MAXSUB regions
! NUMFUN Integer, number of components of the integral.
! VEROLD Real array of dimension (dimens,nrvert)
! Contains the vertices of the given region
! INFOLD
! RINFOL
! Integrand
! Externally declared function for computing all components
! of the integrand at the given evaluation point.
! It must have input parameter X:
! X(1) The x-coordinate of the evaluation point.
! X(2) The y-coordinate of the evaluation point.
! ...
! X(DIMENS) The z-coordinate of the evaluation point.
! and NumFun, the number of components of the integrand.
! It must be compatible with the following interface:
! INTERFACE
! FUNCTION Integrand(NumFun,X)
! USE Precision_Model
! INTEGER :: NumFun
! REAL(kind=stnd) :: X(:)
! REAL(kind=stnd) :: Integrand(NumFun)
! END
! END INTERFACE
!
! Output parameters
! -----------------
! NUM : Integer number of integrand values used to decide subdivision.
! OUTSUB : Integer number of subregions the given region was divided into.
! VERNEW : Real array of dimension (dimens,nrvert,MAXSUB)
! Contains the vertices of the subregions.
! INFNEW : Integer array of dimension (:,MAXSUM)
! Contains additional information for each subregion.
! RINFNE : Real array of dimension (:,MAXSUM)
! Contains additional information for each subregion.
! IFAIL : integer to indicate success or failure
! IFAIL = 0 on normal exit
! IFAIL = 7 if the desired subdivision is not implemented
! for this type of region
! IFAIL = 6 if no subdivsions are implemented for this type
! of region
!
!***ROUTINES CALLED DIVSMP
!***END PROLOGUE DIVIDE
!
! Global variables
INTEGER, INTENT(IN) :: DIMENS,NRVERT,MAXSUB,NUMFUN
INTEGER, DIMENSION(:), INTENT(IN) :: INFOLD
LOGICAL, INTENT(IN) :: UNIFORM_SUBDIV
REAL(kind=stnd), DIMENSION(:,:), INTENT(IN) :: VEROLD
REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: RINFOL
INTEGER, INTENT(OUT) :: NUM,OUTSUB,IFAIL
INTEGER, DIMENSION(:,:), INTENT(OUT) :: INFNEW
REAL(kind=stnd), DIMENSION(:,:,:), INTENT(OUT):: VERNEW
REAL(kind=stnd), DIMENSION(:,:), INTENT(OUT) :: RINFNE
INTERFACE
FUNCTION Integrand(NUMFUN,X) RESULT(Value)
USE Precision_Model
INTEGER, INTENT(IN) :: NUMFUN
REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
END FUNCTION Integrand
END INTERFACE
!
! Local variables
!
! GEOMETRY : specifies the type of region
!
INTEGER :: I,J,GEOMETRY
REAL(kind=stnd) :: VOLUME
REAL(kind=stnd), DIMENSION(3) :: HALF_HEIGHT, HEIGHT
!
!***FIRST EXEUTABLE STATEMENT
IF (MAXSUB == 1) THEN
VERNEW(:,:,1) = VEROLD
INFNEW(:,1) = INFOLD
RINFNE(:,1) = RINFOL
OUTSUB = 1
NUM = 0
IFAIL = 0
RETURN
END IF
!
GEOMETRY = INFOLD(1)
SELECT CASE (GEOMETRY)
CASE (Simplex)
IF (DIMENS == 1 ) THEN
VERNEW(:,:,1) = VEROLD
VERNEW(:,:,2) = VEROLD
VERNEW(:,2,1) = (VEROLD(1,1)+VEROLD(1,2))/2
VERNEW(:,1,2) = VERNEW(:,2,1)
OUTSUB = 2
NUM = 0
!
ELSE IF ((DIMENS == 2) .AND. UNIFORM_SUBDIV ) THEN
VERNEW(:,:,1) = VEROLD
VERNEW(:,:,2) = VEROLD
VERNEW(:,:,3) = VEROLD
VERNEW(:,:,4) = VEROLD
VERNEW(:,2,1) = ( VEROLD(:,1) + VEROLD(:,2) )/2
VERNEW(:,3,1) = ( VEROLD(:,1) + VEROLD(:,3) )/2
VERNEW(:,1,2) = VERNEW(:,2,1)
VERNEW(:,3,2) = ( VEROLD(:,2) + VEROLD(:,3) )/2
VERNEW(:,1,3) = VERNEW(:,3,1)
VERNEW(:,2,3) = VERNEW(:,3,2)
VERNEW(:,1,4) = VERNEW(:,3,2)
VERNEW(:,2,4) = VERNEW(:,2,1)
VERNEW(:,3,4) = VERNEW(:,3,1)
OUTSUB = 4
NUM = 0
!
ELSE IF ((DIMENS == 3) .AND. UNIFORM_SUBDIV ) THEN
VERNEW(:,1,1) = ( VEROLD(:,1) + VEROLD(:,4) )/2
VERNEW(:,2,1) = ( VEROLD(:,1) + VEROLD(:,3) )/2
VERNEW(:,3,1) = ( VEROLD(:,1) + VEROLD(:,2) )/2
VERNEW(:,4,1) = VEROLD(:,1)
VERNEW(:,:,2) = VERNEW(:,:,1)
VERNEW(:,4,2) = ( VEROLD(:,2) + VEROLD(:,4) )/2
VERNEW(:,:,3) = VERNEW(:,:,2)
VERNEW(:,3,3) = ( VEROLD(:,3) + VEROLD(:,4) )/2
VERNEW(:,:,4) = VERNEW(:,:,2)
VERNEW(:,1,4) = ( VEROLD(:,2) + VEROLD(:,3) )/2
VERNEW(:,:,5) = VERNEW(:,:,3)
VERNEW(:,2,5) = VEROLD(:,4)
VERNEW(:,:,6) = VERNEW(:,:,4)
VERNEW(:,3,6) = VERNEW(:,3,3)
VERNEW(:,:,7) = VERNEW(:,:,6)
VERNEW(:,4,7) = VEROLD(:,3)
VERNEW(:,:,8) = VERNEW(:,:,4)
VERNEW(:,2,8) = VEROLD(:,2)
OUTSUB = 8
NUM = 0
!
ELSE IF (.NOT. UNIFORM_SUBDIV ) THEN
CALL DIVSMP( DIMENS, NUMFUN, MAXSUB, VEROLD, Integrand, &
NUM, OUTSUB, VERNEW )
ELSE
IFAIL = 7
RETURN
END IF
!
CASE (Hyperrectangle)
SELECT CASE (DIMENS)
CASE (1)
VERNEW(:,:,1) = VEROLD
VERNEW(:,:,2) = VEROLD
VERNEW(:,2,1) = (VEROLD(1,1)+VEROLD(1,2))/2
VERNEW(:,1,2) = VERNEW(:,2,1)
OUTSUB = 2
NUM = 0
CASE (2)
IF ((MAXSUB == 4) .OR. UNIFORM_SUBDIV) THEN
!
! Divide a parallellogram in 4.
!
VERNEW(:,1,1) = VEROLD(:,1)
VERNEW(:,2,1) = ( VEROLD(:,1) + VEROLD(:,2) )/2
VERNEW(:,3,1) = ( VEROLD(:,1) + VEROLD(:,3) )/2
VERNEW(:,1,2) = ( VEROLD(:,2) + VEROLD(:,3) )/2
VERNEW(:,2,2) = VEROLD(:,2) + ( VEROLD(:,3) - VEROLD(:,1) )/2
VERNEW(:,3,2) = VEROLD(:,3) + ( VEROLD(:,2) - VEROLD(:,1) )/2
VERNEW(:,1,3) = VERNEW(:,2,1)
VERNEW(:,2,3) = VEROLD(:,2)
VERNEW(:,3,3) = VERNEW(:,1,2)
VERNEW(:,1,4) = VERNEW(:,3,1)
VERNEW(:,2,4) = VERNEW(:,1,2)
VERNEW(:,3,4) = VEROLD(:,3)
OUTSUB = 4
NUM = 0
ELSE IF (MAXSUB == 2) THEN
!
! Divide a parallellogram in 2.
!
IF ( INFOLD(4) == 2 ) THEN
! Cut orthogonal to the line through vertices 1-3
VERNEW(:,1,1) = VEROLD(:,1)
VERNEW(:,2,1) = VEROLD(:,2)
VERNEW(:,3,1) = ( VEROLD(:,1) + VEROLD(:,3) )/2
VERNEW(:,1,2) = VERNEW(:,3,1)
VERNEW(:,2,2) = VEROLD(:,2)+( VEROLD(:,3) - VEROLD(:,1) )/2
VERNEW(:,3,2) = VEROLD(:,3)
ELSE
! Cut orthogonal to the line through vertices 1-2
VERNEW(:,1,1) = VEROLD(:,1)
VERNEW(:,2,1) = ( VEROLD(:,1)+VEROLD(:,2) )/2
VERNEW(:,3,1) = VEROLD(:,3)
VERNEW(:,1,2) = VERNEW(:,2,1)
VERNEW(:,3,2) = VEROLD(:,2)
VERNEW(:,2,2) = VEROLD(:,3)+( VEROLD(:,2) - VEROLD(:,1) )/2
END IF
OUTSUB = 2
NUM = 0
ELSE
IFAIL = 7
RETURN
END IF
!
CASE (3)
IF ((MAXSUB == 8) .OR. UNIFORM_SUBDIV) THEN
!
! Divide a 3D-Cube in 8.
!
HALF_HEIGHT = (verold(1:3,4)-verold(1:3,1))/2
vernew(1:3,1,1) = verold(1:3,1)
vernew(1:3,2,1) = (verold(1:3,1)+verold(1:3,2))/2
vernew(1:3,3,1) = (verold(1:3,1)+verold(1:3,3))/2
vernew(1:3,4,1) = (verold(1:3,1)+verold(1:3,4))/2
vernew(1:3,1,2) = vernew(1:3,2,1)
vernew(1:3,2,2) = verold(1:3,2)
vernew(1:3,3,2) = (verold(1:3,2)+verold(1:3,3))/2
vernew(1:3,4,2) = (verold(1:3,2)+verold(1:3,4))/2
vernew(1:3,1,3) = vernew(1:3,3,1)
vernew(1:3,2,3) = vernew(1:3,3,2)
vernew(1:3,3,3) = verold(1:3,3)
vernew(1:3,4,3) = (verold(1:3,3)+verold(1:3,4))/2
vernew(1:3,1,4) = vernew(1:3,3,2)
vernew(1:3,2,4) = verold(1:3,2) + &
(verold(1:3,3)-verold(1:3,1))/2
vernew(1:3,3,4) = verold(1:3,3) + &
(verold(1:3,2)-verold(1:3,1))/2
vernew(1:3,4,4) = vernew(1:3,3,2) + HALF_HEIGHT
vernew(1:3,1,5) = vernew(1:3,4,1)
vernew(1:3,2,5) = vernew(1:3,4,2)
vernew(1:3,3,5) = vernew(1:3,4,3)
vernew(1:3,4,5) = verold(1:3,4)
vernew(1:3,1,6) = vernew(1:3,4,2)
vernew(1:3,2,6) = verold(1:3,2) + HALF_HEIGHT
vernew(1:3,3,6) = vernew(1:3,4,4)
vernew(1:3,4,6) = vernew(1:3,4,2) + HALF_HEIGHT
vernew(1:3,1,7) = vernew(1:3,4,3)
vernew(1:3,2,7) = vernew(1:3,4,4)
vernew(1:3,3,7) = verold(1:3,3) + HALF_HEIGHT
vernew(1:3,4,7) = vernew(1:3,4,3) + HALF_HEIGHT
vernew(1:3,1,8) = vernew(1:3,4,4)
vernew(1:3,2,8) = vernew(1:3,2,4) + HALF_HEIGHT
vernew(1:3,3,8) = vernew(1:3,3,4) + HALF_HEIGHT
vernew(1:3,4,8) = vernew(1:3,4,4) + HALF_HEIGHT
OUTSUB = 8
NUM = 0
ELSE IF (maxsub == 4) THEN
!
! Divide a 3D-Cube in 4.
!
IF (infold(4)/10 == -2) THEN
! Cut orthogonal to the line through vertices 1-2 and vertices 1-3
HEIGHT = (verold(1:3,4) - verold(1:3,1))
vernew(1:3,1,1) = verold(1:3,1)
vernew(1:3,2,1) = (verold(1:3,1)+verold(1:3,2))/2
vernew(1:3,3,1) = (verold(1:3,1)+verold(1:3,3))/2
vernew(1:3,4,1) = verold(1:3,4)
vernew(1:3,1,2) = vernew(1:3,2,1)
vernew(1:3,2,2) = verold(1:3,2)
vernew(1:3,3,2) = (verold(1:3,2)+verold(1:3,3))/2
vernew(1:3,4,2) = vernew(1:3,2,1) + HEIGHT
vernew(1:3,1,3) = vernew(1:3,3,1)
vernew(1:3,2,3) = vernew(1:3,3,2)
vernew(1:3,3,3) = verold(1:3,3)
vernew(1:3,4,3) = vernew(1:3,3,1) + HEIGHT
vernew(1:3,1,4) = vernew(1:3,3,2)
vernew(1:3,2,4) = vernew(1:3,3,1) + &
(verold(1:3,2) - verold(1:3