Fortran数组的C++实现

最近在写CFortranTranslator，一个从Fortran77/Fortran90到C++14的工具，其中涉及到使用C++为Fortran实现一个数组库。
总的来讲，Fortran90的有些语言特性在编写和编译上都让人不是很舒服。比如没有头文件和前置声明，这导致了许多额外的代码和解析工作，比如INTERFACE块因此而生（但其实也是一个比较好的解决方案，和前置声明也差不多）。又比如Fortran兼容老标准的问题，这个写C++的同学也应该深有同感，C++为了保持和C的linkage做了不少擦屁股的事情，像什么POD、函数指针/函数对象啥的。Fortran老标准中允许隐式声明变量，而且可以通过名字推断变量的类型，而且由于COMMON块的存在还要兼容一些奇妙的用法。这使得处理变量声明的工作要延迟到处理函数体时。Fortran90还可以直接根据Attribute specification statements给变量加属性(ISO/IEC 1539 : 1991 ch5.2)，加上隐式声明的情况，实际上类似于INTENT、PARAMETER、DIMENSION这样的语句需要考虑是生成一段新变量声明还是修改老声明，这同样会延迟处理变量声明的工作，还会要求建立符号表。此外各种impied-do结构，虽然可以完全展开成for循环，但是为了保持Fortran源语句的抽象，最好还是做成一个表达式。还有Fortran77里面的nonblock DO construct非常恶心，从语法上完全无法解析了，最偷懒的办法就是先用start condition给flex开洞，再#define YYLEX给bison开洞。不过后来token级的continuation迫使我直接手写词法分析了，果然偷懒还是要不得的。
此外Fortran的传参机制也很特别，类似于宏的形式，基于引用，具体类型和限定要在函数体中才能看到。为了能够兼容语义，我全部使用了右值作为参数限定
关于C++元编程的部分可以参考C++模板编程这篇文章中，关于使用flex/bison进行语法分析的部分我放到了flex和bison使用，其他的部分放在这里。

Fortran数组简介

Fortran语言的实现常常要比C快，其中数组是功不可没的。
在C语言中，数组常常会被decay成指针来处理，而这会在优化时造成pointer aliasing的问题。
一方面C/C++编译器遵循strict aliasing规则，即不同类型的指针绝不会指向同一块内存。具体的说下面的情况被视为相同类型：

1. 被signed/unsigned/const/volatile所修饰
2. 被一个聚合或联合所包含
3. 继承自一个基类
4. 类型char*和void**可以alias任何其他类型
   这样当你写出这样的神奇代码的时候，编译器至少能有个warning

   int a;
   int f(float *b)
   {
       a = 1;
       *b = 0;
       return a;
   }
   int main()
   {
       printf("%d\n", f((float*)&a));
       return 0;
   }

在另一方面，我们还是无法容易得知道两个相同类型的指针是否alias，即指向同一块内存。这时候Pointer aliasing的问题会阻碍编译器进行优化，虽然restrict关键字可以一定程度上减少这个问题。Fortran数组就彻底没有这样的烦恼。

for1array<typename T>

一开始的想法是实现一个for1array<T>，一个一维的数组。可以理解为对C++原生数组的一个包装，使得数组能够自定义上下界，和列优先的存储模式。事实上这种嵌套(nested)数组在实现Fortran的内在函数时存在相当大的麻烦。
例如将初始化序列按列映射到一个嵌套的数组for1array<for1array<...for1array<T>...>> farr中，如果通过通常的递归来做，那么就要自结构最里面从内而外for1array<T>进行递归创建这个数组，但这难于实现。一个可能的解决方案对于嵌套的for1array数组按照通常一样从外到内进行递归遍历，而计算farr[a1][a2]...[an]时对应的一维序列中的位置。对于farr数组中的第k层的指针增加1，实际上相当于一维序列中向后移动size[1] * size[2] * ... size[k-1]，其中size[i]为farr中第i层的大小。
但是对于transpose这样的函数使用嵌套数组实现的代价就相当大了，特别地，Fortran90标准对转置秩大于2的数组并没有规定行为：

