如何在矩阵维度中使用实数数组索引?

How to use Real Array Index in Matrix dimension?

我正在开发一个 Fortran 代码 (.f90),它将在某个时间间隔 (dt1=0.001) 中计算一些矩阵,并且这些矩阵必须在某些时间步长 (dt=0.1) 中进行积分。尽管我在 FORTRAN 77 方面经验丰富,但对 Fortran 90 还是个新手。我无法使矩阵的维数真实(我认为这就是问题所在,我可能错了!)。下面是一个长程序的一部分。我正在尝试(从过去 1 个月开始)不同的方法但没有成功,输出为 NaN。我在 Ubuntu 16.04

中使用 gfortran
PROGRAM MODEL
IMPLICIT NONE
REAL::DT,DT1,DK2
INTEGER::I,J,IL,IM,E(100000),jm,DK1
REAL::XN2(1),XO2(1)
REAL,DIMENSION(100000)::AN,BN,AN1,BN1
REAL,DIMENSION(21)::XA,XB
REAL,DIMENSION(21,21)::WA,WB,WAB,WBA
REAL,DIMENSION(21)::SUMA,SUMAB,SUMB,SUMBA
XN2=1E-7
XO2=1E-8
OPEN(1,FILE='TEST1.dat')
OPEN(2,FILE='TEST2.dat')
OPEN(3,FILE='TEST3.dat')
OPEN(4,FILE='TEST4.dat')
DO I=1,21  
read(1,*)(WA(I,J),J=1,21)
read(2,*)(WB(I,J),J=1,21)
read(3,*)(WAB(I,J),J=1,21)
read(4,*)(WBA(I,J),J=1,21)
!enddo
IM=1
IL=1
AN(IM)=0.0
BN(IM)=0.0
!AN1(IL)=0.0
!BN1(IL)=0.0
JM=1
E(:)=0.0
DT=0.1
DT1=0.01
DK1=1.0
DK2=1.0
DO WHILE(DK1<=10)
!DO I=1,21
DO WHILE(DK2<=100)
  do j=1,21
  CALL XAA(WA,WAB,SUMA,SUMAB,XN2,XO2,AN(IM),BN(IM),XA)
  CALL XBB(WB,WBA,SUMB,SUMBA,XN2,XO2,AN(IM),BN(IM),XB)
  enddo
  AN(IM+1)=AN(IM)+XA(I)*DT1
  BN(IM+1)=BN(IM)+Xb(I)*DT1
  AN1(IL)=AN(IM+1) 
  BN1(IL)=BN(IM+1)
  AN(IM)=AN1(IL)
  BN(IM)=BN1(IL)
  IM=IM+1
  IL=IL+1
  WRITE(1112,203)I,an(Im),bn(Im)
  DK2=DK2+DT1
ENDDO
203     FORMAT(' ',i2,1000000E13.5)
  E(JM)=DK1
  AN(DK1)=AN(IM)
  BN(DK1)=BN(IM)
  WRITE(1111,203)I,AN(DK1),BN(DK1)
  DK1=DK1+DT
  JM=JM+1
  ENDDO
  ENDDO
END PROGRAM MODEL
!!SUBROUTINES XAA

SUBROUTINE XAA(WA,WAB,SUMA,SUMAB,XN2,XO2,AN,BN,XA)
INTEGER::I,J
REAL,DIMENSION(100000)::AN,BN
REAL,intent(out)::XA(21)
REAL::XN2(1),XO2(1)
REAL,DIMENSION(21,21)::WA,WAB
REAL,DIMENSION(21)::SUMA,SUMAB
SUMA(:)=0.0
SUMAB(:)=0.0
DO I=1,21
DO J=1,21
SUMA(I)=SUMA(I)+WA(I,J)*XN2(1)
SUMAB(I)=SUMAB(I)+WAB(I,J)*XO2(1)
ENDDO
XA(I)=1e-5+SUMA(I)*AN(1)+SUMAB(I)*BN(1)
ENDDO
RETURN
END SUBROUTINE XAA
!!SUBROUTINES XBB

