When I compile and link the following piece of code in the relevant program (main), all works perfectly and the executable is formed
successfully. But when running the executable I find that the matrix
"A" in "Ax=b" is singular (only zeros, in fact). I suspect this is due
to the way I have used allocatable arrays (instead of, maybe, automatic
arrays). Actually, all the arrays appearing below have been specified as allocatable arrays. Can anyone please look at my code and help me know
if this is the appropriate use of allocatable arrays (or dynamic memory).
Note: My Newton-Raphson will only terminate based on the number
of iterations (maxiter=2) for now because I'm still testing the
code. Otherwise I also have two other criteria using the infinity
norm of (i) the function values, (ii) the roots.
successfully. But when running the executable I find that the matrix
"A" in "Ax=b" is singular (only zeros, in fact). I suspect this is due
to the way I have used allocatable arrays (instead of, maybe, automatic
arrays). Actually, all the arrays appearing below have been specified as allocatable arrays. Can anyone please look at my code and help me know
if this is the appropriate use of allocatable arrays (or dynamic memory).
Code:
subroutine backEuler_time_step ( )
use sf_1d_mod ! carries definitions of the allocatable arrays
use parameters_mod ! carries parameters definitions.
use new_nr, only : lubksb,ludcmp ! subroutines to solve "Ax=b"
implicit none
integer, dimension(nvars) :: indx
real(8) :: dd
integer :: i , k1, k2, niter
! Set pointers to the current variables
a_cur => a
b_cur => b
x_cur => x
y_cur => y
p_cur => p
q_cur => q
! Store the values of 1D-arrays a,b,x,y in u0 (u0 fixed between Newton iterations)
u0(1, -1:n+1) = a ; u0(2, -1:n+1) = b !u0 = allocatable 2D-array
u0(3, -1:n+1) = x ; u0(4, -1:n+1) = y
!u0(5, -1:n+1) = p ; u0(6, -1:n+1) = q
do niter= 1, maxiter ! Newton loop. mxiter = 2 : testing purposes
u(1, -1:n+1) = a ; u(2, -1:n+1) = b
u(3, -1:n+1) = x ; u(4, -1:n+1) = y
!u(5, -1:n+1) = p ; u(6, -1:n+1) = q
call RHS() !compute right-hand-side functions a_rhs,..,y_rhs
rhs_fun(1, -1:n+1) = a_rhs ; rhs_fun(2, -1:n+1) = b_rhs
rhs_fun(3, -1:n+1) = x_rhs ; rhs_fun(4, -1:n+1) = y_rhs
!rhs_fun(5, -1:n+1) = p_rhs ; rhs_fun(6, -1:n+1) = q_rhs
du = u(:, 1:) - u(:, :n-1)
drhs = rhs_fun(:, 1:) - rhs_fun(:, :n-1)
call jacobians(du, drhs, iter_matrix) ! compute jacobians from the preceding statements; also compute iteration matrix (a 3D-array)
coeff_matrix = iter_matrix ! 3D-arrays dimension(nvars,nvars,0:n)
rhs_vector_tmp = dt*rhs_fun - u ! 2D-arrays, dimension(nvars,0:n)
rhs_vector = rhs_vector_tmp + u0 ! Full step t --> t+dt. this is "b" in "Ax=b"
do i = 0, n ! loop over the r-line i.e. r(0):r(n)
right_vector = rhs_vector(nvars,i) !nvars = 4 (no. of variables)
left_matrix = coeff_matrix(nvars,nvars,i)
write(*,*)
write(*,*) "Printing the iteration/left matrix at i=",i
write(*,"(4(f10.5,x))") ((left_matrix(k1, k2), &
k2=1,nvars), k1=1,nvars)
neg_right_vector = -right_vector
indx(nvars) = 0 ; dd = -1.0
call ludcmp(left_matrix, indx, dd) ! performs lu-decomposition
call lubksb(left_matrix, indx, neg_right_vector)
u(:,i) = u(:,i) + neg_right_vector !
end do ! r-loop ends
call extrapolate() ! calculate a,b,..,y at fictitious points r(-1) & r(n+1)
a = u(1, -1:n+1) ; b = u(2, -1:n+1)
x = u(3, -1:n+1) ; y = u(4, -1:n+1)
!p = u(5, -1:n+1) ; q = u(6, -1:n+1)
! update current values of variables for the next pass thru Newton loop
a_cur => a
b_cur => b
x_cur => x
y_cur => y
!!p_cur => p
!!q_cur => q
end do ! Newton-Raphson loop ends
! Finally set the pointers to the current values of variables
a_cur => a
b_cur => b
x_cur => x
y_cur => y
p_cur => p
q_cur => q
end subroutine backEuler_time_stepNote: My Newton-Raphson will only terminate based on the number
of iterations (maxiter=2) for now because I'm still testing the
code. Otherwise I also have two other criteria using the infinity
norm of (i) the function values, (ii) the roots.