13.13.109 TRANSPOSE (MATRIX)
Description. Transpose an array of rank two.
Class. Transformational function.
Argument. MATRIX may be of any type and must have rank two.

因此实际上对于一个N维/秩数组$X$，假设第$i$维的取值范围是$[0, n_i]$，定义转置$X^T$，满足$X[a_1][a_2]…[a_n]$ = $X^T[a_n][a_{n-1}]…[a_1]$。

farray<typename T, int D>

Fortran90标准中数组的维数，上下界都是确定的，称为显式形状数组(explicit-shape)。此外，除了动态数组(deferred-shape)，即使出现的可调数组(automatic explicit-shape)（类似于C99中的Flex Array）、假定形状数组(assumed-shape)（仅给出维/秩数）和假定大小数组(assumed-size)（最后一维的上下界是不定的）也是作为哑元（形参的）。Fortran77标准更是使用完全静态环境(fully static environment)的运行时（《编译原理及实践》§7），因此都不支持递归调用。因此Fortran多维数组并不需要锯齿数组，也不需要实现类似动态表的自扩，所以内部可以通过一个一维的线性表data来保存数据。相对于C的数组而言，Fortran中数组是按照列优先顺序(column-major order)存储的。在实现数组时，应当考虑空间局部性的问题，将同列的元素尽可能一同存储（注意C++是行优先顺序）。在Fortran数组的内部实现中如果出现嵌套循环、递归，也应当保证内层循环/递归处理同一列，对于多维数组则是遍历最左边下标。

编译期维护数组各维度上下界的尝试

因为要尽可能保证和源代码一致性，所以不能将上界统一改为0开始，在设计时需要使用在对象中维护每一维的上下界。于是自然想到是否可以在编译期进行上下界的相关计算，一个很挫的思路是将每一维的上下界放入模板参数里面，然后通过parameter unpack提取出来。另外还有想法是借助constexpr而不是模板，例如使用constexpr构造函数（注意目前使用的MSVC2015无法编译，原因是尚不支持Extended Constexpr）。但事实上即使能编译，constexpr也不是一个强制性声明，编译器仍然可以选择将其放到运行期完成。其实constexpr也被用来“糊弄”编译器。
StackOverflow上的这篇回答解释了这样做的原因。回答列举了一个希望使用for遍历std::tuple的例子，我们知道C++是静态类型的且std::tuple的参数包中可能包含不同类型，所以for (const auto & x : my_tuple)这样的做法就是在和C++标准过不去；但使用一个index来“遍历”std::tuple中的元素也是不行的，例如for (constexpr i = 0; i < 10; ++i)，程序员完全可以将这个i作为std::get的非类型模板参数来构造出可能属于不同类型的变量。
通过答主的这个例子可以知道一个constexpr的for对当前的C++来说是不切实际的，事实上我们可以采用一些方法可以实现编译期的循环结构，比如变循环为递归，把index放到模板参数里面。对于C++14，可以借助于index_sequence来辅助进行包展开，对于C++17可以去fold一个index_sequence。

实现transpose