SUBROUTINE XBB(WB,WBA,SUMB,SUMBA,XN2,XO2,AN,BN,XB)
INTEGER::I,J
REAL::XN2(1),XO2(1)
REAL,DIMENSION(100000)::AN,BN
REAL,intent(out)::XB(21)
REAL,DIMENSION(21,21)::WB,WBA
REAL,DIMENSION(21)::SUMB,SUMBA
SUMB(:)=0.0
SUMBA(:)=0.0
DO I=1,21
DO J=1,21
SUMB(I)=SUMB(I)+WB(I,J)*XN2(1)
SUMBA(I)=SUMBA(I)+WBA(I,J)*XO2(1)
ENDDO
XB(I)=1e-5+SUMB(I)*AN(1)+SUMBA(I)*BN(1)
ENDDO
RETURN
END SUBROUTINE XBB

您必须 set/increment IM 和 IL(并在使用前为 IL 赋值)! 在 do 循环中没有任何东西,但肯定必须存在。 数组索引必须始终是整数,如 IM 或 IL,也在列表末尾的写入中,将时间转换为例如如果需要,可以通过 NINT(T/DT)+1 发送 IM。 但首先要使程序能够正常工作,请在循环内为 IM 和 IL 设置适当的值。

根据之前聊天中的评论和对代码的更改:

  • 对于这样一个均匀离散化的问题,很容易(也是最好的)对数组索引的所有操作使用整数,并通过乘以步数和时间离散化来恢复时间的真实值间隔。

  • 在这种情况下,我认为最简单的方法是使用两个嵌套循环:一个在 "long" 集成时间步长上的外部循环,以及一个在 "long" 上执行集成的内部循环=57=] 时间步长。在每个长时间步的开始,您从前一个长时间步复制向量的值,并且对于内部循环中的每个短时间步,您添加您需要添加的内容。

  • 在您的情况下,您可以使用数组语法来避免在主程序和子程序中使用如此多的行和列索引。

  • 请务必在程序开始时读取输入矩阵的所有系数,然后再开始使用它们的值。

  • 如果您要对实数值求积分,您很可能要加减不同数量级的实数。为了准确地做到这一点,您应该确保您的实数具有足够高的精度,例如,使用通过 Fortran 95 SELECTED_REAL_KIND 函数获得的 KIND

  • 为您的子例程提供显式接口,例如,将它们设为内部。然后你可以将数组作为隐式形状传递,这简化了声明。

  • 如果可以为您的子程序的每个参数定义一个意图,请提供它们。如果编译器可以访问意图(通过显式接口),它可以进行更多检查和优化。

  • 不再需要通常的 Fortran 77 将工作空间传递给子例程的做法,对于小型工作空间,它不应该以任何明显的方式影响性能。

  • 最好只定义一次参数,而不是让代码充满 21s。

结果代码是:

PROGRAM MODEL

   IMPLICIT NONE

   INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15,300)

   INTEGER, PARAMETER :: N = 21       ! Dimension of arrays.
   INTEGER, PARAMETER :: Nlong = 1000 ! Number of integrated ("long")
                                      ! time steps.
   INTEGER, PARAMETER :: Nshort = 100 ! Number of short time steps
                                      ! within each long time step.

   REAL(kind=DP), PARAMETER :: DT1 = 0.001_DP    ! Small time step
   REAL(kind=DP), PARAMETER :: DT = Nshort * DT1 ! Long time step
