Fortran 函数使用可分配组件重载派生类型之间的乘法

Question

前言

为了存储带状矩阵，其完整对应物可以同时具有从1 以外的索引索引的行和列，我将派生数据类型定义为

TYPE CDS
  REAL, DIMENSION(:,:), ALLOCATABLE :: matrix
  INTEGER, DIMENSION(2) :: lb, ub
  INTEGER :: ld, ud
END TYPE CDS

其中 CDS 代表压缩对角线存储。 鉴于声明TYPE(CDS) :: A，

1. rank-2 组件 matrix 应该包含作为列的实际完整矩阵的对角线（如 here，除了我将对角线存储为列而不是行）。
2. 组件ld和ud应该分别包含上下对角线的数量，即-lbound(A%matrix,2)和+ubound(A%matrix,2)。
3. 2 元素组件lb 和ub 应该包含实际完整矩阵沿两个维度的下限和上限。特别是 lb(1) 和 ub(1) 应该与 lbound(A%matrix,1) 和 lbound(A%matrix,2) 相同。

正如您在第 2 点和第 3 点中看到的，派生类型包含一些冗余信息，但我不在乎，因为它们只是 3 对整数。此外，在我正在编写的代码中，关于实际完整矩阵的边界和带的信息是知道在可以填充矩阵之前。所以我首先将值分配给组件ld、ud、lb和ub，然后我将这些组件用于ALLOCATEmatrix组件（然后我可以正确填充它）。

问题

我必须在这样的稀疏矩阵之间执行矩阵乘法，所以我写了一个FUNCTION 来执行这样的乘积，并用它来重载* 运算符。

目前功能如下，

FUNCTION CDS_mat_x_CDS_mat(A, B)
IMPLICIT NONE
TYPE(CDS), INTENT(IN) :: A, B
TYPE(CDS) :: cds_mat_x_cds_mat

! determine the lower and upper bounds and band of the result based on those of the operands
CDS_mat_x_CDS_mat%lb(1) = A%lb(1)
CDS_mat_x_CDS_mat%ub(1) = A%ub(1)
CDS_mat_x_CDS_mat%lb(2) = B%lb(2)
CDS_mat_x_CDS_mat%ub(2) = B%ub(2)
CDS_mat_x_CDS_mat%ld = A%ld + B%ld
CDS_mat_x_CDS_mat%ud = A%ud + B%ud

! allocate the matrix component
ALLOCATE(CDS_mat_x_CDS_mat%matrix(CDS_mat_x_CDS_mat%lb(1):CDS_mat_x_CDS_mat%ub(1),&
                               & -CDS_mat_x_CDS_mat%ld:+CDS_mat_x_CDS_mat%ud))

! perform the product
:
:

END FUNCTION

这意味着，如果我必须多次执行产品，则分配会在函数中执行多次内部。从性能的角度来看，我认为这并不好。

我请教如何完成带状稀疏矩阵乘以带状稀疏矩阵的任务。我想使用我定义的类型，因为我需要它在边界方面像现在一样通用。但我可以更改执行产品的程序（如有必要，从FUNCTION 更改为SUBROUTINE）。

想法

我可以将过程重写为SUBROUTINE，以便用INTENT(INOUT) 声明CDS_mat_x_CDS_mat，执行matrix 以外的组件的分配以及SUBROUTINE 之外的分配。缺点是我无法重载 * 运算符。

我注意到内部函数 MATMUL 可以对任何 rank-2 操作数进行操作，无论二维的上限和下限。这意味着分配是在函数内部执行的。我想它是有效的（因为它是内在的）。我的函数的不同之处在于它接受任何形状的 rank-2 数组，而我的函数接受具有任何形状的 rank-2 数组组件的派生数据类型对象。

Answer 1

内在函数 MATMUL 等效于自动 (F2008 5.2.2) 结果 - 结果的形状以这样的方式表达，使其成为函数的特征 (F2008 12.3.3) - 形状函数结果在函数的规范部分确定，并且（就实现而言）编译器因此知道如何在实际正确执行函数之前计算函数结果的形状。

因此，没有与内部 MATMUL 函数结果等效的 Fortran 语言 ALLOCATABLE 变量相关联。这与“没有内存分配”不同——编译器可能仍需要在后台分配内存以用于自己的目的——比如临时表达式等。

（我在上面说“等价于”，因为内在过程本质上是特殊的 - 但暂时假设 MATMUL 只是一个用户函数。）

可以通过使用长度类型参数来实现与您的案例类似的自动结果。这是 Fortran 2003 的特性——与引入可分配组件的基本语言标准相同——但它还没有被所有积极维护的编译器实现。

MODULE xyz
  IMPLICIT NONE

  TYPE CDS(lb1, ub1, ld, ud)
    INTEGER, LEN :: lb1, ub1, ld, ud
    REAL :: matrix(lb1:ub1, ld:ud)
    INTEGER :: lb2, ub2
  END TYPE CDS

  INTERFACE OPERATOR(*)
    MODULE PROCEDURE CDS_mat_x_CDS_mat
  END INTERFACE OPERATOR(*)
CONTAINS
  FUNCTION CDS_mat_x_CDS_mat(A, B) RESULT(C)
    TYPE(CDS(*,*,*,*)), INTENT(IN) :: A, B
    TYPE(CDS(A%lb1, A%ub1, A%ld+B%ld, A%ud+B%ud)) :: C

    C%lb2 = B%lb2
    C%ub2 = B%ub2

    ! perform the product.
    ! :

  END FUNCTION CDS_mat_x_CDS_mat
END MODULE xyz

从理论上讲，这为编译器提供了更多优化机会，因为它在调用函数之前对函数结果所需的存储有更深入的了解。这是否真的会带来更好的实际性能取决于编译器的实现和函数引用的性质。