Revision as of 18:03, 5 January 2016

Introduction

MPFUN90 is a Fortran-90 arbitrary precision package (i.e. user specifies the desired precision)written by David H. Bailey et. al from the Lawrence Berkeley National Laboratory. If 64 digits of precision are sufficient for an application then users should use QD instead of MPFUN which would run much faster.

Version Selection

20060123

module purge; module load ldwrapper; module load intel/12.1.3; module load mpfun90/20060123 

20100825

module purge; module load ldwrapper; module load intel/15.0.3; module load mpfun90/intel1503/20100825

Job Submission

To submit a job to the serial queue use:

 sqsub -r 10m --mpp=2G -o ofile.%J ./a.out

Sample Problem

Code: testmp90

The 2006 and 2010 modules are setup slightly differently, so the sample problem will be repeated for each : This example shows how to compile the test program testmp90.f90 using the SHARCNET mpfun build.

20060123

[roberpj@orc-login2:~] module purge; module load ldwrapper; module load intel/12.1.3; module load mpfun90/20060123
[roberpj@orc-login2:~/samples/mpfun90/20060123] cp /opt/sharcnet/mpfun90/20060123/f90/src/testmp90.f90 .
[roberpj@orc-login2:~/samples/mpfun90/20060123] ifort testmp90.f90 $CPPFLAGS $LDFLAGS -lmp
[roberpj@orc-login2:~/samples/mpfun90/20060123] ./a.out
Completed Test  1
Completed Test  2
Completed Test  3
.................
Completed Test 62

20100825

[roberpj@orc-login2:~] module purge; module load ldwrapper; module load intel/15.0.3; module load mpfun90/intel1503/20100825
[roberpj@orc-login2:~/samples/mpfun90/20100825] cp $MPFUN90_ROOT/f90/testmp90.f .
[roberpj@orc-login2:~/samples/mpfun90/20100825] ifort -free testmp90.f $MPFUN90_ROOT/toolkit/mpfun90.o  -I$MPFUN90_ROOT/toolkit
[roberpj@orc-login2:~/samples/mpfun90/20100825] ./a.out
Completed Test  1
Completed Test  2
Completed Test  3
.................
Completed Test 62

factorial_100

Compile and execute the factorial_100_pgm program using mpfun90/20100825 only:

