Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Highlighted
Novice
1,338 Views

Problem with variable that change their values without apparent reason

I would propose a part of my code:INTEGER :: n,a,b
REAL*8 :: c
LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: pop, new_pop
LOGICAL, ALLOCATABLE, DIMENSION(:,:,:;:) :: position
!
OPEN(10, FILE = 'printout.txt')
n = 150
ALLOCATE(pop(1:n,1:n))
ALLOCATE(new_pop(1:n,1:n))
ALLOCATE (position(1:n,1:n,1:n,1:n))
pop(:,:) = .FALSE.
position(:,:,:,:) = .FALSE.
! some assignation to pop matrix
! position should trace the assignation, but it is not important for this example
pop(2,5) = .TRUE.
pop(75,1) = .TRUE. !...ans so on
DO a=1,n
DO b=1,n
IF (a.gt.100.and.b.lt.50) THEN
CYCLE
ELSE
CALL random_number(c)
IF (x.gt.0.5d0) THEN
new_pop(a,b) = .TRUE.
ELSE
new_pop(a,b) = pop(a,b)
END IF
END IF
END DO
END DO
! I would write that
DO a=1,n
WRITE(10,1000) pop(a,:)
END DO

It is old fashion, but it should work.
Inside the two DO cycle "pop" is never on the left side of the assignment.
Anyway, when I try to write, "pop(128,1)" doesn't report the correct value. In fact, it should have the initial value .FALSE., but in debug I found .TRUE.

I fear it can be linked to the memory.
In fact, I use a 32bit compiler on Windows 7 (Intel Fortan 10.x) and I fear Fortran can make some mess with the memory management (these are the only variables and the arraies are very large).
I should set better my compiler?

Regards,
Matteo
 
Labels (2)
0 Kudos
42 Replies
Highlighted
Black Belt Retired Employee
693 Views

It's not possible to diagnose a problem from a fragment of code. Sometimes this symptom can be caused by out-of-bounds array references or mismatched arguments. Do you have all the run-time checks enabled and build with /warn:interface?

Intel Fortran 10.x is VERY old. You probably need to also set /gen-interface with /warn:interface in that version, and even so, the checking there is not as good as in current versions.

It is most likely an error in your code and not the compiler.

0 Kudos
Highlighted
New Contributor II
668 Views

It really helps if you use the code formating item in the message -- tap on the ... and you will find it. 

Fortran is hard enough to read at the best of times, but when placed like FORTRAN it is nie on impossible. 

Highlighted
Novice
644 Views

Sorry: I edited  by my mobile and something went wrong.

INTEGER :: n,a,b
REAL*8 :: c
LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: pop, new_pop
LOGICAL, ALLOCATABLE, DIMENSION(:,:,:;:) :: position
!
OPEN(10, FILE = 'printout.txt')
n = 150
ALLOCATE(pop(1:n,1:n))
ALLOCATE(new_pop(1:n,1:n))
ALLOCATE (position(1:n,1:n,1:n,1:n))
pop(:,:) = .FALSE.
position(:,:,:,:) = .FALSE.
! some assignation to pop matrix
! position should trace the assignation, but it is not important for this example
pop(2,5) = .TRUE.
pop(75,1) = .TRUE. !...ans so on
DO a=1,n
 DO b=1,n
  IF (a.gt.100.and.b.lt.50) THEN
   CYCLE
  ELSE
   CALL random_number(c)
   IF (x.gt.0.5d0) THEN
    new_pop(a,b) = .TRUE.
   ELSE
    new_pop(a,b) = pop(a,b)
   END IF
  END IF
 END DO
END DO
! I would write that
DO a=1,n
 WRITE(10,1000) pop(a,:)
END DO
0 Kudos
Highlighted
Honored Contributor I
640 Views

What is 'x' in your code?

Did you mean 'c .gt. 0.5 ..' instead?
0 Kudos
Highlighted
Novice
634 Views

Hi @Steve_Lionel ,

thank you for your answer.
I found a more recent IFORT. That's the BuildLog:

 	 	
Compiling with Intel(R) Visual Fortran Compiler XE 14.0.6.241 [IA-32]...
ifort /nologo /debug:full /Od /debug-parameters:all /warn:declarations /warn:unused /warn:interfaces /fpe:0 /module:"Debug\\" /object:"Debug\\" /Fd"Debug\vc100.pdb" /traceback /check:all /libs:dll /threads /dbglibs /c /gen-interface /check:noarg_temp_created /Qvc10 /Qlocation,link,"c:\Program Files (x86)\Microsoft Visual Studio 10.0\Intel Fortran\Microsoft Files\VC\\bin" "E:\User\Fortran\SimPop\Simpop\Simpop\main_pop_rev2.f90"
E:\User\Fortran\SimPop\Simpop\Simpop\main_pop_rev2.f90(216): remark #8291: Recommended relationship between field width 'W' and the number of fractional digits 'D' in this edit descriptor is 'W>=D+7'.
 9001    FORMAT(I4, 1X, I6, 1X, E13.7)
---------------------------------^
E:\User\Fortran\SimPop\Simpop\Simpop\main_pop_rev2.f90(215): remark #8290: Recommended relationship between field width 'W' and the number of fractional digits 'D' in this edit descriptor is 'W>=D+3'.
 9000    FORMAT('XXXX = ',I4,' XXXXX = ',I6, ' XX = ' F5.3)
---------------------------------------------------------^

Linking...
Link /OUT:"Debug\Simpop_2.exe" /INCREMENTAL:NO /NOLOGO /MANIFEST /MANIFESTFILE:"E:\User\Fortran\SimPop\Simpop_2\Simpop_2\Debug\Simpop_2.exe.intermediate.manifest" /MANIFESTUAC:"level='asInvoker' uiAccess='false'" /DEBUG /PDB:"E:\User\Fortran\SimPop\Simpop_2\Simpop_2\Debug\Simpop_2.pdb" /SUBSYSTEM:CONSOLE /IMPLIB:"E:\User\Fortran\SimPop\Simpop_2\Simpop_2\Debug\Simpop_2.lib" "Debug\main_pop_rev2.obj"
Link: executing 'link'

