I have a medium-sized (~650 lines) module, which uses Fortran 2003 OOP features, and PGI generates the following error when I try to compile it:
Lowering Error: array upper bound is not a symbol for datatype 176
Lowering Error: array extnt is not a symbol for datatype 176
PGF90-F-0000-Internal compiler error. Errors in Lowering 2 (/home/santos/CAM/gw_fixes/models/atm/cam/src/physics/cam/linear_1d_operators.F90: 202)
PGF90/x86-64 Linux 14.1-0: compilation aborted
Line 202 is simply the “contains” for the module, so it’s not obvious to me where the real issue is. I’m wondering if there’s a way to narrow this down somehow to work around it. There are a lot of things here that I haven’t done with PGI before, so it’s hard to narrow down.
This code is still very much under development, but I think that in its current state it is still valid Fortran; it looks right to me, and compiles with NAG and Intel compilers. Regardless, there should not be an ICE.
Edit: I want to point out that this can compile stand-alone. So I’m hoping that it’s not hard to reproduce the compilation error I had, even without a full program to actually test it.
Here’s the code itself:
module linear_1d_operators
! TODO: evaluate names, more comments, header here, error handling
! (e.g. need to assert array bounds, and that boundary option is
! valid in routines that use it)
! TODO: check to make sure that the input arrays are big enough for second
! order methods
implicit none
private
save
public :: NeighborOp
! Diagonal matrix.
public :: diagonal_operator
! Derivatives on a non-uniform grid.
public :: first_derivative
public :: second_derivative
! Boundary types.
public :: BoundaryType
public :: BoundaryZero
public :: BoundaryFirstOrder
public :: BoundaryExtrap
public :: BoundaryFixedLayer
public :: BoundaryFixedFlux
! Pre-defined boundary condition options for derivatives.
public :: boundary_zero
public :: boundary_first_order
public :: boundary_extrap
public :: boundary_fixed_flux
! Boundary data types.
public :: BoundaryCond
public :: BoundaryNoData
public :: BoundaryData
public :: BoundaryFlux
! Pre-defined boundary condition.
public :: boundary_no_data
integer, parameter :: r8 = selected_real_kind(12)
! NeighborOp represents operators that can work just between nearest
! neighbors, except with some extra logic at the boundaries. The
! implementation is a tridiagonal matrix plus boundary info.
! nlev is the size of the matrix.
! spr, sub, and diag are the super-, sub-, and regular diagonals.
! lbound_term and rbound_term are extra terms to get a second order
! approximation to a value at the boundaries. They represent second-nearest
! neighbor crossterms for the boundary points.
type :: NeighborOp
integer :: ncol
integer :: nlev
real(r8), allocatable :: spr(:,:)
real(r8), allocatable :: sub(:,:)
real(r8), allocatable :: diag(:,:)
real(r8), allocatable :: lbound_term(:)
real(r8), allocatable :: rbound_term(:)
contains
procedure, pass(self) :: apply_2d => apply_neighbor_to_2d
generic, public :: apply => apply_2d
procedure, pass(self) :: deriv_diag => make_neighbor_deriv_diag
procedure, pass(self) :: lmult_as_diag => diagonal_lmult_neighbor
procedure, pass(op1) :: sum_other_neighbor => sum_neighbor_ops
generic, public :: operator(+) => sum_other_neighbor
end type NeighborOp
interface NeighborOp
module procedure new_NeighborOp
end interface
! Boundary conditions for second-order derivative operators.
!
! BoundaryZero means that the operator always sets boundaries to 0.
! BoundaryFirstOrder means a one-sided approximation for the first
! derivative, and a zero boundary condition for higher orders.
! BoundaryExtrap means that a second order approximation will be used,
! even at the boundaries. Boundary points do this by using their next-
! nearest neighbor to "extrapolate".
! BoundaryFixedLayer means that there's an extra layer outside of the given
! grid, and with data that will be given separately from the mutable
! data when the operator is actually applied.
! BoundaryFixedFlux is intended to provide a fixed-flux condition for
! typical advection/diffusion operators. It tweaks the edge condition
! to work on an input current rather than a value.
!
! At this time, it is not possible to request circular boundary conditions.
type, abstract :: BoundaryType
contains
procedure(make_boundary), deferred :: make_left
procedure(make_boundary), deferred :: make_right
end type BoundaryType
abstract interface
subroutine deriv_seed(del_p, del_m, sub, spr)
import :: r8
real(r8), intent(in) :: del_p(:)
real(r8), intent(in) :: del_m(:)
real(r8), intent(out) :: sub(:)
real(r8), intent(out) :: spr(:)
end subroutine deriv_seed
subroutine make_boundary(self, grid_spacing, seed, &
term1, term2)
import :: BoundaryType, r8, deriv_seed
class(BoundaryType) :: self
real(r8), intent(in) :: grid_spacing(:,:)
procedure(deriv_seed) :: seed
real(r8), intent(out) :: term1(:)
real(r8), intent(out) :: term2(:)
end subroutine make_boundary
end interface
type, extends(BoundaryType) :: BoundaryZero
contains
procedure :: make_left => make_either_boundary_zero
procedure :: make_right => make_either_boundary_zero
end type BoundaryZero
type, extends(BoundaryType) :: BoundaryFirstOrder
contains
procedure :: make_left => make_left_boundary_1st_ord
procedure :: make_right => make_right_boundary_1st_ord
end type BoundaryFirstOrder
type, extends(BoundaryType) :: BoundaryExtrap
contains
procedure :: make_left => make_left_boundary_extrap
procedure :: make_right => make_right_boundary_extrap
end type BoundaryExtrap
type, extends(BoundaryType) :: BoundaryFixedLayer
real(r8), allocatable :: del_edge(:)
contains
procedure :: make_left => make_left_boundary_fixed
procedure :: make_right => make_right_boundary_fixed
end type BoundaryFixedLayer
type, extends(BoundaryType) :: BoundaryFixedFlux
contains
procedure :: make_left => make_left_boundary_fixed_flux
procedure :: make_right => make_right_boundary_fixed_flux
end type BoundaryFixedFlux
type(BoundaryZero), target :: boundary_zero = BoundaryZero()
type(BoundaryFirstOrder), target :: boundary_first_order = BoundaryFirstOrder()
type(BoundaryExtrap), target :: boundary_extrap = BoundaryExtrap()
type(BoundaryFixedFlux), target :: boundary_fixed_flux = BoundaryFixedFlux()
! Data for boundary conditions themselves.
type, abstract :: BoundaryCond
contains
procedure(apply_boundary), deferred :: apply_left
procedure(apply_boundary), deferred :: apply_right
end type BoundaryCond
abstract interface
function apply_boundary(self, op, array) result(output_edge)
import :: BoundaryCond, NeighborOp, r8
class(BoundaryCond), intent(in) :: self
class(NeighborOp), intent(in) :: op
real(r8), intent(in) :: array(:,:)
real(r8) :: output_edge(size(array, 1))
end function apply_boundary
end interface
type, extends(BoundaryCond) :: BoundaryNoData
contains
procedure :: apply_left => apply_left_bndry_extrap
procedure :: apply_right => apply_right_bndry_extrap
end type BoundaryNoData
type, extends(BoundaryCond) :: BoundaryData
real(r8), allocatable :: data_edge(:)
contains
procedure :: apply_left => apply_left_bndry_data
procedure :: apply_right => apply_right_bndry_data
end type BoundaryData
type, extends(BoundaryCond) :: BoundaryFlux
real(r8), allocatable :: delta_edge(:)
contains
procedure :: apply_left => apply_bndry_flux
procedure :: apply_right => apply_bndry_flux
end type BoundaryFlux
type(BoundaryNoData), target :: boundary_no_data = BoundaryNoData()
! Constructor for BoundaryFlux; this is provided to help the user remember
! that the time/space resolution has to be accounted for in constructing
! this condition.
interface BoundaryFlux
module procedure new_BoundaryFlux
end interface
contains
function diagonal_operator(diag) result(op)
real(r8), intent(in) :: diag(:,:)
type(NeighborOp) :: op
op = NeighborOp(size(diag, 1), size(diag, 2))
op%spr = 0._r8
op%sub = 0._r8
op%diag = diag
op%lbound_term = 0._r8
op%rbound_term = 0._r8
end function diagonal_operator
subroutine first_derivative_seed(del_p, del_m, sub, spr)
real(r8), intent(in) :: del_p(:)
real(r8), intent(in) :: del_m(:)
real(r8), intent(out) :: sub(:)
real(r8), intent(out) :: spr(:)
real(r8) :: del_sum(size(del_p))
del_sum = del_p + del_m
sub = - del_p / (del_m*del_sum)
spr = del_m / (del_p*del_sum)
end subroutine first_derivative_seed
function first_derivative(grid_spacing, l_bndry, r_bndry) result(d)
real(r8), intent(in) :: grid_spacing(:,:)
class(BoundaryType), intent(in), optional :: l_bndry, r_bndry
type(NeighborOp) :: d
d = deriv_op_from_seed(grid_spacing, first_derivative_seed, &
l_bndry, r_bndry)
end function first_derivative
subroutine second_derivative_seed(del_p, del_m, sub, spr)
real(r8), intent(in) :: del_p(:)
real(r8), intent(in) :: del_m(:)
real(r8), intent(out) :: sub(:)
real(r8), intent(out) :: spr(:)
real(r8) :: del_sum(size(del_p))
del_sum = del_p + del_m
sub = 2._r8 / (del_m*del_sum)
spr = 2._r8 / (del_p*del_sum)
end subroutine second_derivative_seed
function second_derivative(grid_spacing, l_bndry, r_bndry) result(d)
real(r8), intent(in) :: grid_spacing(:,:)
class(BoundaryType), intent(in), optional :: l_bndry, r_bndry
type(NeighborOp) :: d
d = deriv_op_from_seed(grid_spacing, second_derivative_seed, &
l_bndry, r_bndry)
end function second_derivative
function deriv_op_from_seed(grid_spacing, seed, l_bndry, r_bndry) result(d)
real(r8), intent(in) :: grid_spacing(:,:)
procedure(deriv_seed) :: seed
class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry
type(NeighborOp) :: d
class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc
real(r8) :: del_p(size(grid_spacing,1))
real(r8) :: del_m(size(grid_spacing,1))
integer :: k
if (present(l_bndry)) then
l_bndry_loc => l_bndry
else
l_bndry_loc => boundary_extrap
end if
if (present(r_bndry)) then
r_bndry_loc => r_bndry
else
r_bndry_loc => boundary_extrap
end if
! Number of grid points is one greater than the spacing.
d = NeighborOp(size(grid_spacing, 1), size(grid_spacing, 2) + 1)
! Left boundary condition.
call l_bndry_loc%make_left(grid_spacing, seed, &
d%lbound_term, d%spr(:,1))
do k = 2, d%nlev-1
del_p = grid_spacing(:,k)
del_m = grid_spacing(:,k-1)
call seed(del_p, del_m, d%sub(:,k-1), d%spr(:,k))
end do
! Right boundary condition.
call r_bndry_loc%make_right(grid_spacing, seed, &
d%sub(:,d%nlev-1), d%rbound_term)
! Above, we found all off-diagonals. Now get the diagonal.
call d%deriv_diag()
end function deriv_op_from_seed
subroutine make_either_boundary_zero(self, grid_spacing, seed, &
term1, term2)
class(BoundaryZero) :: self
real(r8), intent(in) :: grid_spacing(:,:)
procedure(deriv_seed) :: seed
real(r8), intent(out) :: term1(:)
real(r8), intent(out) :: term2(:)
term1 = 0._r8
term2 = 0._r8
end subroutine make_either_boundary_zero
subroutine make_left_boundary_1st_ord(self, grid_spacing, seed, &
term1, term2)
class(BoundaryFirstOrder) :: self
real(r8), intent(in) :: grid_spacing(:,:)
procedure(deriv_seed) :: seed
real(r8), intent(out) :: term1(:)
real(r8), intent(out) :: term2(:)
real(r8) :: del_p(size(term1)), del_m(size(term1))
! To calculate to first order, just use a really huge del_m (i.e.
! pretend that there's a point so far away it doesn't matter).
del_p = grid_spacing(:,1)
del_m = del_p * 4._r8 / epsilon(1._r8)
call seed(del_p, del_m, term1, term2)
end subroutine make_left_boundary_1st_ord
subroutine make_right_boundary_1st_ord(self, grid_spacing, seed, &
term1, term2)
class(BoundaryFirstOrder) :: self
real(r8), intent(in) :: grid_spacing(:,:)
procedure(deriv_seed) :: seed
real(r8), intent(out) :: term1(:)
real(r8), intent(out) :: term2(:)
real(r8) :: del_p(size(term1)), del_m(size(term1))
! Use huge del_p, analogous to how left boundary works.
del_m = grid_spacing(:,size(grid_spacing, 2))
del_p = del_m * 4._r8 / epsilon(1._r8)
call seed(del_p, del_m, term1, term2)
end subroutine make_right_boundary_1st_ord
subroutine make_left_boundary_extrap(self, grid_spacing, seed, &
term1, term2)
class(BoundaryExtrap) :: self
real(r8), intent(in) :: grid_spacing(:,:)
procedure(deriv_seed) :: seed
real(r8), intent(out) :: term1(:)
real(r8), intent(out) :: term2(:)
real(r8) :: del_p(size(term1)), del_m(size(term1))
! To extrapolate from the boundary, use distance from the nearest
! neighbor (as usual) and the second nearest neighbor (with a negative
! sign, since we are using two points on the same side).
del_p = grid_spacing(:,1)
del_m = - (grid_spacing(:,1) + grid_spacing(:,2))
call seed(del_p, del_m, term1, term2)
end subroutine make_left_boundary_extrap
subroutine make_right_boundary_extrap(self, grid_spacing, seed, &
term1, term2)
class(BoundaryExtrap) :: self
real(r8), intent(in) :: grid_spacing(:,:)
procedure(deriv_seed) :: seed
real(r8), intent(out) :: term1(:)
real(r8), intent(out) :: term2(:)
real(r8) :: del_p(size(term1)), del_m(size(term1))
integer :: n_space
n_space = size(grid_spacing, 2)
! Same strategy as left boundary, but reversed.
del_p = - (grid_spacing(:,n_space - 1) + grid_spacing(:,n_space))
del_m = grid_spacing(:,n_space)
call seed(del_p, del_m, term1, term2)
end subroutine make_right_boundary_extrap
subroutine make_left_boundary_fixed(self, grid_spacing, seed, &
term1, term2)
class(BoundaryFixedLayer) :: self
real(r8), intent(in) :: grid_spacing(:,:)
procedure(deriv_seed) :: seed
real(r8), intent(out) :: term1(:)
real(r8), intent(out) :: term2(:)
real(r8) :: del_p(size(term1)), del_m(size(term1))
! Use edge value to extend the grid.
del_p = grid_spacing(:,1)
del_m = self%del_edge
call seed(del_p, del_m, term1, term2)
end subroutine make_left_boundary_fixed
subroutine make_right_boundary_fixed(self, grid_spacing, seed, &
term1, term2)
class(BoundaryFixedLayer) :: self
real(r8), intent(in) :: grid_spacing(:,:)
procedure(deriv_seed) :: seed
real(r8), intent(out) :: term1(:)
real(r8), intent(out) :: term2(:)
real(r8) :: del_p(size(term1)), del_m(size(term1))
! Use edge value to extend the grid.
del_p = self%del_edge
del_m = grid_spacing(:,size(grid_spacing, 2))
call seed(del_p, del_m, term1, term2)
end subroutine make_right_boundary_fixed
subroutine make_left_boundary_fixed_flux(self, grid_spacing, seed, &
term1, term2)
class(BoundaryFixedFlux) :: self
real(r8), intent(in) :: grid_spacing(:,:)
procedure(deriv_seed) :: seed
real(r8), intent(out) :: term1(:)
real(r8), intent(out) :: term2(:)
real(r8) :: del_p(size(term1)), del_m(size(term1))
! Use edge value to extend the grid.
del_p = grid_spacing(:,1)
del_m = del_p
call seed(del_p, del_m, term1, term2)
term1 = 0._r8
end subroutine make_left_boundary_fixed_flux
subroutine make_right_boundary_fixed_flux(self, grid_spacing, seed, &
term1, term2)
class(BoundaryFixedFlux) :: self
real(r8), intent(in) :: grid_spacing(:,:)
procedure(deriv_seed) :: seed
real(r8), intent(out) :: term1(:)
real(r8), intent(out) :: term2(:)
real(r8) :: del_p(size(term1)), del_m(size(term1))
! Use edge value to extend the grid.
del_p = grid_spacing(:,size(grid_spacing, 2))
del_m = del_p
call seed(del_p, del_m, term1, term2)
term2 = 0._r8
end subroutine make_right_boundary_fixed_flux
! Constructor for NeighborOp; this just sets the size and allocates
! arrays.
type(NeighborOp) function new_NeighborOp(ncol, nlev)
integer, intent(in) :: ncol, nlev
new_NeighborOp%ncol = ncol
new_NeighborOp%nlev = nlev
allocate(new_NeighborOp%spr(ncol,nlev-1), &
new_NeighborOp%sub(ncol,nlev-1), &
new_NeighborOp%diag(ncol,nlev), &
new_NeighborOp%lbound_term(ncol), &
new_NeighborOp%rbound_term(ncol))
end function new_NeighborOp
function apply_left_bndry_extrap(self, op, array) result(delta_edge)
class(BoundaryNoData), intent(in) :: self
class(NeighborOp), intent(in) :: op
real(r8), intent(in) :: array(:,:)
real(r8) :: delta_edge(size(array, 1))
delta_edge = op%lbound_term*array(:,3)
end function apply_left_bndry_extrap
function apply_right_bndry_extrap(self, op, array) result(delta_edge)
class(BoundaryNoData), intent(in) :: self
class(NeighborOp), intent(in) :: op
real(r8), intent(in) :: array(:,:)
real(r8) :: delta_edge(size(array, 1))
delta_edge = op%rbound_term*array(:,op%nlev-2)
end function apply_right_bndry_extrap
function apply_left_bndry_data(self, op, array) result(delta_edge)
class(BoundaryData), intent(in) :: self
class(NeighborOp), intent(in) :: op
real(r8), intent(in) :: array(:,:)
real(r8) :: delta_edge(size(array, 1))
delta_edge = op%lbound_term*self%data_edge
end function apply_left_bndry_data
function apply_right_bndry_data(self, op, array) result(delta_edge)
class(BoundaryData), intent(in) :: self
class(NeighborOp), intent(in) :: op
real(r8), intent(in) :: array(:,:)
real(r8) :: delta_edge(size(array, 1))
delta_edge = op%rbound_term*self%data_edge
end function apply_right_bndry_data
function new_BoundaryFlux(flux, dt, spacing) result(new_flux)
real(r8), intent(in) :: flux(:)
real(r8), intent(in) :: dt
real(r8), intent(in) :: spacing(:)
type(BoundaryFlux) :: new_flux
new_flux%delta_edge = flux*dt/spacing
end function new_BoundaryFlux
function apply_bndry_flux(self, op, array) result(delta_edge)
class(BoundaryFlux), intent(in) :: self
class(NeighborOp), intent(in) :: op
real(r8), intent(in) :: array(:,:)
real(r8) :: delta_edge(size(array, 1))
delta_edge = self%delta_edge
end function apply_bndry_flux
function apply_neighbor_to_2d(self, array, l_cond, r_cond) result(output)
class(NeighborOp), intent(in) :: self
real(r8), intent(in) :: array(:,:)
class(BoundaryCond), target, intent(in), optional :: l_cond, r_cond
real(r8) :: output(size(array, 1), size(array, 2))
class(BoundaryCond), pointer :: l_cond_loc, r_cond_loc
integer :: k
if (present(l_cond)) then
l_cond_loc => l_cond
else
l_cond_loc => boundary_no_data
end if
if (present(r_cond)) then
r_cond_loc => r_cond
else
r_cond_loc => boundary_no_data
end if
! Left boundary.
output(:,1) = self%diag(:,1)*array(:,1) + &
self%spr(:,1)*array(:,2) + &
l_cond_loc%apply_left(self, array)
do k = 2, self%nlev-1
output(:,k) = &
self%sub(:,k-1)*array(:,k-1) + &
self%diag(:,k)*array(:,k ) + &
self%spr(:,k)*array(:,k+1)
end do
! Right boundary.
output(:,self%nlev) = &
self%sub(:,self%nlev-1)*array(:,self%nlev-1) + &
self%diag(:,self%nlev)*array(:,self%nlev) + &
r_cond_loc%apply_right(self, array)
end function apply_neighbor_to_2d
! Fill in the diagonal for a NeighborOp for a derivative operator, where
! the off diagonal elements are already filled in.
subroutine make_neighbor_deriv_diag(self)
class(NeighborOp), intent(inout) :: self
! If a derivative operator operates on a constant function, it must
! return 0 everywhere. To force this, make sure that all rows add to
! zero in the matrix.
self%diag(:,:self%nlev-1) = - self%spr
self%diag(:,self%nlev) = - self%rbound_term
self%diag(:,1) = self%diag(:,1) - self%lbound_term
self%diag(:,2:) = self%diag(:,2:) - self%sub
end subroutine make_neighbor_deriv_diag
! Sum two NeighborOp objects into a new one; this is just the addition of
! all the entries.
function sum_neighbor_ops(op1, op2) result(new_op)
class(NeighborOp), intent(in) :: op1, op2
type(NeighborOp) :: new_op
! TODO: Add assert that inputs are the same size.
new_op = NeighborOp(op1%ncol, op1%nlev)
new_op%spr = op1%spr + op2%spr
new_op%sub = op1%sub + op2%sub
new_op%diag = op1%diag + op2%diag
new_op%lbound_term = op1%lbound_term + op2%lbound_term
new_op%rbound_term = op1%rbound_term + op2%rbound_term
end function sum_neighbor_ops
! Multiply in a rank 1 array as if it contained the entries of a diagonal
! matrix being multiplied from the left.
subroutine diagonal_lmult_neighbor(self, diag_array)
class(NeighborOp), intent(inout) :: self
real(r8), intent(in) :: diag_array(:,:)
self%spr = self%spr * diag_array(:,:self%nlev-1)
self%sub = self%sub * diag_array(:,2:)
self%diag = self%diag * diag_array(:,:)
self%lbound_term = self%lbound_term * diag_array(1,1)
self%rbound_term = self%rbound_term * diag_array(1,self%nlev)
end subroutine diagonal_lmult_neighbor
end module linear_1d_operators