How to use Real Array Index in Matrix dimension?


How to use Real Array Index in Matrix dimension?



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


dt1=0.001


dt=0.1


NaN


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





Comments are not for extended discussion; this conversation has been moved to chat.
– Samuel Liew♦
Jun 21 at 2:57




2 Answers
2



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.


KIND


SELECTED_REAL_KIND



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.


21



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



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.





You mean I can write like below? IM=1 IL=1
– user50695
Jun 20 at 5: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 at 5:58






By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Comments

Post a Comment

Popular posts from this blog

paramiko-expect timeout is happening after executing the command

paramiko-expect timeout is happening after executing the command Can someone help me why timeout is happening after executing the command .I am trying to SSH to a machine and execute a command.I want to store the result of the executed command . I have referred paramiko-timeout but i didn't get my expected result. import traceback try: import paramiko except ImportError: raise Exception("Please install paramiko , pip install paramiko") import pexpect from paramiko_expect import SSHClientInteraction def main(): HOSTNAME = "minion5.net" USERNAME = "root" PASSWORD = "help@432" PROMPT = "root@*:~$s+" try: client = paramiko.SSHClient() client.load_system_host_keys() client.set_missing_host_key_policy(paramiko.AutoAddPolicy()) client.connect(hostname=HOSTNAME, username=USERNAME, password=PASSWORD) with SSHClientInteraction(client, timeout=20, display=True) as interact:
Read more

how to run turtle graphics in Colaboratory

how to run turtle graphics in Colaboratory Given that we can do other types of drawing/visualization in Colaboratory, I would think running turtle graphics would be pretty simple (but not so simple I've figured it out!) I can import the turtle library, but I get the following error when any draw command occurs: TclError: no display name and no $DISPLAY environment variable Thanks in advance. 1 Answer 1 I think P5.js is your best bet. I test it in this notebook. https://colab.research.google.com/drive/1TqF6RObllyslad3KrO4SiAuegRjuJDCF By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.
Read more

Export result set on Dbeaver to CSV

Export result set on Dbeaver to CSV Normally I use Dbeaver for windows and always export my result set like this: run my query --> select the result --> export the result set --> select export to clipboard --> done This step by step puts my result set in my clipboard and I can paste it wherever I want to work with it. The problem is that now I am using dbeaver for mac and this guide is not working. I can go on until the moment that I select my result set like in the image below: exporting data set But once I go further in the process, in the last step I get: no query Note that in "source" it was suppose to show the query that originated the result set, but instead it says just "select. As a result it does't select my result or anything (besides being "successful"). Normally my query would show up there automatically and I couldn't find any option that corrects this problem in the menus. Please help me! Thanks. Aft
Read more