Embedding manifest...
mt.exe /nologo /outputresource:"E:\User\Fortran\SimPop\Simpop_2\Simpop_2\Debug\Simpop_2.exe;#1" /manifest "E:\User\Fortran\SimPop\Simpop_2\Simpop_2\Debug\Simpop_2.exe.intermediate.manifest"

Simpop_2 - 0 error(s), 0 warning(s)

I have some warnings, but it doesn't seem serious problem.
As you can see, my code it is built in 32 bit (on a 64 bit machine). 
From Windows Task Manager, I found the executable file requires about 1.5 Gb of memory.
In fact, if I launch that in Release modality, it reports an error about "insufficient virtual memory".Do you think this size can be compatible with my architecture?
Anyway the same problem occurs reducing the arrays dimension.

Furthermore, I found the same problem in my Watch windows for other array, before calling any subroutine.
Both arrays are allocatable and, since they are in the main program, I never deallocate them.
The STAT during the ALLOCATE reports 0 for each allocation.
I don't know if it can help.
I don't know if I have some problem with my installation.

Thank you very much for your support

0 Kudos
Highlighted
Novice
636 Views

Hi @FortranFan .
yes, it simulates a coin toss.
Sorry for my mistake, the real code is correct 

0 Kudos
Highlighted
Black Belt Retired Employee
637 Views

This is now a very different problem than you first described. I agree the warnings you show can be ignored. 

1.5GB is approaching the limit on a 32-bit system. I assume that it was an ALLOCATE that got the insufficient virtual memory error. How much did you try to allocate?

Do you have some large COMMON blocks? It's hard to get to 1.5GB memory use from an ordinary program.

0 Kudos
Highlighted
Black Belt
632 Views

The array position(150,150,150,150) alone would take up about 2 GB or allocated memory, even if that array is not used in the rest of the code. Try removing that array from the code, if that is permissible.

0 Kudos
Highlighted
Novice
607 Views

Hi @Steve_Lionel  and @mecej4 ,

you are right.
I have just re-ordered my code. Now position has the same role and size (720000,6) with integer value (17.3 Mb ?)
Apart a better memory management, the issue is still present.

I putted 2 conditional breakpoints, for example:

  1. when a = 31 and b = 4
    (so, I monitor when the program draws this couple of value. Anyway pop is never on left side of the assignment, but to be sure...)
  2. when pop(31,4) = .TRUE. (the original value at the beginning of the program is .FALSE.)

The first breakpoint is not triggered during the whole program, but the second one just after the first cycle.
I should debug again, but at the moment I don't have any path to follow.

0 Kudos
Highlighted
Black Belt
596 Views

As I suspected, one of the allocate statements is failing when you build and run for IA32. Here is modified code, with STAT= clauses added to the ALLOCATE statements.

program matteo
implicit none
INTEGER :: n,a,b,astat
REAL*8 :: c
LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: pop, new_pop
LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: position
!
n = 150
ALLOCATE(pop(1:n,1:n),stat=astat)
if(astat /= 0)stop 'Error allocating POP'
ALLOCATE(new_pop(1:n,1:n),stat=astat)
if(astat /= 0)stop 'Error allocating NEW_POP'
ALLOCATE (position(1:n,1:n,1:n,1:n),stat=astat)
if(astat /= 0)stop 'Error allocating POSITION'
pop(:,:) = .FALSE.
position(:,:,:,:) = .FALSE.
! some assignation to pop matrix
! position should trace the assignation, but it is not important for this example
pop(2,5) = .TRUE.
pop(75,1) = .TRUE. !...ans so on
DO a=1,n
 DO b=1,n
  IF (a.gt.100.and.b.lt.50) THEN
   CYCLE
  ELSE
   CALL random_number(c)
   IF (c.gt.0.5d0) THEN
    new_pop(a,b) = .TRUE.
   ELSE
    new_pop(a,b) = pop(a,b)
   END IF
  END IF
 END DO
END DO
print *,pop(31,4)
end program

When I compile for IA32 and run, the program prints "Error allocating POSITION" and stops.

When an allocation fails, or an array subscript is out of bounds, etc., from that point onward anything that you see in the debugger can be wrong; if you trust what the debugger shows you, you (and we) can spend a lot of time on wild goose chases. 

When you post again to this forum, please post actual code in a code box, make sure that the code is complete (all details and files needed to compile, link and run are provided) and avoid stating misleading conclusions.

 

0 Kudos
Highlighted
Novice
571 Views

Dear @mecej4 ,

