gfortran ---problem with random umber
rudra
here is the full code that is creating the problem!!! while running
with ifort, the configuration is created as wished, but not in
gfortran....
!===========================================
! This is the main driver routine
! of the reverse monte carlo program
program rmonte
!===========================================
implicit none
integer::i,j
integer:: seed,tr,rmcs
integer::ntyp,nap,ntot,it1,it2,atnm
real(selected_real_kind(8))::rmin,rmax,boxl,dr
real(selected_real_kind(8))::rr,rdf
real(selected_real_kind(8)),dimension(100)::x,y,z
real(selected_real_kind(8)),dimension(100):: gr,rrl
real(selected_real_kind(8)):: dx,dy,dz,dist,rcut
character(2)::elmnt
!
real::t1,t2
t1=secnds(t1)
!
seed=12345
BOXL=20
RCUT=1.5
DR=5.0
RMCS=1000
ntot=100
tr=1
write(*,*) &
"===============Creating Initial Configuration=============="
write(*,&
'(1x,"Minimum distance between any two particle is",1x,f6.4)'),
&
rcut
!-----------------------------------------------
! INITIALISE POSITION
!-----------------------------------------------
10 do i=1,ntot
x(i)=boxl*(ran(seed)-0.0)
y(i)=boxl*(ran(seed)-0.0)
z(i)=boxl*(ran(seed)-0.0)
call pbc(x,y,z,boxl)
end do
!- - - - - - - - - - - - - - - - - - - - - - - -
! KEEPING MINIMUM DIST. BETWEEN ANY TWO
! ATOM GREATER THEN RCUT
!- - - - - - - - - - - - - - - - - - - - - - - -
do i=1,ntot-1
do j=i+1,ntot
dx=x(i)-x(j)
dy=y(i)-y(j)
dz=z(i)-z(j)
dist=dsqrt(dx*dx+dy*dy+dz*dz)
tr=tr+1
! write(*,*) tr
if (dist<rcut)go to 10
end do
end do
write(*,'(1x,"Success after",1x,i0,1x,"try")')tr
call mindist(x,y,z)
write(*,*) &
"==========================================================="
write(*,*) ""
!-----------------------------------------------
! INITIAL CONFIGURATION DONE
!-----------------------------------------------
t2=secnds(t1)
write(*,'("Elapsed Time =",1x,f10.4)'), t2
end
!===========================================
! Subroutine to apply periodic
! boundary value
!===========================================
subroutine pbc(xx,yy,zz,boxl)
implicit none
real(selected_real_kind(4))::xx,yy,zz,boxl
if(xx.gt. boxl)xx=xx-boxl
if(xx.lt.-boxl)xx=xx+boxl
if(yy.gt. boxl)yy=yy-boxl
if(yy.lt.-boxl)yy=yy+boxl
if(zz.gt. boxl)zz=zz-boxl
if(zz.lt.-boxl)zz=zz+boxl
return
end
!==============================================
! Subroutine to calculate minimum
! distance between two atom
subroutine mindist(x,y,z)
!==============================================
implicit none
integer::i,j,p1,p2
real(selected_real_kind(8)),dimension(100):: x,y,z
real(selected_real_kind(8)):: dist,dx,dy,dz,dmini,dmin
write(*, &
'(" Calculating minimum distance between atoms....",$)')
dx=x(1)-x(2)
dy=y(1)-y(2)
dz=z(1)-z(2)
dmini=dsqrt(dx*dx+dy*dy+dz*dz)
do i=1,99
do j=i+1,100
dx=x(i)-x(j)
dy=y(i)-y(j)
dz=z(i)-z(j)
dmin=dsqrt(dx*dx+dy*dy+dz*dz)
if (dmini>dmin)then
dmini=dmin
p1=i
p2=j
end if
end do
end do
write(*,*) "DONE"
write &
(*,'(1x,"Shortest distance between to atom is",2x,f8.6)')dmini
write(*,'(1x,"between particle pair",1x,i0," - ",i0)') p1,p2
return
end
while running with ifort, its prompting:
$ ./a.out
===============Creating Initial Configuration==============
Minimum distance between any two particle is 1.5000
Success after 1680724 try
Calculating minimum distance between atoms.... DONE
Shortest distance between to atom is 1.626997
between particle pair 40 - 92
===========================================================
Elapsed Time = 0.0623
but while running with gfortran, its never success...why?
Steven G. Kargl
rudra <bnrj....@gmail.com> writes:
> dear friends,
> here is the full code that is creating the problem!!! while running
> with ifort, the configuration is created as wished, but not in
> gfortran....
>
> seed=12345
> 10 do i=1,ntot
> x(i)=boxl*(ran(seed)-0.0)
> y(i)=boxl*(ran(seed)-0.0)
> z(i)=boxl*(ran(seed)-0.0)
> call pbc(x,y,z,boxl)
> end do
> but while running with gfortran, its never success...why?
Read the f*cking manual!!!!!!
I've already pointed you at the docs twice, and I gave you
test code.
RAN() does not do want you think.
Think about the following.
10 do i=1,ntot
x(i)=boxl*(3.14159-0.0)
y(i)=boxl*(3.14159-0.0)
z(i)=boxl*(3.14159-0.0)
call pbc(x,y,z,boxl)
end do
Ron Ford
6.160 RAN — Real pseudo-random number
Description:
For compatibility with HP FORTRAN 77/iX, the RAN intrinsic is provided as
an
alias for RAND. See Section 6.161 [RAND], page 121 for complete
documentation.
Standard: GNU extension
Class: Function
See also: Section 6.161 [RAND], page 121, Section 6.162 [RANDOM NUMBER],
That's what the manual says.
--
Steven G. Kargl
And, you obvious followed the link. You can't possibly be that stupid!
_Description_:
`RAND(FLAG)' returns a pseudo-random number from a uniform
distribution between 0 and 1. If FLAG is 0, the next number in the
current sequence is returned; if FLAG is 1, the generator is
restarted by `CALL SRAND(0)'; if FLAG has any other value, it is
used as a new seed with `SRAND'.
Ron Ford
"Steven G. Kargl" <ka...@troutmask.apl.washington.edu> wrote in message
news:g0j4va$8ad$1...@gnus01.u.washington.edu...
No, Steve, I did not follow the link, but I'm a dumb union guy.
>
> _Description_:
> `RAND(FLAG)' returns a pseudo-random number from a uniform
> distribution between 0 and 1. If FLAG is 0, the next number in the
> current sequence is returned; if FLAG is 1, the generator is
> restarted by `CALL SRAND(0)'; if FLAG has any other value, it is
> used as a new seed with `SRAND'.
This seems much more imformative than the entry in gfortran.pdf.
--
"We wouldn't have a police department that had to show a profit every
year - or a fire department that had to show a profit every year. We
shouldn't do that with our health care system, either."
~~ Michael Moore
steven....@gmail.com
> "Steven G. Kargl" <ka...@troutmask.apl.washington.edu> wrote in messagenews:g0j4va$8ad$1...@gnus01.u.washington.edu...
>
>
>
> > Ron Ford <r...@nowhere.net> writes:
> >> On Fri, 16 May 2008 04:52:17 +0000 (UTC), Steven G. Kargl wrote:
>
> >>> In article
> >>> rudra <bnrj.ru...@gmail.com> writes:
> >>>> dear friends,
> >>>> here is the full code that is creating the problem!!! while running
> >>>> with ifort, the configuration is created as wished, but not in
> >>>> gfortran....
>
> >>>> seed=12345
>
> >>>> 10 do i=1,ntot
> >>>> x(i)=boxl*(ran(seed)-0.0)
> >>>> y(i)=boxl*(ran(seed)-0.0)
> >>>> z(i)=boxl*(ran(seed)-0.0)
> >>>> call pbc(x,y,z,boxl)
> >>>> end do
>
> >>>> but while running with gfortran, its never success...why?
>
> >>> Read the f*cking manual!!!!!!
>
> >>> I've already pointed you at the docs twice, and I gave you
> >>> test code.
>
> >>> RAN() does not do want you think.
>
> >>> Think about the following.
>
> >>> 10 do i=1,ntot
> >>> x(i)=boxl*(3.14159-0.0)
> >>> y(i)=boxl*(3.14159-0.0)
> >>> z(i)=boxl*(3.14159-0.0)
> >>> call pbc(x,y,z,boxl)
> >>> end do
>
> >> Description:
> >> For compatibility with HP FORTRAN 77/iX, the RAN intrinsic is provided as
> >> an
> >> alias for RAND. See Section 6.161 [RAND], page 121 for complete
> >> documentation.
> >> Standard: GNU extension
> >> Class: Function
> >> See also: Section 6.161 [RAND], page 121, Section 6.162 [RANDOM NUMBER],
>
> >> That's what the manual says.
>
> > And, you obvious followed the link. You can't possibly be that stupid!
>
> No, Steve, I did not follow the link, but I'm a dumb union guy.
>
>
>
> > _Description_:
> > `RAND(FLAG)' returns a pseudo-random number from a uniform
> > distribution between 0 and 1. If FLAG is 0, the next number in the
> > current sequence is returned; if FLAG is 1, the generator is
> > restarted by `CALL SRAND(0)'; if FLAG has any other value, it is
> > used as a new seed with `SRAND'.
>
> This seems much more imformative than the entry in gfortran.pdf.
This *is* an entry in gfortran.pdf. Under RAND. Geez.
Gr.
Steven
FX
The entry Steve quoted *is* in your doc (whether HTML of PDF). The RAN
entry, which you quoted, says "see RAND for complete documentation". Now,
the RAND entry is just next to that one, and even if it weren't both HTML
and PDF docs have plenty of hyperlinks to get you from one place to
another: in the case here, clicking on *RAND* in "see RAND for complete
documentation" would have scrolled down to the RAND entry.
When I look at "colour" in my New Oxford American Dictionary, it says
"Bristish spelling of COLOR". And I go to COLOR (and they also have
hyperlinks) if I want more.
--
FX
e p chandler
How did you get this program to compile?
1. The source code looks like fixed format, but it is not. After
correcting one line break it does compile as free format.
2. G95 does not recognize ran(). As others have noted, RAN() is not
defined by the language standard. Its behavior is whatever the
"vendor" says it is. After replacing ran() by rand(), the program
still will not compile.
3. Using "magic" numbers like 4 or 8 as arguments to
selected_real_kind is even less portable than the old fashioned
real*4, real*8 notation.
4. Your actual arguments to pbc() are real(8) but the routine is
written to accept real(4).
5. pbc() is written to accept simple variables (scalars) as arguments
but your actual arguments are the names of arrays.
So even if it does compile and run you are probably passing it
x(1),y(1),z(1) not x(i),y(i),z(i).
6. Adding or subtracting BOX just ONCE does not guarantee that x,y,z
have absoulute value in range!
7. Most probably you are mis-applying ran(). If ran() is really an
alias for RAND(), then you really do not want to use it at all.
8. [Rhetorical] How is it possible for such a short program to have SO
many problems????
9. I think we should sue ALL compiler vendors for NOT writing this
warning in 72 point type on EACH and EVERY compiler box:
"WARNING. Just because your program compiles with our compiler does
not insure that it is either correct or standard conforming."
:-).
-- e
Richard Maine
> 3. Using "magic" numbers like 4 or 8 as arguments to
> selected_real_kind is even less portable than the old fashioned
> real*4, real*8 notation.
That one I disagree with. Nothing wrong with such numbers as arguments
to selected_real_kind. That's the number of decimal digits you want and
is pretty much the normal use of selected_real kind. I normally put the
selected_real_kind invocation in the definition of a named constant in a
module instead of using the whole long-winded thing in each type
definition, but otherwise literal constants are exactly what I tend to
use with it.
I suspect that the particular numbers used (4 and 8) are confusing you
into thinking that they are kind values instead of arguments to
selected_real_kind. Of course, I'd posit that, considering the number of
other errors here, and that 4 is a pretty small number of decimal digits
to ask for, the OP quite likely has this same confusion and actually
meant to use the 4 and 8 as kind values, which is indeed non-portable
and non-advisable.
Most of your other comments I agree with.
--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain
e p chandler
> e p chandler <e...@juno.com> wrote:
>
> > 3. Using "magic" numbers like 4 or 8 as arguments to
> > selected_real_kind is even less portable than the old fashioned
> > real*4, real*8 notation.
>
> That one I disagree with. Nothing wrong with such numbers as arguments
> to selected_real_kind. That's the number of decimal digits you want and
> is pretty much the normal use of selected_real kind. I normally put the
> selected_real_kind invocation in the definition of a named constant in a
> module instead of using the whole long-winded thing in each type
> definition, but otherwise literal constants are exactly what I tend to
> use with it.
>
> I suspect that the particular numbers used (4 and 8) are confusing you
> into thinking that they are kind values instead of arguments to
> selected_real_kind. Of course, I'd posit that, considering the number of
> other errors here, and that 4 is a pretty small number of decimal digits
> to ask for, the OP quite likely has this same confusion and actually
> meant to use the 4 and 8 as kind values, which is indeed non-portable
> and non-advisable.
Yes indeed. The 4 and 8 fooled me. I misread the meaning.
The error message when compiling with g95 only made it worse for me:
C:\Users\epc\temp>g95 x.f95
In file x.f95:43
call pbc(x,y,z,boxl)
1
In file x.f95:81
subroutine pbc(xx,yy,zz,boxl)
2
Warning (155): Inconsistent types (REAL(8)/REAL(4)) in actual argument
lists at
(1) and (2)
OOPS.
- e
e p chandler
[Please excuse reply to self.]
The declarations also helped reel me in. :-(((.
real(selected_real_kind(8)),dimension(100)::x,y,z
real(selected_real_kind(8)),dimension(100):: gr,rrl
1. I think if the OP had written
integer,parameter :: dp = selected_real_kind(8)
real(dp),dimension(100)::x,y,z
real(dp),dimension(100)::gr,rrl
I _would_ have let it go.
2. Since I'm swimming against the current today, I'll also express my
opinion that using dimension as the OP did is not such a bad idea.
IMO it highlights a common structure and makes the code more compact.
-- e
rudra
> e p chandler <e...@juno.com> wrote:
>but your actual arguments are the names of arrays.
> So even if it does compile and run you are probably passing it
>x(1),y(1),z(1) not x(i),y(i),z(i).
Since i am calling the pbc subroutine inside the do loop, i expect it
to check if the coordinate is greater then the length
plz correct me if i am wrong!!
thnx for pointing out other mistakes!!
and your 9th point will be really helpful for beginners like us
Richard Maine
> On May 16, 8:07 pm, nos...@see.signature (Richard Maine) wrote:
> > e p chandler <e...@juno.com> wrote:
> >5. pbc() is written to accept simple variables (scalars) as arguments.
> >but your actual arguments are the names of arrays.
> > So even if it does compile and run you are probably passing it
> >x(1),y(1),z(1) not x(i),y(i),z(i).
>
> Since i am calling the pbc subroutine inside the do loop, i expect it
> to check if the coordinate is greater then the length
> plz correct me if i am wrong!!
That's wrong on two points of substance (I'll not count the terminology
issues).
First, array bounds checking is not guaranteed to happen at all. That is
a debugging option provided by some (well, most) compilers to help find
programming errors. But it is not part of the language, and they are
programming errors.
Second, and more fundamental, I suggest you actually read what ep wrote
and look at the code. Your code didn't have *ANY* index values - in
bounds or not. No, just because you are doing something in a loop, that
doesn't somehow magically transform every array reference into an array
element reference. If you want x(i), you actually have to type x(i).
e p chandler
> here is the full code that is creating the problem!!! while running
> with ifort, the configuration is created as wished, but not in
> gfortran....
>
[rest of code snipped]
> while running with ifort, its prompting:
>
> $ ./a.out
> ===============Creating Initial Configuration==============
> Minimum distance between any two particle is 1.5000
> Success after 1680724 try
> Calculating minimum distance between atoms.... DONE
> Shortest distance between to atom is 1.626997
> between particle pair 40 - 92
> ===========================================================
>
> Elapsed Time = 0.0623
>
> but while running with gfortran, its never success...why?
OK. Time to deal with the _substantive_ issues in your program.
Random number routines may be built-in to your compiler's library, but
they vary from compiler to compiler. Even routines with the same
interface may output different results.
Examine the following session (most recent MinGW g95 & gfortran under
Win XP):
C:\temp>type x1.f95
implicit none
integer seed
integer i
real r
seed = 12345
do i=1,10
r = ran(seed)
print *,i,seed,r
end do
end
C:\temp>g95 x1.f95
In file x1.f95:9
r = ran(seed)
1
Error: Function 'ran' at (1) has no implicit type
G95 does not recognize RAN() but gfortran does.
C:\temp>gfortran x1.f95
C:\temp>a
1 12345 9.66165066E-02
2 12345 9.66165066E-02
3 12345 9.66165066E-02
4 12345 9.66165066E-02
5 12345 9.66165066E-02
6 12345 9.66165066E-02
7 12345 9.66165066E-02
8 12345 9.66165066E-02
9 12345 9.66165066E-02
10 12345 9.66165066E-02
Note that the generator returns the exact same value on every draw.
1. You need to read the manual.
2. Even after reading the manual, it is worth running a small test
program as given above to see what the RNG *actually* does!
RAN() may be an alias for RAND()
C:\temp>type x2.f95
implicit none
integer seed
integer i
real r
seed = 12345
do i=1,10
r = rand(seed)
print *,i,seed,r
end do
end
C:\temp>g95 x2.f95
C:\temp>a
1 12345 0.28407806
2 12345 0.28407806
3 12345 0.28407806
4 12345 0.28407806
5 12345 0.28407806
6 12345 0.28407806
7 12345 0.28407806
8 12345 0.28407806
9 12345 0.28407806
10 12345 0.28407806
G95 outputs different results.
C:\temp>gfortran x2.f95
C:\temp>a
1 12345 9.66165066E-02
2 12345 9.66165066E-02
3 12345 9.66165066E-02
4 12345 9.66165066E-02
5 12345 9.66165066E-02
6 12345 9.66165066E-02
7 12345 9.66165066E-02
8 12345 9.66165066E-02
9 12345 9.66165066E-02
10 12345 9.66165066E-02
It looks like they are the same.
1. Again note that the output stays the same from draw to draw.
Read the manual to see how to draw a stream of random numbers.
2. Later on in your program you call a routine that supposedly
restricts the absolute value of x,y and z to box1. WHY???
Most routines that generate uniform random numbers (like RAN or RAND
or whatever) return a real number approximately uniformly distributed
between 0 and 1. Usually 1 is excluded but 0 is included. (But not
always.). A simple linear transform A + B*X allows you to change the
scale and location of your random variable to something else. For more
complicated distributions you MAY need to generate some values and
REJECT part of them. You don't need to do that here.
Suppose you wanted to generate values between -box1 and box1. If X is
uniform on (0,1) then 2*X - 1 is uniform on (-1,1). Multiply by box1
and voila.
-- e
Gib Bogle
simple cases writing output to see what the code is actually doing.
Even after more than 20 years of coding Fortran I still make silly
mistakes that are (mostly) easily caught when I look at intermediate
printout.
Thomas Koenig
RANDOM_NUMBER instead of RAND or RAN. The pseudo-random number
generator is of much better quality.
If you want a good seed, you can use the example routine in
the documentation of the RANDOM_SEED subroutine.
Ron Ford
>> This seems much more i[n]formative than the entry in gfortran.pdf.
Fair enough. I think it's a bit of an uphill battle with the documentation
on the gnu project. As I haven't contributed anything to the project, I
don't feel justified to criticize. If someone told me to #define RAN RAND,
I'd just get C flashbacks. I commented on it because it seemed like OP had
an issue that wasn't addressed by the manual, and was being emphatically
told to the contrary.
For whatever computer time I've had today, I've been trying to get the
upper hand on gdb. Am I correct that I must compile a whole bunch of files
in order for gdb.exe to exist on my machine? I made the download from
http://www.mingw.org/download.shtml .
This is how gdb.c reads:
#include "defs.h"
#include "main.h"
#include "gdb_string.h"
#include "interps.h"
int
main (int argc, char **argv)
{
struct captured_main_args args;
memset (&args, 0, sizeof args);
args.argc = argc;
args.argv = argv;
args.use_windows = 0;
args.interpreter_p = INTERP_CONSOLE;
return gdb_main (&args);
}
Many of the instructions in readme.txt seem to be for a differing command
prompt than winxp's.
--
Ron Ford
FX
> upper hand on gdb. Am I correct that I must compile a whole bunch of
> files in order for gdb.exe to exist on my machine?
No, you should download the executable. It comes with the binaries I
distribute at http://gcc.gnu.org/wiki/GFortranBinaries, for example.
--
FX
e p chandler
> here is the full code that is creating the problem!!! while running
> with ifort, the configuration is created as wished, but not in
> gfortran....
[snip]
> real::t1,t2
On G95, secnds() subtracts t1 from the elapsed time since midnight, so
t1 should be set to 0 on the first call.
t1=0
> t1=secnds(t1)
Otherwise later on
t2=secnds(t1) may return wild (NaN) values.
> dist=dsqrt(dx*dx+dy*dy+dz*dz) [and later uses of dsqrt]
This does work, but only when selected_real_kind(8) is the kind of
double precision.
Fortran 77 and later offer generic functions. In this case it is
sqrt(). So it is best to use the generics.
1. You might scale back to single precision.
2. What happens if the default integer is 64 bits?
- e
Ron Ford
Well, goody gumdrops. The download and install are painless. I'm a little
surprised that set >before.txt and after.txt were identical *and* I can
have the gfortran and gdb commands understood in a place other than the
bin. --version yields 4.4.0 experimental, so I don't need to argue.
I really like the tool I use to analyze the difference between two files,
in particular files that have a lot of extraneous, computer-generated data.
It's the soln to K&R 7.6. clc-wiki now has the solns to K&R material
reviewed by Heathfield, so it's rigorous. It's a way to deal with K&R
material without having to deal with an unfriendly ng. For interested
parties (it looks like there's 4 lines with usenet fold-over):
// kr76.c
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define BUFF_SIZE 1000
/* uses fgets, removes the '\n' at the end of the string if it exists */
char *safegets(char *buffer, int length, FILE *file)
{
char *ptr;
int len;
if (buffer != NULL)
{
ptr = fgets(buffer, length, file);
if (ptr != NULL)
{
len = strlen(buffer);
if (len > 0)
{
if (buffer[len - 1] == '\n')
{
buffer[len - 1] = '\0';
}
}
}
return ptr;
}
return NULL;
}
int main(int argc, char *argv[])
{
FILE *leftFile, *rightFile;
char buff1[BUFF_SIZE], buff2[BUFF_SIZE];
char *ptr1, *ptr2;
unsigned long lineNum = 0;
if (argc < 3)
{
fprintf(stderr, "Usage : 7_6 <path to file> <path to
file>\n");
return 0;
}
if (!(leftFile = fopen(argv[1], "r")))
{
fprintf(stderr, "Couldn't open %s for reading\n",
argv[1]);
return 0;
}
if (!(rightFile = fopen(argv[2], "r")))
{
fprintf(stderr, "Couldn't open %s for reading\n",
argv[2]);
fclose(leftFile); /* RJH 10 Jul 2000 */
return 0;
}
/* read through each file, line by line */
ptr1 = safegets(buff1, BUFF_SIZE, leftFile);
ptr2 = safegets(buff2, BUFF_SIZE, rightFile);
++lineNum;
/* stop when either we've exhausted either file's data */
while (ptr1 != NULL && ptr2 != NULL)
{
/* compare the two lines */
if (strcmp(buff1, buff2) != 0)
{
printf("Difference:\n");
printf("%lu\t\"%s\" != \"%s\"\n", lineNum, buff1,
buff2);
goto CleanUp;
}
ptr1 = safegets(buff1, BUFF_SIZE, leftFile);
ptr2 = safegets(buff2, BUFF_SIZE, rightFile);
++lineNum;
}
/*
* if one of the files ended prematurely, it definitely
* isn't equivalent to the other
*/
if (ptr1 != NULL && ptr2 == NULL)
{
printf("Difference:\n");
printf("%lu\t\"%s\" != \"EOF\"\n", lineNum, buff1);
}
else if (ptr1 == NULL && ptr2 != NULL)
{
printf("Difference:\n");
printf("%lu\t\"EOF\" != \"%s\"\n", lineNum, buff2);
}
else
{
printf("No differences\n");
}
CleanUp:
fclose(leftFile);
fclose(rightFile);
return EXIT_SUCCESS;
}
// gcc -o kr kr76.c
// kr before.txt after.txt 2>text41.txt >text42.txt
--
Ron Ford