-2

I am developing a Fortran code (.f90) which shall calculate some matrix in some time interval (dt1=0.001) and these matrices have to be integrated in some time steps (dt=0.1). Though I am experienced in FORTRAN 77, new to Fortran 90. I am unable to make dimension of matrix real (I think that is the problem, I may be wrong!). Below is the part of a long program. I am trying (from last 1 month) different methods but no success, output are NaN. I am using gfortran in Ubuntu 16.04 

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
user50695
  • 3
  • 2

2 Answers2

0

You have to set/increment IM and IL (and assign somewhere a value to IL prior using it)! Within the do loops there is nothing of this but certainly must be there. The array indices have always to be the integers like IM or IL also in the write at the end of the listing, convert time to e.g. IM by NINT(T/DT)+1 if required. But first for the program to be able to work properly set proper values to IM and IL inside the loops.

Michael
  • 42
  • 3
  • You mean I can write like below? IM=1 IL=1 – user50695 Jun 20 '18 at 05:52
  • You mean I can write like ? IM=1 IL=1 ....AN1(I,IL)=AN(I,IM+1) BN1(I,IL)=BN(I,IM+1) AN(I,IM)=AN1(I,IL) BN(I,IM)=BN1(I,IL) with IL=IL+1 and IM=IM+1 (or IL=IL+NINT(DK2/DT)+1 and IM=IM+NINT(DK2/DT)+1) inside inner loop and keep same as earlier in outer loop. – user50695 Jun 20 '18 at 05:58
0

Further to previous comments in the chat and changes to the code:

  • For such a uniformly discretized problem, it is very easy (and best) to use integers for all the manipulations of array indices, and recover real values of time by multiplying numbers of steps and time discretization intervals.

  • In this case, I think it is easiest if you use two nested loops: an external loop over "long" integrated time steps, and an internal one where you perform the integration over "short" time steps. At the beginning of each long time step you copy the value of your vectors from the previous long time step, and for each short time step in the inner loop you add what you need to add to them.

  • In your case, you can use array syntax to avoid so many row and column indices, both in the main program and in the subroutines.

  • Please do read all the coefficients of your input matrices at the beginning of your program, before starting to use their values.

  • If you are going to integrate real values, you are likely going to add or subtract real numbers with different orders of magnitude. In order to do this accurately, you should ensure a high enough precision in your reals, e.g., using a KIND obtained via the Fortran 95 SELECTED_REAL_KIND function.

  • Provide explicit interfaces for your subroutines, e.g., by making them internal. Then you can pass arrays as implicit-shape, which simplifies the declarations.

  • If an intent can be defined for each of the arguments of your subprograms, please provide them. If the compiler has access to intents (via explicit interfaces) it can do more checks and more optimisations.

  • The usual Fortran 77 practice of passing workspaces to subroutines is no longer needed, and for small workspaces it should not affect performance in any noticeable way.

  • It is good practice to define parameters only once, rather than having code filled with 21s.

The resulting code is:

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
jme52
  • 1,123
  • 9
  • 18