! file = factorial_100.f90

      program factorial_100_pgm

      use mpmodule
      implicit none
      integer, parameter :: Nmx = 100
      integer            :: AllocateStatus
      integer            :: J
      type(mp_real) zero, one

      type(mp_real), dimension (:), allocatable        :: fac

      call mpinit(800)
      call mpsetprec(800)
      call mpsetoutputprec(700)

      zero='0.'
      one ='1.'

      allocate(fac(0:Nmx),STAT=AllocateStatus)
      IF (AllocateStatus /= 0) STOP "Not enough memory for fac "

      fac(0)=one

      do J=1,Nmx
        fac(J)=fac(J-1)*J
      end do

      J=4
      write(6,1001) J
 1001 format(/"fac(",i3.3,") =",$)
      call mpwrite(6,fac(J))

      J=100
      write(6,1001) J
      call mpwrite(6,fac(J))

      write(6,1005)
 1005 format(/"DONE"//)

      end

Compile

[roberpj@orc-login2:~/samples/mpfun90/factorial_100_pgm] ifort factorial_100.f90 $MPFUN90_ROOT/toolkit/mpmod90.o $MPFUN90_ROOT/toolkit/mpfun90.o -I$MPFUN90_ROOT/toolkit 

Execute

[roberpj@orc-login2:~/samples/mpfun90/factorial_100_pgm] ./a.out
fac(004) =10 ^         1 x  2.4,
fac(100) =10 ^       157 x  9.3326215443944152681699238856266700490715968264381621468592
 963895217599993229915608941463976156518286253697920827223758251185210916864,
DONE

General Notes

OpenMP with MPFUN and QD

The Multi Precision packages MPFUN90 and QD use "user defined data type" for their variables. Since the OpenMP reduction clause cannot be used for "user defined data types" to run MPFUN90 or QD with OpenMP the reduction must be done manually in a PARALLEL OpenMP region.

Detailed documentation on how to use OpenMP with MPFUN90 and QD are presented in the following document:

https://www.sharcnet.ca/help/index.php/OpenMP_Reduction_for_User_Defined_Data_Types

Handling large arrays in an MPFUN program

Declaring large arrays in an MPFUN program often leads to compile errors like:

 "relocation truncated to fit r_x86_64_pc32 against .bss"

Even with the flag '-mcmodel large' in the compile command the error will persist, because there is a limit to the amount of data you can have allocated in static arrays for x86_64 architecture. This data segment is limited to 2GB. The data in COMMON plus any other statically declared arrays PLUS your code must fit in 2GB.

In order to overcome this problem use allocatable arrays declared in a module together with a short subroutine which allocates the arrays.

We illustrate these concepts by starting with a program 'betax.f90' which work fine for values of parameter 'Nsz' less than 37. Thus, if you compile it with

      integer, parameter :: Nsz=36

it works, but for any higher values it fails with error message:

 "relocation truncated to fit: R_X86_64_PC32 against symbol '...'"

To overcome this problem, we first copy program betax.f90 into a new file, say, betay.f90, and write a module called 'mp_data' and a subroutine called 'allocate_arrays' at the top of the new file. We move all the declarations of large arrays into the module, and declare them as 'allocatable'.

In the subroutine 'allocate_arrays' we do the allocations and in the main program we add the statements:

     USE mp_data
     ...
     call allocate_arrays
     ...

Below you will find the listings for the files 'betax.f90', 'betay.f90' and the 'HIST' file containing the commands to compile these files and submit jobs:

! file betax.f90

      program BETAX

      use mpmodule
      implicit none
!     integer, parameter :: Nsz=12   ! compilation fails for 37
      integer, parameter :: Nsz=37

      INTEGER :: I,J,K,L

      type(mp_real) zero, one

      type(mp_real) fac(0:Nsz),fci(0:Nsz)
      type(mp_real) BETA2(Nsz,Nsz)
      type(mp_real) BETA3(Nsz,Nsz,Nsz)
      type(mp_real) BETA4(Nsz,Nsz,Nsz,Nsz)


      call mpinit(800)
      call mpsetprec(800)
      call mpsetoutputprec(700)

      zero='0.'
      one ='1.'

      fac(0)=one
      fci(0)=one

      do j=1,Nsz
        fac(j)=fac(j-1)*j
        fci(j)=one/fac(j)
      end do

      I=Nsz/2
      J=I-1
      K=J-1
      L=K-1

      write(6,1000) I,J,K,L
 1000 format(/"I,J,K,L = ",4i3)
 
      write(6,1001) I
 1001 format(/"fac(",i3.3,") =",$)
      call mpwrite(6,fac(I))

      BETA2(I,J)      =fac(I)*fac(J)*fci(I+J)
      BETA3(I,J,K)    =fac(I)*fac(J)*fac(K)*fci(I+J)*fci(J+K)
      BETA4(I,J,K,L)  =fac(I)*fac(J)*fac(K)*fac(L)*                     &
     &                 fci(I+J)*fci(J+K)*fci(K+L)

      write(6,1002) I,J
 1002 format(/"BETA2(",i3.3,",",i3.3,") =")
      call mpwrite(6,BETA2(I,J))

      write(6,1003) I,J,K
 1003 format(/"BETA3(",i3.3,",",i3.3,",",i3.3,") =")
      call mpwrite(6,BETA3(I,J,K))

      write(6,1004) I,J,K,L
 1004 format(/"BETA4(",i3.3,",",i3.3,",",i3.3,",",i3.3,") =")
      call mpwrite(6,BETA4(I,J,K,L))

      write(6,1005)
 1005 format(/"DONE"//)

      stop
      end
 

! file betay.f90

      module mp_data

      use mpmodule
      use mpfunmod

      implicit none

      integer, parameter :: Nsz =48

      integer            :: AllocateStatus

!     type(mp_real) zero, one
!
!     type(mp_real) fac(0:Nsz),fci(0:Nsz)
!     type(mp_real) BETA2(Nsz,Nsz)
!     type(mp_real) BETA3(Nsz,Nsz,Nsz)
!     type(mp_real) BETA4(Nsz,Nsz,Nsz,Nsz)

      type(mp_real) zero, one

      type(mp_real), dimension (:), allocatable        :: fac,fci
      type(mp_real), dimension (:,:), allocatable      :: BETA2
      type(mp_real), dimension (:,:,:), allocatable    :: BETA3
      type(mp_real), dimension (:,:,:,:), allocatable  :: BETA4

      end module mp_data

      subroutine allocate_arrays

      USE mp_data
      implicit none

      allocate(fac(0:Nsz),STAT=AllocateStatus)
      IF (AllocateStatus /= 0) STOP "Not enough memory for fac "

      allocate(fci(0:Nsz),STAT=AllocateStatus)
      IF (AllocateStatus /= 0) STOP "Not enough memory for fci "

      allocate(BETA2(Nsz,Nsz),STAT=AllocateStatus)
      IF (AllocateStatus /= 0) STOP "Not enough memory for BETA2"

      allocate(BETA3(Nsz,Nsz,Nsz),STAT=AllocateStatus)
      IF (AllocateStatus /= 0) STOP "Not enough memory for BETA3"

      allocate(BETA4(Nsz,Nsz,Nsz,Nsz),STAT=AllocateStatus)
      IF (AllocateStatus /= 0) STOP "Not enough memory for BETA4"

      print *,"Done with allocations"

      return
      END

! ----------------------------------------------------------------------

      program BETAY
      USE mp_data

      implicit none

      INTEGER :: I,J,K,L

      call mpinit(800)
      call mpsetprec(800)
      call mpsetoutputprec(700)

      zero='0.'
      one ='1.'

      call allocate_arrays

      fac(0)=one
      fci(0)=one

      do j=1,Nsz
        fac(j)=fac(j-1)*j
        fci(j)=one/fac(j)
      end do

      I=Nsz/2
      J=I-1
      K=J-1
      L=K-1

      write(6,1000) I,J,K,L
 1000 format(/"I,J,K,L = ",4i3)

      write(6,1001) I
 1001 format(/"fac(",i3.3,") =",$)
      call mpwrite(6,fac(I))

      BETA2(I,J)      =fac(I)*fac(J)*fci(I+J)
      BETA3(I,J,K)    =fac(I)*fac(J)*fac(K)*fci(I+J)*fci(J+K)
      BETA4(I,J,K,L)  =fac(I)*fac(J)*fac(K)*fac(L)*                     &
     &                 fci(I+J)*fci(J+K)*fci(K+L)

      write(6,1002) I,J
 1002 format(/"BETA2(",i3.3,",",i3.3,") =")
      call mpwrite(6,BETA2(I,J))

      write(6,1003) I,J,K
 1003 format(/"BETA3(",i3.3,",",i3.3,",",i3.3,") =")
      call mpwrite(6,BETA3(I,J,K))

      write(6,1004) I,J,K,L
 1004 format(/"BETA4(",i3.3,",",i3.3,",",i3.3,",",i3.3,") =")
      call mpwrite(6,BETA4(I,J,K,L))

      write(6,1005)
 1005 format(/"DONE"//)

      stop
      end

# file HIST
#
# Set Nsz=12

        vi betax.f90

and set Nsz=12 in file betax.f90

# Compile (copy/paste the three lines)

ifort betax.f90                                   \
-I/opt/sharcnet/mpfun/6.1.23/f90/include          \
-L/opt/sharcnet/mpfun/6.1.23/f90/lib -lmp >& AOUT1

# Execute code (it should work)

        ./a.out

# Set Nsz=40

        vi betax.f90

and set Nsz=40 in file betax.f90

# Re-compile and look at file AOUT1

ifort betax.f90                                   \
-I/opt/sharcnet/mpfun/6.1.23/f90/include          \
-L/opt/sharcnet/mpfun/6.1.23/f90/lib -lmp >& AOUT1


============================================================================

# File betay.f90 was copied from betax.f90 and modified as described earlier.

# Compile betay.f90

ifort betay.f90 -mcmodel medium -shared-intel \
-I/opt/sharcnet/mpfun/6.1.23/f90/include          \
-L/opt/sharcnet/mpfun/6.1.23/f90/lib -lmp >& AOUT2

============================================================================

# The executable might require more than 2.0G of memory: 

# This will pruduce an error: 

sqsub -r 5m  --mpp=4.0G  -o OUTPUTFILE4 ./a.out

# This should work:

sqsub -r 5m  --mpp=8.0G  -o OUTPUTFILE8 ./a.out

============================================================================

MPFUN's TOOLKIT

This toolkit provides and extensive example of how the mpfun package can be used. Read the header of files

mathinit.f90 
    and 
mathtool.f90 

found under

/opt/sharcnet/mpfun90/version/[flavor]/toolkit 

for an explanation what the package does. The Makefile in this directory shows how the mpfun libraries were linked to compile this sample program pair. The following steps shows the first initialization step followed by actually running the mathtool program to perform some examples such as performing an integration and displaying the value of pie. Note that the displayed zero cputimes are due to a bug in the mpfun current build which will be fixed in future updates on the clusters.

Initializing Mathtool

To run Mathtool you must first generate following two input files

const.dat 
  and 
quadts.dat 

which can be accomplished by running following command:

[roberpj@bal34:~/samples/mpfun] /opt/sharcnet/mpfun90/version/[flavor]/toolkit/mathinit
 mathinit: start
 const complete
 file const.dat written
 cpu time = 0.E+0
initqts: Tanh-sinh quadrature initialization
 0 18632
 1000 18632
 2000 18632
 3000 18632
 4000 18632
 5000 18632
 6000 18632
 7000 18632
 8000 18632
initqts: Table spaced used =    8177
 initqts complete
 file quadts.dat written
 cpu time = 0.E+0
 total cpu time = 0.E+0

[roberpj@bal34:~/samples/mpfun] ls
const.dat  quadts.dat

Running Mathtool

Use the files created in the previous section to run mathtool as follows:

[roberpj@bal34:~/samples/mpfun] /opt/sharcnet/mpfun90/version/[flavor]/toolkit/mathtool
Welcome to the Experimental Mathematician's Toolkit
Initializing...

Current settings:
Debug level =   2
Primary precision level =   100 digits
Secondary precision level =   200 digits
Primary epsilon = 10^  -100
Secondary epsilon = 10^  -200
PSLQ bound = 100
PSLQ level =   1
Quadrature level =   6
Quadrature type =   3  (Tanh-Sinh)
-- snip --

TYPE: integrate[1/(1+x^2), {x, 0, 1}]   (hit enter)
quadts: Iteration  1 of  6; est error = 10^    0; approx value =
       7.853157298876894309110693190063938396203990510665917998384243e-1
quadts: Iteration  2 of  6; est error = 10^    0; approx value =
       7.853981605125528448099840689980111720633592558858121242731082e-1
quadts: Iteration  3 of  6; est error = 10^  -17; approx value =
       7.853981633974483082028533923930157568609007576776270726821597e-1
quadts: Iteration  4 of  6; est error = 10^  -36; approx value =
       7.853981633974483096156608458198757190586422063114153956005278e-1
quadts: Iteration  5 of  6; est error = 10^  -71; approx value =
       7.853981633974483096156608458198757210492923498437764552437361e-1
quadts: Iteration  6 of  6; est error = 10^ -101; approx value =
       7.853981633974483096156608458198757210492923498437764552437361e-1
Result[  1] =
     7.85398163397448309615660845819875721049292349843776455243736148076954e-1
CPU time =      0.0000

TYPE: functions  (hit enter)
Abs              Arccos           Arcsin           Arctan          
Arctan2          Bessel           Besselexp        Binomial        
Cos              Erf              Exp              Factorial       
Gamma            Integrate        Log              Max             
Min              Polyroot         Pslq             Result          
Sin              Sqrt             Sum              Table           
Tan              Zeta             Zetap            Zetaz           

TYPE: variables  (hit enter)
E                Log2             Log10            Pi              
Catalan          Eulergamma       Infinity         Arg1            
Arg2             Arg3             Arg4             Arg5            
Arg6             Arg7             Arg8             Arg9            
Debug            Digits           Digits2          Epsilon         
Epsilon2         Pslqbound        Pslqlevel        Quadlevel       
Quadtype         x                

TYPE: Pi  (hit enter)
Result[  2] =
      3.14159265358979323846264338327950288419716939937510582097494459230781e0
CPU time =      0.0000

More Mathtool commands

To start mathtool (assuming you have issued the command /opt/sharcnet/mpfun/current/toolkit/mathinit and have these two files: const.dat HIST quadts.dat - you type:

/opt/sharcnet/mpfun/current/toolkit/mathtool

Then try following commands:

integrate[2*x,{x,0,1}]
sqrt[1600]
help  Binomial
Binomial[3,2]
hypo[x,y] = sqrt[x^2+y^2]
z=hypo[3,4]
z
exit

This is what you should get from above commands:

integrate[2*x,{x,0,1}]
quadts: Iteration  1 of  6; est error = 10^    0; approx  value =
        1.000003359570811224742398159504692881030640064061042259417294e0
quadts: Iteration  2 of  6; est error = 10^    0; approx value =
         1.000000000000036584185557611216071731517571134591404403031033e0
quadts: Iteration  3 of  6; est error = 10^  -27; approx  value =
        1.000000000000000000000000000001046972220803162709182580101627e0
quadts: Iteration  4 of  6; est error = 10^  -60; approx value =
         1.000000000000000000000000000000000000000000000000000000000000e0
quadts: Iteration  5 of  6; est error = 10^ -101; approx  value =
      9.999999999999999999999999999999999999999999999999999999999999e-1
Result[  1] =
    9.99999999999999999999999999999999999999999999999999999999999999999999e-1
CPU time =      0.0000

sqrt[1600]
Result[  2] =
      4.00000000000000000000000000000000000000000000000000000000000000000000e1
CPU time =      0.0000

help  Binomial
Binomial[m,n] computes the binomial coefficient of integers (m,n).

Binomial[3,2]
Result[  3] =
     3.00000000000000000000000000000000000000000000000000000000000000000000e0
CPU time =      0.0000

hypo[x,y] = sqrt[x^2+y^2]
parse: new function name = hypo
CPU time =      0.0000

z=hypo[3,4]
parse: new variable name = z
CPU time =      0.0000

z
Result[  4] =
      5.00000000000000000000000000000000000000000000000000000000000000000000e0
CPU time =      0.0000

exit

References

o High-Precision Software Directory
http://crd.lbl.gov/~dhbailey/mpdist/index.html