- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Did you mean 'c .gt. 0.5 ..' instead?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi @FortranFan .
yes, it simulates a coin toss.
Sorry for my mistake, the real code is correct
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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:
- 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...) - 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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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 programWhen 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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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:
- BuildLog, so you can understand how I set my compiler
- The original F90 file. It is the file which is required to compile the code
- Res.txt is the most important output file
- JPEG post processing of res.txt
- 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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
@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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
@Matteo , you misunderstood mecej4's comment.
For example, you have:
IF (pop(a,b).ne. .TRUE.) THENThis is simply wrong. Do not use .eq. or .ne. on LOGICAL values. You could write:
IF (pop(a,b).neqv. .TRUE.) THENbut better would be simply:
IF (.not. pop(a,b)) THENI discuss this issue in Doctor Fortran in “To .EQV. or to .NEQV., that is the question”, or “It’s only LOGICAL”
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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 IFApart 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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
@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_newThe code works in this way:
- 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 - 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 [...] - 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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I wanted to correct a couple of typos in my most recent post, but I am unable to do so.
The Forum software shows only the initial post if I am logged in. I can see my response that I want to edit if I am logged out, but then I cannot edit it. Catch 22, anyone?
Here are the corrections.
1. Change
":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."
to
"Before adding a line, however, the writer has to check all the preceding entries (as many as 150 X 150 X 32, which are scanned using a sequential search), and this makes the program run very slowly."
2. Change
"for the computer run a simulation"
to
"for the computer to run a simulation"
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page