Fortran runtime-error 'memory allocation failed' occurs while I'm only asking for about 100 real(8)'s of memory.
I'm writing some code in Fortran. I used several (about ten) allocatable real(8) arrays. When I allocate them, my command line shows me 'Memory allocation failed'. But I'm only allocating their sizes to about ten. I've briefly checked my computer and I don't think it's in a condition which cannot handle 100 float-point numbers.
The code is something like this:
module m
contains
function foo (a, b, c) result (...)
implicit none
real(8), allocatable, intent(inout) :: a(:), b(:), c(:)
allocate(a(0:10)) ! Actual size include an integer n passed to here. I'm pretty sure n is no greater then 10
allocate(b(1:10))
allocate(c(0:9))
! Some calculation
a_variable = bar(...)
endfunction foo
function bar (...) result (x)
implicit none
real(8), allocatable :: x(:)
allocate(x(8)) ! This is line X
! Some calculation
endfunction bar
endmodule m
program testm
use m
implicit none
! Some code that calls foo
endprogram testm
When I compile and run the code, the command line tells me 'Memory allocation failed' and an error message from the operating system (Windows 10) 'Not enough memory resources to complete this command'. That's really weird since I'm only asking for about 100 float-point numbers of memory.
Something even more puzzling happens when I try to ensure that I didn't pass unexpectedly a big number to n.
When I add a print*, n to line X, the exception stays and the indicated line remains the same line X, which is now the print code!
I'm just completely puzzled what kind of bug could this be. Since I'm a beginner in Fortran, I'm not sure if I've made any undefined behaviors here.
I'm using g95 to compile a .f90 file on Windows 10.
Here's the complete code:
module SplineInterp
! MODULE COMMENT
! Cubic Spline Interpolation Module
! Does Cubic Spline Interpolation with natural boundary conditions
! All calculation is made in real(8). Parameter rk = 8 is defined for real kind specifications.
! In all places where an unexpected zero on the denominator appears, the program will stop with a message,
! usually "The matrix is singular!". A number n is considered to be zero if dabs(n) < e_min, e_min is a
! pre-defined parameter.
integer, parameter :: rk = 8
real(rk), parameter :: e_min = 1d-5
contains
function CH_SOLVE (a, b, c, f, n) result (x)
implicit none
! FUNCTION COMMENT
! Thomas Chasing Method, solves a tri-diagonal-matrix linear equation set Ax=f
! First does LU-decomposition and then solves the upper- and lower-triangular equations
! ARGUMENTS:
! a(2:n), b(1:n), c(1:n-1) :: elements of the coefficient matrix A. See lecture notes for convention
! f(1:n) :: the inhomogeneous term column-vector in the equation
! n :: size of coefficient matrix
! RETURNS:
! x :: solution to the equation
real(rk), allocatable, intent(inout) :: a(:)
real(rk), allocatable, intent(inout) :: b(:)
real(rk), allocatable, intent(in) :: c(:)
real(rk), allocatable, intent(inout) :: f(:)
integer, intent(in) :: n
! Iteration variable
integer :: i
real(rk), allocatable :: x(:)
print*, n
allocate(x(n))
! LU-decomposition
! because L has only one non-trivial diagonal line, we can multiply L^-1 to f during the process
! See lecture notes for formulae used
do i = 2, n
if (dabs(b(i - 1)) < e_min) then
stop "The matrix is singular!"
endif
! Calculate elements of L
a(i) = a(i) / b(i - 1)
! Calculate elements of U
b(i) = b(i) - a(i) * c(i - 1)
! Multiply L^-1 to f
f(i) = f(i) - a(i) * f(i - 1)
enddo
if (dabs(b(n)) < e_min) then
stop "The matrix is singular!"
endif
! Solve Ux=(L^-1 f) using front-iteration method
! See lecture notes for formulae used
x(n) = f(n) / b(n)
do i = n - 1, 1, -1
if (dabs(b(i)) < e_min) then
stop "The matrix is singular!"
endif
x(i) = (f(i) - c(i) * x(i + 1)) / b(i)
enddo
endfunction CH_SOLVE
function interp (x, y, n, xi) result (yi)
implicit none
! FUNCTION COMMENT
! Interpolates the values of y at points xi based on the given data points (x, y)
! Using Cubic-Spline Interpolation, solving the intermediate equation set with Thomas Method
! ARGUMENTS:
! x, y :: Known data points, should be one-dimensional arrays of same shape (0:n)
! n :: size of given datas points
! xi :: The points at which the interpolation is to be evaluated (assumed to be in increasing order)
! RETURNS:
! yi :: Interpolation results, an array with the same size of xi
real(rk), allocatable, intent(in) :: x(:), y(:), xi(:)
integer, intent(in) :: n
real(rk), allocatable :: yi(:)
! Intermediate variables
! h(0:n-1) :: length of intervals
! M(0:n) :: to be solved, second-order-derivatives of spline function at data points
! mu(2:n-1), l(1:n-2), diag2(1:n-1) :: coefficients in the tri-diagonal-matrix
! d(1:n-1) :: non-homogeneous term in the linear equation set
real(rk), allocatable :: h(:), M(:), mu(:), l(:), diag2(:), d(:)
! Iteration variable
integer :: i, j
print*, n
allocate(yi(size(xi)))
allocate(h(0 : n-1))
allocate(M(0 : n))
allocate(mu(2 : n-1))
allocate(l(1 : n-2))
allocate(diag2(1 : n-1))
allocate(d(1 : n-1))
! Construction of linear equation set for M
! Calculation of h, m & d, see lecture notes for formulae used
h(1) = x(2) - x(1)
do i = 2, n-1
h(i) = x(i + 1) - x(i)
! Used a trick here, calculate d(i) in two parts to save calculation time
d(i+1) = - (y(i + 1) - y(i)) / h(i)
d(i) = (d(i) - d(i + 1)) * 6 / (h(i - 1) + h(i))
mu(i) = h(i - 1) / (h(i - 1) + h(i))
if (i <= n-2) then
l(i) = h(i) / (h(i - 1) + h(i))
endif
enddo
! Do Thomas Method solving
diag2 = 2 ! Diagonal elements of coefficient matrix are all 2
M(1 : n - 1) = CH_SOLVE(mu, diag2, l, d, n-1)
! Boundary conditions
M(0) = 0
M(n) = 0
! Do evaluation
! Here j represents the index of interval where the current xi points is
j = 0
! Interpolation point is invalid if it's smaller than the smallest x
if (xi(1) < x(0) - e_min) then
stop "Invalid interpolation point!"
endif
do i = 1, size(xi)
! Locate the interval where xi(i) is
do while (xi(i) > x(j + 1) + e_min)
j = j + 1
! Invalid interpolation point if xi(i) > x(n)
if (j >= n) then
stop "Invalid interpolation point!"
endif
enddo
! See lecture notes for this formula
yi(i) = (x(j+1)-xi(i))**3*M(j)/(6*h(j))&
+(xi(i)-x(j))**3*M(j+1)/(6*h(j))&
+(y(j)-M(j)*h(j)**2/6)*(x(j+1)-xi(i))/h(j)&
+(y(j+1)-M(j+1)*h(j)**2/6)*(xi(i)-x(j))/h(j)
enddo
endfunction interp
endmodule SplineInterp
program DoSplineInterp
use SplineInterp
implicit none
! PROGRAM COMMENT
! A program to use to above algorithm to interpole the profile of an airplane wing
real(rk), allocatable :: x(:), y(:), xi(:), yi(:)
! Iteration variable
integer :: i
integer, parameter :: n = 9
allocate(x(0:9))
allocate(y(0:9))
allocate(xi(151))
allocate(yi(151))
x = (/0d0, 3d0, 5d0, 7d0, 9d0, 11d0, 12d0, 13d0, 14d0, 15d0/)
y = (/0d0, 1.2d0, 1.7d0, 2d0, 2.1d0, 2d0, 1.8d0, 1.2d0, 1d0, 1.6d0/)
! xi are evenly-distributed points with interval 0.1 between 0 and 15
do i = 1, 151
xi(i) = (i - 1) * 0.1
enddo
yi = interp(x, y, n, xi)
! Write results to file
open (1, file = "Spline.txt", status = "new")
write (1, *) "x:"
write (1, *) x
write (1, *) "y:"
write (1, *) y
write (1, *) "xi:"
write (1, *) xi
write (1, *) "yi:"
write (1, *) yi
close(1)
endprogram DoSplineInterp