现在实现transpose函数，transpose函数主要是通过变换$n_i$和对应的一维序列，满足上面式子。
对于原矩阵$X[a_1][a_2]…[a_n]$，可以映射到$X$对应的一维序列的第$\sum_{i=1}^{N}{(a_i \times \prod_{j = i + 1}}^{N}{n_j})$项，令$\prod_{N + 1}^{N}{} = 1 $；
对于转置后的矩阵$X^T[b_1][b_2]…[b_n]$，令其每一维大小为$n’_i $，可得$n’_i = n_{N + 1 - i}$，可以映射到$X^T$对应的一维序列的第$\sum_{i=1}^{N}{(b_i \times \prod_{j = i + 1}}^{N}{n’_j})$项，其中令$\prod_{1}^{0}{} = 1 $。
根据$X[a_1][a_2]…[a_n]$ = $X^T[a_n][a_{n-1}]…[a_1]$，带入$b_i = a_{N + 1 - i}$，可以得到映射规则$X: \sum_{i=1}^{N}{(a_i \times \prod_{j = i + 1}}^{N}{n_j})$ $\Rightarrow$ $X^T: \sum_{i=1}^{N}{(b_i \times \prod_{j = i + 1}}^{N}{n’_j})$ = $\sum_{i=1}^{N}{(a_{N + 1 - i} \times \prod_{j = i + 1}}^{N}{n_{N + 1 - j}})$。
以上是对于C-style的数组讨论的，对于Fortran-style的数组，将$\Sigma$的上下界改成$0 .. i - 1$
此外，对于2维非方阵矩阵转置有$O(1)$空间的方法，主要是按照序号从小到大遍历矩阵中的每个元素，希望能够找到以它开始的环。容易看出一个环没有被遍历过当前仅当当前序号是环中最小的序号。

slice数组片段

fortran中的slice是by reference的，对slice的修改也是对原数组的修改，于是要实现一个forslice函数，其返回的slice应当是数组片段的引用。对于这种情况，首先可以考虑使用std::vector<std::shared_ptr<T>> data来实现内部的线性表data，但是这样会产生不必要的空间开销。因为虽然每一个slice是不一样的，没办法对std::shared_ptr<std::vector<T>> data，但是每个slice的data部分是作为一个整体的。所以还是使用的裸指针RAII，对于从forslice函数构造的farray<T>，在构造时加一个tag，析构时不释放指针即可，类似view的行为。考虑到Fortran77标准本身的完全静态环境（Fortran90可能是基于栈的），所以不会出现悬挂指针导致的AV错误。
此外fortran中还允许使用(/ /)这样的Delimiters来具体指定位置组成数组片段，例如A((/1,2/))表示选取A中的第1和第2个元素。要实现这个功能只需要给slice_info做一些修改，加一个迭代器即可，考虑到这个功能并不常用，所以并没有实现。

数组的秩

此外fortran中有使用到数组的秩作为参数（通常参数名为dim）的函数，要注意秩是从1开始计算的。

farray<typename T>

fortran有一个slice操作，它返回的维数是不固定的，例如a(:, 1, 1:2)返回一个2维的数组，但是a(:, 1, 1:1)返回一个三维的数组，虽然这两个数组所承载的数据是一样的。这给实现带来了困难，因为返回类型是根据slice_info的初始化列表来决定的，为了追求语法尽可能简单，贴近fortran源码，对slice_info这样的工具类模板开洞是不太好的，因此决定将维数参数也移出模板参数。
因此将farray<T>的数据分成两块，第一块是有关形状的元数据，通过reset_array来设置；另一块是数据本身，通过reset_value来设置。在数组声明（构造函数）时，可以确定数组形状，可以顺带赋值；在数组（operator=）赋值时只能赋值，而不使用目标右值的形状。
再一次说明，这样带来了额外的空间开销lb用来存储上界，sz用来存储大小。由于fortran数组是按列存储，越往后的下边周期越长，因此预先计算delta也就是第d维是i时数组的长度，用来辅助运算。可以看到是把整个数组的运算都放到了运行期，而这些东西完全是编译器决定的。显然我们可以在把fortran编译到C++的阶段处理这些东西，但是要不这让生成的代码因为带上这些附加的数据显得非常难看，要不就会损失掉fortran原来的语义，例如将fortran的数组整合成上界始终是0，按行存储的数组。所以唯一的办法是利用元编程，让C++编译器在模板阶段处理这些东西。

fortran数组的内在函数

处理dim参数

minloc、maxloc、all、any等函数可以传dim参数。这个dim参数实际上就是数组的秩(rank)。fortran中数组的秩类似于C++中数组的维数(dimension)，但是秩是从1开始算的（可是参数叫dim呢）。
在fortran90标准中并没有签名为RESULT = MAXLOC(ARRAY, DIM [, MASK]) ，从gfortran的文档中了解到这个函数时从95标准后新增加的。在实现这个函数的时候产生了困惑，对于高维数组，这个maxloc函数到底返回什么呢？gfortran文档中只给出：

