0

我正在尝试计算一个相当复杂的函数，比如 func() - 涉及 fortran 中几个二维数组的几个加法、减法、乘法、除法和三角函数。计算是大规模并行的，因为每个 func() 独立于其行和列位置。每个矩阵的大小都是千兆字节，大约有十几个作为参数。

我想利用英特尔 MKL 函数（调用 --mkl-parallel），特别是 VML 函数来加、减、除等。我的问题是：如何渲染复杂的函数表达式，例如，

例如： func(x,y,z) = x*y+cos(z*xx) 其中 x,y,z 是几个 GB 的二维数组

在 VML 函数方面，但使用更熟悉的二元运算符。您会看到我的问题原则上需要将所有二进制运算符（例如“+”和“*”）转换为以 ?vadd(x,y) 为参数的二进制函数。当然，这对于大型表达式来说会非常麻烦且难看。有没有办法重载二进制算术运算符，例如“+”、“-”，以便在 fortran 中优先使用 MKL/VML 版本。一个例子会很好！谢谢！

4

2 回答 2

1

我知道这个答案有点离题。

由于所有操作都是按元素进行的，并且您的操作很简单，因此这func()可能是一项内存带宽受限的任务。在这种情况下，使用 VML 可能不是最大化性能的好选择。

假设您的每个数组大小为 10GB，如下所示的 uisng VML 将需要至少 9 x 10GB 读取和 5 x 10GB 写入。

func(...) {
    tmp1=x*z
    tmp1=tmp1-x;
    tmp1=cos(tmp1);
    tmp2=x*y;
    return tmp1+tmp2;
}

其中所有操作都为 2d 数组重载。

相反，您可能会发现以下方法具有更少的内存访问（3 x 10GB 读取和 1 x 10GB 写入）因此可能更快（伪代码）。

$omp parallel for
for i in 1 to m
    for j in 1 to n
        result(i,j)= x(i,j)*y(i,j)+cos(z(i,j)*x(i,j)-x(i,j));
    end
end
于 2013-10-23T11:51:20.870 回答
0

我开发了一个小例子来展示两个向量的相加。由于我不再安装 MKL，因此我使用了SAXPYBLAS 的命令。原理应该是一样的。

首先，您使用适当的定义定义一个模块。在我的情况下，这将是在我的数据类型中保存一个真实数组的赋值（这只是一个方便的函数，因为您也可以直接访问该array变量）和加法的定义。两者都是+运算符和=赋值的新重载。

在程序中，我定义了三个字段。其中两个被分配随机数，然后相加得到第三个字段。然后前两个字段存储在我的特殊变量中，并且此添加的结果存储在此类型的第三个变量中。

最后，通过直接访问数组来比较结果。请注意，从自定义数据类型到相同数据类型的分配已经定义（例如ffield3 = ffield1已经定义。）

我的模块：

MODULE fasttype

    IMPLICIT NONE

    PRIVATE

    PUBLIC :: OPERATOR(+), ASSIGNMENT(=)

    TYPE,PUBLIC :: fastreal
        REAL,DIMENSION(:),ALLOCATABLE :: array
    END TYPE

    INTERFACE OPERATOR(+)
        MODULE PROCEDURE fast_add
    END INTERFACE

    INTERFACE ASSIGNMENT(=)
        MODULE PROCEDURE fast_assign
    END INTERFACE

CONTAINS

    FUNCTION fast_add(fr1, fr2) RESULT(fr3)
        TYPE(FASTREAL), INTENT(IN) :: fr1, fr2
        TYPE(FASTREAL) :: fr3
        INTEGER :: L

        L = SIZE(fr2%array)
        fr3 = fr2       

        CALL SAXPY(L, 1., fr1%array, 1, fr3%array, 1)
    END FUNCTION

    SUBROUTINE fast_assign(fr1, r2)
        TYPE(FASTREAL), INTENT(OUT) :: fr1
        REAL, DIMENSION(:), INTENT(IN) :: r2
        INTEGER :: L

        IF (.NOT. ALLOCATED(fr1%array)) THEN
            L = SIZE(r2)
            ALLOCATE(fr1%array(L))
        END IF

        fr1%array = r2
    END SUBROUTINE
END MODULE

我的程序：

PROGRAM main
    USE fasttype
    IMPLICIT NONE

    REAL, DIMENSION(:), ALLOCATABLE :: field1, field2, field3
    TYPE(fastreal) :: ffield1, ffield2, ffield3

    ALLOCATE(field1(10),field2(10),field3(10))
    CALL RANDOM_NUMBER(field1)
    CALL RANDOM_NUMBER(field2)

    field3 = field1 + field2

    ffield1 = field1
    ffield2 = field2

    ffield3 = ffield1 + ffield2

    WRITE(*,*) field3 == ffield3%array
END PROGRAM
于 2013-10-23T13:26:46.113 回答