thank you for your suggestion that I am going to follow.
Maybe you have lost my answer, where I reduced the memory usage of the position, despite the fact that the program is slower now.
Anyway, since it is a program for a no-profit educational purpose, I can share the whole code in this forum. That's the original code running on my workstation.

      PROGRAM Sim_pop
         IMPLICIT NONE
         INTEGER :: n,j,a,b, k, aa, bb
         REAL*8 :: ar, br
         LOGICAL, DIMENSION(:,:), ALLOCATABLE :: pop,  pop_new
         INTEGER(kind=2), DIMENSION(:,:), ALLOCATABLE :: R0, immuni, degenza
         LOGICAL, DIMENSION(:,:), ALLOCATABLE :: guarito
         INTEGER, DIMENSION(:,:), ALLOCATABLE :: visite
         INTEGER :: ivisit, ivisit_old, nvis
         INTEGER :: days, t0, startc
         INTEGER, DIMENSION(1:2) :: r
         REAL*8, DIMENSION(1:2) :: p
         REAL*8:: cont_rate, r0_print
         REAL*8 :: c, raggio, angolo, std
         INTEGER :: uscite, usc, quarantena, immuno_day
         INTEGER :: contagiati, infetti
         REAL*8 :: abitanti
         INTEGER, DIMENSION(1:2) :: z
         REAL*8, DIMENSION(1:2) :: rr
         INTEGER, ALLOCATABLE, DIMENSION(:) :: x,y
         INTEGER :: xmin, xmax, ymin, ymax, day_free
         LOGICAL :: separa
         INTEGER(kind=1), ALLOCATABLE, DIMENSION(:,:) :: reg0
         INTEGER(kind=1), ALLOCATABLE, DIMENSION(:) :: vec
         INTEGER(kind=1), DIMENSION(1:9) :: error
         REAL :: tic, toc
         !
         CALL CPU_TIME(tic)
         abitanti = 22500.d0
         n = idnint(dsqrt(abitanti))
         z = [4,4]
         day_free = 100
         days = 60
         immuno_day = 15
         cont_rate = 0.01d0
         usc = 15
         std = DBLE(n)*dsqrt(2.d0)/2.d0
         t0 = 4
         quarantena = 15
         CALL random_seed(PUT = [23,42])
!
         OPEN(10, FILE ='res.txt')
         OPEN(11, FILE ='r0.txt')
         OPEN(12, FILE ='summary.txt')
         OPEN(13, FILE ='track.txt')
         OPEN(45, FILE ='debug46.txt')
!
         ALLOCATE(pop(1:n, 1:n), STAT = error(1))
         ALLOCATE(R0(1:n,1:n), STAT = error(2))
         ALLOCATE(pop_new(1:n, 1:n), STAT = error(3))
         ALLOCATE(degenza(1:n, 1:n), STAT = error(4))
         ALLOCATE(visite(1:(n*n*(usc+1)*2),1:7), STAT = error(5))
         ALLOCATE(guarito(1:n, 1:n), STAT = error(6))
         ALLOCATE(immuni(1:n, 1:n), STAT = error(7)) 
         ALLOCATE(reg0(1:t0,1:2), STAT = error(8))     
         ALLOCATE(vec(1:n), STAT = error(9))
         IF (sum(error).eq.0) THEN
             WRITE(*,*) 'Not Error in allocation'
         ELSE
             WRITE(*,*) 'Error'
             DO j=1,9
                 WRITE(*,*) error(J)
             END DO
             STOP
         END IF
         nvis = SIZE(visite, DIM=1)
         startc = 0
         pop(:,:) = .FALSE.
         R0(:,:) = 0
         pop_new(:,:) = 0
         degenza(:,:) = -1
         immuni(:,:) = 0
         guarito(:,:) = .FALSE.
         separa = .FALSE.
!
         DO j=1,2
             IF (Z(j).gt.1) THEN
                separa =.TRUE.
                IF (j.eq.1) THEN
                    ALLOCATE(x(1:z(j)))
                    CALL linspace(DBLE(1),DBLE(n),z(j),x)
                ELSE
                    ALLOCATE(y(1:z(j)))
                    CALL linspace(DBLE(1),DBLE(n),z(j),y)
                END IF
             END IF
         END DO
