FORTRAN - Correct values are calculated in simple precision, but erroneous solutions in double precision? -
i have following code, regarding newton's method calculating zeros in multivariable vectorial functions:
program newton_method_r_n_r_n_1_numeric use inverse_matrices implicit none integer :: i, j !index real (kind = 8) :: x(0:501,3) !vectors in newton method real (kind = 8) :: a(3,3) !jacobi's matrix each vector integer , parameter :: n = 3 !dimension of jacobi's matriz real (kind = 8) :: inv(3,3) !inverse jacobi's matrix real (kind = 8), parameter :: eps = 1d-8 !near-zero value real (kind = 8) :: det !determinant x(0,:) = (/2.d0,2.d0,2.d0/) = 0, 500 = h(x(i,:)) call lu_factorization (n, a) call determinant (n, a, det) if (abs(det)<eps) exit end if call inverse_matrix (n, a, inv) x(i + 1,:) = x(i,:) - matmul(inv,g(x(i,:))) end write(*,*) ' 0 of function is: ', x(i,:) read(*,*) contains function g(t) implicit none real (kind = 8), intent(in) :: t(3) real (kind = 8) :: g(3) g(1) = t(1)**2 - t(2)**2 + 1.0d0 g(2) = 2.0d0*t(1)*t(2) g(3) = t(3)**2 - 2.0d0 end function g function g_1(t_1) implicit none real (kind = 8), intent(in) :: t_1(3) real (kind = 8) :: g_1 g_1 = t_1(1)**2 - t_1(2)**2 + 1.0d0 end function g_1 function g_2(t_2) implicit none real (kind = 8), intent(in) :: t_2(3) real (kind = 8) :: g_2 g_2 = (2.0d0)*t_2(1)*t_2(2) end function g_2 function g_3(t_3) implicit none real (kind = 8), intent(in) :: t_3(3) real (kind = 8) :: g_3 g_3 = t_3(3)**2 - 2.0d0 end function g_3 function h(q) implicit none real (kind = 8), intent(in) :: q(3) real (kind = 8) :: h(3,3) real (kind = 8), parameter :: dx = 1d-4 real (kind = 8) :: a(3) real (kind = 8) :: b(3) real (kind = 8) :: c(3) a(1) = dx a(2) = 0.d0 a(3) = 0.d0 b(1) = 0.d0 b(2) = dx b(3) = 0.d0 c(1) = 0.d0 c(2) = 0.d0 c(3) = dx h(1,1) = ((g_1(q + a)) - g_1(q)) / dx h(1,2) = ((g_1(q + b)) - g_1(q)) / dx h(1,3) = ((g_1(q + c)) - g_1(q)) / dx h(2,1) = ((g_2(q + a)) - g_2(q)) / dx h(2,2) = ((g_2(q + b)) - g_2(q)) / dx h(2,3) = ((g_2(q + c)) - g_2(q)) / dx h(3,1) = ((g_3(q + a)) - g_3(q)) / dx h(3,2) = ((g_3(q + b)) - g_3(q)) / dx h(3,3) = ((g_3(q + c)) - g_3(q)) / dx end function h end program in simple precision solutions correctly displayed. curious thing when make program write vectors each step, final solution comes correctly, eliminating write statement loop, final solution comes not number. double precision seems messing results, in simple precision works properly, , cannot find error might be. perhaps writing functions wrong in double precision?
any appreciated.
this module, if helpful:
module matrices_inversas_2 implicit none contains subroutine factorizacion_lu (n, a) implicit none integer, intent(in) :: n !dimensiÓn del sistema (matriz de coef.) real (kind = 8), intent(inout) :: a(n,n) !matriz de coeficientes del sistema integer :: i, j, k !Índices real (kind = 8) :: v(n), w(n) !vectores auxiliares !factorizaciÓn lu k = 1, n !columnas = 1, (k - 1) v(i) = a(k,i) end j = k, n = 1, (k - 1) w(i) = a(i,j) end a(k,j) = a(k,j) - dot_product(v,w) end !filas = 1, (k - 1) w(i) = a(i,k) end = (k + 1), n j = 1, (k - 1) v(j) = a(i,j) end a(i,k) = (a(i,k) - dot_product(v,w)) / a(k,k) end end end subroutine factorizacion_lu subroutine sustitucion_directa (n, a, b, c, v) implicit none integer, intent(in) :: n !dimensiÓn del sistema real (kind = 8), intent(inout) :: a(n,n) !matrices del sistema real (kind = 8), intent (in) :: b(n) !matriz de tÉrminos independientes real (kind = 8), intent(out) :: c(n) !matriz de incÓgnitas real (kind = 8), intent(out) :: v(n) !vector auxiliar integer :: i, j !Índices en las matrices real (kind = 8) :: s !sumatorio = 1, n v(i) = a(i,i) a(i,i) = 1.d0 end !sustituciÓn directa c(1) = b(1) / a(1,1) = 2, n s = 0.d0 j = 1, (i - 1) s = s + a(i,j)*c(j) end c(i) = (b(i) - s) / a(i,i) end end subroutine sustitucion_directa subroutine sustitucion_regresiva (n, a, c, x, v) implicit none integer, intent(in) :: n !dimensiÓn del sistema real (kind = 8), intent(inout) :: a(n,n) !matriz de coeficientes real (kind = 8), intent(in) :: c(n) !matriz de tÉrminos independientes real (kind = 8), intent(out) :: x(n) !matriz de incÓgnitas real (kind = 8), intent(in) :: v(n) !vector auxiliar integer :: i,j !Índices en las matrices real (kind = 8) :: s !sumatorio en sustituciÓn regresiva = 1, n a(i,i) = v(i) end !sustitcuiÓn regresiva x(n) = c(n) / a(n,n) = (n - 1), 1, -1 s = 0.d0 j = (i + 1), n s = s + a(i,j) * x(j) end x(i) = (c(i) - s) / a(i,i) end end subroutine sustitucion_regresiva subroutine matriz_inversa (n, a, inv) implicit none integer, intent(in) :: n !dimensiÓn del sistema real (kind = 8), intent(inout) :: a(n,n) !matriz de la que hallar la inversa real (kind = 8), intent(out) :: inv(n,n) !inversa integer :: i, j !Índices real (kind = 8) :: b(n) !matriz de tÉrminos independientes real (kind = 8) :: c(n) !matriz intermedia de incÓgnitas real (kind = 8) :: x(n) !matriz de incÓgnitas real (kind = 8) :: v(n) !vector auxiliar = 1, n b = 0.d0 b(i) = 1.d0 call sustitucion_directa (n, a, b, c, v) call sustitucion_regresiva (n, a, c, x, v) j = 1, n inv(j,i) = x(j) end end end subroutine matriz_inversa subroutine determinante (n, a, det) integer, intent(in) :: n !dimensiÓn del sistema real (kind = 8), intent(inout) :: a(n,n)!matriz de entrada real(kind = 8), intent(out) :: det !valor del determinante integer :: !Índice det = 1.d0 = 1, n det = det * a(i,i) end end subroutine determinante end module
you have provided module matrices_inversas_2 containing apparently-spanish named procedures, program uses module called inverse_matrices, still not have fully-compilable set of code work with.
you have not specified consider 'correctly displayed'.
i made obvious changes module code in order procedure names match , compiled using nagfor -kind=byte -u -c=all -c=undefined -gline. running, get
mod.f90: [nag fortran compiler normal termination] test.f90: warning: test.f90, line 93: unused local variable j [nag fortran compiler normal termination, 1 warning] loading... runtime error: mod.f90, line 27: reference undefined variable v program terminated fatal error mod.f90, line 27: error occurred in matrices_inversas_2:factorizacion_lu test.f90, line 18: called newton_method_r_n_r_n_1_numeric abort (core dumped) you attempting dot_product of 2 size-n arrays (v , w) seems code first (k-1) elements initialized. if replace both of dot_product calls dot_product(v(1:(k-1)),w(1:(k-1))) program runs completion , displays
a 0 of function is: 0.0000000000000000 1.0000000000000000 1.4142135623730951 
Comments
Post a Comment