我为 ARPACK 例程 znaupd 编写了 Fortran 2003 包装器，基本上将示例驱动例程 zndrv1 翻译成具有自动数组的现代 Fortran 2003 语言。我初始化传递给 znaupd 例程的每个变量。只要我在同一程序运行期间不多次调用 zndrv1 方法，一切都会完美运行。它第一次计算正确的解决方案，但第二次导致“未收敛”错误并以零特征对退出。使用编译器标志 -frecursive 编译 ARPACK 并不能解决问题

我的 fortran 2003 子程序是：

子程序 zndrv1(mat,gs_vec,gs_val)

complex(R_P), dimension(:,:), intent(in) :: mat
complex(R_P), dimension(:), intent(inout) :: gs_vec
real(R_P), intent(out) :: gs_val

integer, parameter :: nev = 1, ncv = 10
integer :: iparam(11), ipntr(14)
logical :: select(ncv)
complex(R_P) :: d(ncv), v(size(mat,1),ncv), workd(3*size(mat,1)), workev(3*ncv), resid(size(mat,1))
complex(R_P) :: workl(3*ncv*ncv+5*ncv)
real(R_P) :: rwork(ncv)

character(len=1) :: bmat = 'I'
character(len=2) :: which = 'SR'
integer, parameter :: ishfts = 1, maxitr = 300, mode = 1
integer :: ido = 0, lworkl = 3 * ncv**2+5*ncv, info = 1, n, ldv, ierr = 0
complex(R_P) :: sigma = (0.0,0.0_R_P)
real(R_P) :: tol = 0.0_R_P
logical :: rvec = .true.

iparam = 0
ipntr = 0
select  = .true.
d = (0.0,0.0_R_P)
v = (0.0,0.0_R_P)
workd = (0.0,0.0_R_P)
workev = (0.0,0.0_R_P)
workl = (0.0,0.0_R_P)
rwork = 0.0_R_P
resid = gs_vec
n = size(mat,1)
ldv = n
iparam(1) = ishfts
iparam(3) = maxitr
iparam(7) = mode

do while (ido == 0 .or. ido == 1 .or. ido == -1)

     call znaupd( ido, bmat, n, which, nev, tol, resid, ncv,&
            v, ldv, iparam, ipntr, workd, workl, lworkl,&
            rwork,info )

     if (ido .eq. -1 .or. ido .eq. 1) then

        call av(workd(ipntr(1):(ipntr(1)+n-1)), workd(ipntr(2):(ipntr(2)+n-1)),mat)

     end if
end do

if ( info .lt. 0 ) then

    print *, ' '
    print *, ' Error with _naupd, info = ', info
    print *, ' Check the documentation of _naupd'
    print *, ' '

else

    rvec = .true.

    call zneupd (rvec, 'A', select, d, v, ldv, sigma,&
         workev, bmat, n, which, nev, tol, resid, ncv,&
         v, ldv, iparam, ipntr, workd, workl, lworkl, &
         rwork, ierr)

    if ( ierr .ne. 0) then

        print *, ' '
        print *, ' Error with _neupd, info = ', ierr
        print *, ' Check the documentation of _neupd. '
        print *, ' '

    end if

    gs_vec = v(:,1)
    gs_val = d(1)

!c %-------------------------------------------% !c | 打印附加收敛信息。| ！C ％ - - - - - - - - - - - - - - - - - - - - - -％

     if ( info .eq. 1) then
         print *, ' '
         print *, ' Maximum number of iterations reached.'
         print *, ' '
     else if ( info .eq. 3) then
         print *, ' '
         print *, ' No shifts could be applied during implicit Arnoldi update, try increasing NCV.'
         print *, ' '
     end if

     print *, ' '
     print *, '_NDRV1'
     print *, '====== '
     print *, ' '
     print *, ' Size of the matrix is ', n
     print *, ' The number of Ritz values requested is ', nev
     print *, ' The number of Arnoldi vectors generated (NCV) is ', ncv
     print *, ' What portion of the spectrum: ', which
     print *, ' The number of converged Ritz values is ', iparam(5)
     print *, ' The number of Implicit Arnoldi update iterations taken is ', iparam(3)
     print *, ' The number of OP*x is ', iparam(9)
     print *, ' The convergence criterion is ', tol
     print *, ' '

  end if



contains

    subroutine av ( v, w, mat)

        complex(R_P), dimension(:), intent(in) :: v
        complex(R_P), dimension(:), intent(out) :: w
        complex(R_P), dimension(:,:), intent(in) :: mat

        w = matmul(mat,v)

    end subroutine av

end subroutine zndrv1