!
         DO WHILE (startc.lt.t0)
            CALL random_number(p)
            ar = (1.d0 + p(1) * (DBLE(n)))
            br = (1.d0 + p(2) * (DBLE(n)))
            a = floor(ar)
            b = floor(br)
            IF (pop(a,b).ne. .TRUE.) THEN
               pop(a,b) = .TRUE.
               degenza(a,b) = -1
               startc = startc + 1
               reg0(startc,:) = [a,b]
            END IF
         END DO
            DO k=1,n
                vec = 0
                WHERE (pop(k,:).eq. .TRUE.)
                    vec = 1
                END WHERE
                WRITE(10,1000) vec
                WRITE(11,1001) R0(k,:)
            END DO
         DO aa=1,n
            pop_new(aa,:) = pop(aa,:)
         END DO
         DO j=1,days
            infetti = 0
            DO a=1,n
                DO b=1,n
                    IF (pop(a,b)) THEN
                        infetti = infetti + 1
                    END IF
                END DO
            END DO
            contagiati = SUM(SUM(R0, DIM =1), DIM = 1)
            IF (infetti.ne.0) THEN
                r0_print = DBLE(contagiati)/DBLE(infetti)
            ELSE
                r0_print = 0.d0
            END IF
            WRITE(*,9000) j-1, infetti, r0_print
            WRITE(12,9001) j-1, infetti, r0_print
            R0(:,:) = 0
            visite(:,:) = 0
            ivisit = 0
            DO a=1,n
               DO b=1,n
                  uscite = 0
                  IF ((separa).and.(j.le.day_free)) THEN
                    CALL confina(a,z(1),x,xmin,xmax)   
                    CALL confina(b,z(2),y,ymin,ymax)  
                  ELSE
                    xmin = 1
                    xmax = n
                    ymin = 1
                    ymax = n
                  END IF
                  DO WHILE (uscite.le.usc)
                     CALL random_norm(raggio)
                     raggio = dabs(raggio) * std * 0.38532d0
                     CALL random_number(c)
                     angolo = (0.d0 + c * (360.d0))
                     rr(1) = (DBLE(a) + raggio*dcosd(angolo))
                     rr(2) = (DBLE(b) + raggio*dsind(angolo))
                     r(1)  = idnint(rr(1))
                     r(2)  = idnint(rr(2))
                     IF (a.eq.r(1).and.b.eq.r(2)) THEN
                         CYCLE
                     END IF
                     ! Nessuno può uscire dalla sua "cella" dentro pop
                     IF (r(1).ge.xmin.and.r(1).le.xmax) THEN
                        IF (r(2).ge.ymin.and.r(2).le.ymax) THEN
                           ivisit_old = ivisit
                           CALL cercavisite (j,[a,b,r(1),r(2)],[pop(a,b),pop(r(1),r(2))], ivisit, nvis, visite)
                           IF (ivisit.gt.ivisit_old) THEN
                               uscite = uscite + 1
                               IF (((pop(a,b).eq. .FALSE.).and.(pop(r(1),r(2)).eq.  .FALSE.)).or. &
                                  ((pop(a,b).eq. .TRUE.).and.(pop(r(1),r(2)).eq. .TRUE.))) THEN
                                CYCLE
                               ELSE
                                  CALL random_number(c)
                                  IF (c.le.cont_rate) THEN
                                  ! Calcolo R0 e propago virus
                                     IF ((pop(a,b).eq. .FALSE.).and.(immuni(a,b).eq.0)) THEN
                                            R0(r(1),r(2)) = R0(r(1),r(2)) + 1
                                            pop_new(a,b) = .TRUE.
                                     ELSEIF ((pop(r(1), r(2)).eq. .FALSE.).and.(immuni(r(1), r(2)).eq.0)) THEN
                                            R0(a,b) = R0(a,b) + 1
                                            pop_new(r(1),r(2)) = .TRUE.
                                     END IF
                                  END IF
                               END IF
                           END IF
                        END IF
                     END IF
                  END DO
               END DO
            END DO
            DO aa=1,n
                pop(aa,:) = pop_new(aa,:)
            END DO
            ! verifico le guarigioni
            DO aa=1,n
               DO bb=1,n
                  IF (pop(aa,bb)) THEN
                     degenza(aa,bb) = degenza(aa,bb) + 1
                  END IF
                  IF (guarito(aa,bb)) THEN
                    immuni(aa,bb) = immuni(aa,bb) -1
                  END IF
                  IF (immuni(aa,bb).le.0) THEN
                    guarito(aa,bb) = .FALSE.
                    immuni(aa,bb) = 0
                  END IF
                  IF ((pop(aa,bb)).and.(degenza(aa,bb).ge.quarantena)) THEN
                     degenza(aa,bb) = -1
                     pop(aa,bb) = .FALSE.
                     guarito(aa,bb) = .TRUE.
                     immuni(aa,bb) = immuno_day
                  END IF
               END DO
            END DO
            DO aa=1,n
                pop_new(aa,:) = pop(aa,:)
            END DO
            DO k=1,n
                vec = 0
                WHERE (pop(k,:).eq. .TRUE.)
                    vec = 1
                END WHERE
                WRITE(10,1000) vec
                IF (j.eq.46) THEN
                    WRITE(45,1000) vec
                END IF
                WRITE(11,1001) R0(k,:)
            END DO
            DO k=1,nvis
                WRITE(13,1001) visite(k,:)
            END DO
         END DO
         WRITE(*,*) 'End of simulation. Diseases = ',infetti
         CLOSE(10)
         CLOSE(11)
         CLOSE(45)
         CALL CPU_TIME(toc)
         print '("Time = ",f6.3," minutes.")',(tic-toc)/60.d0
!
 1000    FORMAT(1000(1X, I1))
 1001   FORMAT(1000(1X, I6))         
 9000    FORMAT('Day = ',I4,' Diseases = ',I6, ' R0 = ' F5.3)
 9001    FORMAT(I4, 1X, I6, 1X, E13.7)
         READ(*,*)
         STOP
      END PROGRAM
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Random_norm
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      SUBROUTINE random_norm(c)
         IMPLICIT NONE
         ReAL*8, INTENT(OUT) :: c
         REAL*8, DIMENSION(1:2) :: u
         REAL*8 :: pi = 3.14d0
!
         CALL random_number(u)
         C = dsqrt(-2.d0*dlog(u(1)))*dcos(2.d0*pi*u(2))
!
         RETURN
    END SUBROUTINE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Linpsace
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
      SUBROUTINE linspace(a,b,n,yi)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: n
      REAL*8, INTENT(IN) :: a,b
      REAL*8, DIMENSION(1:n) :: y
      INTEGER, INTENT(OUT), DIMENSION(1:n) :: yi
      INTEGER :: j
      REAL*8 :: p
!
      p = (b-a)/DBLE(n-1)
      DO j=1,n
        y(j) = a + p*DBLE(j-1)
        yi(j) = idnint(y(j))
      END DO
!
      RETURN
      END SUBROUTINE  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Confina
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       
      SUBROUTINE confina(a,n,x,xmin,xmax)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: a, n
      INTEGER, DIMENSION(1:n), INTENT(IN) :: x
      INTEGER, INTENT(OUT) :: xmin, xmax
      INTEGER :: j
!
      DO j=1,(n-1)
        IF(a.ge.x(j).and.a.le.x(j+1)) THEN
            xmin = x(j)
            xmax = x(j+1)
            EXIT
        END IF
      END DO
!
    RETURN
    END SUBROUTINE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! cercavisite
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
SUBROUTINE cercavisite (day,v,inf, ind_visite, n, visite)
IMPLICIT NONE
INTEGER, INTENT(IN) :: day
INTEGER, DIMENSION(1:4), INTENT(IN) :: v
LOGICAL, DIMENSION(1:2), INTENT(IN) :: inf
INTEGER, INTENT(INOUT) :: ind_visite
INTEGER, INTENT(INOUT), DIMENSION(1:n, 1:7) :: visite
INTEGER, INTENT(IN) :: n
!
LOGICAL :: succ
INTEGER :: j, k, somma
INTEGER, DIMENSION(1:2) :: i
!
succ = .TRUE.
DO j=1,ind_visite
    somma = 0
    DO k=1,4
        IF (v(k).eq.visite(j,k)) THEN
            somma = somma + 1
        END IF
    END DO
    IF (somma.eq.4) THEN
        succ = .FALSE.
        EXIT
    END IF
