Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

  • Congratulations strongm on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

Appropriate use of allocatable arrays

Status
Not open for further replies.

dibas

Programmer
Nov 23, 2011
24
ZA
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).


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_step

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.
 
To help you, a glimpse at your type declarations, apparently in sf_1d_mod, would be needed.



Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.
 
Just found some more time to have a look at your code.

There is nothing unusual in how you build your matrices u and u0 (they are the 'A' in Ax = b, right ?), provided matching dimensions.

Where are the numerical values of a, b, x, y, p and q set ? Did you check that they found their ways into your arrays before this subroutine is called ?

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.
 
Here are the type declarations in "sf_1d_mod" below.
The numerical values of a, b, x, y, p and q are computed in some other
subroutine before they are "pointed" to in this "backEuler_time_step"
one. Even though I haven't done a specific check, I believe they do find
their ways into the arrays (u & u0) because I checked the the jacobian
computations where the array u is also used and found them to be ok.

The matrix "A" in "Ax=b" is not necessarily the u & u0 themselves, but
sections of some array (coeff_matrix(nvars, nvars, 0:n)) that is formed
from them. So, for each 0<=i<=n, "A" is an nvars-by-nvars (2D) array,
which is a section of "coeff_matrix" for a fixed i. In the code, the "A"
is the array "left_matrix".

Code:
module sf_1d_mod
  
   use real_type_mod

   integer  :: n , index

   real(kind=wp), dimension(:), allocatable :: r, alpha, beta, &
                                             a_rhs, b_rhs, x_rhs, &
                                             y_rhs, p_rhs, q_rhs

   real(kind=wp), dimension(:), allocatable :: a_schw, b_schw, x_schw, &
                                              y_schw, p_schw, q_schw

                                                                           real(kind=wp), dimension(:), allocatable, target :: a, b, x, y, p, q, &
                                                       a_tmp, b_tmp, &
                                                       x_tmp, &
                                                       y_tmp, p_tmp, q_tmp

! The pointers used to select the active variables.

   real(kind=wp), dimension(:), pointer  :: a_cur, b_cur, x_cur, &
                                            y_cur, p_cur, q_cur


! Analysis variables: Ricci scalar, constraints, mass energy, expansion
! light cone variables

! The grid spacing (dr), the time (t) and the timestep (dt)

   real(kind=wp) :: dr, t, dt, dt_cur

  real(kind=wp), dimension(:, :, :), allocatable :: coeff_matrix, &
                                     iter_matrix
  real(kind=wp), dimension(:, :), allocatable :: rhs_fun, u, u0
  real(kind=wp), dimension(:, :), allocatable :: drhs, du
  real(kind=wp), dimension(:,:), allocatable  :: left_matrix
  real(kind=wp), dimension(:, :), allocatable :: rhs_vector_tmp, &
                                                 rhs_vector
  real(kind=wp), dimension(:), allocatable   :: right_vector, &
                                       neg_right_vector
end module sf_1d_mod
 
So, let's see if I got you right:

Your problem is, that a matrix which is compiled of elements form u, which itself is compiled from a, b etc., does show zeros only. At one point in time you found that u holds all the values you expect it to do which are not zeros.

If this matrix in question is your left_matrix then look at this piece of your code:

Code:
         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)

In first line you assign one single value to all of left matrix, nvars is 4 so for i = 0 this assigns coeff_matrix (4, 4, 0) to all elements of left_matrix. If this is value is 0, then all of coeff_matrix will be zeroed. In the write statements follows you write out all of left matrix, which consists of the same value all over.

And btw, it would be worthwhile to printout iter_matrix as it is returned from jacobians().

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.
 
Status
Not open for further replies.

Part and Inventory Search

Sponsor

Back
Top