CODEWITHSUNDEEP
arraysfortranrankmismatch
Noco
Noco

Reputation: 70

Rank mismatch in array reference at (1) (1/2)

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

Upvotes: 0

Views: 863

Answers (1)

Noco
Noco

Reputation: 70

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

Upvotes: 2

Related Questions