END DO
IF (succ) THEN
    DO j=1,2
        IF (inf(j)) THEN
            i(j) = 1
        ELSE 
            i(j) = 0
        END IF
    END DO
    ind_visite = ind_visite + 1
    visite(ind_visite,:) = [v,i,day]
    ind_visite = ind_visite + 1
    visite(ind_visite,:) = [v(3),v(4),v(1),v(2),i(2),i(1), day]
END IF
RETURN
END SUBROUTINE

  For your convenience I attached a zip file with following data:

  1. BuildLog, so you can understand how I set my compiler
  2. The original F90 file. It is the file which is required to compile the code
  3. Res.txt is the most important output file
  4. JPEG post processing of res.txt
  5. Summary.txt that is printed at your screen, too and they show something is wrong in my code.

Pay attention: the code is slower now, but it need only few Mb in RAM.
Anyway, a track.txt is very large (about 2 Gb). It could be interesting only in a debug phase.
If you comment line 224-226, this file will be not written.

Thank you for your support.

 

0 Kudos
Highlighted
Black Belt
564 Views

Matteo,

You have now posted a larger and substantially different code. Apart from the size, the question now is: "In what way is the compiler is treating this larger code that, to you, is wrong? Is there still a variable that is "changing in value without apparent reason" -- as you said earlier? If so, where do you see it?

You have to provide the details, since most of us neither know nor may be particularly interested in the application domain.

In the larger code, you have several places where you compare logical variables using the arithmetic comparison operator .EQ., for example, pop(a,b).eq. .FALSE. . See the cautionary notes in Equivalence-Versus-Equality  For related reasons, it is incorrect to set the logical variable pop_new to the integer value 0, as is done at Line 70.

I must admit that the JPG files that you provided and the numbers in the output files mean nothing to me.

0 Kudos
Highlighted
Novice
557 Views

@mecej4 ,

the previous code is just an example, while this is the real original code.
I understood my misuse of equivalence .vs. equality and I am going to re-arrange this part.
You are right, line 90 should be corrected with pop_new(:,:) = .FALSE.

Informatics aspects

pop is initialized at line 96 via random process.
A do-loop from line 133 to line 185 never change these values.
Anyway, a breakpoint at line 186 will show some variation in pop.

Context

If you are interesting in the frame, I can give additional information.
Pop array is "divided" in 9 square sectors, where x and y give their boundaries.
Any pop(a,b) can move only inside its sector, computed any time by confina subroutine.
Pop(a,b) can "visit" pop(r(1),r(2)) and can modify that and it can be modified, where pop_new reports the new status (not pop!). The xmin, xmax, ymin and ymax computed by confina, constraint pop(a,b) to not exit from its own sector.

The jpeg reports in red the .TRUE. (1) of pop matrix at any time j and in blue the FALSE (0).
Since pop has only 4 true values, at the end of the program I should see 5 all blue sector in jpeg.

I hope it is clearer what I have in mind.
Regards.

0 Kudos
Highlighted
Valued Contributor III
551 Views

Mecej4 is correct there are many non standard things in your code. It is a small code and correcting should be quick and easy. It may not fix any problem you have but it is a good basis to start with a conforming code. Compile with /stand.

I will also note your loops are all very inefficient at you loop  array(i,:) everywhere so you jump around memory as you are never working on adjacent memory locations, array(:,i).

Things like 

 

 ! a bit insane
           DO aa=1,n
                pop(aa,:) = pop_new(aa,:)
            END DO
! why not 
            DO aa=1,n
                pop(:,aa) = pop_new(:,aa)
            END DO
! or better
             pop = pop_new 

 

Highlighted
Black Belt Retired Employee
532 Views

@Matteo , you misunderstood mecej4's comment.

For example, you have:

IF (pop(a,b).ne. .TRUE.) THEN

This is simply wrong. Do not use .eq. or .ne. on LOGICAL values. You could write:

IF (pop(a,b).neqv. .TRUE.) THEN

but better would be simply:

IF (.not. pop(a,b)) THEN

I discuss this issue in Doctor Fortran in “To .EQV. or to .NEQV., that is the question”, or “It’s only LOGICAL”

 

Highlighted
Novice
517 Views

Thank you all.
I am always glad to learn something.

@Steve_Lionel 
I understood the two remarks, the bad pop_new initialization and the wrong use of IF argument for logical variable.
Now it should be correct.

@andrew_4619 
you are right too. vec(:,j) is faster than vec(j,:).
As I said, this is only a prototype code for no-profit activity.
It will be slow at the beginning, but when I will be sure it works, I am going to optimize that (or someone else will do it).

In any case that's the new code that it is running on my workstation.
In attachment you can find the F90 file if you would like to run it.
When I will get the results, I will share with you.

