arrays - Fortran 算法如何使用这些数组？

Question

我正在用 Fortran90 编写一些子例程来执行一些数值计算。但是，作为其中的一部分，我需要包含一些来自 netlib 模板库的代码，这些代码是用 Fortran77 编写的。我很难让他们工作 - 特别是理解数组的用法。

例如，我需要使用一个名为 GMRES 的子程序。以下是论点:

  SUBROUTINE GMRES( N, B, X, RESTRT, WORK, LDW, WORK2, LDW2, 
 $                  ITER, RESID, MATVEC, PSOLVE, INFO )

*     .. Scalar Arguments ..
  INTEGER            N, RESTRT, LDW, LDW2, ITER, INFO
  DOUBLE PRECISION   RESID
*     ..
*     .. Array Arguments ..
  DOUBLE PRECISION   B( * ), X( * ), WORK( * ), WORK2( * )
*     ..
*     .. Function Arguments ..
*
  EXTERNAL           MATVEC, PSOLVE
*
*
*  Arguments
*  =========
*
*  N       (input) INTEGER. 
*          On entry, the dimension of the matrix.
*          Unchanged on exit.
* 
*  B       (input) DOUBLE PRECISION array, dimension N.
*          On entry, right hand side vector B.
*          Unchanged on exit.
*
*  X       (input/output) DOUBLE PRECISION array, dimension N.
*          On input, the initial guess; on exit, the iterated       solution.
*
*  RESTRT  (input) INTEGER
*          Restart parameter, <= N. This parameter controls the amount
*          of memory required for matrix WORK2.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (LDW,5).
*          Note that if the initial guess is the zero vector, then 
*          storing the initial residual is not necessary.
*
*  LDW     (input) INTEGER
*          The leading dimension of the array WORK. LDW >= max(1,N).
*
*  WORK2   (workspace) DOUBLE PRECISION array, dimension     (LDW2,2*RESTRT+2).
*          This workspace is used for constructing and storing the
*          upper Hessenberg matrix. The two extra columns are used to
*          store the Givens rotation matrices.
*
*  LDW2    (input) INTEGER
*          The leading dimension of the array WORK2.
*          LDW2 >= max(1,RESTRT).
*
*  ITER    (input/output) INTEGER
*          On input, the maximum iterations to be performed.
*          On output, actual number of iterations performed.
*
*  RESID   (input/output) DOUBLE PRECISION
*          On input, the allowable error tolerance.
*          On ouput, the norm of the residual vector if solution
*          approximated to tolerance, otherwise reset to input
*          tolerance.
*
*  INFO    (output) INTEGER
*          =  0:  successful exit
*          =  1:  maximum number of iterations performed;
*                 convergence not achieved.
*
*  MATVEC  (external subroutine)
*          The user must provide a subroutine to perform the
*          matrix-vector product A*x = y.
*          Vector x must remain unchanged. The solution is
*          over-written on vector y.
*
*          The call is:
*
*             CALL MATVEC( X, Y )

注意数组是如何定义的 WORK( * ), WORK2( * )。所以在我看来，这些是任意长度的一维数组。但是在参数描述中，似乎暗示它们是二维数组，矩阵，维度为 WORK(LDW, 5)。那么它们是一维的还是二维的？

另外，在GMRES算法中，这些WORK数组是这样使用的:

CALL MATVEC(SCLR1, WORK(NDX1), SCLR2, WORK(NDX2))

因此，如果 WORK 数组是二维的，上面的内容不会给出秩不匹配吗？访问一个只有一个索引的二维数组是什么意思？我应该将 WORK 数组定义为 2D 还是 1D？

编辑

GMRES 例程需要提供一个 matvec 例程。在 GMRES 代码中，它被称为

CALL MATVEC(SCLR1, X, SCLR2, WORK(NDX2))

也喜欢

CALL MATVEC(SCLR1, WORK(NDX1), SCLR2, WORK(NDX2))

我尝试提供的子程序 MATVEC 如下所示:

subroutine matvec(alpha, x, beta, y)

    real(dp), intent(in) :: alpha, beta
    real(dp), dimension(*), intent(in) :: x
    real(dp), dimension(*), intent(inout) :: y
    real(dp), dimension(*,*) :: A
    integer :: n

    call make_Jac(n,A)
    call dgemv('notranspose', n, n, alpha, A, n, x, 1, beta, y, 1)

end subroutine matvec

make_Jac 返回我正在处理的问题的矩阵及其维度 n。 blas 例程 dgemv 然后处理矩阵向量积。

Answer 1

WORK( * ) 声明一个假设大小的数组，它可以是二维的。参见 here .

如果您向子例程提供一维数组，编译器不会报错，但可能会发生奇怪的事情(包括段错误)。

最好使用符合规范的数组。