diff --git a/GNU_Linux/gfortran_case/prec_real128/LibDualzn128/diff_mod.mod b/GNU_Linux/gfortran_case/prec_real128/LibDualzn128/diff_mod.mod deleted file mode 100644 index a4e4faf..0000000 Binary files a/GNU_Linux/gfortran_case/prec_real128/LibDualzn128/diff_mod.mod and /dev/null differ diff --git a/GNU_Linux/gfortran_case/prec_real128/LibDualzn128/dualzn_mod.mod b/GNU_Linux/gfortran_case/prec_real128/LibDualzn128/dualzn_mod.mod deleted file mode 100644 index f69e479..0000000 Binary files a/GNU_Linux/gfortran_case/prec_real128/LibDualzn128/dualzn_mod.mod and /dev/null differ diff --git a/GNU_Linux/gfortran_case/prec_real128/LibDualzn128/libdualzn.a b/GNU_Linux/gfortran_case/prec_real128/LibDualzn128/libdualzn.a deleted file mode 100644 index 2c1b0c6..0000000 Binary files a/GNU_Linux/gfortran_case/prec_real128/LibDualzn128/libdualzn.a and /dev/null differ diff --git a/GNU_Linux/gfortran_case/prec_real128/LibDualzn128/precision_mod.mod b/GNU_Linux/gfortran_case/prec_real128/LibDualzn128/precision_mod.mod deleted file mode 100644 index 470c0c1..0000000 Binary files a/GNU_Linux/gfortran_case/prec_real128/LibDualzn128/precision_mod.mod and /dev/null differ diff --git a/GNU_Linux/gfortran_case/prec_real128/examples/create_exec.sh b/GNU_Linux/gfortran_case/prec_real128/examples/create_exec.sh deleted file mode 100755 index 4fb2327..0000000 --- a/GNU_Linux/gfortran_case/prec_real128/examples/create_exec.sh +++ /dev/null @@ -1,13 +0,0 @@ -#!/bin/bash - -#Check if at least one argument (the source file) is provided -if [ "$#" -lt 1 ]; then - echo "Usage: $0 [output_executable_name]" - exit 1 -fi - -#Use the provided source file and output executable name or default to "a.out" -SOURCE_FILE=$1 -OUTPUT_NAME=${2:-a.out} - -gfortran -I../LibDualzn128 -o "$OUTPUT_NAME" "$SOURCE_FILE" -L../LibDualzn128 -ldualzn diff --git a/GNU_Linux/gfortran_case/prec_real128/examples/ex1.f90 b/GNU_Linux/gfortran_case/prec_real128/examples/ex1.f90 deleted file mode 100644 index 66e94b2..0000000 --- a/GNU_Linux/gfortran_case/prec_real128/examples/ex1.f90 +++ /dev/null @@ -1,37 +0,0 @@ -!gfortran -I../LibDualzn128 -o e1 ex1.f90 -L../LibDualzn128 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - - call set_order(5) !before anything we set the 'order' to work with - - z0 = (1.1_prec,2.2_prec) - r = 0 !we initialize the dual number to 0 - - !we set the 0-th and 1-th components. If dual numbers are used to - !calculate D^n f(z0) then r must be of the form r = r0 + 1*eps_1 - r%f(0) = z0 - r%f(1) = 1 - - fval = sin(r)**log(r*r) - - !writing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - print*,fval%f(k) - end do - - deallocate(r%f) - deallocate(fval%f) - -end program main - - - - diff --git a/GNU_Linux/gfortran_case/prec_real128/examples/ex2.f90 b/GNU_Linux/gfortran_case/prec_real128/examples/ex2.f90 deleted file mode 100644 index c935727..0000000 --- a/GNU_Linux/gfortran_case/prec_real128/examples/ex2.f90 +++ /dev/null @@ -1,53 +0,0 @@ -!gfortran -I../LibDualzn128 -o e2 ex2.f90 -L../LibDualzn128 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - real :: t1,t2 - - call set_order(5) !we set the order to work with - - !since a dualzn numbers is an allocatable entity, do not forget to - !initialize it - r = 0 !<--- initializing r to 0 - r%f(0) = (1.1_prec,0.0_prec) - r%f(1) = 1 !since we want to differentiate, r = r0 +1*eps_1 - !all the other components are 0 as r was initialized to 0 - - call cpu_time(t1) - fval = ftest(r) - call cpu_time(t2) - - !Computing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - write(*,"(i0,a,f0.1,a,e17.10)") k,"-th derivative at x = ", & - real(r%f(0)),":",real(fval%f(k)) - end do - - print*,"elapsed time (s):",t2-t1 - - deallocate(r%f) - deallocate(fval%f) - -contains - function ftest(x) result(fr) - type(dualzn), intent(in) :: x - type(dualzn) :: fr - integer :: k - - !nested function f(x) = sin(x) * exp(-x^2), f(f(...(f(x))...)) - !applied 1000 times - fr = sin(x)*exp(-x*x) - do k=1, 1000-1 - fr = sin(fr)*exp(-fr*fr) - end do - end function ftest -end program main - - diff --git a/GNU_Linux/gfortran_case/prec_real128/examples/ex3.f90 b/GNU_Linux/gfortran_case/prec_real128/examples/ex3.f90 deleted file mode 100644 index 041011e..0000000 --- a/GNU_Linux/gfortran_case/prec_real128/examples/ex3.f90 +++ /dev/null @@ -1,90 +0,0 @@ -!gfortran -I../LibDualzn128 -o ExDO ex3.f90 -L../LibDualzn128 -ldualzn - -!module with example of functions -module function_mod - use dualzn_mod - implicit none - private - - public :: fstest, fvectest - -contains - !Example of scalar function f = sin(x*y*z) + cos(x*y*z) - function fstest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn) :: fr - type(dualzn) :: x,y,z - - x = r(1); y = r(2); z = r(3) - fr = sin(x*y*z) + cos(x*y*z) - end function fstest - - !Example of vector function f = [f1,f2,f3] - !f = fvectest(r) is a function f:D^m ---> Dn - function fvectest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn), allocatable, dimension(:) :: fr - type(dualzn) :: f1,f2,f3 - type(dualzn) :: x,y,z,w - - x = r(1); y = r(2); z = r(3); w = r(4) - - f1 = sin(x*y*z*w) - f2 = cos(x*y*z*w)*sqrt(w/y - x/z) - f3 = sin(log(x*y*z*w)) - - allocate(fr(3)) - fr = [f1,f2,f3] - end function fvectest -end module function_mod - -!----------------------------------------------------------------------- -!main program -program main - use precision_mod - use dualzn_mod - use diff_mod - use function_mod - implicit none - - integer, parameter :: nf =3, mq = 4 - complex(prec), parameter :: ii = (0,1) - complex(prec), dimension(mq) :: q, vec - complex(prec), dimension(nf,mq) :: Jmat - complex(prec), dimension(nf) :: JV - complex(prec), dimension(3) :: GV - complex(prec), dimension(3,3) :: Hmat - integer :: i - - vec = [1.0_prec,2.0_prec,3.0_prec,4.0_prec] - q = vec/10.0_prec + ii - - print*,"Jv using matmul" - Jmat = Jacobian(fvectest, q , nf) - JV = matmul(Jmat,vec) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"Jv using vector directional derivative" - JV = d1fvector(fvectest,vec,q,nf) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"---Hessian matrix---" - Hmat = Hessian(fstest,q(1:3)) - do i=1,3 - write(*,"(A,i0)") "row:",i - write(*,*) Hmat(i,:) - end do - write(*,*) - - print*,"---Gradient---" - GV = gradient(fstest,q(1:3)) - do i=1,3 - write(*,*) GV(i) - end do -end program main diff --git a/GNU_Linux/gfortran_case/prec_real64/LibDualzn64/diff_mod.mod b/GNU_Linux/gfortran_case/prec_real64/LibDualzn64/diff_mod.mod deleted file mode 100644 index 3610a69..0000000 Binary files a/GNU_Linux/gfortran_case/prec_real64/LibDualzn64/diff_mod.mod and /dev/null differ diff --git a/GNU_Linux/gfortran_case/prec_real64/LibDualzn64/dualzn_mod.mod b/GNU_Linux/gfortran_case/prec_real64/LibDualzn64/dualzn_mod.mod deleted file mode 100644 index a7500d3..0000000 Binary files a/GNU_Linux/gfortran_case/prec_real64/LibDualzn64/dualzn_mod.mod and /dev/null differ diff --git a/GNU_Linux/gfortran_case/prec_real64/LibDualzn64/libdualzn.a b/GNU_Linux/gfortran_case/prec_real64/LibDualzn64/libdualzn.a deleted file mode 100644 index 5e339c0..0000000 Binary files a/GNU_Linux/gfortran_case/prec_real64/LibDualzn64/libdualzn.a and /dev/null differ diff --git a/GNU_Linux/gfortran_case/prec_real64/LibDualzn64/precision_mod.mod b/GNU_Linux/gfortran_case/prec_real64/LibDualzn64/precision_mod.mod deleted file mode 100644 index a6ca63f..0000000 Binary files a/GNU_Linux/gfortran_case/prec_real64/LibDualzn64/precision_mod.mod and /dev/null differ diff --git a/GNU_Linux/gfortran_case/prec_real64/examples/create_exec.sh b/GNU_Linux/gfortran_case/prec_real64/examples/create_exec.sh deleted file mode 100755 index d80997e..0000000 --- a/GNU_Linux/gfortran_case/prec_real64/examples/create_exec.sh +++ /dev/null @@ -1,13 +0,0 @@ -#!/bin/bash - -#Check if at least one argument (the source file) is provided -if [ "$#" -lt 1 ]; then - echo "Usage: $0 [output_executable_name]" - exit 1 -fi - -#Use the provided source file and output executable name or default to "a.out" -SOURCE_FILE=$1 -OUTPUT_NAME=${2:-a.out} - -gfortran -I../LibDualzn64 -o "$OUTPUT_NAME" "$SOURCE_FILE" -L../LibDualzn64 -ldualzn diff --git a/GNU_Linux/gfortran_case/prec_real64/examples/ex1.f90 b/GNU_Linux/gfortran_case/prec_real64/examples/ex1.f90 deleted file mode 100644 index dbe9810..0000000 --- a/GNU_Linux/gfortran_case/prec_real64/examples/ex1.f90 +++ /dev/null @@ -1,35 +0,0 @@ -!gfortran -I../LibDualzn64 -o e1 ex1.f90 -L../LibDualzn64 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - - call set_order(5) !before anything we set the 'order' to work with - - z0 = (1.1_prec,2.2_prec) - r = 0 !we initialize the dual number to 0 - - !we set the 0-th and 1-th components. If dual numbers are used to - !calculate D^n f(z0) then r must be of the form r = r0 + 1*eps_1 - r%f(0) = z0 - r%f(1) = 1 - - fval = sin(r)**log(r*r) - - !writing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - print*,fval%f(k) - end do - deallocate(r%f) - deallocate(fval%f) -end program main - - - - diff --git a/GNU_Linux/gfortran_case/prec_real64/examples/ex2.f90 b/GNU_Linux/gfortran_case/prec_real64/examples/ex2.f90 deleted file mode 100644 index 2d86ba2..0000000 --- a/GNU_Linux/gfortran_case/prec_real64/examples/ex2.f90 +++ /dev/null @@ -1,53 +0,0 @@ -!gfortran -I../LibDualzn64 -o e2 ex2.f90 -L../LibDualzn64 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - real :: t1,t2 - - call set_order(5) !we set the order to work with - - !since a dualzn numbers is an allocatable entity, do not forget to - !initialize it - r = 0 !<--- initializing r to 0 - r%f(0) = (1.1_prec,0.0_prec) - r%f(1) = 1 !since we want to differentiate, r = r0 +1*eps_1 - !all the other components are 0 as r was initialized to 0 - - call cpu_time(t1) - fval = ftest(r) - call cpu_time(t2) - - !Computing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - write(*,"(i0,a,f0.1,a,e17.10)") k,"-th derivative at x = ", & - real(r%f(0)),":",real(fval%f(k)) - end do - - print*,"elapsed time (s):",t2-t1 - - deallocate(r%f) - deallocate(fval%f) - -contains - function ftest(x) result(fr) - type(dualzn), intent(in) :: x - type(dualzn) :: fr - integer :: k - - !nested function f(x) = sin(x) * exp(-x^2), f(f(...(f(x))...)) - !applied 1000 times - fr = sin(x)*exp(-x*x) - do k=1, 1000-1 - fr = sin(fr)*exp(-fr*fr) - end do - end function ftest -end program main - - diff --git a/GNU_Linux/gfortran_case/prec_real64/examples/ex3.f90 b/GNU_Linux/gfortran_case/prec_real64/examples/ex3.f90 deleted file mode 100644 index 73fdc70..0000000 --- a/GNU_Linux/gfortran_case/prec_real64/examples/ex3.f90 +++ /dev/null @@ -1,90 +0,0 @@ -!gfortran -I../LibDualzn64 -o e3 ex3.f90 -L../LibDualzn64 -ldualzn - -!module with example of functions -module function_mod - use dualzn_mod - implicit none - private - - public :: fstest, fvectest - -contains - !Example of scalar function f = sin(x*y*z) + cos(x*y*z) - function fstest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn) :: fr - type(dualzn) :: x,y,z - - x = r(1); y = r(2); z = r(3) - fr = sin(x*y*z) + cos(x*y*z) - end function fstest - - !Example of vector function f = [f1,f2,f3] - !f = fvectest(r) is a function f:D^m ---> Dn - function fvectest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn), allocatable, dimension(:) :: fr - type(dualzn) :: f1,f2,f3 - type(dualzn) :: x,y,z,w - - x = r(1); y = r(2); z = r(3); w = r(4) - - f1 = sin(x*y*z*w) - f2 = cos(x*y*z*w)*sqrt(w/y - x/z) - f3 = sin(log(x*y*z*w)) - - allocate(fr(3)) - fr = [f1,f2,f3] - end function fvectest -end module function_mod - -!----------------------------------------------------------------------- -!main program -program main - use precision_mod - use dualzn_mod - use diff_mod - use function_mod - implicit none - - integer, parameter :: nf =3, mq = 4 - complex(prec), parameter :: ii = (0,1) - complex(prec), dimension(mq) :: q, vec - complex(prec), dimension(nf,mq) :: Jmat - complex(prec), dimension(nf) :: JV - complex(prec), dimension(3) :: GV - complex(prec), dimension(3,3) :: Hmat - integer :: i - - vec = [1.0_prec,2.0_prec,3.0_prec,4.0_prec] - q = vec/10.0_prec + ii - - print*,"Jv using matmul" - Jmat = Jacobian(fvectest, q , nf) - JV = matmul(Jmat,vec) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"Jv using vector directional derivative" - JV = d1fvector(fvectest,vec,q,nf) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"---Hessian matrix---" - Hmat = Hessian(fstest,q(1:3)) - do i=1,3 - write(*,"(A,i0)") "row:",i - write(*,*) Hmat(i,:) - end do - write(*,*) - - print*,"---Gradient---" - GV = gradient(fstest,q(1:3)) - do i=1,3 - write(*,*) GV(i) - end do -end program main diff --git a/GNU_Linux/ifx_case/prec_real128/LibDualzn128/diff_mod.mod b/GNU_Linux/ifx_case/prec_real128/LibDualzn128/diff_mod.mod deleted file mode 100644 index fd1b25e..0000000 Binary files a/GNU_Linux/ifx_case/prec_real128/LibDualzn128/diff_mod.mod and /dev/null differ diff --git a/GNU_Linux/ifx_case/prec_real128/LibDualzn128/dualzn_mod.mod b/GNU_Linux/ifx_case/prec_real128/LibDualzn128/dualzn_mod.mod deleted file mode 100644 index 7962cce..0000000 Binary files a/GNU_Linux/ifx_case/prec_real128/LibDualzn128/dualzn_mod.mod and /dev/null differ diff --git a/GNU_Linux/ifx_case/prec_real128/LibDualzn128/libdualzn.a b/GNU_Linux/ifx_case/prec_real128/LibDualzn128/libdualzn.a deleted file mode 100644 index e8780af..0000000 Binary files a/GNU_Linux/ifx_case/prec_real128/LibDualzn128/libdualzn.a and /dev/null differ diff --git a/GNU_Linux/ifx_case/prec_real128/LibDualzn128/precision_mod.mod b/GNU_Linux/ifx_case/prec_real128/LibDualzn128/precision_mod.mod deleted file mode 100644 index 92d26db..0000000 Binary files a/GNU_Linux/ifx_case/prec_real128/LibDualzn128/precision_mod.mod and /dev/null differ diff --git a/GNU_Linux/ifx_case/prec_real128/examples/create_exec.sh b/GNU_Linux/ifx_case/prec_real128/examples/create_exec.sh deleted file mode 100755 index a322a79..0000000 --- a/GNU_Linux/ifx_case/prec_real128/examples/create_exec.sh +++ /dev/null @@ -1,13 +0,0 @@ -#!/bin/bash - -#Check if at least one argument (the source file) is provided -if [ "$#" -lt 1 ]; then - echo "Usage: $0 [output_executable_name]" - exit 1 -fi - -#Use the provided source file and output executable name or default to "a.out" -SOURCE_FILE=$1 -OUTPUT_NAME=${2:-a.out} - -ifx -I../LibDualzn128 -o "$OUTPUT_NAME" "$SOURCE_FILE" -L../LibDualzn128 -ldualzn diff --git a/GNU_Linux/ifx_case/prec_real128/examples/ex1.f90 b/GNU_Linux/ifx_case/prec_real128/examples/ex1.f90 deleted file mode 100644 index e46e527..0000000 --- a/GNU_Linux/ifx_case/prec_real128/examples/ex1.f90 +++ /dev/null @@ -1,36 +0,0 @@ -!ifx -I../LibDualzn128 -o e1 ex1.f90 -L../LibDualzn128 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - - call set_order(5) !before anything we set the 'order' to work with - - z0 = (1.1_prec,2.2_prec) - r = 0 !we initialize the dual number to 0 - - !we set the 0-th and 1-th components. If dual numbers are used to - !calculate D^n f(z0) then r must be of the form r = r0 + 1*eps_1 - r%f(0) = z0 - r%f(1) = 1 - - fval = sin(r)**log(r*r) - - !writing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - print*,fval%f(k) - end do - - deallocate(r%f) - deallocate(fval%f) -end program main - - - - diff --git a/GNU_Linux/ifx_case/prec_real128/examples/ex2.f90 b/GNU_Linux/ifx_case/prec_real128/examples/ex2.f90 deleted file mode 100644 index d5fce47..0000000 --- a/GNU_Linux/ifx_case/prec_real128/examples/ex2.f90 +++ /dev/null @@ -1,52 +0,0 @@ -!ifx -I../LibDualzn128 -o e2 ex2.f90 -L../LibDualzn128 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - real :: t1,t2 - - call set_order(5) !we set the order to work with - - !since a dualzn numbers is an allocatable entity, do not forget to - !initialize it - r = 0 !<--- initializing r to 0 - r%f(0) = (1.1_prec,0.0_prec) - r%f(1) = 1 !since we want to differentiate, r = r0 +1*eps_1 - !all the other components are 0 as r was initialized to 0 - - call cpu_time(t1) - fval = ftest(r) - call cpu_time(t2) - - !Computing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - write(*,"(i0,a,f0.1,a,e17.10)") k,"-th derivative at x = ", & - real(r%f(0)),":",real(fval%f(k)) - end do - - print*,"elapsed time (s):",t2-t1 - deallocate(r%f) - deallocate(fval%f) - -contains - function ftest(x) result(fr) - type(dualzn), intent(in) :: x - type(dualzn) :: fr - integer :: k - - !nested function f(x) = sin(x) * exp(-x^2), f(f(...(f(x))...)) - !applied 1000 times - fr = sin(x)*exp(-x*x) - do k=1, 1000-1 - fr = sin(fr)*exp(-fr*fr) - end do - end function ftest -end program main - - diff --git a/GNU_Linux/ifx_case/prec_real128/examples/ex3.f90 b/GNU_Linux/ifx_case/prec_real128/examples/ex3.f90 deleted file mode 100644 index d1fe232..0000000 --- a/GNU_Linux/ifx_case/prec_real128/examples/ex3.f90 +++ /dev/null @@ -1,90 +0,0 @@ -!ifx -I../LibDualzn128 -o e3 ex3.f90 -L../LibDualzn128 -ldualzn - -!module with example of functions -module function_mod - use dualzn_mod - implicit none - private - - public :: fstest, fvectest - -contains - !Example of scalar function f = sin(x*y*z) + cos(x*y*z) - function fstest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn) :: fr - type(dualzn) :: x,y,z - - x = r(1); y = r(2); z = r(3) - fr = sin(x*y*z) + cos(x*y*z) - end function fstest - - !Example of vector function f = [f1,f2,f3] - !f = fvectest(r) is a function f:D^m ---> Dn - function fvectest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn), allocatable, dimension(:) :: fr - type(dualzn) :: f1,f2,f3 - type(dualzn) :: x,y,z,w - - x = r(1); y = r(2); z = r(3); w = r(4) - - f1 = sin(x*y*z*w) - f2 = cos(x*y*z*w)*sqrt(w/y - x/z) - f3 = sin(log(x*y*z*w)) - - allocate(fr(3)) - fr = [f1,f2,f3] - end function fvectest -end module function_mod - -!----------------------------------------------------------------------- -!main program -program main - use precision_mod - use dualzn_mod - use diff_mod - use function_mod - implicit none - - integer, parameter :: nf =3, mq = 4 - complex(prec), parameter :: ii = (0,1) - complex(prec), dimension(mq) :: q, vec - complex(prec), dimension(nf,mq) :: Jmat - complex(prec), dimension(nf) :: JV - complex(prec), dimension(3) :: GV - complex(prec), dimension(3,3) :: Hmat - integer :: i - - vec = [1.0_prec,2.0_prec,3.0_prec,4.0_prec] - q = vec/10.0_prec + ii - - print*,"Jv using matmul" - Jmat = Jacobian(fvectest, q , nf) - JV = matmul(Jmat,vec) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"Jv using vector directional derivative" - JV = d1fvector(fvectest,vec,q,nf) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"---Hessian matrix---" - Hmat = Hessian(fstest,q(1:3)) - do i=1,3 - write(*,"(A,i0)") "row:",i - write(*,*) Hmat(i,:) - end do - write(*,*) - - print*,"---Gradient---" - GV = gradient(fstest,q(1:3)) - do i=1,3 - write(*,*) GV(i) - end do -end program main diff --git a/GNU_Linux/ifx_case/prec_real64/LibDualzn64/diff_mod.mod b/GNU_Linux/ifx_case/prec_real64/LibDualzn64/diff_mod.mod deleted file mode 100644 index c141158..0000000 Binary files a/GNU_Linux/ifx_case/prec_real64/LibDualzn64/diff_mod.mod and /dev/null differ diff --git a/GNU_Linux/ifx_case/prec_real64/LibDualzn64/dualzn_mod.mod b/GNU_Linux/ifx_case/prec_real64/LibDualzn64/dualzn_mod.mod deleted file mode 100644 index 149bf37..0000000 Binary files a/GNU_Linux/ifx_case/prec_real64/LibDualzn64/dualzn_mod.mod and /dev/null differ diff --git a/GNU_Linux/ifx_case/prec_real64/LibDualzn64/libdualzn.a b/GNU_Linux/ifx_case/prec_real64/LibDualzn64/libdualzn.a deleted file mode 100644 index 1305204..0000000 Binary files a/GNU_Linux/ifx_case/prec_real64/LibDualzn64/libdualzn.a and /dev/null differ diff --git a/GNU_Linux/ifx_case/prec_real64/LibDualzn64/precision_mod.mod b/GNU_Linux/ifx_case/prec_real64/LibDualzn64/precision_mod.mod deleted file mode 100644 index 8038966..0000000 Binary files a/GNU_Linux/ifx_case/prec_real64/LibDualzn64/precision_mod.mod and /dev/null differ diff --git a/GNU_Linux/ifx_case/prec_real64/examples/create_exec.sh b/GNU_Linux/ifx_case/prec_real64/examples/create_exec.sh deleted file mode 100755 index 13113b7..0000000 --- a/GNU_Linux/ifx_case/prec_real64/examples/create_exec.sh +++ /dev/null @@ -1,13 +0,0 @@ -#!/bin/bash - -#Check if at least one argument (the source file) is provided -if [ "$#" -lt 1 ]; then - echo "Usage: $0 [output_executable_name]" - exit 1 -fi - -#Use the provided source file and output executable name or default to "a.out" -SOURCE_FILE=$1 -OUTPUT_NAME=${2:-a.out} - -ifx -I../LibDualzn64 -o "$OUTPUT_NAME" "$SOURCE_FILE" -L../LibDualzn64 -ldualzn diff --git a/GNU_Linux/ifx_case/prec_real64/examples/ex1.f90 b/GNU_Linux/ifx_case/prec_real64/examples/ex1.f90 deleted file mode 100644 index 55005f6..0000000 --- a/GNU_Linux/ifx_case/prec_real64/examples/ex1.f90 +++ /dev/null @@ -1,37 +0,0 @@ -!ifx -I../LibDualzn64 -o e1 ex1.f90 -L../LibDualzn64 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - - call set_order(5) !before anything we set the 'order' to work with - - z0 = (1.1_prec,2.2_prec) - r = 0 !we initialize the dual number to 0 - - !we set the 0-th and 1-th components. If dual numbers are used to - !calculate D^n f(z0) then r must be of the form r = r0 + 1*eps_1 - r%f(0) = z0 - r%f(1) = 1 - - fval = sin(r)**log(r*r) - - !writing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - print*,fval%f(k) - end do - - - deallocate(r%f) - deallocate(fval%f) -end program main - - - - diff --git a/GNU_Linux/ifx_case/prec_real64/examples/ex2.f90 b/GNU_Linux/ifx_case/prec_real64/examples/ex2.f90 deleted file mode 100644 index 41e2177..0000000 --- a/GNU_Linux/ifx_case/prec_real64/examples/ex2.f90 +++ /dev/null @@ -1,52 +0,0 @@ -!ifx -I../LibDualzn64 -o e2 ex2.f90 -L../LibDualzn64 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - real :: t1,t2 - - call set_order(15) !we set the order to work with - - !since a dualzn numbers is an allocatable entity, do not forget to - !initialize it - r = 0 !<--- initializing r to 0 - r%f(0) = (1.1_prec,0.0_prec) - r%f(1) = 1 !since we want to differentiate, r = r0 +1*eps_1 - !all the other components are 0 as r was initialized to 0 - - call cpu_time(t1) - fval = ftest(r) - call cpu_time(t2) - - !Computing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - write(*,"(i0,a,f0.1,a,e17.10)") k,"-th derivative at x = ", & - real(r%f(0)),":",real(fval%f(k)) - end do - - print*,"elapsed time (s):",t2-t1 - deallocate(r%f) - deallocate(fval%f) - -contains - function ftest(x) result(fr) - type(dualzn), intent(in) :: x - type(dualzn) :: fr - integer :: k - - !nested function f(x) = sin(x) * exp(-x^2), f(f(...(f(x))...)) - !applied 1000 times - fr = sin(x)*exp(-x*x) - do k=1, 1000-1 - fr = sin(fr)*exp(-fr*fr) - end do - end function ftest -end program main - - diff --git a/README.md b/README.md index d504aed..5855847 100644 --- a/README.md +++ b/README.md @@ -2,37 +2,18 @@ **DNAOAD** (Dual Number Arbitrary Order Automatic Differentiation) is a Fortran implementation of dual numbers for arbitrary order automatic differentiation. -## ๐Ÿงช Running Examples +## ๐Ÿงช Building & Running Examples -To run the examples, navigate to the `examples` directory and use a terminal. For example: +To build and run the examples we recommend installing the fortran package manger ([fpm](https://fpm.fortran-lang.org/index.html)) first. Then the examples can by executed by simply running, for example -```sh -sh create_exec.sh ex1.f90 +```bash +fpm run --example ex1 ``` -On **Windows**, use: - -```bat -create_exec.bat ex1.f90 -``` - -> **Note:** Alternatively (and this is the recommended method), you can manually compile and run the source files. - -## โš™๏ธ Manual Compilation - -The examples in the `SourcesF` folder can be compiled directly using the **Intel Fortran** compiler: - -```sh -ifx ex3.f90 -``` - -If you're using **gfortran**, use the `SourcesF_gfortran` folder instead. +from the command line. ## ๐Ÿ“– More Information -For further details, see **Section 4** of the article: +For further details, see the article: ๐Ÿ“Ž [https://doi.org/10.48550/arXiv.2501.04159](https://doi.org/10.48550/arXiv.2501.04159) - - - diff --git a/SourcesF/diff_mod.f90 b/SourcesF/diff_mod.f90 deleted file mode 100644 index 58883ed..0000000 --- a/SourcesF/diff_mod.f90 +++ /dev/null @@ -1,178 +0,0 @@ -module diff_mod - use precision_mod !Depending on how to link, it will set prec to be - !real64 or real128 - use dualzn_mod - implicit none - - private - public :: d1fscalar, d2fscalar, d1fvector - public :: Jacobian, Hessian, gradient - - !Abstract interfaces for functions passed as arguments to other - !functions. - ! - !Interface for a scalar dual function f: D^m --> D (similar to - !f: R^m --> R) - abstract interface - function fsdual(xd) result(frsd) - use dualzn_mod - type(dualzn), intent(in), dimension(:) :: xd - type(dualzn) :: frsd - end function fsdual - end interface - - !Interface for a vector dual function f: D^m --> D^n - interface - function fvecdual(xd) result(frd) - use dualzn_mod - type(dualzn), intent(in), dimension(:) :: xd - type(dualzn), allocatable, dimension(:) :: frd - end function fvecdual - end interface - - interface d2fscalar - module procedure d2fscalarvv - module procedure d2fscalaruv - end interface d2fscalar - !--------------------------------------------------------------------- - -contains - !Hessian operator - function Hessian(fsd,qcmplx) result(Hmat) - procedure(fsdual) :: fsd - complex(prec), intent(in), dimension(:) :: qcmplx - complex(prec), dimension(size(qcmplx),size(qcmplx)) :: Hmat - - complex(prec), dimension(size(qcmplx)) :: ei, ej - integer :: i,j,m - - m = size(qcmplx) - - Hmat = 0 - do i=1,m - ei = 0 - ei(i) = 1 - Hmat(i,i) = d2fscalarvv(fsd,ei,qcmplx) - end do - - do i=1, m - ei = 0 - ei(i) = 1 - do j=i+1,m - ej = 0 - ej(j) = 1 - Hmat(i,j) = 0.5_prec * (d2fscalar(fsd,ei + ej,qcmplx) - Hmat(i,i) & - Hmat(j,j)) - Hmat(j,i) = Hmat(i,j) - end do - end do - end function Hessian - - !d2fscalaruv(f,u,v,q) gives the second-order directional derivative - !of the scalar function f: D^m --> D, where f = f(q), along vectors u - !and v, evaluated at point q - function d2fscalaruv(fsd,x,y,q) result(fr) - procedure(fsdual) :: fsd - complex(prec), intent(in), dimension(:) :: x, y, q - complex(prec) :: fr - - if(all(x == y)) then - fr = d2fscalarvv(fsd,x,q) - else - fr = 0.5_prec*(d2fscalarvv(fsd,x + y,q) - d2fscalarvv(fsd,x,q) -& - d2fscalarvv(fsd,y,q)) - end if - end function d2fscalaruv - - !Second-order directional derivative of a scalar function along vector - !v, evaluated at point q - function d2fscalarvv(fsd,v,q) result(fr) - procedure(fsdual) :: fsd - complex(prec), intent(in), dimension(:) :: v, q - complex(prec) :: fr - type(dualzn) :: eps1 - integer :: original_order - - original_order = get_order() - call set_order(2) - eps1 = 0 - eps1%f(1) = 1 - - fr = f_part(fsd(q + eps1*v),2) - call set_order(original_order) - end function d2fscalarvv - - !Jacobian operator: To optimize efficiency, we include the parameter - !'n', representing the dimension of 'fvecd'. - function Jacobian(fvecd,qcmplx,n) result(Jmat) - procedure(fvecdual) :: fvecd - complex(prec), intent(in), dimension(:) :: qcmplx - integer, intent(in) :: n - complex(prec), dimension(n,size(qcmplx)) :: Jmat - - complex(prec), dimension(size(qcmplx),n) :: TrpJmat - complex(prec), dimension(size(qcmplx)) :: ei - integer :: i - - do i = 1,size(qcmplx) - ei = 0 - ei(i) = 1 - TrpJmat(i,:) = d1fvector(fvecd,ei,qcmplx,n) - end do - - Jmat = transpose(TrpJmat) - end function Jacobian - - !First-order directional derivative of a vector function along vector - !v, evaluated at point q. To optimize efficiency, we include the - !parameter 'n', representing the dimension of 'fvecd'. - ! - function d1fvector(fvecd,v,q,n) result(fr) - procedure(fvecdual) :: fvecd - complex(prec), intent(in), dimension(:) :: v, q - integer, intent(in) :: n - complex(prec), dimension(n) :: fr - type(dualzn) :: eps1 - integer :: original_order - - original_order = get_order() - call set_order(1) - eps1 = 0 - eps1%f(1) = 1 - - fr = f_part(fvecd(q + eps1*v),1) - call set_order(original_order) - end function d1fvector - - function gradient(fsd,q) result(fr) - procedure(fsdual) :: fsd - complex(prec), intent(in), dimension(:) :: q - complex(prec), dimension(size(q)) :: fr - complex(prec), dimension(size(q)) :: ei - integer :: i - - do i = 1,size(q) - ei = 0 - ei(i) = 1 - fr(i) = d1fscalar(fsd,ei,q) - end do - end function gradient - - !First-order directional derivative of a scalar function along vector - !v, evaluated at point q - function d1fscalar(fsd,v,q) result(fr) - procedure(fsdual) :: fsd - complex(prec), intent(in), dimension(:) :: v, q - complex(prec) :: fr - type(dualzn) :: eps1 - integer :: original_order - - original_order = get_order() - call set_order(1) - eps1 = 0 - eps1%f(1) = 1 - - fr = f_part(fsd(q + eps1*v),1) - call set_order(original_order) - end function d1fscalar -end module diff_mod diff --git a/SourcesF/ex1.f90 b/SourcesF/ex1.f90 deleted file mode 100644 index 9e60eaf..0000000 --- a/SourcesF/ex1.f90 +++ /dev/null @@ -1,42 +0,0 @@ -include './precision_mod.f90' -include './dualzn_mod.f90' - -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k, local_order - - - call set_order(5) !before anything we set the 'order' to work with - - local_order = get_order() - - z0 = (1.1_prec,2.2_prec) - r = 0 !we initialize the dual number to 0 - - !we set the 0-th and 1-th components. If dual numbers are used to - !calculate D^n f(z0) then r must be of the form r = r0 + 1*eps_1 - r%f(0) = z0 - r%f(1) = 1 - - fval = sin(r)**log(r*r) - - !writing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0, local_order - print*,fval%f(k) - end do - - deallocate(r%f) - deallocate(fval%f) - -end program main - - - - diff --git a/SourcesF/ex2.f90 b/SourcesF/ex2.f90 deleted file mode 100644 index c14b7b7..0000000 --- a/SourcesF/ex2.f90 +++ /dev/null @@ -1,54 +0,0 @@ -include './precision_mod.f90' -include './dualzn_mod.f90' - -program main - use precision_mod - use dualzn_mod - implicit none - - type(dualzn) :: r, fval - integer :: k, local_order - real :: t1,t2 - - - call set_order(15) !we set the order to work with - local_order = get_order() - - !since a dualzn numbers is an allocatable entity, do not forget to - !initialize it - r = 0 !<--- initializing r to 0 - r%f(0) = (1.1_prec,0.0_prec) - r%f(1) = 1 !since we want to differentiate, r = r0 +1*eps_1 - !all the other components are 0 as r was initialized to 0 - - call cpu_time(t1) - fval = ftest(r) - call cpu_time(t2) - - !Computing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0, local_order - write(*,"(i0,a,f0.1,a,e17.10)") k,"-th derivative at x = ", & - real(r%f(0)),":",real(fval%f(k)) - end do - - print*,"elapsed time (s):",t2-t1 - - deallocate(r%f) - deallocate(fval%f) - -contains - function ftest(x) result(fr) - type(dualzn), intent(in) :: x - type(dualzn) :: fr - integer :: k - - !nested function f(x) = sin(x) * exp(-x^2), f(f(...(f(x))...)) - !applied 1000 times - fr = sin(x)*exp(-x*x) - do k=1, 1000-1 - fr = sin(fr)*exp(-fr*fr) - end do - end function ftest -end program main diff --git a/SourcesF/ex3.f90 b/SourcesF/ex3.f90 deleted file mode 100644 index 36d87d5..0000000 --- a/SourcesF/ex3.f90 +++ /dev/null @@ -1,92 +0,0 @@ -include './precision_mod.f90' -include './dualzn_mod.f90' -include './diff_mod.f90' - -!module with example of functions -module function_mod - use dualzn_mod - implicit none - private - - public :: fstest, fvectest - -contains - !Example of scalar function f = sin(x*y*z) + cos(x*y*z) - function fstest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn) :: fr - type(dualzn) :: x,y,z - - x = r(1); y = r(2); z = r(3) - fr = sin(x*y*z) + cos(x*y*z) - end function fstest - - !Example of vector function f = [f1,f2,f3] - !f = fvectest(r) is a function f:D^m ---> Dn - function fvectest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn), allocatable, dimension(:) :: fr - type(dualzn) :: f1,f2,f3 - type(dualzn) :: x,y,z,w - - x = r(1); y = r(2); z = r(3); w = r(4) - - f1 = sin(x*y*z*w) - f2 = cos(x*y*z*w)*sqrt(w/y - x/z) - f3 = sin(log(x*y*z*w)) - - allocate(fr(3)) - fr = [f1,f2,f3] - end function fvectest -end module function_mod - -!----------------------------------------------------------------------- -!main program -program main - use precision_mod - use dualzn_mod - use diff_mod - use function_mod - implicit none - - integer, parameter :: nf =3, mq = 4 - complex(prec), parameter :: ii = (0,1) - complex(prec), dimension(mq) :: q, vec - complex(prec), dimension(nf,mq) :: Jmat - complex(prec), dimension(nf) :: JV - complex(prec), dimension(3) :: GV - complex(prec), dimension(3,3) :: Hmat - integer :: i - - vec = [1.0_prec,2.0_prec,3.0_prec,4.0_prec] - q = vec/10.0_prec + ii - - print*,"Jv using matmul" - Jmat = Jacobian(fvectest, q , nf) - JV = matmul(Jmat,vec) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"Jv using vector directional derivative" - JV = d1fvector(fvectest,vec,q,nf) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"---Hessian matrix---" - Hmat = Hessian(fstest,q(1:3)) - do i=1,3 - write(*,"(A,i0)") "row:",i - write(*,*) Hmat(i,:) - end do - write(*,*) - - print*,"---Gradient---" - GV = gradient(fstest,q(1:3)) - do i=1,3 - write(*,*) GV(i) - end do -end program main diff --git a/SourcesF_gfortran/Readme.txt b/SourcesF_gfortran/Readme.txt deleted file mode 100644 index 3d09cc5..0000000 --- a/SourcesF_gfortran/Readme.txt +++ /dev/null @@ -1,2 +0,0 @@ -The addition of dualzn numbers is implemented as an elemental operation, which allows expressions like xvec + yvec for dualzn vectors. However, the gfortran compiler produces errors when such expressions are passed directly as arguments to a function (e.g., f(x + y)). In such cases, gfortran requires the result to be stored in an intermediate variable before being passed. To avoid this issue, a dedicated function for adding dual vectors is defined in the diff_mod.f90 file. - diff --git a/SourcesF_gfortran/dualzn_mod.f90 b/SourcesF_gfortran/dualzn_mod.f90 deleted file mode 100644 index b2259b5..0000000 --- a/SourcesF_gfortran/dualzn_mod.f90 +++ /dev/null @@ -1,1816 +0,0 @@ -!Dual numbers to order n of complex components -!This can be used to compute first, second,... nth order derivatives -!F. Pe~nu~nuri -!UADY, Merida Yucatan Mexico -!2024 -! - -module dualzn_mod - use precision_mod - implicit none - - private - !--------------------------------------------------------------------- - !Some Module variables - !Default order, can be modified with set_order - integer, protected :: order = 1 - real(prec), public, parameter :: Pi = 4.0_prec*atan(1.0_prec) - !--------------------------------------------------------------------- - - !dual number definition to any order - type, public :: dualzn - complex(prec), allocatable, dimension(:) :: f - end type dualzn - - public :: set_order, get_order, initialize_dualzn, f_part - public :: binomial, BellY, Dnd - public :: Mset_fpart - public :: xto_dzn - public :: itodn, realtodn, cmplxtodn - - public :: inv, sin, cos, tan, exp, log, sqrt, asin, acos, atan, asinh - public :: acosh, atanh, sinh, cosh, tanh, absx, atan2 - public :: conjg, sum, product, matmul - - public :: assignment (=) - public :: operator(*), operator(/), operator(+), operator(-) - public :: operator(**), operator(==), operator(/=) - - !equal assignment - interface assignment (=) - module procedure igualz16_dzn !dualzn <--- complex16 - module procedure igualz8_dzn !dualzn <--- complex8 - module procedure igualz4_dzn !dualzn <--- complex4 - module procedure igualr16_dzn !dualzn <--- real16 - module procedure igualr8_dzn !dualzn <--- real8 - module procedure igualr4_dzn !dualzn <--- real4 - module procedure iguali8_dzn !dualzn <--- integer8 - module procedure iguali4_dzn !dualzn <--- integer4 - end interface assignment (=) - !--------------------------------------------------------------------- - - !functions to change from integer, real and complex numbers to dual - ! can use xto_dzn instead - interface itodn - module procedure i4todn - module procedure i8todn - end interface itodn - - interface realtodn - module procedure r4todn - module procedure r8todn - module procedure r16todn - end interface realtodn - - interface cmplxtodn - module procedure c4todn - module procedure c8todn - module procedure c16todn - end interface cmplxtodn - !--------------------------------------------------------------------- - - !Logical equal operator - interface operator (==) - module procedure eq_dzn - end interface operator (==) - !--------------------------------------------------------------------- - - !Logical not equal operator - interface operator (/=) - module procedure noteq_dzn - end interface operator (/=) - !--------------------------------------------------------------------- - - !overloaded operators - interface operator(*) - module procedure timesdX - module procedure timesXd - end interface operator(*) - - interface operator(/) - module procedure divdX - module procedure divXd - end interface operator(/) - - interface operator (+) - module procedure masd !unary - module procedure sumadX - module procedure sumaXd - end interface operator (+) - - interface operator (-) - module procedure menosd !unary - module procedure restadX - module procedure restaXd - end interface operator (-) - - interface operator(**) - module procedure powerdX - module procedure powerXd - end interface operator(**) - !--------------------------------------------------------------------- - !some matrix and vector operations - interface sum - module procedure sumR2dzn !sum(dual_rank2,dir) - module procedure sumR20dzn !sum(dual_rank2) - module procedure sumR1dzn !sum(dual_rank1) - end interface sum - !--------------------------------------------------------------------- - - !product for dual - interface product - module procedure prodR2dzn !product(dual_rank2,dir) - module procedure prodR20dzn !product(dual_rank2) - module procedure prodR1dzn !product(dual_rank1) - end interface product - !--------------------------------------------------------------------- - - interface matmul - module procedure MtimesdX - module procedure Mtimesc128d - module procedure Mtimesc64d - module procedure Mtimesc32d - module procedure Mtimesr128d - module procedure Mtimesr64d - module procedure Mtimesr32d - module procedure Mtimesi64d - module procedure Mtimesi32d - end interface matmul - !--------------------------------------------------------------------- - - !overloaded functions - interface conjg - module procedure conjg_dzn - end interface conjg - - interface atan2 - module procedure atan2d_ - end interface atan2 - - interface tanh - module procedure tanhd - end interface tanh - - interface cosh - module procedure coshd - end interface cosh - - interface sinh - module procedure sinhd - end interface sinh - - interface atanh - module procedure atanhd - end interface atanh - - interface acosh - module procedure acoshd - end interface acosh - - interface asinh - module procedure asinhd - end interface asinh - - interface atan - module procedure atand_ - ! Fortran 2008 standard for 2 argument atan: - module procedure atan2d_ - end interface atan - - interface acos - module procedure acosd_ - end interface acos - - interface asin - module procedure asind_ - end interface asin - - interface log - module procedure logd - end interface log - - interface exp - module procedure expd - end interface exp - - interface sqrt - module procedure sqrtd - end interface sqrt - - interface sin - module procedure sind_ - end interface sin - - interface cos - module procedure cosd_ - end interface cos - - interface tan - module procedure tand_ - end interface tan - !--------------------------------------------------------------------- - - !interface to define a function of complex variable returning a dual - !variable - !f:z --> dual - abstract interface - pure function funzdual(z_val) result(f_result) - use precision_mod - import :: dualzn - complex(prec), intent(in) :: z_val - type(dualzn) :: f_result - end function funzdual - end interface - !--------------------------------------------------------------------- - - !Functions -contains - !to set the order for dual numbers - subroutine set_order(new_order) - integer, intent(in) :: new_order - if(new_order<1) then - stop "order must be greater than 0" - end if - order = new_order - end subroutine set_order - - !check the order for dual numbers - pure function get_order() result(order_value) - integer :: order_value - order_value = order - end function get_order - !--------------------------------------------------------------------- - - !Assignment, equal operator - !dualzn <--- complex16 - elemental subroutine igualz16_dzn(A, z) - type(dualzn), intent(out) :: A - complex(real128), intent(in) :: z - - A = cmplxtodn(z) - end subroutine igualz16_dzn - - !dualzn <--- complex8 - elemental subroutine igualz8_dzn(A, z) - type(dualzn), intent(out) :: A - complex(real64), intent(in) :: z - - A = cmplxtodn(z) - end subroutine igualz8_dzn - - !dualzn <--- complex4 - elemental subroutine igualz4_dzn(A, z) - type(dualzn), intent(out) :: A - complex(real32), intent(in) :: z - - A = cmplxtodn(z) - end subroutine igualz4_dzn - !--------------------------------------------------------------------- - - !dualzn <--- real16 - elemental subroutine igualr16_dzn(A, x) - type(dualzn), intent(out) :: A - real(real128), intent(in) :: x - - A = realtodn(x) - end subroutine igualr16_dzn - - !dualzn <--- real8 - elemental subroutine igualr8_dzn(A, x) - type(dualzn), intent(out) :: A - real(real64), intent(in) :: x - - A = realtodn(x) - end subroutine igualr8_dzn - - !dualzn <--- real4 - elemental subroutine igualr4_dzn(A, x) - type(dualzn), intent(out) :: A - real(real32), intent(in) :: x - - A = realtodn(x) - end subroutine igualr4_dzn - !--------------------------------------------------------------------- - - !dualzn <--- integer8 - elemental subroutine iguali8_dzn(A, ix) - type(dualzn), intent(out) :: A - integer(int64), intent(in) :: ix - - A = itodn(ix) - end subroutine iguali8_dzn - - !dualzn <--- integer4 - elemental subroutine iguali4_dzn(A, ix) - type(dualzn), intent(out) :: A - integer(int32), intent(in) :: ix - - A = itodn(ix) - end subroutine iguali4_dzn - !--------------------------------------------------------------------- - - !to initialize a dual number to zero - elemental subroutine initialize_dualzn(zdn) - type(dualzn), intent(out) :: zdn - - allocate(zdn%f(0:order)) !notice, it will be order+1 components - zdn%f=(0.0_prec,0.0_prec) - end subroutine initialize_dualzn - !--------------------------------------------------------------------- - - !to extract the f%(k) part of a dual number or array of duals - elemental function f_part(x,k) result(fr) - type(dualzn), intent(in) :: x - integer, intent(in) :: k - complex(prec) :: fr - - fr = x%f(k) - end function f_part - !--------------------------------------------------------------------- - - !Logical equal operator - elemental function eq_dzn(lhs, rhs) result(fr) - type (dualzn), intent(in) :: lhs, rhs - logical :: fr - logical :: eqfk - integer :: k - - fr = .true. - do k=0,order - eqfk = lhs%f(k) == rhs%f(k) - if(.not. eqfk) then - fr = .false. - exit - end if - end do - end function eq_dzn - !--------------------------------------------------------------------- - !Logical not equal operator - elemental function noteq_dzn(lhs, rhs) result(f_res) - type (dualzn), intent(in) :: lhs, rhs - logical :: f_res - - f_res = .not.(lhs == rhs) - end function noteq_dzn - !--------------------------------------------------------------------- - - ! Converts a value (integer, real, complex, dualzn) `X` to a dualzn - ! number. - elemental function xto_dzn(X) result(fr) - class(*), intent(in) :: X - type(dualzn) :: fr - - call initialize_dualzn(fr) - - select type (X) - type is (dualzn) - fr = X - return - type is (complex(kind=real128)) - fr%f(0) = cmplx(X, kind=prec) - type is (complex(kind=real64)) - fr%f(0) = cmplx(X, kind=prec) - type is (complex(kind=real32)) - fr%f(0) = cmplx(X, kind=prec) - type is (real(kind=real128)) - fr%f(0) = cmplx(X, 0.0_prec, kind=prec) - type is (real(kind=real64)) - fr%f(0) = cmplx(X, 0.0_prec, kind=prec) - type is (real(kind=real32)) - fr%f(0) = cmplx(X, 0.0_prec, kind=prec) - type is (integer(kind=int64)) - fr%f(0) = cmplx(X, 0.0_prec, kind=prec) - type is (integer(kind=int32)) - fr%f(0) = cmplx(X, 0.0_prec, kind=prec) - class default - continue - end select - end function xto_dzn - - !the beolow functions are private, use xto_dzn instead - elemental function c16todn(x) result(fr) - complex(real128), intent(in) :: x - type(dualzn) :: fr - - call initialize_dualzn(fr) - fr%f(0)=cmplx(x,kind=prec) - end function c16todn - - elemental function c8todn(x) result(fr) - complex(real64), intent(in) :: x - type(dualzn) :: fr - - call initialize_dualzn(fr) - fr%f(0)=cmplx(x,kind=prec) - end function c8todn - - elemental function c4todn(x) result(fr) - complex(real32), intent(in) :: x - type(dualzn) :: fr - - call initialize_dualzn(fr) - fr%f(0)=cmplx(x,kind=prec) - end function c4todn - !--------------------------------------------------------------------- - - elemental function r16todn(x) result(fr) - real(real128), intent(in) :: x - type(dualzn) :: fr - - call initialize_dualzn(fr) - fr%f(0)=real(x,kind=prec) - end function r16todn - - elemental function r8todn(x) result(fr) - real(real64), intent(in) :: x - type(dualzn) :: fr - - call initialize_dualzn(fr) - fr%f(0)=real(x,kind=prec) - end function r8todn - - elemental function r4todn(x) result(fr) - real(real32), intent(in) :: x - type(dualzn) :: fr - - call initialize_dualzn(fr) - fr%f(0)=real(x,kind=prec) - end function r4todn - !--------------------------------------------------------------------- - - elemental function i4todn(ix) result(fr) - integer(int32), intent(in) :: ix - type(dualzn) :: fr - - call initialize_dualzn(fr) - fr%f(0)=real(ix,kind=prec) - end function i4todn - - elemental function i8todn(ix) result(fr) - integer(int64), intent(in) :: ix - type(dualzn) :: fr - - call initialize_dualzn(fr) - fr%f(0)=real(ix,kind=prec) - end function i8todn - !--------------------------------------------------------------------- - - !dual + class - elemental function sumadX(B,X) result(fr) - class(*), intent(in) :: X - type(dualzn), intent(in) :: B - type(dualzn) :: fr - type(dualzn) :: Xd - - select type (X) - type is (dualzn) - Xd = X - type is(complex(kind=real128)) - Xd = cmplxtodn(X) - type is (complex(kind=real64)) - Xd =cmplxtodn(X) - type is (complex(kind=real32)) - Xd =cmplxtodn(X) - type is (real(kind=real128)) - Xd = realtodn(X) - type is (real(kind=real64)) - Xd = realtodn(X) - type is (real(kind=real32)) - Xd = realtodn(X) - type is (integer(kind=int64)) - Xd = itodn(X) - type is (integer(kind=int32)) - Xd = itodn(X) - end select - fr = sumad(B,Xd) - end function sumadX - !--------------------------------------------------------------------- - - !class + dual - elemental function sumaXd(XX,BB) result(fr) - class(*), intent(in) :: XX - type(dualzn), intent(in) :: BB - type(dualzn) :: fr - type(dualzn) :: Xd - - select type (XX) - type is (dualzn) - Xd = XX - type is(complex(kind=real128)) - Xd = cmplxtodn(XX) - type is (complex(kind=real64)) - Xd =cmplxtodn(XX) - type is (complex(kind=real32)) - Xd =cmplxtodn(XX) - type is (real(kind=real128)) - Xd = realtodn(XX) - type is (real(kind=real64)) - Xd = realtodn(XX) - type is (real(kind=real32)) - Xd = realtodn(XX) - type is (integer(kind=int64)) - Xd = itodn(XX) - type is (integer(kind=int32)) - Xd = itodn(XX) - end select - fr = sumad(Xd, BB) - end function sumaXd - !--------------------------------------------------------------------- - - !A+B - elemental function sumad(A,B) result(fr) - type(dualzn), intent(in) :: A,B - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - do k=0,order - fr%f(k) = A%f(k) + B%f(k) - end do - end function sumad - !--------------------------------------------------------------------- - - !+dualzn (unary) - elemental function masd(A) result(fr) - type(dualzn), intent(in) :: A - type(dualzn) :: fr - - fr = A - end function masd - !--------------------------------------------------------------------- - - !dual - class - elemental function restadX(B,X) result(fr) - class(*), intent(in) :: X - type(dualzn), intent(in) :: B - type(dualzn) :: fr - type(dualzn) :: Xd - - select type (X) - type is (dualzn) - Xd = X - type is(complex(kind=real128)) - Xd = cmplxtodn(X) - type is (complex(kind=real64)) - Xd =cmplxtodn(X) - type is (complex(kind=real32)) - Xd =cmplxtodn(X) - type is (real(kind=real128)) - Xd = realtodn(X) - type is (real(kind=real64)) - Xd = realtodn(X) - type is (real(kind=real32)) - Xd = realtodn(X) - type is (integer(kind=int64)) - Xd = itodn(X) - type is (integer(kind=int32)) - Xd = itodn(X) - end select - fr = restad(B,Xd) - end function restadX - !--------------------------------------------------------------------- - - !class - dual - elemental function restaXd(XX,BB) result(fr) - class(*), intent(in) :: XX - type(dualzn), intent(in) :: BB - type(dualzn) :: fr - type(dualzn) :: Xd - - select type (XX) - type is (dualzn) - Xd = XX - type is(complex(kind=real128)) - Xd = cmplxtodn(XX) - type is (complex(kind=real64)) - Xd =cmplxtodn(XX) - type is (complex(kind=real32)) - Xd =cmplxtodn(XX) - type is (real(kind=real128)) - Xd = realtodn(XX) - type is (real(kind=real64)) - Xd = realtodn(XX) - type is (real(kind=real32)) - Xd = realtodn(XX) - type is (integer(kind=int64)) - Xd = itodn(XX) - type is (integer(kind=int32)) - Xd = itodn(XX) - end select - fr = restad(Xd, BB) - end function restaXd - !--------------------------------------------------------------------- - - !-dualzn (unary) - elemental function menosd(A) result(fr) - type(dualzn), intent(in) :: A - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - do k=0,order - fr%f(k) = -A%f(k) - end do - end function menosd - - !A-B - elemental function restad(A,B) result(fr) - type(dualzn), intent(in) :: A,B - type(dualzn) :: fr - - fr = -B+A - end function restad - !--------------------------------------------------------------------- - - !dual**class - elemental function powerdX(B, X) result(fr) - class(*), intent(in) :: X - type(dualzn), intent(in) :: B - type(dualzn) :: fr - - select type (X) - type is (dualzn) - fr = powerd(B, X) - type is (integer(kind=int64)) - fr = power_dint8(B, X) - type is (integer(kind=int32)) - fr = power_dint(B, X) - class default - fr = powerd(B, xto_dzn(X)) - end select - end function powerdX - !--------------------------------------------------------------------- - - !class**dual (X**B) - elemental function powerXd(XX,BB) result(fr) - class(*), intent(in) :: XX - type(dualzn), intent(in) :: BB - type(dualzn) :: fr - - select type (XX) - type is (dualzn) - fr = powerd(XX, BB) - class default - fr = powerd(xto_dzn(XX),BB) - end select - end function powerXd -!--------------------------------------------------------------------- - - function MtimesdX(A,X) result(fr) - type(dualzn), intent(in), dimension(:,:) :: A - class(*), intent(in), dimension(:,:) :: X - type(dualzn), dimension(size(A,1),size(X,2)) :: fr - type(dualzn), dimension(size(X,1),size(X,2)) :: Xd - - select type (X) - type is (dualzn) - Xd = X - type is(complex(kind=real128)) - Xd = X - type is(complex(kind=real64)) - Xd = X - type is(complex(kind=real32)) - Xd = X - type is(real(kind=real128)) - Xd = X - type is(real(kind=real64)) - Xd = X - type is(real(kind=real32)) - Xd = X - type is(integer(kind=int64)) - Xd = X - type is(integer(kind=int32)) - Xd = X - end select - fr = Mtimesd(A,Xd) - !fr = Mtimesd_LbnzRule(A,Xd) - end function MtimesdX - !--------------------------------------------------------------------- - - !cmplx*dual - function Mtimesc128d(X,B) result(fr) - complex(real128), intent(in), dimension(:,:) :: X - type(dualzn), intent(in), dimension(:,:) :: B - type(dualzn), dimension(size(X,1),size(B,2)) :: fr - type(dualzn), dimension(size(X,1),size(X,2)) :: Xd - - Xd = X - - fr = Mtimesd(Xd,B) - !fr = Mtimesd_LbnzRule(Xd,B) - end function Mtimesc128d - - - function Mtimesc64d(X,B) result(fr) - complex(real64), intent(in), dimension(:,:) :: X - type(dualzn), intent(in), dimension(:,:) :: B - type(dualzn), dimension(size(X,1),size(B,2)) :: fr - type(dualzn), dimension(size(X,1),size(X,2)) :: Xd - - Xd = X - - fr = Mtimesd(Xd,B) - !fr = Mtimesd_LbnzRule(Xd,B) - end function Mtimesc64d - - function Mtimesc32d(X,B) result(fr) - complex(real32), intent(in), dimension(:,:) :: X - type(dualzn), intent(in), dimension(:,:) :: B - type(dualzn), dimension(size(X,1),size(B,2)) :: fr - type(dualzn), dimension(size(X,1),size(X,2)) :: Xd - - Xd = X - - fr = Mtimesd(Xd,B) - !fr = Mtimesd_LbnzRule(Xd,B) - end function Mtimesc32d - !--------------------------------------------------------------------- - - !real*dual - function Mtimesr128d(X,B) result(fr) - real(real128), intent(in), dimension(:,:) :: X - type(dualzn), intent(in), dimension(:,:) :: B - type(dualzn), dimension(size(X,1),size(B,2)) :: fr - type(dualzn), dimension(size(X,1),size(X,2)) :: Xd - - Xd = X - - fr = Mtimesd(Xd,B) - !fr = Mtimesd_LbnzRule(Xd,B) - end function Mtimesr128d - - function Mtimesr64d(X,B) result(fr) - real(real64), intent(in), dimension(:,:) :: X - type(dualzn), intent(in), dimension(:,:) :: B - type(dualzn), dimension(size(X,1),size(B,2)) :: fr - type(dualzn), dimension(size(X,1),size(X,2)) :: Xd - - Xd = X - - fr = Mtimesd(Xd,B) - !fr = Mtimesd_LbnzRule(Xd,B) - end function Mtimesr64d - - function Mtimesr32d(X,B) result(fr) - real(real32), intent(in), dimension(:,:) :: X - type(dualzn), intent(in), dimension(:,:) :: B - type(dualzn), dimension(size(X,1),size(B,2)) :: fr - type(dualzn), dimension(size(X,1),size(X,2)) :: Xd - - Xd = X - - fr = Mtimesd(Xd,B) - !fr = Mtimesd_LbnzRule(Xd,B) - end function Mtimesr32d - !--------------------------------------------------------------------- - - function Mtimesi64d(X,B) result(fr) - integer(real64), intent(in), dimension(:,:) :: X - type(dualzn), intent(in), dimension(:,:) :: B - type(dualzn), dimension(size(X,1),size(B,2)) :: fr - type(dualzn), dimension(size(X,1),size(X,2)) :: Xd - - Xd = X - - fr = Mtimesd(Xd,B) - !fr = Mtimesd_LbnzRule(Xd,B) - end function Mtimesi64d - - function Mtimesi32d(X,B) result(fr) - integer(real32), intent(in), dimension(:,:) :: X - type(dualzn), intent(in), dimension(:,:) :: B - type(dualzn), dimension(size(X,1),size(B,2)) :: fr - type(dualzn), dimension(size(X,1),size(X,2)) :: Xd - - Xd = X - - fr = Mtimesd(Xd,B) - !fr = Mtimesd_LbnzRule(Xd,B) - end function Mtimesi32d - !--------------------------------------------------------------------- - - !direct implementation of dual matrix multiplication - function Mtimesd(A,B) result(fr) - type(dualzn), intent(in), dimension(:,:) :: A, B - type(dualzn), dimension(size(A,1),size(B,2)) :: fr - integer :: i, j, k, m, n, p, q - - m = size(A,1) - n = size(A,2) - p = size(B,1) - q = size(B,2) - - if (n /= p) then - print*,"Error: Matrix dimensions do not align for multiplication" - return - end if - - call initialize_dualzn(fr) - - do i = 1, m - do j = 1, q - do k = 1, n - fr(i,j) = fr(i,j) + A(i,k)*B(k,j) - end do - end do - end do - end function Mtimesd - !................... - - - !to set the dual-k component to cm in matrix A - subroutine Mset_fpart(k,cm,A) - integer, intent(in) :: k - complex(prec),intent(in), dimension(:,:) :: cm - type(dualzn), intent(inout), dimension(size(cm,1),size(cm,1)) :: A - integer :: i,j - - do i=1,size(cm,1) - do j=1,size(cm,2) - A(i,j)%f(k) = cm(i,j) - end do - end do - end subroutine Mset_fpart - !--------------------------------------------------------------------- - - !A**B - elemental function powerd(A,B) result(fr) - type(dualzn), intent(in) :: A, B - type(dualzn) :: fr - type(dualzn) :: intdual - integer :: iaux - logical :: BIQ - - iaux = int(real(B%f(0))) - intdual = itodn(iaux) - - BIQ = intdual == B - - if(BIQ) then - fr = A**iaux - else - fr = expd(B*logd(A)) - end if - end function powerd - - !A**n (n integer) - elemental function power_dint(A,n) result(fr) - type(dualzn), intent(in) :: A - integer, intent(in) :: n - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - fr%f=(0.0_prec,0.0_prec) - - if(A==itodn(1)) then - fr=A - return - elseif(A==itodn(0) .and. n>0)then - fr=A - return - end if - - !0**0 ---> 1 - if(n==0) then - fr%f(0) = (1.0_prec,0.0_prec) - elseif(n>=1) then - fr%f(0) = (1.0_prec,0.0_prec) - do k=1,n - fr = fr*A - end do - elseif(n<0) then - fr%f(0) = (1.0_prec,0.0_prec) - do k=1,-n - fr = fr*A - end do - fr = inv(fr) - end if - end function power_dint - !--------------------------------------------------------------------- - - !A**n (n integer) - elemental function power_dint8(A,n) result(fr) - type(dualzn), intent(in) :: A - integer(int64), intent(in) :: n - type(dualzn) :: fr - integer(int64) :: k - - allocate(fr%f(0:order)) - fr%f=(0.0_prec,0.0_prec) - - if(A==itodn(1)) then - fr=A - return - elseif(A==itodn(0) .and. n>0)then - fr=A - return - end if - - !0**0 ---> 1 - if(n==0) then - fr%f(0) = (1.0_prec,0.0_prec) - elseif(n>=1) then - fr%f(0) = (1.0_prec,0.0_prec) - do k=1,n - fr = fr*A - end do - elseif(n<0) then - fr%f(0) = (1.0_prec,0.0_prec) - do k=1,-n - fr = fr*A - end do - fr = inv(fr) - end if - end function power_dint8 - !--------------------------------------------------------------------- - - - !chain rule, D^n (f(g)) - recursive pure function Dnd(fc,gdual,n) result(dnfc) - procedure(funzdual) :: fc - type(dualzn), intent(in) :: gdual - integer, intent(in) :: n - complex(prec) :: dnfc - type(dualzn) :: fvd - complex(prec) :: g0, suma - complex(prec), allocatable, dimension(:) :: xvg - integer :: k, j - - g0 = gdual%f(0) - fvd = fc(g0) - if(n==0) then - dnfc = fvd%f(0) - else - suma = 0 - do k=1,n - allocate(xvg(1:n-k+1)) - do j=1,n-k+1 - xvg(j) = gdual%f(j) - end do - suma = suma + fvd%f(k)*BellY(n,k,xvg) - deallocate(xvg) - end do - dnfc = suma - end if - end function Dnd - !--------------------------------------------------------------------- - - !dual/class - elemental function divdX(B,X) result(fr) - class(*), intent(in) :: X - type(dualzn), intent(in) :: B - type(dualzn) :: fr - type(dualzn) :: Xd - - select type (X) - type is (dualzn) - Xd = X - type is(complex(kind=real128)) - Xd = cmplxtodn(X) - type is (complex(kind=real64)) - Xd =cmplxtodn(X) - type is (complex(kind=real32)) - Xd =cmplxtodn(X) - type is (real(kind=real128)) - Xd = realtodn(X) - type is (real(kind=real64)) - Xd = realtodn(X) - type is (real(kind=real32)) - Xd = realtodn(X) - type is (integer(kind=int64)) - Xd = itodn(X) - type is (integer(kind=int32)) - Xd = itodn(X) - end select - fr = divd(B,Xd) - end function divdX - !--------------------------------------------------------------------- - - !class/dual - elemental function divXd(XX,BB) result(fr) - class(*), intent(in) :: XX - type(dualzn), intent(in) :: BB - type(dualzn) :: fr - type(dualzn) :: Xd - - select type (XX) - type is (dualzn) - Xd = XX - type is(complex(kind=real128)) - Xd = cmplxtodn(XX) - type is (complex(kind=real64)) - Xd =cmplxtodn(XX) - type is (complex(kind=real32)) - Xd =cmplxtodn(XX) - type is (real(kind=real128)) - Xd = realtodn(XX) - type is (real(kind=real64)) - Xd = realtodn(XX) - type is (real(kind=real32)) - Xd = realtodn(XX) - type is (integer(kind=int64)) - Xd = itodn(XX) - type is (integer(kind=int32)) - Xd = itodn(XX) - end select - fr = divd(Xd, BB) - end function divXd - !--------------------------------------------------------------------- - - !A/B - elemental function divd(A,B) result(fr) - type(dualzn), intent(in) :: A, B - type(dualzn) :: fr - - fr = A*inv(B) - end function divd - !--------------------------------------------------------------------- - - !dual*class - elemental function timesdX(B,X) result(fr) - class(*), intent(in) :: X - type(dualzn), intent(in) :: B - type(dualzn) :: fr - type(dualzn) :: Xd - - select type (X) - type is (dualzn) - Xd = X - type is(complex(kind=real128)) - Xd = cmplxtodn(X) - type is (complex(kind=real64)) - Xd =cmplxtodn(X) - type is (complex(kind=real32)) - Xd =cmplxtodn(X) - type is (real(kind=real128)) - Xd = realtodn(X) - type is (real(kind=real64)) - Xd = realtodn(X) - type is (real(kind=real32)) - Xd = realtodn(X) - type is (integer(kind=int64)) - Xd = itodn(X) - type is (integer(kind=int32)) - Xd = itodn(X) - end select - fr = timesd(B,Xd) - end function timesdX - !--------------------------------------------------------------------- - - !class*dual - elemental function timesXd(XX,BB) result(fr) - class(*), intent(in) :: XX - type(dualzn), intent(in) :: BB - type(dualzn) :: fr - type(dualzn) :: Xd - - select type (XX) - type is (dualzn) - Xd = XX - type is(complex(kind=real128)) - Xd = cmplxtodn(XX) - type is (complex(kind=real64)) - Xd =cmplxtodn(XX) - type is (complex(kind=real32)) - Xd =cmplxtodn(XX) - type is (real(kind=real128)) - Xd = realtodn(XX) - type is (real(kind=real64)) - Xd = realtodn(XX) - type is (real(kind=real32)) - Xd = realtodn(XX) - type is (integer(kind=int64)) - Xd = itodn(XX) - type is (integer(kind=int32)) - Xd = itodn(XX) - end select - fr = timesd(Xd, BB) - end function timesXd - !--------------------------------------------------------------------- - - !A*B - elemental function timesd(A,B) result(fr) - type(dualzn), intent(in) :: A, B - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - do k=0,order - fr%f(k)=timesdzn(A,B,k) - end do - end function timesd - - pure function timesdzn(A,B,k) result(fr) - type(dualzn), intent(in) :: A, B - integer, intent(in) :: k - complex(prec) :: fr - integer :: i - - fr=0.0_prec - do i=0,k - fr=fr+binomial(k,i)*A%f(i)*B%f(k-i) - end do - end function timesdzn - !--------------------------------------------------------------------- - - !inverse function - elemental function inv(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = Dnd(invzdn,g,k) - end do - end function inv - - pure function invzdn(z) result(fr) - complex(prec), intent(in) :: z - type(dualzn) :: fr - integer :: k,signo - - allocate(fr%f(0:order)) - - signo=1 - do k=0,order - fr%f(k)=signo*gamma(real(k+1,kind=prec))/(z**(k+1)) - signo=-signo - end do - end function invzdn - !--------------------------------------------------------------------- - - !acos function - elemental function acosd_(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = Dnd(acoszdn,g,k) - end do - end function acosd_ - - pure function acoszdn(z) result(fr) - complex(prec), intent(in) :: z - type(dualzn) :: fr - integer :: k, i - complex(prec) :: sumterm, sqz - type(dualzn) :: sqzdn, zdn - - allocate(fr%f(0:order)) - fr%f(0) = acos(z) - if (order == 0) return - - fr%f(1) = -1.0_prec/sqrt(1.0_prec - z**2) - if (order == 1) return - - call initialize_dualzn(zdn) - zdn%f(0) = z - zdn%f(1) = 1.0_prec - - sqz = sqrt(1.0_prec-z*z) - sqzdn = sqrt(realtodn(1.0_prec) - zdn*zdn) - do k=2,order - sumterm = 0.0_prec - do i=1,k-1 - sumterm = sumterm + binomial(k-1,i)*sqzdn%f(i)*fr%f(k-i)/sqz - end do - fr%f(k) = -sumterm - end do - end function acoszdn - !--------------------------------------------------------------------- - - !asin function - elemental function asind_(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = Dnd(asinzdn,g,k) - end do - end function asind_ - - pure function asinzdn(z) result(fr) - complex(prec), intent(in) :: z - type(dualzn) :: fr - integer :: k, i - complex(prec) :: sumterm, sqz - type(dualzn) :: sqzdn, zdn - - allocate(fr%f(0:order)) - fr%f(0) = asin(z) - if (order == 0) return - - fr%f(1) = 1.0_prec/sqrt(1.0_prec - z**2) - if (order == 1) return - - call initialize_dualzn(zdn) - zdn%f(0) = z - zdn%f(1) = 1.0_prec - - sqz = sqrt(1.0_prec-z*z) - sqzdn = sqrt(realtodn(1.0_prec) - zdn*zdn) - do k=2,order - sumterm = 0.0_prec - do i=1,k-1 - sumterm = sumterm + binomial(k-1,i)*sqzdn%f(i)*fr%f(k-i)/sqz - end do - fr%f(k) = -sumterm - end do - end function asinzdn - !--------------------------------------------------------------------- - - !atan function - elemental function atand_(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = Dnd(atanzdn,g,k) - end do - end function atand_ - - pure function atanzdn(z) result(fr) - complex(prec), intent(in) :: z - type(dualzn) :: fr - integer :: k, i - complex(prec) :: sumterm, den - type(dualzn) :: auxdn - - allocate(fr%f(0:order)) - fr%f(0) = atan(z) - - if (order == 0) return - - fr%f(1) = 1.0_prec/(1.0_prec + z**2) - if (order == 1) return - - den = 1.0_prec + z*z - - call initialize_dualzn(auxdn) - auxdn%f(0) = den - auxdn%f(1) = 2.0_prec * z - auxdn%f(2) = 2.0_prec - - do k=2,order - sumterm = 0.0_prec - do i=1,k-1 - sumterm = sumterm + binomial(k-1,i)*auxdn%f(i)*fr%f(k-i)/den - end do - fr%f(k) = -sumterm - end do - end function atanzdn - !--------------------------------------------------------------------- - - !asinh function - elemental function asinhd(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = Dnd(asinhzdn,g,k) - end do - end function asinhd - - pure function asinhzdn(z) result(fr) - complex(prec), intent(in) :: z - type(dualzn) :: fr - integer :: k, i - complex(prec) :: sumterm, den - type(dualzn) :: auxdn, zdn - - allocate(fr%f(0:order)) - fr%f(0) = asinh(z) - if (order == 0) return - - fr%f(1) = 1.0_prec/sqrt(1.0_prec + z**2) - if (order == 1) return - - den = sqrt(1.0_prec + z*z) - - call initialize_dualzn(zdn) - zdn%f(0) = z - zdn%f(1) = 1.0_prec - - auxdn = sqrt(realtodn(1.0_prec) + zdn*zdn) - - do k=2,order - sumterm = 0.0_prec - do i=1,k-1 - sumterm = sumterm + binomial(k-1,i)*auxdn%f(i)*fr%f(k-i)/den - end do - fr%f(k) = -sumterm - end do - end function asinhzdn - !--------------------------------------------------------------------- - - !acosh function - elemental function acoshd(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = Dnd(acoshzdn,g,k) - end do - end function acoshd - - pure function acoshzdn(z) result(fr) - complex(prec), intent(in) :: z - type(dualzn) :: fr - integer :: k, i - complex(prec) :: sumterm, den - type(dualzn) :: auxdn, zdn, onedn - - allocate(fr%f(0:order)) - fr%f(0) = acosh(z) - if (order == 0) return - - !do not 'simplify' they are complex - den = sqrt(z - 1.0_prec) * sqrt(1.0_prec + z) - fr%f(1) = 1.0_prec/den - if (order == 1) return - - call initialize_dualzn(zdn) - zdn%f(0) = z - zdn%f(1) = 1.0_prec - - onedn = realtodn(1.0_prec) - auxdn = sqrt(zdn - onedn) * sqrt(zdn + onedn) - do k=2,order - sumterm = 0.0_prec - do i=1,k-1 - sumterm = sumterm + binomial(k-1,i)*auxdn%f(i)*fr%f(k-i)/den - end do - fr%f(k) = -sumterm - end do - end function acoshzdn - !--------------------------------------------------------------------- - - !atanh function - elemental function atanhd(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = Dnd(atanhzdn,g,k) - end do - end function atanhd - - pure function atanhzdn(z) result(fr) - complex(prec), intent(in) :: z - type(dualzn) :: fr - integer :: k, i - complex(prec) :: sumterm, den - type(dualzn) :: auxdn - - allocate(fr%f(0:order)) - fr%f(0) = atanh(z) - - if (order == 0) return - - !do not 'simplify' they are complex - den = 1.0_prec - z*z - fr%f(1) = 1.0_prec/den - if (order == 1) return - - call initialize_dualzn(auxdn) - auxdn%f(0) = den - auxdn%f(1) = -2.0_prec * z - auxdn%f(2) = -2.0_prec - - do k=2,order - sumterm = 0.0_prec - do i=1,k-1 - sumterm = sumterm + binomial(k-1,i)*auxdn%f(i)*fr%f(k-i)/den - end do - fr%f(k) = -sumterm - end do - end function atanhzdn - !--------------------------------------------------------------------- - - !atan2d function - elemental function atan2d_(y,x) result(fr) - type(dualzn), intent(in) :: y, x - type(dualzn) :: fr - complex(prec) :: x0, y0 - - fr = atan(y/x) - x0 = x%f(0) - y0 = y%f(0) - fr%f(0) = atan2_z(y0,x0) - end function atan2d_ - - !Atan2 for complex arguments - elemental function atan2_z(zy,zx) result (f_res) - complex(prec), intent(in) :: zy, zx - complex(prec) :: f_res - complex(prec), parameter :: ii = cmplx(0,1,prec) - complex(prec) :: num, den, divnd, t1, t2 - real(prec) :: x1, x2, y1, y2 - complex(prec) :: x1c, x2c, y1c, y2c - - x1 = real(zx,kind=prec) - x2 = aimag(zx) - - y1 = real(zy,kind=prec) - y2 = aimag(zy) - - x1c = x1 - x2c = x2 - y1c = y1 - y2c = y2 - - num = x1c + ii*x2c + ii*(y1c + ii*y2c) - den = sqrt((x1c + ii*x2c)**2 + (y1c + ii*y2c)**2) - divnd = num/den; - t1 = atan2(aimag(divnd),real(divnd,kind=prec)) - t2 = ii*log(sqrt((x2c + y1c)**2 + (x1c - y2c)**2)/((2.0_prec*x1c*& - x2c + 2.0_prec*y1c*y2c)**2 + (x1c**2 - x2c**2 + y1c**2 - & - y2c**2)**2)**0.25_prec) - - f_res = t1 - t2 - end function atan2_z - !--------------------------------------------------------------------- - - !sinh function - elemental function sinhd(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - fr = (exp(g)-exp(-g)) - do k=0,order - fr%f(k) = 0.5_prec*fr%f(k) - end do - end function sinhd - !--------------------------------------------------------------------- - - !cosh function - elemental function coshd(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - fr = (exp(g)+exp(-g)) - do k=0,order - fr%f(k) = 0.5_prec*fr%f(k) - end do - end function coshd - !--------------------------------------------------------------------- - - !tanh function - elemental function tanhd(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - - fr = (exp(g) - exp(-g))/(exp(g) + exp(-g)) - end function tanhd - !--------------------------------------------------------------------- - - !absx = sqrt(z*z) is not sqrt(z*conjg(z)) - elemental function absx(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - complex(prec) :: g0 - integer :: k - - g0 = g%f(0) - - fr = g - do k=0,order - fr%f(k) = g%f(k) * g0/sqrt(g0*g0) - end do - end function absx - !--------------------------------------------------------------------- - - !conjg - !notice tat the conjugation operation is not differentiable. In the - !below definitions we mean (df)* not d(f*) - elemental function conjg_dzn(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - call initialize_dualzn(fr) - do k=0, order - fr%f(k) = conjg(g%f(k)) - end do - end function conjg_dzn - !-------------------------------------------------------------------- - - !sum(AR2,dir) for rank2 arrays - !the result is given as array of rank 1 - function sumR2dzn(AR2,dir) result(fr) - type(dualzn), intent(in), dimension(:,:) :: AR2 - integer, intent(in) :: dir - type(dualzn), dimension(size(AR2,2/dir)) :: fr - complex(prec), dimension(size(AR2,2/dir)) :: sAk - complex(prec), dimension(size(AR2,1),size(AR2,2)) :: Ak - integer :: k,j, dimfr - - call initialize_dualzn(fr) - dimfr = size(AR2,2/dir) - - do k=0,order - Ak = f_part(AR2,k) - sAk = sum(Ak,dir) - do j=1,dimfr - fr(j)%f(k) = sAk(j) - end do - end do - end function sumR2dzn - !--------------------------------------------------------------------- - - !sum for Rank 2 array - function sumR20dzn(AR2) result(fr) - type(dualzn), intent(in), dimension(:,:) :: AR2 - type(dualzn) :: fr - complex(prec), dimension(size(AR2,1),size(AR2,2)) :: Azaux - integer :: k - - call initialize_dualzn(fr) - do k=0,order - Azaux = f_part(AR2,k) - fr%f(k) = sum(Azaux) - end do - end function sumR20dzn - !--------------------------------------------------------------------- - - !sum for Rank 1 array - function sumR1dzn(AR1) result(fr) - type(dualzn), intent(in), dimension(:) :: AR1 - type(dualzn) :: fr - complex(prec), dimension(size(AR1)) :: Azaux - integer :: k - - call initialize_dualzn(fr) - do k=0,order - Azaux = f_part(AR1,k) - fr%f(k) = sum(Azaux) - end do - end function sumR1dzn - !--------------------------------------------------------------------- - - !product(AR2,dir) for Rank 2 array - function prodR2dzn(AR2,dir) result(fr) - type(dualzn), intent(in), dimension(:,:) :: AR2 - integer, intent(in) :: dir - type(dualzn), dimension(size(AR2,2/dir)) :: fr - type(dualzn), dimension(size(AR2,dir)) :: vkdir - integer :: k - - if(dir==1) then - do k = 1, size(AR2,2) - vkdir = AR2(:,k) - fr(k) = prodR1dzn(vkdir) - end do - else if(dir==2) then - do k = 1, size(AR2,1) - vkdir = AR2(k,:) - fr(k) = prodR1dzn(vkdir) - end do - else - stop 'use 1 (2) to collapse rows (columns) in product function' - end if - end function prodR2dzn - - !product for Rank 2 array - function prodR20dzn(AR2) result(fr) - type(dualzn), intent(in), dimension(:,:) :: AR2 - type(dualzn) :: fr - integer :: k, m - - m=size(AR2,1) - fr = 1.0_prec - do k=1,m - fr = fr*prodR1dzn(AR2(k,:)) - end do - end function prodR20dzn - - !product for Rank 1 array - function prodR1dzn(x) result(fr) - type(dualzn), intent(in), dimension(:) :: x - type(dualzn) :: fr - integer :: k - - fr = 1.0_prec - do k=1,size(x) - fr = fr*x(k) - end do - end function prodR1dzn - !--------------------------------------------------------------------- - - !log function - elemental function logd(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = Dnd(logzdn,g,k) - end do - end function logd - - pure function logzdn(z) result(fr) - complex(prec), intent(in) :: z - type(dualzn) :: fr - integer :: k, signo - - allocate(fr%f(0:order)) - fr%f(0) = log(z) - - signo = 1 - do k=1,order - fr%f(k) = signo*gamma(real(k,kind=prec))/(z**k) - signo = -signo - end do - end function logzdn - !--------------------------------------------------------------------- - - !exp function - elemental function expd(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = Dnd(expzdn,g,k) - end do - end function expd - - pure function expzdn(z) result(fr) - complex(prec), intent(in) :: z - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = exp(z) - end do - end function expzdn - !--------------------------------------------------------------------- - - !sin function - elemental function sind_(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = Dnd(sinzdn,g,k) - end do - end function sind_ - - pure function sinzdn(z) result(fr) - complex(prec), intent(in) :: z - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = sin(z + k*Pi/2) - end do - end function sinzdn - !--------------------------------------------------------------------- - - !cos function - elemental function cosd_(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = Dnd(coszdn,g,k) - end do - end function cosd_ - - pure function coszdn(z) result(fr) - complex(prec), intent(in) :: z - type(dualzn) :: fr - integer :: k - - allocate(fr%f(0:order)) - - do k=0,order - fr%f(k) = cos(z + k*Pi/2) - end do - end function coszdn - !--------------------------------------------------------------------- - - !tan function - elemental function tand_(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - - fr = sind_(g)*inv(cosd_(g)) - end function tand_ - !--------------------------------------------------------------------- - - !sqrt function - elemental function sqrtd(g) result(fr) - type(dualzn), intent(in) :: g - type(dualzn) :: fr - - fr = g**(0.5_prec) - end function sqrtd - !--------------------------------------------------------------------- - - !Combinations - !m!/((m-n)! * n!) - pure function binomial(m, n) result(binom) - integer, intent(in) :: m, n - real(prec) :: binom - integer :: j - - if (n == 0 .or. n == m) then - binom = 1.0 - else - binom = 1.0 - do j = 1, n - binom = binom*(m-j+1)/j - end do - endif - end function binomial - - pure function BellY(n, k, x) result(result_value) - integer, intent(in) :: n, k - complex(prec), dimension(:), intent(in) :: x - complex(prec) :: result_value - complex(prec), dimension(n+1,k+1) :: dp - complex(prec) :: sum_val - complex(prec), dimension(:), allocatable :: newx - integer :: nn, kk, ii, LX, nx - - dp = 0.0 - dp(1, 1) = 1.0 - - do nn = 1, n - dp(nn+1, 1) = 0.0 - end do - - do kk = 1, k - dp(1, kk+1) = 0.0 - end do - - !special cases - if (n == 0 .and. k == 0) then - result_value = 1.0 - return - elseif (size(x) == 0) then - result_value = 0.0 - return - end if - - !Main loop to compute BellY[n, k, x] - do nn = 1, n - do kk = 1, k - LX = size(x) - nx = max(nn - kk + 1, LX) - if (nx > 0) then - allocate(newx(nx)) - newx = 0.0 - newx(1:LX) = x - sum_val = 0.0 - do ii = 0, nn - kk - sum_val = sum_val + binomial(nn - 1, ii) * & - newx(ii + 1)*dp(nn - ii, kk) - end do - dp(nn + 1, kk + 1) = sum_val - deallocate(newx) - end if - end do - end do - - result_value = dp(n + 1, k + 1) - end function BellY -end module dualzn_mod diff --git a/SourcesF_gfortran/ex3.f90 b/SourcesF_gfortran/ex3.f90 deleted file mode 100644 index 36d87d5..0000000 --- a/SourcesF_gfortran/ex3.f90 +++ /dev/null @@ -1,92 +0,0 @@ -include './precision_mod.f90' -include './dualzn_mod.f90' -include './diff_mod.f90' - -!module with example of functions -module function_mod - use dualzn_mod - implicit none - private - - public :: fstest, fvectest - -contains - !Example of scalar function f = sin(x*y*z) + cos(x*y*z) - function fstest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn) :: fr - type(dualzn) :: x,y,z - - x = r(1); y = r(2); z = r(3) - fr = sin(x*y*z) + cos(x*y*z) - end function fstest - - !Example of vector function f = [f1,f2,f3] - !f = fvectest(r) is a function f:D^m ---> Dn - function fvectest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn), allocatable, dimension(:) :: fr - type(dualzn) :: f1,f2,f3 - type(dualzn) :: x,y,z,w - - x = r(1); y = r(2); z = r(3); w = r(4) - - f1 = sin(x*y*z*w) - f2 = cos(x*y*z*w)*sqrt(w/y - x/z) - f3 = sin(log(x*y*z*w)) - - allocate(fr(3)) - fr = [f1,f2,f3] - end function fvectest -end module function_mod - -!----------------------------------------------------------------------- -!main program -program main - use precision_mod - use dualzn_mod - use diff_mod - use function_mod - implicit none - - integer, parameter :: nf =3, mq = 4 - complex(prec), parameter :: ii = (0,1) - complex(prec), dimension(mq) :: q, vec - complex(prec), dimension(nf,mq) :: Jmat - complex(prec), dimension(nf) :: JV - complex(prec), dimension(3) :: GV - complex(prec), dimension(3,3) :: Hmat - integer :: i - - vec = [1.0_prec,2.0_prec,3.0_prec,4.0_prec] - q = vec/10.0_prec + ii - - print*,"Jv using matmul" - Jmat = Jacobian(fvectest, q , nf) - JV = matmul(Jmat,vec) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"Jv using vector directional derivative" - JV = d1fvector(fvectest,vec,q,nf) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"---Hessian matrix---" - Hmat = Hessian(fstest,q(1:3)) - do i=1,3 - write(*,"(A,i0)") "row:",i - write(*,*) Hmat(i,:) - end do - write(*,*) - - print*,"---Gradient---" - GV = gradient(fstest,q(1:3)) - do i=1,3 - write(*,*) GV(i) - end do -end program main diff --git a/SourcesF_gfortran/precision_mod.f90 b/SourcesF_gfortran/precision_mod.f90 deleted file mode 100644 index c7a6948..0000000 --- a/SourcesF_gfortran/precision_mod.f90 +++ /dev/null @@ -1,7 +0,0 @@ -module precision_mod - use iso_fortran_env, only: int32, int64, real32, real64, real128 - implicit none - - integer, parameter :: prec = real64 - !integer, parameter :: prec = real128 -end module precision_mod diff --git a/Windows/gfortranCase/prec_real128/LibDualzn128/diff_mod.mod b/Windows/gfortranCase/prec_real128/LibDualzn128/diff_mod.mod deleted file mode 100644 index 45e332b..0000000 Binary files a/Windows/gfortranCase/prec_real128/LibDualzn128/diff_mod.mod and /dev/null differ diff --git a/Windows/gfortranCase/prec_real128/LibDualzn128/dualzn_mod.mod b/Windows/gfortranCase/prec_real128/LibDualzn128/dualzn_mod.mod deleted file mode 100644 index 19fad01..0000000 Binary files a/Windows/gfortranCase/prec_real128/LibDualzn128/dualzn_mod.mod and /dev/null differ diff --git a/Windows/gfortranCase/prec_real128/LibDualzn128/libdualzn.a b/Windows/gfortranCase/prec_real128/LibDualzn128/libdualzn.a deleted file mode 100644 index 9a1534d..0000000 Binary files a/Windows/gfortranCase/prec_real128/LibDualzn128/libdualzn.a and /dev/null differ diff --git a/Windows/gfortranCase/prec_real128/LibDualzn128/precision_mod.mod b/Windows/gfortranCase/prec_real128/LibDualzn128/precision_mod.mod deleted file mode 100644 index 6647765..0000000 Binary files a/Windows/gfortranCase/prec_real128/LibDualzn128/precision_mod.mod and /dev/null differ diff --git a/Windows/gfortranCase/prec_real128/examples/create_exec.bat b/Windows/gfortranCase/prec_real128/examples/create_exec.bat deleted file mode 100644 index a34e894..0000000 --- a/Windows/gfortranCase/prec_real128/examples/create_exec.bat +++ /dev/null @@ -1,20 +0,0 @@ -@echo off - -:: Check if at least one argument (the source file) is provided -if "%~1"=="" ( - echo Usage: %~nx0 ^ [output_executable_name] - exit /b 1 -) - -:: Use the provided source file and output executable name or default to "a.exe" -set "SOURCE_FILE=%~1" -set "OUTPUT_NAME=%~2" - -if "%OUTPUT_NAME%"=="" ( - set "OUTPUT_NAME=a.exe" -) - -:: Run the compilation command -gfortran -I../LibDualzn128 -o "%OUTPUT_NAME%" "%SOURCE_FILE%" -L../LibDualzn128 -ldualzn - - diff --git a/Windows/gfortranCase/prec_real128/examples/ex1.f90 b/Windows/gfortranCase/prec_real128/examples/ex1.f90 deleted file mode 100644 index 7d9f062..0000000 --- a/Windows/gfortranCase/prec_real128/examples/ex1.f90 +++ /dev/null @@ -1,35 +0,0 @@ -!gfortran -I../LibDualzn128 -o e1 ex1.f90 -L../LibDualzn128 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - - call set_order(5) !before anything we set the 'order' to work with - - z0 = (1.1_prec,2.2_prec) - r = 0 !we initialize the dual number to 0 - - !we set the 0-th and 1-th components. If dual numbers are used to - !calculate D^n f(z0) then r must be of the form r = r0 + 1*eps_1 - r%f(0) = z0 - r%f(1) = 1 - - fval = sin(r)**log(r*r) - - !writing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - print*,fval%f(k) - end do - deallocate(r%f) - deallocate(fval%f) -end program main - - - - diff --git a/Windows/gfortranCase/prec_real128/examples/ex2.f90 b/Windows/gfortranCase/prec_real128/examples/ex2.f90 deleted file mode 100644 index 2cc3abe..0000000 --- a/Windows/gfortranCase/prec_real128/examples/ex2.f90 +++ /dev/null @@ -1,52 +0,0 @@ -!gfortran -I../LibDualzn128 -o e2 ex2.f90 -L../LibDualzn128 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - real :: t1,t2 - - call set_order(5) !we set the order to work with - - !since a dualzn numbers is an allocatable entity, do not forget to - !initialize it - r = 0 !<--- initializing r to 0 - r%f(0) = (1.1_prec,0.0_prec) - r%f(1) = 1 !since we want to differentiate, r = r0 +1*eps_1 - !all the other components are 0 as r was initialized to 0 - - call cpu_time(t1) - fval = ftest(r) - call cpu_time(t2) - - !Computing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - write(*,"(i0,a,f0.1,a,e17.10)") k,"-th derivative at x = ", & - real(r%f(0)),":",real(fval%f(k)) - end do - - print*,"elapsed time (s):",t2-t1 - deallocate(r%f) - deallocate(fval%f) - -contains - function ftest(x) result(fr) - type(dualzn), intent(in) :: x - type(dualzn) :: fr - integer :: k - - !nested function f(x) = sin(x) * exp(-x^2), f(f(...(f(x))...)) - !applied 1000 times - fr = sin(x)*exp(-x*x) - do k=1, 1000-1 - fr = sin(fr)*exp(-fr*fr) - end do - end function ftest -end program main - - diff --git a/Windows/gfortranCase/prec_real128/examples/ex3.f90 b/Windows/gfortranCase/prec_real128/examples/ex3.f90 deleted file mode 100644 index 041011e..0000000 --- a/Windows/gfortranCase/prec_real128/examples/ex3.f90 +++ /dev/null @@ -1,90 +0,0 @@ -!gfortran -I../LibDualzn128 -o ExDO ex3.f90 -L../LibDualzn128 -ldualzn - -!module with example of functions -module function_mod - use dualzn_mod - implicit none - private - - public :: fstest, fvectest - -contains - !Example of scalar function f = sin(x*y*z) + cos(x*y*z) - function fstest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn) :: fr - type(dualzn) :: x,y,z - - x = r(1); y = r(2); z = r(3) - fr = sin(x*y*z) + cos(x*y*z) - end function fstest - - !Example of vector function f = [f1,f2,f3] - !f = fvectest(r) is a function f:D^m ---> Dn - function fvectest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn), allocatable, dimension(:) :: fr - type(dualzn) :: f1,f2,f3 - type(dualzn) :: x,y,z,w - - x = r(1); y = r(2); z = r(3); w = r(4) - - f1 = sin(x*y*z*w) - f2 = cos(x*y*z*w)*sqrt(w/y - x/z) - f3 = sin(log(x*y*z*w)) - - allocate(fr(3)) - fr = [f1,f2,f3] - end function fvectest -end module function_mod - -!----------------------------------------------------------------------- -!main program -program main - use precision_mod - use dualzn_mod - use diff_mod - use function_mod - implicit none - - integer, parameter :: nf =3, mq = 4 - complex(prec), parameter :: ii = (0,1) - complex(prec), dimension(mq) :: q, vec - complex(prec), dimension(nf,mq) :: Jmat - complex(prec), dimension(nf) :: JV - complex(prec), dimension(3) :: GV - complex(prec), dimension(3,3) :: Hmat - integer :: i - - vec = [1.0_prec,2.0_prec,3.0_prec,4.0_prec] - q = vec/10.0_prec + ii - - print*,"Jv using matmul" - Jmat = Jacobian(fvectest, q , nf) - JV = matmul(Jmat,vec) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"Jv using vector directional derivative" - JV = d1fvector(fvectest,vec,q,nf) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"---Hessian matrix---" - Hmat = Hessian(fstest,q(1:3)) - do i=1,3 - write(*,"(A,i0)") "row:",i - write(*,*) Hmat(i,:) - end do - write(*,*) - - print*,"---Gradient---" - GV = gradient(fstest,q(1:3)) - do i=1,3 - write(*,*) GV(i) - end do -end program main diff --git a/Windows/gfortranCase/prec_real64/LibDualzn64/diff_mod.mod b/Windows/gfortranCase/prec_real64/LibDualzn64/diff_mod.mod deleted file mode 100644 index 94fa76d..0000000 Binary files a/Windows/gfortranCase/prec_real64/LibDualzn64/diff_mod.mod and /dev/null differ diff --git a/Windows/gfortranCase/prec_real64/LibDualzn64/dualzn_mod.mod b/Windows/gfortranCase/prec_real64/LibDualzn64/dualzn_mod.mod deleted file mode 100644 index 978612f..0000000 Binary files a/Windows/gfortranCase/prec_real64/LibDualzn64/dualzn_mod.mod and /dev/null differ diff --git a/Windows/gfortranCase/prec_real64/LibDualzn64/libdualzn.a b/Windows/gfortranCase/prec_real64/LibDualzn64/libdualzn.a deleted file mode 100644 index 7d8eb02..0000000 Binary files a/Windows/gfortranCase/prec_real64/LibDualzn64/libdualzn.a and /dev/null differ diff --git a/Windows/gfortranCase/prec_real64/LibDualzn64/precision_mod.mod b/Windows/gfortranCase/prec_real64/LibDualzn64/precision_mod.mod deleted file mode 100644 index 7709459..0000000 Binary files a/Windows/gfortranCase/prec_real64/LibDualzn64/precision_mod.mod and /dev/null differ diff --git a/Windows/gfortranCase/prec_real64/examples/create_exec.bat b/Windows/gfortranCase/prec_real64/examples/create_exec.bat deleted file mode 100644 index 2e134a4..0000000 --- a/Windows/gfortranCase/prec_real64/examples/create_exec.bat +++ /dev/null @@ -1,20 +0,0 @@ -@echo off - -:: Check if at least one argument (the source file) is provided -if "%~1"=="" ( - echo Usage: %~nx0 ^ [output_executable_name] - exit /b 1 -) - -:: Use the provided source file and output executable name or default to "a.exe" -set "SOURCE_FILE=%~1" -set "OUTPUT_NAME=%~2" - -if "%OUTPUT_NAME%"=="" ( - set "OUTPUT_NAME=a.exe" -) - -:: Run the compilation command -gfortran -I../LibDualzn64 -o "%OUTPUT_NAME%" "%SOURCE_FILE%" -L../LibDualzn64 -ldualzn - - diff --git a/Windows/gfortranCase/prec_real64/examples/ex1.f90 b/Windows/gfortranCase/prec_real64/examples/ex1.f90 deleted file mode 100644 index dbe9810..0000000 --- a/Windows/gfortranCase/prec_real64/examples/ex1.f90 +++ /dev/null @@ -1,35 +0,0 @@ -!gfortran -I../LibDualzn64 -o e1 ex1.f90 -L../LibDualzn64 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - - call set_order(5) !before anything we set the 'order' to work with - - z0 = (1.1_prec,2.2_prec) - r = 0 !we initialize the dual number to 0 - - !we set the 0-th and 1-th components. If dual numbers are used to - !calculate D^n f(z0) then r must be of the form r = r0 + 1*eps_1 - r%f(0) = z0 - r%f(1) = 1 - - fval = sin(r)**log(r*r) - - !writing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - print*,fval%f(k) - end do - deallocate(r%f) - deallocate(fval%f) -end program main - - - - diff --git a/Windows/gfortranCase/prec_real64/examples/ex2.f90 b/Windows/gfortranCase/prec_real64/examples/ex2.f90 deleted file mode 100644 index ccee209..0000000 --- a/Windows/gfortranCase/prec_real64/examples/ex2.f90 +++ /dev/null @@ -1,52 +0,0 @@ -!gfortran -I../LibDualzn64 -o e2 ex2.f90 -L../LibDualzn64 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - real :: t1,t2 - - call set_order(5) !we set the order to work with - - !since a dualzn numbers is an allocatable entity, do not forget to - !initialize it - r = 0 !<--- initializing r to 0 - r%f(0) = (1.1_prec,0.0_prec) - r%f(1) = 1 !since we want to differentiate, r = r0 +1*eps_1 - !all the other components are 0 as r was initialized to 0 - - call cpu_time(t1) - fval = ftest(r) - call cpu_time(t2) - - !Computing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - write(*,"(i0,a,f0.1,a,e17.10)") k,"-th derivative at x = ", & - real(r%f(0)),":",real(fval%f(k)) - end do - - print*,"elapsed time (s):",t2-t1 - deallocate(r%f) - deallocate(fval%f) - -contains - function ftest(x) result(fr) - type(dualzn), intent(in) :: x - type(dualzn) :: fr - integer :: k - - !nested function f(x) = sin(x) * exp(-x^2), f(f(...(f(x))...)) - !applied 1000 times - fr = sin(x)*exp(-x*x) - do k=1, 1000-1 - fr = sin(fr)*exp(-fr*fr) - end do - end function ftest -end program main - - diff --git a/Windows/gfortranCase/prec_real64/examples/ex3.f90 b/Windows/gfortranCase/prec_real64/examples/ex3.f90 deleted file mode 100644 index 73fdc70..0000000 --- a/Windows/gfortranCase/prec_real64/examples/ex3.f90 +++ /dev/null @@ -1,90 +0,0 @@ -!gfortran -I../LibDualzn64 -o e3 ex3.f90 -L../LibDualzn64 -ldualzn - -!module with example of functions -module function_mod - use dualzn_mod - implicit none - private - - public :: fstest, fvectest - -contains - !Example of scalar function f = sin(x*y*z) + cos(x*y*z) - function fstest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn) :: fr - type(dualzn) :: x,y,z - - x = r(1); y = r(2); z = r(3) - fr = sin(x*y*z) + cos(x*y*z) - end function fstest - - !Example of vector function f = [f1,f2,f3] - !f = fvectest(r) is a function f:D^m ---> Dn - function fvectest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn), allocatable, dimension(:) :: fr - type(dualzn) :: f1,f2,f3 - type(dualzn) :: x,y,z,w - - x = r(1); y = r(2); z = r(3); w = r(4) - - f1 = sin(x*y*z*w) - f2 = cos(x*y*z*w)*sqrt(w/y - x/z) - f3 = sin(log(x*y*z*w)) - - allocate(fr(3)) - fr = [f1,f2,f3] - end function fvectest -end module function_mod - -!----------------------------------------------------------------------- -!main program -program main - use precision_mod - use dualzn_mod - use diff_mod - use function_mod - implicit none - - integer, parameter :: nf =3, mq = 4 - complex(prec), parameter :: ii = (0,1) - complex(prec), dimension(mq) :: q, vec - complex(prec), dimension(nf,mq) :: Jmat - complex(prec), dimension(nf) :: JV - complex(prec), dimension(3) :: GV - complex(prec), dimension(3,3) :: Hmat - integer :: i - - vec = [1.0_prec,2.0_prec,3.0_prec,4.0_prec] - q = vec/10.0_prec + ii - - print*,"Jv using matmul" - Jmat = Jacobian(fvectest, q , nf) - JV = matmul(Jmat,vec) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"Jv using vector directional derivative" - JV = d1fvector(fvectest,vec,q,nf) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"---Hessian matrix---" - Hmat = Hessian(fstest,q(1:3)) - do i=1,3 - write(*,"(A,i0)") "row:",i - write(*,*) Hmat(i,:) - end do - write(*,*) - - print*,"---Gradient---" - GV = gradient(fstest,q(1:3)) - do i=1,3 - write(*,*) GV(i) - end do -end program main diff --git a/Windows/ifxCase/prec_real128/LibDualzn128/diff_mod.mod b/Windows/ifxCase/prec_real128/LibDualzn128/diff_mod.mod deleted file mode 100644 index 54e9ddd..0000000 Binary files a/Windows/ifxCase/prec_real128/LibDualzn128/diff_mod.mod and /dev/null differ diff --git a/Windows/ifxCase/prec_real128/LibDualzn128/dualzn_mod.mod b/Windows/ifxCase/prec_real128/LibDualzn128/dualzn_mod.mod deleted file mode 100644 index b51167e..0000000 Binary files a/Windows/ifxCase/prec_real128/LibDualzn128/dualzn_mod.mod and /dev/null differ diff --git a/Windows/ifxCase/prec_real128/LibDualzn128/libdualzn.lib b/Windows/ifxCase/prec_real128/LibDualzn128/libdualzn.lib deleted file mode 100644 index 8cbc9e2..0000000 Binary files a/Windows/ifxCase/prec_real128/LibDualzn128/libdualzn.lib and /dev/null differ diff --git a/Windows/ifxCase/prec_real128/LibDualzn128/precision_mod.mod b/Windows/ifxCase/prec_real128/LibDualzn128/precision_mod.mod deleted file mode 100644 index db5873b..0000000 Binary files a/Windows/ifxCase/prec_real128/LibDualzn128/precision_mod.mod and /dev/null differ diff --git a/Windows/ifxCase/prec_real128/examples/create_exec.bat b/Windows/ifxCase/prec_real128/examples/create_exec.bat deleted file mode 100644 index b738a5f..0000000 --- a/Windows/ifxCase/prec_real128/examples/create_exec.bat +++ /dev/null @@ -1,19 +0,0 @@ -@echo off - -:: Check if at least one argument (the source file) is provided -if "%~1"=="" ( - echo Usage: %~nx0 ^ [output_executable_name] - exit /b 1 -) - -:: Use the provided source file and output executable name or default to "a.exe" -set "SOURCE_FILE=%~1" -set "OUTPUT_NAME=%~2" - -if "%OUTPUT_NAME%"=="" ( - set "OUTPUT_NAME=a.exe" -) - -:: Run the compilation command -ifx -o "%OUTPUT_NAME%" "%SOURCE_FILE%" -I../LibDualzn128 ../LibDualzn128/libdualzn.lib - diff --git a/Windows/ifxCase/prec_real128/examples/ex1.f90 b/Windows/ifxCase/prec_real128/examples/ex1.f90 deleted file mode 100644 index b8bbb11..0000000 --- a/Windows/ifxCase/prec_real128/examples/ex1.f90 +++ /dev/null @@ -1,37 +0,0 @@ -!ifx -I../LibDualzn128 -o e1 ex1.f90 -L../LibDualzn128 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - - call set_order(5) !before anything we set the 'order' to work with - - z0 = (1.1_prec,2.2_prec) - r = 0 !we initialize the dual number to 0 - - !we set the 0-th and 1-th components. If dual numbers are used to - !calculate D^n f(z0) then r must be of the form r = r0 + 1*eps_1 - r%f(0) = z0 - r%f(1) = 1 - - fval = sin(r)**log(r*r) - - !writing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - print*,fval%f(k) - end do - - deallocate(r%f) - deallocate(fval%f) - -end program main - - - - diff --git a/Windows/ifxCase/prec_real128/examples/ex2.f90 b/Windows/ifxCase/prec_real128/examples/ex2.f90 deleted file mode 100644 index 00296b7..0000000 --- a/Windows/ifxCase/prec_real128/examples/ex2.f90 +++ /dev/null @@ -1,52 +0,0 @@ -!ifx -I../LibDualzn128 -o e2 ex2.f90 -L../LibDualzn128 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - real :: t1,t2 - - call set_order(5) !we set the order to work with - - !since a dualzn numbers is an allocatable entity, do not forget to - !initialize it - r = 0 !<--- initializing r to 0 - r%f(0) = (1.1_prec,0.0_prec) - r%f(1) = 1 !since we want to differentiate, r = r0 +1*eps_1 - !all the other components are 0 as r was initialized to 0 - - call cpu_time(t1) - fval = ftest(r) - call cpu_time(t2) - - !Computing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - write(*,"(i0,a,f0.1,a,e17.10)") k,"-th derivative at x = ", & - real(r%f(0)),":",real(fval%f(k)) - end do - - print*,"elapsed time (s):",t2-t1 - deallocate(r%f) - deallocate(fval%f) - -contains - function ftest(x) result(fr) - type(dualzn), intent(in) :: x - type(dualzn) :: fr - integer :: k - - !nested function f(x) = sin(x) * exp(-x^2), f(f(...(f(x))...)) - !applied 1000 times - fr = sin(x)*exp(-x*x) - do k=1, 1000-1 - fr = sin(fr)*exp(-fr*fr) - end do - end function ftest -end program main - - diff --git a/Windows/ifxCase/prec_real128/examples/ex3.f90 b/Windows/ifxCase/prec_real128/examples/ex3.f90 deleted file mode 100644 index d1fe232..0000000 --- a/Windows/ifxCase/prec_real128/examples/ex3.f90 +++ /dev/null @@ -1,90 +0,0 @@ -!ifx -I../LibDualzn128 -o e3 ex3.f90 -L../LibDualzn128 -ldualzn - -!module with example of functions -module function_mod - use dualzn_mod - implicit none - private - - public :: fstest, fvectest - -contains - !Example of scalar function f = sin(x*y*z) + cos(x*y*z) - function fstest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn) :: fr - type(dualzn) :: x,y,z - - x = r(1); y = r(2); z = r(3) - fr = sin(x*y*z) + cos(x*y*z) - end function fstest - - !Example of vector function f = [f1,f2,f3] - !f = fvectest(r) is a function f:D^m ---> Dn - function fvectest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn), allocatable, dimension(:) :: fr - type(dualzn) :: f1,f2,f3 - type(dualzn) :: x,y,z,w - - x = r(1); y = r(2); z = r(3); w = r(4) - - f1 = sin(x*y*z*w) - f2 = cos(x*y*z*w)*sqrt(w/y - x/z) - f3 = sin(log(x*y*z*w)) - - allocate(fr(3)) - fr = [f1,f2,f3] - end function fvectest -end module function_mod - -!----------------------------------------------------------------------- -!main program -program main - use precision_mod - use dualzn_mod - use diff_mod - use function_mod - implicit none - - integer, parameter :: nf =3, mq = 4 - complex(prec), parameter :: ii = (0,1) - complex(prec), dimension(mq) :: q, vec - complex(prec), dimension(nf,mq) :: Jmat - complex(prec), dimension(nf) :: JV - complex(prec), dimension(3) :: GV - complex(prec), dimension(3,3) :: Hmat - integer :: i - - vec = [1.0_prec,2.0_prec,3.0_prec,4.0_prec] - q = vec/10.0_prec + ii - - print*,"Jv using matmul" - Jmat = Jacobian(fvectest, q , nf) - JV = matmul(Jmat,vec) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"Jv using vector directional derivative" - JV = d1fvector(fvectest,vec,q,nf) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"---Hessian matrix---" - Hmat = Hessian(fstest,q(1:3)) - do i=1,3 - write(*,"(A,i0)") "row:",i - write(*,*) Hmat(i,:) - end do - write(*,*) - - print*,"---Gradient---" - GV = gradient(fstest,q(1:3)) - do i=1,3 - write(*,*) GV(i) - end do -end program main diff --git a/Windows/ifxCase/prec_real64/LibDualzn64/diff_mod.mod b/Windows/ifxCase/prec_real64/LibDualzn64/diff_mod.mod deleted file mode 100644 index 1b5dd5f..0000000 Binary files a/Windows/ifxCase/prec_real64/LibDualzn64/diff_mod.mod and /dev/null differ diff --git a/Windows/ifxCase/prec_real64/LibDualzn64/dualzn_mod.mod b/Windows/ifxCase/prec_real64/LibDualzn64/dualzn_mod.mod deleted file mode 100644 index ffd946a..0000000 Binary files a/Windows/ifxCase/prec_real64/LibDualzn64/dualzn_mod.mod and /dev/null differ diff --git a/Windows/ifxCase/prec_real64/LibDualzn64/libdualzn.lib b/Windows/ifxCase/prec_real64/LibDualzn64/libdualzn.lib deleted file mode 100644 index 3ad53e2..0000000 Binary files a/Windows/ifxCase/prec_real64/LibDualzn64/libdualzn.lib and /dev/null differ diff --git a/Windows/ifxCase/prec_real64/LibDualzn64/precision_mod.mod b/Windows/ifxCase/prec_real64/LibDualzn64/precision_mod.mod deleted file mode 100644 index e220496..0000000 Binary files a/Windows/ifxCase/prec_real64/LibDualzn64/precision_mod.mod and /dev/null differ diff --git a/Windows/ifxCase/prec_real64/examples/create_exec.bat b/Windows/ifxCase/prec_real64/examples/create_exec.bat deleted file mode 100644 index e8ee318..0000000 --- a/Windows/ifxCase/prec_real64/examples/create_exec.bat +++ /dev/null @@ -1,19 +0,0 @@ -@echo off - -:: Check if at least one argument (the source file) is provided -if "%~1"=="" ( - echo Usage: %~nx0 ^ [output_executable_name] - exit /b 1 -) - -:: Use the provided source file and output executable name or default to "a.exe" -set "SOURCE_FILE=%~1" -set "OUTPUT_NAME=%~2" - -if "%OUTPUT_NAME%"=="" ( - set "OUTPUT_NAME=a.exe" -) - -:: Run the compilation command -ifx -o "%OUTPUT_NAME%" "%SOURCE_FILE%" -I../LibDualzn64 ../LibDualzn64/libdualzn.lib - diff --git a/Windows/ifxCase/prec_real64/examples/ex1.f90 b/Windows/ifxCase/prec_real64/examples/ex1.f90 deleted file mode 100644 index 92d406b..0000000 --- a/Windows/ifxCase/prec_real64/examples/ex1.f90 +++ /dev/null @@ -1,35 +0,0 @@ -!ifx -I../LibDualzn64 -o e1 ex1.f90 -L../LibDualzn64 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - - call set_order(5) !before anything we set the 'order' to work with - - z0 = (1.1_prec,2.2_prec) - r = 0 !we initialize the dual number to 0 - - !we set the 0-th and 1-th components. If dual numbers are used to - !calculate D^n f(z0) then r must be of the form r = r0 + 1*eps_1 - r%f(0) = z0 - r%f(1) = 1 - - fval = sin(r)**log(r*r) - - !writing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - print*,fval%f(k) - end do - deallocate(r%f) - deallocate(fval%f) -end program main - - - - diff --git a/Windows/ifxCase/prec_real64/examples/ex2.f90 b/Windows/ifxCase/prec_real64/examples/ex2.f90 deleted file mode 100644 index 058fc3f..0000000 --- a/Windows/ifxCase/prec_real64/examples/ex2.f90 +++ /dev/null @@ -1,52 +0,0 @@ -!ifx -I../LibDualzn64 -o e2 ex2.f90 -L../LibDualzn64 -ldualzn -program main - use precision_mod - use dualzn_mod - implicit none - - complex(prec) :: z0 - type(dualzn) :: r, fval - integer :: k - real :: t1,t2 - - call set_order(5) !we set the order to work with - - !since a dualzn numbers is an allocatable entity, do not forget to - !initialize it - r = 0 !<--- initializing r to 0 - r%f(0) = (1.1_prec,0.0_prec) - r%f(1) = 1 !since we want to differentiate, r = r0 +1*eps_1 - !all the other components are 0 as r was initialized to 0 - - call cpu_time(t1) - fval = ftest(r) - call cpu_time(t2) - - !Computing the derivatives, from the 0th derivative up to the - !order-th derivative. - print*,"derivatives" - do k=0,order - write(*,"(i0,a,f0.1,a,e17.10)") k,"-th derivative at x = ", & - real(r%f(0)),":",real(fval%f(k)) - end do - - print*,"elapsed time (s):",t2-t1 - deallocate(r%f) - deallocate(fval%f) - -contains - function ftest(x) result(fr) - type(dualzn), intent(in) :: x - type(dualzn) :: fr - integer :: k - - !nested function f(x) = sin(x) * exp(-x^2), f(f(...(f(x))...)) - !applied 1000 times - fr = sin(x)*exp(-x*x) - do k=1, 1000-1 - fr = sin(fr)*exp(-fr*fr) - end do - end function ftest -end program main - - diff --git a/Windows/ifxCase/prec_real64/examples/ex3.f90 b/Windows/ifxCase/prec_real64/examples/ex3.f90 deleted file mode 100644 index 09882b2..0000000 --- a/Windows/ifxCase/prec_real64/examples/ex3.f90 +++ /dev/null @@ -1,90 +0,0 @@ -!ifx -I../LibDualzn64 -o e3 ex3.f90 -L../LibDualzn64 -ldualzn - -!module with example of functions -module function_mod - use dualzn_mod - implicit none - private - - public :: fstest, fvectest - -contains - !Example of scalar function f = sin(x*y*z) + cos(x*y*z) - function fstest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn) :: fr - type(dualzn) :: x,y,z - - x = r(1); y = r(2); z = r(3) - fr = sin(x*y*z) + cos(x*y*z) - end function fstest - - !Example of vector function f = [f1,f2,f3] - !f = fvectest(r) is a function f:D^m ---> Dn - function fvectest(r) result(fr) - type(dualzn), intent(in), dimension(:) :: r - type(dualzn), allocatable, dimension(:) :: fr - type(dualzn) :: f1,f2,f3 - type(dualzn) :: x,y,z,w - - x = r(1); y = r(2); z = r(3); w = r(4) - - f1 = sin(x*y*z*w) - f2 = cos(x*y*z*w)*sqrt(w/y - x/z) - f3 = sin(log(x*y*z*w)) - - allocate(fr(3)) - fr = [f1,f2,f3] - end function fvectest -end module function_mod - -!----------------------------------------------------------------------- -!main program -program main - use precision_mod - use dualzn_mod - use diff_mod - use function_mod - implicit none - - integer, parameter :: nf =3, mq = 4 - complex(prec), parameter :: ii = (0,1) - complex(prec), dimension(mq) :: q, vec - complex(prec), dimension(nf,mq) :: Jmat - complex(prec), dimension(nf) :: JV - complex(prec), dimension(3) :: GV - complex(prec), dimension(3,3) :: Hmat - integer :: i - - vec = [1.0_prec,2.0_prec,3.0_prec,4.0_prec] - q = vec/10.0_prec + ii - - print*,"Jv using matmul" - Jmat = Jacobian(fvectest, q , nf) - JV = matmul(Jmat,vec) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"Jv using vector directional derivative" - JV = d1fvector(fvectest,vec,q,nf) - do i=1,nf - write(*,*) JV(i) - end do - write(*,*) - - print*,"---Hessian matrix---" - Hmat = Hessian(fstest,q(1:3)) - do i=1,3 - write(*,"(A,i0)") "row:",i - write(*,*) Hmat(i,:) - end do - write(*,*) - - print*,"---Gradient---" - GV = gradient(fstest,q(1:3)) - do i=1,3 - write(*,*) GV(i) - end do -end program main diff --git a/SourcesF_gfortran/ex1.f90 b/example/ex1.f90 similarity index 85% rename from SourcesF_gfortran/ex1.f90 rename to example/ex1.f90 index 09704db..d09f22f 100644 --- a/SourcesF_gfortran/ex1.f90 +++ b/example/ex1.f90 @@ -1,6 +1,3 @@ -include './precision_mod.f90' -include './dualzn_mod.f90' - program main use precision_mod use dualzn_mod @@ -12,14 +9,14 @@ program main call set_order(5) !before anything we set the 'order' to work with local_order = get_order() - + z0 = (1.1_prec,2.2_prec) r = 0 !we initialize the dual number to 0 !we set the 0-th and 1-th components. If dual numbers are used to - !calculate D^n f(z0) then r must be of the form r = r0 + 1*eps_1 + !calculate D^n f(z0) then r must be of the form r = r0 + 1*eps_1 r%f(0) = z0 - r%f(1) = 1 + r%f(1) = 1 fval = sin(r)**log(r*r) @@ -31,10 +28,6 @@ program main end do deallocate(r%f) - deallocate(fval%f) + deallocate(fval%f) end program main - - - - diff --git a/SourcesF_gfortran/ex2.f90 b/example/ex2.f90 similarity index 84% rename from SourcesF_gfortran/ex2.f90 rename to example/ex2.f90 index a9ae322..348788c 100644 --- a/SourcesF_gfortran/ex2.f90 +++ b/example/ex2.f90 @@ -1,6 +1,3 @@ -include './precision_mod.f90' -include './dualzn_mod.f90' - program main use precision_mod use dualzn_mod @@ -15,15 +12,15 @@ program main !since a dualzn numbers is an allocatable entity, do not forget to !initialize it - r = 0 !<--- initializing r to 0 + r = 0 !<--- initializing r to 0 r%f(0) = (1.1_prec,0.0_prec) r%f(1) = 1 !since we want to differentiate, r = r0 +1*eps_1 - !all the other components are 0 as r was initialized to 0 + !all the other components are 0 as r was initialized to 0 call cpu_time(t1) fval = ftest(r) call cpu_time(t2) - + !Computing the derivatives, from the 0th derivative up to the !order-th derivative. print*,"derivatives" @@ -32,7 +29,7 @@ program main real(r%f(0)),":",real(fval%f(k)) end do - print*,"elapsed time (s):",t2-t1 + print*,"elapsed time (s):",t2-t1 deallocate(r%f) deallocate(fval%f) @@ -51,5 +48,3 @@ function ftest(x) result(fr) end do end function ftest end program main - - diff --git a/GNU_Linux/ifx_case/prec_real64/examples/ex3.f90 b/example/ex3.f90 similarity index 92% rename from GNU_Linux/ifx_case/prec_real64/examples/ex3.f90 rename to example/ex3.f90 index 09882b2..98f3c59 100644 --- a/GNU_Linux/ifx_case/prec_real64/examples/ex3.f90 +++ b/example/ex3.f90 @@ -1,5 +1,3 @@ -!ifx -I../LibDualzn64 -o e3 ex3.f90 -L../LibDualzn64 -ldualzn - !module with example of functions module function_mod use dualzn_mod @@ -20,7 +18,7 @@ function fstest(r) result(fr) end function fstest !Example of vector function f = [f1,f2,f3] - !f = fvectest(r) is a function f:D^m ---> Dn + !f = fvectest(r) is a function f:D^m ---> Dn function fvectest(r) result(fr) type(dualzn), intent(in), dimension(:) :: r type(dualzn), allocatable, dimension(:) :: fr @@ -55,25 +53,25 @@ program main complex(prec), dimension(3) :: GV complex(prec), dimension(3,3) :: Hmat integer :: i - + vec = [1.0_prec,2.0_prec,3.0_prec,4.0_prec] q = vec/10.0_prec + ii print*,"Jv using matmul" Jmat = Jacobian(fvectest, q , nf) - JV = matmul(Jmat,vec) + JV = matmul(Jmat,vec) do i=1,nf write(*,*) JV(i) end do write(*,*) - + print*,"Jv using vector directional derivative" JV = d1fvector(fvectest,vec,q,nf) do i=1,nf write(*,*) JV(i) end do write(*,*) - + print*,"---Hessian matrix---" Hmat = Hessian(fstest,q(1:3)) do i=1,3 @@ -81,10 +79,10 @@ program main write(*,*) Hmat(i,:) end do write(*,*) - + print*,"---Gradient---" GV = gradient(fstest,q(1:3)) do i=1,3 write(*,*) GV(i) - end do + end do end program main diff --git a/fpm.toml b/fpm.toml new file mode 100644 index 0000000..3c82fb1 --- /dev/null +++ b/fpm.toml @@ -0,0 +1,21 @@ +name = "dnaoad" +version = "0.1.0" +license = "license" +author = "jaosch" +maintainer = "jasper.schommartz@gmail.com" +copyright = "Copyright 2025, jaosch" +[build] +auto-executables = true +auto-tests = true +auto-examples = true +module-naming = false +[install] +library = false +[fortran] +implicit-typing = false +implicit-external = false +source-form = "free" + +[dev-dependencies] +test-drive.git = "https://github.com/fortran-lang/test-drive" +test-drive.tag = "v0.5.0" diff --git a/SourcesF_gfortran/diff_mod.f90 b/src/diff_mod.f90 similarity index 99% rename from SourcesF_gfortran/diff_mod.f90 rename to src/diff_mod.f90 index 7a0a45c..d6c5639 100644 --- a/SourcesF_gfortran/diff_mod.f90 +++ b/src/diff_mod.f90 @@ -62,7 +62,7 @@ function Hessian(fsd,qcmplx) result(Hmat) ej = 0 ej(j) = 1 Hmat(i,j) = 0.5_prec * (d2fscalar(fsd,ei + ej,qcmplx) - Hmat(i,i) & - Hmat(j,j)) + - Hmat(j,j)) Hmat(j,i) = Hmat(i,j) end do end do diff --git a/SourcesF/dualzn_mod.f90 b/src/dualzn_mod.f90 similarity index 100% rename from SourcesF/dualzn_mod.f90 rename to src/dualzn_mod.f90 diff --git a/SourcesF/precision_mod.f90 b/src/precision_mod.f90 similarity index 100% rename from SourcesF/precision_mod.f90 rename to src/precision_mod.f90 diff --git a/test/test_dualzn_mod.f90 b/test/test_dualzn_mod.f90 new file mode 100644 index 0000000..b2cf1af --- /dev/null +++ b/test/test_dualzn_mod.f90 @@ -0,0 +1,37 @@ +module test_dualzn_mod + use testdrive, only: error_type, unittest_type, new_unittest, check + use dualzn_mod + use precision_mod, only: dp => prec + implicit none + private + + public :: collect_dualzn_mod + +contains + subroutine collect_dualzn_mod(testsuite) + type(unittest_type), allocatable, intent(out) :: testsuite(:) + testsuite = [new_unittest("add", test_add)] + end subroutine collect_dualzn_mod + + subroutine test_add(error) + type(error_type), allocatable, intent(out) :: error + + integer, parameter :: order = 2 + type(dualzn) :: x, y + + call set_order(order) + call initialize_dualzn(x) + call initialize_dualzn(y) + + x%f(0) = 2.0_dp + x%f(1) = 1.0_dp + + ! Add two dual numbers + y = x + x + call check(error, y%f(0), (4.0_dp, 0.0_dp)) + call check(error, y%f(1), (2.0_dp, 0.0_dp)) + call check(error, y%f(2), (0.0_dp, 0.0_dp)) + + end subroutine test_add + +end module test_dualzn_mod diff --git a/test/tester.f90 b/test/tester.f90 new file mode 100644 index 0000000..99a67e8 --- /dev/null +++ b/test/tester.f90 @@ -0,0 +1,16 @@ +program tester + use, intrinsic :: iso_fortran_env, only: error_unit + use testdrive, only: run_testsuite + use test_dualzn_mod, only: collect_dualzn_mod + implicit none + + integer :: stat + stat = 0 + + call run_testsuite(collect_dualzn_mod, error_unit, stat) + + if (stat > 0) then + write(error_unit, '(i0, 1x, "test(s) failed!")') stat + end if + +end program tester