Regards.

 

      PROGRAM Sim_pop
         IMPLICIT NONE
         INTEGER :: n,j,a,b, k, aa, bb
         REAL*8 :: ar, br
         LOGICAL, DIMENSION(:,:), ALLOCATABLE :: pop,  pop_new
         INTEGER(kind=2), DIMENSION(:,:), ALLOCATABLE :: R0, immuni, degenza
         LOGICAL, DIMENSION(:,:), ALLOCATABLE :: guarito
         INTEGER, DIMENSION(:,:), ALLOCATABLE :: visite
         INTEGER :: ivisit, ivisit_old, nvis
         INTEGER :: days, t0, startc
         INTEGER, DIMENSION(1:2) :: r
         REAL*8, DIMENSION(1:2) :: p
         REAL*8:: cont_rate, r0_print
         REAL*8 :: c, raggio, angolo, std
         INTEGER :: uscite, usc, quarantena, immuno_day
         INTEGER :: contagiati, infetti
         REAL*8 :: abitanti
         INTEGER, DIMENSION(1:2) :: z
         REAL*8, DIMENSION(1:2) :: rr
         INTEGER, ALLOCATABLE, DIMENSION(:) :: x,y
         INTEGER :: xmin, xmax, ymin, ymax, day_free
         LOGICAL :: separa
         INTEGER(kind=1), ALLOCATABLE, DIMENSION(:,:) :: reg0
         INTEGER(kind=1), ALLOCATABLE, DIMENSION(:) :: vec
         INTEGER(kind=1), DIMENSION(1:9) :: error
         REAL :: tic, toc
         !
         CALL CPU_TIME(tic)
         abitanti = 22500.d0
         n = idnint(dsqrt(abitanti))
         z = [4,4]
         day_free = 100
         days = 60
         immuno_day = 15
         cont_rate = 0.01d0
         usc = 15
         std = DBLE(n)*dsqrt(2.d0)/2.d0
         t0 = 4
         quarantena = 15
         CALL random_seed(PUT = [23,42])
!
         OPEN(10, FILE ='res.txt')
         OPEN(11, FILE ='r0.txt')
         OPEN(12, FILE ='summary.txt')
         OPEN(13, FILE ='track.txt')
         OPEN(45, FILE ='debug46.txt')
!
         ALLOCATE(pop(1:n, 1:n), STAT = error(1))
         ALLOCATE(R0(1:n,1:n), STAT = error(2))
         ALLOCATE(pop_new(1:n, 1:n), STAT = error(3))
         ALLOCATE(degenza(1:n, 1:n), STAT = error(4))
         ALLOCATE(visite(1:(n*n*(usc+1)*2),1:7), STAT = error(5))
         ALLOCATE(guarito(1:n, 1:n), STAT = error(6))
         ALLOCATE(immuni(1:n, 1:n), STAT = error(7)) 
         ALLOCATE(reg0(1:t0,1:2), STAT = error(8))     
         ALLOCATE(vec(1:n), STAT = error(9))
         IF (sum(error).eq.0) THEN
             WRITE(*,*) 'Not Error in allocation'
         ELSE
             WRITE(*,*) 'Error'
             DO j=1,9
                 WRITE(*,*) error(J)
             END DO
             STOP
         END IF
         nvis = SIZE(visite, DIM=1)
         startc = 0
         pop(:,:) = .FALSE.
         R0(:,:) = 0
         pop_new(:,:) = .TRUE.
         degenza(:,:) = -1
         immuni(:,:) = 0
         guarito(:,:) = .FALSE.
         separa = .FALSE.
!
         DO j=1,2
             IF (Z(j).gt.1) THEN
                separa =.TRUE.
                IF (j.eq.1) THEN
                    ALLOCATE(x(1:z(j)))
                    CALL linspace(DBLE(1),DBLE(n),z(j),x)
                ELSE
                    ALLOCATE(y(1:z(j)))
                    CALL linspace(DBLE(1),DBLE(n),z(j),y)
                END IF
             END IF
         END DO
