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

Question

前言

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

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

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

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

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

问题

我必须在这些稀疏矩阵之间执行矩阵乘法，所以我写了一个 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 以便声明 CDS_mat_x_CDS_mat 和 INTENT(INOUT) 分配 matrix 以外的组件以及分配, 在 SUBROUTINE 之外。缺点是我无法重载 * 运算符。

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

Answer 1

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

因此，没有与内在 MATMUL 函数结果的等效项相关联的 Fortran 语言 ALLOCATABLE 变量。这与 "no memory allocation" 不同——编译器可能仍然需要在幕后为自己的目的分配内存——用于临时表达式等。

（我在上面说 "equivalent of"，因为内在过程本质上是特殊的 - 但暂时假装 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

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

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

Fortran function to overload multiplication between derived types with allocatable components

fortran

operator-overloading

sparse-matrix

matrix-multiplication

derived-types