!  REAL(kind=DP), PARAMETER :: T = DT * Nlong    ! Total time

   REAL(kind=DP), DIMENSION(N,N)::WA,WB,WAB,WBA  ! Input arrays.
   REAL(kind=DP) :: XN2(1), XO2(1)               ! Other coefficients. 

   REAL(kind=DP), DIMENSION(N,0:Nlong) :: AN, BN ! Output vectors.

   REAL(kind=DP), DIMENSION(N) :: XA, XB ! Vectors used in integrations.

   INTEGER :: I, J  ! Indices for rows and columns of arrays.
   INTEGER :: IM    ! Index for long time steps.
   INTEGER :: IL    ! Index for short time steps.

   REAL(kind=DP) :: time_now

   ! Set coefficient values.
   XN2=1E-7_DP
   XO2=1E-8_DP

   ! Open input files,
   OPEN(1,FILE='TEST1.dat')
   OPEN(2,FILE='TEST2.dat')
   OPEN(3,FILE='TEST3.dat')
   OPEN(4,FILE='TEST4.dat')
   ! read input arrays
   DO I=1,N
      read(1,*)(WA(I,J),J=1,N)
      read(2,*)(WB(I,J),J=1,N)
      read(3,*)(WAB(I,J),J=1,N)
      read(4,*)(WBA(I,J),J=1,N)
   END DO
   ! and close input files.
   DO I = 1, 4
      CLOSE(I)
   END DO

   ! Initialize output vectors for time=0.
   AN(:,0) = 0.0_DP
   BN(:,0) = 0.0_DP

   ! Loop over long time steps.
   DO IM = 1, Nlong

      ! Initialize vectors from final values at previous time step.
      AN(:,IM) = BN(:,IM-1)
      BN(:,IM) = BN(:,IM-1)

      ! Loop over short time steps.
      DO IL = 1, Nshort

         time_now = (IM-1) * DT + IL * DT1

         CALL XAA(WA,WAB,XN2,XO2,AN(:,IM),BN(:,IM),XA)
         CALL XBB(WB,WBA,XN2,XO2,AN(:,IM),BN(:,IM),XB)

         ! Update vectors.
         AN(:,IM) = AN(:,IM) + XA * DT1
         BN(:,IM) = BN(:,IM) + XB * DT1

         ! Print their values.
         WRITE(1112,'(F12.5,2(1X,i7),42(1X,E13.5))') &
              time_now, IM, IL, AN(:,IM), BN(:,IM)

      END DO

      ! Print integrated values.
      WRITE(1111,'(F12.5,1X,i7,42(1X,E13.5))') &
           time_now, IM, AN(:,IM), BN(:,IM)

   END DO

CONTAINS

   !!SUBROUTINES XAA
   SUBROUTINE XAA(WA,WAB,XN2,XO2,AN,BN,XA)
      ! Arguments
      REAL(kind=DP), DIMENSION(:,:) :: WA,WAB
      REAL(kind=DP), DIMENSION(:), INTENT(IN) :: AN,BN
      REAL(kind=DP), INTENT(OUT) :: XA(:)
      REAL(kind=DP), INTENT(IN) :: XN2(1),XO2(1)

      ! Local variables.
      REAL(kind=DP), DIMENSION(size(XA)) :: SUMA,SUMAB
      INTEGER::I

      DO I=1,N
         SUMA(I)=SUM(WA(I,:))*XN2(1)
         SUMAB(I)=SUM(WAB(I,:))*XO2(1)
      END DO
      XA(:)=1e-5_DP+SUMA(:)*AN(:)+SUMAB(:)*BN(:)

   END SUBROUTINE XAA

   !!SUBROUTINES XBB
   SUBROUTINE XBB(WB,WBA,XN2,XO2,AN,BN,XB)
      ! Arguments
      REAL(kind=DP), DIMENSION(:,:), INTENT(IN) :: WB,WBA
      REAL(kind=DP), INTENT(IN) ::XN2(1),XO2(1)
      REAL(kind=DP), DIMENSION(:), INTENT(IN) ::AN,BN
      REAL(kind=DP), INTENT(OUT) :: XB(:)

      ! Local variables.
      INTEGER::I
      REAL(kind=DP), DIMENSION(size(XB)) :: SUMB,SUMBA

      DO I=1,N
         SUMB(I)=SUM(WB(I,:))*XN2(1)
         SUMBA(I)=SUM(WBA(I,:))*XO2(1)
      ENDDO
      XB(:)=1e-5_DP+SUMB(:)*AN(:)+SUMBA(:)*BN(:)

   END SUBROUTINE XBB

END PROGRAM MODEL