making sure -DFORCE VECTORIZATION always works fine in SPECFEM3D_Cartesian using the CONTIGUOUS Fortran2008 statement

[komatits]

The bug below is probably not fixed yet, but we should check with Daniel Peter because he worked on this a bit in August 2013:

The vectorized version of the code (which is compiled when flag -DFORCE VECTORIZATION is used in flags.guess) can fail on some machines because it replaces multiple indices such as (i,j,k) with a 1D version (ijk,1,1). This works fine only if the array is contiguous in memory. This is always the case for static arrays (as used for instance in SPECFEM3D_GLOBE) but the Fortran norm says that it is NOT necessarily the case when dynamic memory allocation is used. Many compilers will use a continuous memory chunk, but some will not and then the vectorized code of FORCE VECTORIZATION will fail.

In Fortran 2008 this is easy to fix by using the "CONTIGUOUS" option in the call to "allocate()", but some compilers may not support it yet. (many compilers still use F95 or F2003, and only parts of F2008).

Another (less elegant) option would be to allocate a 1D array in the allocate statement in the main program, use it in subroutine calls, and then declare and use the arrays as 3D inside all the routines.

[komatits]

This only applies to SPECFEM3D_Cartesian; SPECFEM3D_GLOBE is fine (because it uses static arrays in the solver), and there is no FORCE VECTORIZATION in SPECFEM2D, which is thus fine as well.

[komatits]

The solution would indeed be to use the "CONTIGUOUS" option of Fortran2008, once it is supported by all compilers we routinely use. From https://software.intel.com/sites/products/documentation/doclib/stdxe/2013/composerxe/compiler/fortran-mac/GUID-1C133644-3883-4146-A6BB-1B948B29EDC7.htm :

"CONTIGUOUS

Statement and Attribute: Specifies that the target of a pointer or an assumed-sized array is contiguous.

The CONTIGUOUS attribute can be specified in a type declaration statement or an CONTIGUOUS statement, and takes one of the following forms:

Type Declaration Statement:

type, [att-ls,] CONTIGUOUS [, att-ls] :: object [, object] ...

Statement:

CONTIGUOUS [::] object [, object] ...

type

Is a data type specifier.

att-ls

Is an optional list of attribute specifiers.

object

Is an assumed-shape array or an array pointer. Description

This attribute explicitly indicates that an assumed-shape array is contiguous or that a pointer will only be associated with a contiguous object.

An entity can be contiguous even if CONTIGUOUS is not specified. An object is contiguous if it is one of the following:

An object with the CONTIGUOUS attribute

A nonpointer whole array that is not assumed-shape

An assumed-shape array that is argument associated with an array that is contiguous

An array allocated by an ALLOCATE statement

An pointer associated with a contiguous target

A nonzero-sized array section in which the following is true:

    Its base object is contiguous.

    It does not have a vector subscript.

    The elements of the section, in array element order, are a subset of the base object elements that are consecutive in array element order.

    If the array is of type character and a substring-range appears, the substring-range specifies all of the characters of the parent-string.

    Only its final reference to a structure component, if any, has nonzero rank

    It is not the real or imaginary part of an array of type complex.

An object is not contiguous if it is an array subobject, and all of the following are true:

The object has two or more elements.

The elements of the object in array element order are not consecutive in the elements of the base object.

The object is not of type character with length zero.

The object is not of a derived type that has no ultimate components other than zero-sized arrays and characters with length zero.

The CONTIGUOUS attribute can make it easier to enable optimizations that rely on the memory layout of an object occupying a contiguous block of memory. Examples

The following examples show valid CONTIGUOUS statements:

REAL, CONTIGUOUS, DIMENSION(:,:) :: A

REAL, POINTER, CONTIGUOUS :: MY_POINTER(:)

" (end of citation from https://software.intel.com/sites/products/documentation/doclib/stdxe/2013/composerxe/compiler/fortran-mac/GUID-1C133644-3883-4146-A6BB-1B948B29EDC7.htm )

[komatits]

To be safe I have turned -DFORCE VECTORIZATION off by default.

[komatits]

We could probably just get rid of FORCE_VECTORIZATION support in the code, since modern compilers can do that automatically if the loops are well written. However, since there is currently full support for FORCE_VECTORIZATION in acoustic (fluid) elements, it could turn out to be useful in its current state in the case of mostly acoustic simulations (e.g. in ocean acoustics).

[danielpeter]

with PR #310, FORCE_VECTORIZATION is implemented in both fluid and elastic domains. it will again be used by default configuration since it improves code performance.

if we still encounter discrepancies between solutions on different systems then we'll turn it off again in the configuration scripts.

[komatits]

It seems Daniel @danielpeter got rid of FORCE_VECTORIZATION; good idea; thus closing this issue.

[Ram1708]

What are the computer system specification need to run the SPECFEM3D.