if the DIM argument is supplied, determines the locations of the maximum element along each row of the array in the DIM direction.
If DIM is present, the result is an array with a rank one less than the rank of ARRAY, and a size corresponding to the size of ARRAY with the DIM dimension removed.

后来我在这里找到了一个例程：

integer :: i6(6) = (/-14,3,0,-2,19,1/)
integer :: i23(2,3) = reshape((/-14,3,0,-2,19,1/),shape(i23))
write(*,'(2i4)') i23         ! writes -14   3
							!          0  -2
							!         19   1
write(*,*) maxloc(i6)        ! writes 5
write(*,*) maxloc(i23)       ! writes 1  3
write(*,*) maxloc(i23,dim=1) ! writes 2  1  1
write(*,*) maxloc(i23,dim=2) ! writes 3  1
write(*,*) maxloc(i23,dim=1,mask=(i23 < 10))
							! writes 2  1  2

maxloc(a, dim=1)这个函数返回的是[lb[dim], sz[dim]]构成的dim-1维数组，令剩下来的维数取遍所有组合，对于每一种取组合，给出取得最大值时候对应原数组第dim维的下标，即$ \underset{i_{dim}}{\arg\min}(b(i_{dim})), \forall b = a(i_1, …, i_{dim - 1}, :, i_{dim + 1}, …, i_n)$
这个例子中，列是第一维，行是第二维。
于是对于此类函数可以抽象成reduce函数，对于不带dim参数的重载版本，这是对于原数列的reduce操作；对于带dim的版本，这是对于一个rank-1的数组中的每个元素（是一个一维数组）进行reduce操作。

处理mask参数

minloc、maxloc、all、any等函数还可以传mask参数。这个参数是一个逻辑表达式，实际上类似于一个谓词，并且这个谓词需要使用当前作用域内的变量。至此为止用实现一个[&](){/* expressions directly translated from fortran */}的lambda表达式即可。比较麻烦的是这lambda的函数体内部直接使用数组名，而谓词是作用于数组的元素上的。并且在整个函数调用的可见范围内，无法判断lambda函数体中的变量究竟是数组还是标量，例如下面的情况：

integer,dimension(8)::N1, N2
integer::N3
print *, maxval(N2, mask = N1 < 5) ! N1 is array
print *, maxval(N2, mask = N3 < 5) ! N3 is scalar

由于编译器实际上并没有建立符号表，在处理maxval(N2, mask = N1 < 5)这个调用时并不能知道N1，N2，N3的类型与定义位置。所以也就无法将数组改为对应的元素形式了。
有两种方法解决这个问题：
第一种，将这个谓词作为表达式来计算，最后返回一个布尔值数组。这需要重载运算符和fortran的内在函数。
第二种，是实现一个模板函数formap(F f, T x)和重载版本formap(F f, farray<T> x)，然后对于所有的变量x，用formap(predicate, x)包一下就好了

fortran数组中的运算符重载

fortran中数组可以对同样长度的数组和标量进行算术运算和比较运算等，所以需要实现一系列重载函数

数组声明

数组声明注意点比较多。首先给出fortran中声明(declare)和定义(define)的区别

The term declaration refers to the specification of attributes for various program entities. Often this involves
specifying the data type of a named data object or specifying the shape of a named array object.
The term definition is used in two ways. First, when a data object is given a valid value during program
execution, it is said to become defined. This is often accomplished by execution of an assignment statement or
input statement. Under certain circumstances, a variable does not have a predictable value and is said to be
undefined. Section 14 describes the ways in which variables may become defined and undefined. The second
use of the term definition refers to the declaration of derived types and procedures.

简而言之，声明(declaration)，指的是设置变量的type或者attribute，而定义(definition)是给变量初始化。

下面考虑数组的声明方式。首先fortran77中可以使用下面的方式声明