!
         DO WHILE (startc.lt.t0)
            CALL random_number(p)
            ar = (1.d0 + p(1) * (DBLE(n)))
            br = (1.d0 + p(2) * (DBLE(n)))
            a = floor(ar)
            b = floor(br)
            IF (.not. pop(a,b)) THEN
               pop(a,b) = .TRUE.
               degenza(a,b) = -1
               startc = startc + 1
               reg0(startc,:) = [a,b]
            END IF
         END DO
            DO k=1,n
                vec = 0
                WHERE (pop(k,:))
                    vec = 1
                END WHERE
                WRITE(10,1000) vec
                WRITE(11,1001) R0(k,:)
            END DO
         pop_new = pop
         DO j=1,days
            infetti = 0
            DO a=1,n
                DO b=1,n
                    IF (pop(a,b)) THEN
                        infetti = infetti + 1
                    END IF
                END DO
            END DO
            contagiati = SUM(SUM(R0, DIM =1), DIM = 1)
            IF (infetti.ne.0) THEN
                r0_print = DBLE(contagiati)/DBLE(infetti)
            ELSE
                r0_print = 0.d0
            END IF
            WRITE(*,9000) j-1, infetti, r0_print
            WRITE(12,9001) j-1, infetti, r0_print
            R0(:,:) = 0
            visite(:,:) = 0
            ivisit = 0
            DO a=1,n
               DO b=1,n
                  uscite = 0
                  IF ((separa).and.(j.le.day_free)) THEN
                    CALL confina(a,z(1),x,xmin,xmax)   
                    CALL confina(b,z(2),y,ymin,ymax)  
                  ELSE
                    xmin = 1
                    xmax = n
                    ymin = 1
                    ymax = n
                  END IF
                  DO WHILE (uscite.le.usc)
                     CALL random_norm(raggio)
                     raggio = dabs(raggio) * std * 0.38532d0
                     CALL random_number(c)
                     angolo = (0.d0 + c * (360.d0))
                     rr(1) = (DBLE(a) + raggio*dcosd(angolo))
                     rr(2) = (DBLE(b) + raggio*dsind(angolo))
                     r(1)  = idnint(rr(1))
                     r(2)  = idnint(rr(2))
                     IF (a.eq.r(1).and.b.eq.r(2)) THEN
                         CYCLE
                     END IF
                     ! Nessuno può uscire dalla sua "cella" dentro pop
                     IF (r(1).ge.xmin.and.r(1).le.xmax) THEN
                        IF (r(2).ge.ymin.and.r(2).le.ymax) THEN
                           ivisit_old = ivisit
                           CALL cercavisite (j,[a,b,r(1),r(2)],[pop(a,b),pop(r(1),r(2))], ivisit, nvis, visite)
                           IF (ivisit.gt.ivisit_old) THEN
                               uscite = uscite + 1
                               IF (((.not. pop(a,b)).and.(.not. pop(r(1),r(2)))).or. &
                                  ((pop(a,b)).and.(pop(r(1),r(2))))) THEN
                                CYCLE
                               ELSE
                                  CALL random_number(c)
                                  IF (c.le.cont_rate) THEN
                                  ! Calcolo R0 e propago virus
                                     IF ((.not. pop(a,b)).and.(immuni(a,b).eq.0)) THEN
                                            R0(r(1),r(2)) = R0(r(1),r(2)) + 1
                                            pop_new(a,b) = .TRUE.
                                     ELSEIF ((.not. pop(r(1), r(2))).and.(immuni(r(1), r(2)).eq.0)) THEN
                                            R0(a,b) = R0(a,b) + 1
                                            pop_new(r(1),r(2)) = .TRUE.
                                     END IF
                                  END IF
                               END IF
                           END IF
                        END IF
                     END IF
                  END DO
               END DO
            END DO
            pop = pop_new
            ! verifico le guarigioni
            DO aa=1,n
               DO bb=1,n
                  IF (pop(aa,bb)) THEN
                     degenza(aa,bb) = degenza(aa,bb) + 1
                  END IF
                  IF (guarito(aa,bb)) THEN
                    immuni(aa,bb) = immuni(aa,bb) -1
                  END IF
                  IF (immuni(aa,bb).le.0) THEN
                    guarito(aa,bb) = .FALSE.
                    immuni(aa,bb) = 0
                  END IF
                  IF ((pop(aa,bb)).and.(degenza(aa,bb).ge.quarantena)) THEN
                     degenza(aa,bb) = -1
                     pop(aa,bb) = .FALSE.
                     guarito(aa,bb) = .TRUE.
                     immuni(aa,bb) = immuno_day
                  END IF
               END DO
            END DO
            pop_new = pop
            DO k=1,n
                vec = 0
                WHERE (pop(k,:))
                    vec = 1
                END WHERE
                WRITE(10,1000) vec
                IF (j.eq.46) THEN
                    WRITE(45,1000) vec
                END IF
                WRITE(11,1001) R0(k,:)
            END DO
            DO k=1,nvis
                WRITE(13,1001) visite(k,:)
            END DO
         END DO
         WRITE(*,*) 'End of simulation. Diseases = ',infetti
         CLOSE(10)
         CLOSE(11)
         CLOSE(45)
         CALL CPU_TIME(toc)
         print '("Time = ",f12.3," minutes.")',(tic-toc)/60.d0
!
 1000    FORMAT(1000(1X, I1))
 1001   FORMAT(1000(1X, I6))         
 9000    FORMAT('Day = ',I4,' Diseases = ',I6, ' R0 = ' F5.3)
 9001    FORMAT(I4, 1X, I6, 1X, E13.7)
         READ(*,*)
         STOP
      END PROGRAM
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Random_norm
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      SUBROUTINE random_norm(c)
         IMPLICIT NONE
         ReAL*8, INTENT(OUT) :: c
         REAL*8, DIMENSION(1:2) :: u
         REAL*8 :: pi = 3.14d0
!
         CALL random_number(u)
         C = dsqrt(-2.d0*dlog(u(1)))*dcos(2.d0*pi*u(2))
!
         RETURN
    END SUBROUTINE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Linpsace
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
      SUBROUTINE linspace(a,b,n,yi)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: n
      REAL*8, INTENT(IN) :: a,b
      REAL*8, DIMENSION(1:n) :: y
      INTEGER, INTENT(OUT), DIMENSION(1:n) :: yi
      INTEGER :: j
      REAL*8 :: p
!
      p = (b-a)/DBLE(n-1)
      DO j=1,n
        y(j) = a + p*DBLE(j-1)
        yi(j) = idnint(y(j))
      END DO
!
      RETURN
      END SUBROUTINE  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Confina
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       
      SUBROUTINE confina(a,n,x,xmin,xmax)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: a, n
      INTEGER, DIMENSION(1:n), INTENT(IN) :: x
      INTEGER, INTENT(OUT) :: xmin, xmax
      INTEGER :: j
!
      DO j=1,(n-1)
        IF(a.ge.x(j).and.a.le.x(j+1)) THEN
            xmin = x(j)
            xmax = x(j+1)
            EXIT
        END IF
      END DO
!
    RETURN
    END SUBROUTINE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! cercavisite
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
SUBROUTINE cercavisite (day,v,inf, ind_visite, n, visite)
IMPLICIT NONE
INTEGER, INTENT(IN) :: day
INTEGER, DIMENSION(1:4), INTENT(IN) :: v
LOGICAL, DIMENSION(1:2), INTENT(IN) :: inf
INTEGER, INTENT(INOUT) :: ind_visite
INTEGER, INTENT(INOUT), DIMENSION(1:n, 1:7) :: visite
INTEGER, INTENT(IN) :: n
!
LOGICAL :: succ
INTEGER :: j, k, somma
INTEGER, DIMENSION(1:2) :: i
!
succ = .TRUE.
DO j=1,ind_visite
    somma = 0
    DO k=1,4
        IF (v(k).eq.visite(j,k)) THEN
            somma = somma + 1
        END IF
    END DO
    IF (somma.eq.4) THEN
        succ = .FALSE.
        EXIT
    END IF
END DO
IF (succ) THEN
    DO j=1,2
        IF (inf(j)) THEN
            i(j) = 1
        ELSE 
            i(j) = 0
        END IF
    END DO
    ind_visite = ind_visite + 1
    visite(ind_visite,:) = [v,i,day]
    ind_visite = ind_visite + 1
    visite(ind_visite,:) = [v(3),v(4),v(1),v(2),i(2),i(1), day]
