I am new to Fortran and basically just playing around at university. My code represents a more dimensional Newton Algorithm to solve residual problems. I put all functions and subroutines in one file to exclude the possibility to get errors for wrong inclusion. My problem is that I get the following error in line 52:
J = DFuncDx(x)
Error: Rank mismatch in array reference at (1) (1/2)
I know that there must be a problem with the dimensions of the passed argument, but I cannot figure out where that might be. In my understanding the variable x is declared as a 3 dimensional array and the function expects 3 dimensions.
I use "gfortran -o nDNewtonMain nDNewtonMain.f90 -l lapack" to compile on Linux.
Additionally, I get a "REAL array index at (1)" warning whenever I call the functions Funct or DFuncDx.
program nDNewton
IMPLICIT NONE
integer, parameter :: N = 3
double precision, dimension(N) :: x
integer, parameter :: NMAXIT = 10
double precision, parameter :: MAXERR = 1.d-8
integer :: Nit, i
double precision :: ERR
do i=1,N
x(i) = 0.d0
end do
Nit = 0; ERR = 1.d0
call myNewton( x, Nit, NMAXIT, MAXERR, ERR, N)
write(*,*) "Solution: x="
do i=1,N
write(*,*) x(i)
end do
write(*,*) "Found after ", Nit, " iterations. Error: ", ERR
end program nDNewton
subroutine myNewton(x, Nit, NMAXIT, MAXERR, ERR, N)
IMPLICIT NONE
integer, intent(in) :: N
double precision, intent(inout), dimension(N) :: x
integer, intent(in) :: NMAXIT
integer, intent(inout) :: Nit
double precision, intent(in) :: MAXERR
double precision, intent(inout) :: ERR
double precision, dimension(N, N) :: J
integer, dimension(N) :: ipiv
integer :: info, i
double precision, dimension(N) :: F
double precision, dimension(N) :: Funct
double precision, dimension(N, N) :: DFuncDx
Nit = 0
F = Funct(x)
ERR = norm2(F)
do while (Nit < NMAXIT .AND. ERR > abs(MAXERR))
J = DFuncDx(x)
call DGESV( 3, 1, J, 3, ipiv, -F, 3, info )
x = x + F
Nit = Nit + 1
F = Funct(x)
ERR = norm2(F)
do i=1, 3
write(*,*) x(i)
end do
write(*,*) "Error: ", ERR,", Nit: ", Nit
end do
end subroutine myNewton
function Funct(x)
IMPLICIT NONE
double precision, dimension(3), intent(in) :: x
double precision, dimension(3, 3) :: ASys
double precision, dimension(3) :: a, e1
double precision, dimension(3) :: Funct
double precision :: exp_part
ASys = reshape( (/ 1.d0,2.d0/3.d0,-2.d0,0.d0,7.d0,6.d0/7.d0,-1.d0/3.d0,0.d0,1.d0 /), (/ 3,3/) )
ASys = transpose(ASys)
a = (/ 5.0, 3.0, 2.0 /)
e1 = (/ 1.0, 0.0, 0.0 /)
Funct = matmul(ASys, x) + a + exp_part(x(1), x(2), x(3)) * e1
end function
function DFuncDx(x)
IMPLICIT NONE
double precision, dimension(3), intent(in) :: x
double precision, dimension(3, 3) :: DFDx
double precision :: val
double precision, dimension(3, 3) :: DFuncDx
double precision :: exp_part
val = exp_part(x(1), x(2), x(3))
DFDx = reshape( (/ 1.d0 + val, 2.d0/3.d0 + 7.d-1*val, -2.d0 + 5.d-1*val, &
0.d0, 7.d0, 6.d0/7.d0, -1.d0/3.d0, 0.d0, 1.d0 /), (/ 3,3/) )
DFuncDx = transpose(DFDx)
end function
double precision function exp_part(x1, x2, x3)
IMPLICIT NONE
double precision, intent(in) :: x1, x2, x3
double precision, dimension(3) :: b
b = (/ 1.0, 0.7, 0.5 /)
exp_part = exp(b(1)*x1 + b(2)*x2 + b(3)*x3)
end function exp_part
Thanks to High Performance Mark I was able to figure out a solution. Maybe this will help another user in the future:
module A
contains
subroutine myNewton(x, Nit, NMAXIT, MAXERR, ERR, N)
IMPLICIT NONE
integer, intent(in) :: N
double precision, intent(inout), dimension(N) :: x
integer, intent(in) :: NMAXIT
integer, intent(inout) :: Nit
double precision, intent(in) :: MAXERR
double precision, intent(inout) :: ERR
double precision, dimension(N, N) :: J
integer, dimension(N) :: ipiv
integer :: info, i
double precision, dimension(N) :: F, B
Nit = 0
F = Funct(x)
ERR = norm2(F)
do while (Nit < NMAXIT .AND. ERR > abs(MAXERR))
J = DFuncDx(x)
B = -F
call DGESV( 3, 1, J, 3, ipiv, B, 3, info )
x = x + B
Nit = Nit + 1
F = Funct(x)
ERR = norm2(F)
do i=1, 3
write(*,*) x(i)
end do
write(*,*) "Error: ", ERR,", Nit: ", Nit
end do
end subroutine myNewton
function Funct(x) result(res)
IMPLICIT NONE
double precision, dimension(3), intent(in) :: x
double precision, dimension(3, 3) :: ASys
double precision, dimension(3) :: a, e1
double precision, dimension(3) :: res
ASys = reshape( (/ 1.d0,2.d0/3.d0,-2.d0,0.d0,7.d0,6.d0/7.d0,-1.d0/3.d0,0.d0,1.d0 /), (/ 3,3/) )
ASys = transpose(ASys)
a = (/ 5.0, 3.0, 2.0 /)
e1 = (/ 1.0, 0.0, 0.0 /)
res = matmul(ASys, x) + a + exp_part(x(1), x(2), x(3)) * e1
end function
function DFuncDx(x) result(res)
IMPLICIT NONE
double precision, dimension(3), intent(in) :: x
double precision, dimension(3, 3) :: DFDx
double precision :: val
double precision, dimension(3, 3) :: res
val = exp_part(x(1), x(2), x(3))
DFDx = reshape( (/ 1.d0 + val, 2.d0/3.d0 + 7.d-1*val, -2.d0 + 5.d-1*val, &
0.d0, 7.d0, 6.d0/7.d0, -1.d0/3.d0, 0.d0, 1.d0 /), (/ 3,3/) )
!DFDx
res = transpose(DFDx)
end function
double precision function exp_part(x1, x2, x3)
IMPLICIT NONE
double precision, intent(in) :: x1, x2, x3
double precision, dimension(3) :: b
b = (/ 1.0, 0.7, 0.5 /)
exp_part = exp(b(1)*x1 + b(2)*x2 + b(3)*x3)
end function exp_part
end module
program nDNewton
use A
IMPLICIT NONE
integer, parameter :: N = 3
double precision, dimension(N) :: x
integer, parameter :: NMAXIT = 10
double precision, parameter :: MAXERR = 1.d-8
integer :: Nit, i
double precision :: ERR
do i=1,N
x(i) = 0.d0
end do
Nit = 0; ERR = 1.d0
call myNewton( x, Nit, NMAXIT, MAXERR, ERR, N)
write(*,*) "Solution: x="
do i=1,N
write(*,*) x(i)
end do
write(*,*) "Found after ", Nit, " iterations. Error: ", ERR
x = Funct(x)
write(*,*) "Residium Test: r="
do i=1,N
write(*,*) x(i)
end do
end program nDNewton