REAL        A(10,2,3)             ! 类型说明（常用）
DIMENSION   A(10,2,3)             ! DIMENSI0N语句（常用）
! COMMON

fortran90中可以使用下面的方式声明

REAL，DIMENSION(2,5):: D          ! 类型说明中的DIMENSION属性（最常用）
REAL，ALLOCATABLE::    E(:,:,:)   ! 类型说明中的ALLOCATABLE属性
REAL，POINTER::        F(:,:)     ! 类型说明中的POINTER属性
POINTER     C(:,:,:)              ! POINTER语句
ALLOCATABLE B(:,:)                ! ALLOCATABLE语句

根据fortran90标准section5

R501 type-declaration-stmt is type-spec [ [ , attr-spec ] … :: ] entity-decl-list
entity-decl is object-name [ ( array-spec ) ] [ * char-length ] [ = initialization-expr ]
or function-name] [ * char-length ]

其中array-spec就是dimension()语句括号中的东西，表示数组的形状。
注意：

1. 声明中并不一定要出现double colon separator::，但是当有初始化语句initialization-expr时，必须要有这个分隔符
   这在语法分析时可能出现语义冲突导致二义性文法，原因是integer、real等既可以作为类型限定符，也可以作为函数名，所以实际上要去掉callable_head规则，将它统一合并入callable规则
2. char-length仅允许给字符串指定长度，由于实现上使用的是std::string，所以可以直接去掉这个，我们只需要对变量声明中出现的*号进行特判处理即可。

考虑到可以通过REAL R(10) = (/1,2,3,4,5,6,7,8,9,0/)声明数组，所以不能在扫描完type-spec和attr-spec就确定类型，而要根据每一个entity-decl后面是否出现( array-spec )来决定。

integer::Z
dimension:: Z(10)
! dimension(10)::Z 注意这样的语句是不行的

这样的说明也是允许的，考虑到fortran还有隐式声明，事实上第一行还可以直接去掉。所以说即使整个声明语句读下来，仍然不能决定具体的类型是什么：

5.1.2.4 DIMENSION attribute
The DIMENSION attribute specifies that entities whose names are declared in this statement are arrays. The
rank or the rank and shape are specified by the array-spec, if there is one, in the entity-decl, or by the array-spec
in the DIMENSION attribute otherwise. An array-spec in an entity-decl specifies either the rank or the rank and
shape for a single array and overrides the array-spec in the DIMENSION attribute. If the DIMENSION attribute
is omitted, an array-spec must be specified in the entity-decl to declare an array in this statement.
为了解决这个问题，变量的声明必须要在整个变量声明语句块结束之后才能确定（fortran的命令语句必须在所有声明语句后面）。所以只有在program或者function_decl规则中才能对其中的suite规则中的变量名生成定义和初始化语句。
这里还要注意一点，变量的隐式声明和隐式类型(implicit type indicated by the first letter)还不一样，后者指的是在声明变量的时候不指出变量的类型，这时候根据变量名的第一个字母来决定是实型还是整型。当然隐式声明和隐式类型是可以同时存在的。

数组初始化和赋值

使用数组构造器

Fortran中的数组构造器(array constructor)我认为是一个很方便实用的功能，根据Fortran90标准，数组构造器使用下面的语法

array-constructor is (/ ac-value-list /)
R432 ac-value is expr or ac-implied-do
R433 ac-implied-do is ( ac-value-list , ac-implied-do-control )
R434 ac-implied-do-control is ac-do-variable = scalar-int-expr , scalar-int-expr [ , scalar-int-expr ]
R435 ac-do-variable is scalar-int-variable

数组构造器始终是一个一维的数组，不能使用数组构造器直接初始化高维数组，而应该使用reshape函数。又注意数组构造器在赋值和构造时都可能用到，所以比较方便的方法是对(/ /)直接生成一维的farray实体而不是将它直接生成代码到构造函数中作为一个brace-init-list的参数：

An array constructor is defined as a sequence of specified scalar values and is interpreted as a rank-one array
whose element values are those specified in the sequence.