END IF
RETURN
END SUBROUTINE

 

 

0 Kudos
Highlighted
Black Belt
490 Views

Matteo,

In a previous post in this thread you wrote:

Pop array is "divided" in 9 square sectors, where x and y give their boundaries.
Any pop(a,b) can move only inside its sector, computed any time by confina subroutine. Pop(a,b) can "visit" pop(r(1),r(2)) and can modify that and it can be modified, where pop_new reports the new status (not pop!). The xmin, xmax, ymin and ymax computed by confina, constraint pop(a,b) to not exit from its own sector.

The jpeg reports in red the .TRUE. (1) of pop matrix at any time j and in blue the FALSE (0). Since pop has only 4 true values, at the end of the program I should see 5 all blue sector in jpeg.

Are you implying that the number of elements of POP that are true cannot change? If so, the code does not satisfy that requirement. Specifically, you have 

IF ((.not. pop(a,b)).and.(immuni(a,b).eq.0)) THEN
  R0(r(1),r(2)) = R0(r(1),r(2)) + 1
  pop_new(a,b) = .TRUE.
ELSEIF ((.not. pop(r(1), r(2))).and.(immuni(r(1), r(2)).eq.0)) THEN
  R0(a,b) = R0(a,b) + 1
  pop_new(r(1),r(2)) = .TRUE.
END IF

Either of the clauses in the IF-THEN-ELSEIF-ENDIF construct will cause the number of nodes where pop_new is .TRUE. to increase. Perhaps you need to reset the  old locations to .FALSE., as in:

IF ((.not. pop(a,b)).and.(immuni(a,b).eq.0)) THEN
  R0(r(1),r(2)) = R0(r(1),r(2)) + 1
  pop_new(a,b) = .TRUE.
  pop_new(r(1),r(2)) = .FALSE. 
ELSEIF ((.not. pop(r(1), r(2))).and.(immuni(r(1), r(2)).eq.0)) THEN
  R0(a,b) = R0(a,b) + 1
  pop_new(r(1),r(2)) = .TRUE.
  pop_new(a,b) = .FALSE.
END IF

 Apart from this issue, I have not seen a single piece of evidence that POP(:,:) is getting any values changed inside the DO a = .. and DO b = .. loops.

0 Kudos
Highlighted
Novice
476 Views

@mecej4 ,

the number of TRUE elements pop change in this way.
Any j iteration some pop_new can change in TRUE and, at the end:

pop = pop_new

The code works in this way:

  1. any pop(a,b) can keep in touch with a casual pop(r(1), r(2)) 
    CALL random_norm(raggio)
                         raggio = dabs(raggio) * std * 0.38532d0
                         CALL random_number(c)
                         angolo = (0.d0 + c * (360.d0))
                         rr(1) = (DBLE(a) + raggio*dcosd(angolo))
                         rr(2) = (DBLE(b) + raggio*dsind(angolo))
                         r(1)  = idnint(rr(1))
                         r(2)  = idnint(rr(2))
                         IF (a.eq.r(1).and.b.eq.r(2)) THEN
                             CYCLE
                         END IF
  2. The code check that pop(r(1), r(2))  is in the same pop(a,b) "sector" computed by confina
    IF (r(1).ge.xmin.and.r(1).le.xmax) THEN
                            IF (r(2).ge.ymin.and.r(2).le.ymax) THEN
    [...]​
  3. A random system, if pop(a,b) or pop(r(1), r(2))  are .TRUE., both pop(a,b) or pop(r(1), r(2)) can switch from .FALSE. to .TRUE. but not viceversa (in the part of code you quoted).

So, pop_new is reset to pop values only at the end of the j iteration.

pop = pop_new
! other operations on pop
pop_new = pop

 

So, you didn't find any evidence that pop change during a e b DO. In fact it should work in this way!
But, if you run this code in the debug with breakpoints positioned where I suggested, you will find some pop values flip from .FALSE. to .TRUE.
In fact, the pop matrix, written in res.txt, shows each sector (they are 9) are infected by the original 4 pop TRUE values. The confina subroutine (read the previous bullets) avoids the pop values can go out from their own sector.

By the way, this night the code ran with the modification that I have just reported.
The res.txt files is identical to the previous one.

0 Kudos
Highlighted
Black Belt
497 Views

During the last couple of days, I thought about why Matteo's code ran slowly, and I have found a solution.

The code performs a Monte-Carlo simulation of disease propagation. Points in a 150 X 150 grid are traversed in a random walk with some restricting rules. The simulation spans 60 days worth of disease spreading. To advance the positions of the "infected" entities, it is necessary to answer the question, "Have I been in this location before?" The method presently used is to keep adding to a large, unsorted array of one record for each location visited, somewhat similar to writing a new line in a registration book with name, phone number and address. For each grid point, up to 32 records may be entered.

Before adding a line, however, the writer has to check all the preceding entries (as many as 150 X 150 X 32) are scanned using a sequential search, and this makes the program run very slowly. As the days progress, the haystack grows taller, and we have to find thousands of needles in it and many times over. I let the program run only for a 5-day simulation, and that took 30 minutes on a PC with an i5-8400 CPU. Linear extrapolation gives an estimate of 6 hours for the computer run a simulation for the full 60 days.

Here is the fix. Instead of keeping a single huge register of visitation records, keep a 150 X 150 array of one page registers, with each page capable of holding about 70 visit records. When searching, we only search in the page for that particular square of the grid. I have attached a modified program that I constructed with this idea. It runs in less than a minute of CPU time for the full simulation.

Permit me to drive home an important conclusion: Optimize the algorithm before you optimize the code.

0 Kudos