因为farray的构造函数用来确定数组的形状（顺带赋值是可选的），而赋值放到farray<T>::operator=去做。在实现时因为farray<T>::operator=仅仅是复制内部的线性表data，并不改变形状，可以在不同维数之间给数组赋值。

implied-do

在上面的定义中ac-implied-do是个特别有意思的东西，与之类似的还有io-implied-do（用于read等io语句）和data-implied-do（用于data语句），称为隐式Do循环(Implied Do Loop)。
标准中对循环的特性做出了这样的规定

If an ac-value is an ac-implied-do, it is expanded to form an ac-value sequence under the control of the ac-do-variable, as in the DO construct (8.1.4.4).

以ac-implied-do举例，隐式Do循环的的语法通常为

ac-value ::= expr | ac-implied-do
ac-implied-do ::= ( ac-value-list , scalar-int-expr , scalar-int-expr [ , scalar-int-expr ] )

在转换中为了保留隐式DO循环的结构，所以我们尽可能避免将这个Implied Do展开为C++中的for循环语句。
第一个思路是把整个嵌套的Implied Do深度优先到底，得到最里面的表达式放入一个[&](int do_variable1, int do_variable2, ...){...}里面，然后借助于farray里面的map函数喂给这个lambda块各层do-variable的值。不过由于ac-value-list可能包含不止一个表达式，所以这个方法实际上麻烦。
于是有了第二个思路，就是按部就班一层一层处理。但这种方法仍然是要在Implied Do的上下文中处理，因为内层Do的循环过程需要依赖于外层Do对应的do-variable的当前取值。

使用DATA

data语句可以用来初始化数组，可以采用下面的语法

data a1/1.,2.,3.,4./, a2/2.,3.,4./
data b/3*1.,4*2.,3*3./
data c(1),c(2),c(3),c(4)/1.,2.,3.,4./
data (d(j),j=1,4)/1.,2.,3.,4./

其中值得注意的是第二条语句中的*并不是二元算术运算符，而是表示该元素的重复次数

使用WHERE

生成AST

会碰到这样的问题，例如翻译上面的数组生成器，可以生成一个`NT_ARRAYBUILDER`节点，所以`NT_ARRAYBUILDER`也是一个`NT_EXPRESSION`节点，而数组生成器是可以出现的等号右边构成表达式，于是也可以作为赋值式归约成`NT_EXPRESSION`节点。所以问题存在于直接将`NT_ARRAYBUILDER`立即归约成表达式`NT_EXPRESSION`节点还是当`NT_ARRAYBUILDER`和新的串构成例如赋值式时再归约成表达式。 对于第一种方法，会在AST的深度方向产生较多的**不确定深度**的节点，例如在`NT_FUNCTION_ARRAY`会被归约成`NT_EXPRESSION`。出现这种问题主要是我在列表方面有点偷懒，只设置了一个`paramtable`，一般非终结符必须要归约成`NT_EXPRESSION`才行，此外，只要列表中含有slice，所有元素都被promote成slice；进一步地，只要列表中含有键值对，所有的元素都被promote成键值对。如果修改语法，可能也要修改语义处理程序。其次是因为对于所有的运算符（除了赋值），我都将其设置为一个`NT_EXPRESSION`节点 对于第二种方法，容易在单独使用`NT_ARRAYBUILDER`的时候出现问题，例如使用`if(X == NT_EXPRESSION)`筛选`NT_EXPRESSION`节点时会因为`NT_ARRAYBUILDER != NT_EXPRESSION`造成被误筛掉

总结

Fortran数组相对于C++数组还是有很大的不同的，正因为如此，Fortran数组具有比C++要稍好的（数值计算）性能（如果能够正确使用的话）。例如Fortran的数组模型非常适合编译器针对处理器做矢量优化，而到了C++中数组往往会退化成一个指针，这可能会妨碍编译器了解内存布局，从而进行优化。此外就指针本身，Fortran也没有C++的Pointer Alias的开销。