PGI Folks,
I have an…oddity I need help explaining. I’m trying to figure out a nice, portable way to get random numbers with different compilers. Obviously, RANDOM_NUMBER is the way to go. Now, the code I’m using will also need a static seed, so I’ll need to use RANDOM_SEED as well. Following a chunk of code I found via Google at Lahey, I made this program:
program randomTest
implicit none
real*8 :: rn(10)
integer :: seed_size
integer, allocatable :: seed(:)
integer :: i
call random_seed()
call random_seed(size=seed_size)
write (*,*) 'seed_size for RANDOM_SEED is: '
write (*,*) seed_size
allocate(seed(seed_size))
call random_seed(get=seed)
write (*,*) 'System generated RANDOM_SEED is: '
write(*,*) seed
seed=12
!seed = [ (i,i=1,seed_size) ]
call random_seed(put=seed)
call random_seed(get=seed)
write (*,*) 'Our specified RANDOM_SEED is: '
write(*,*) seed
deallocate(seed)
call random_number(rn)
write (*,*) rn
end program randomTest
When I compile and run this code I get a…curiosity:
$ pgfortran randomtest.F90 && ./a.out
seed_size for RANDOM_SEED is:
34
System generated RANDOM_SEED is:
1120731 7216334 1854449 4125410 8077957 55766
6275042 7511160 8086732 6595401 175655 3844873
2093857 1818440 6712803 7103074 5115164 839058
4556604 467731 284012 2784412 2530688 7365308
6706998 3320922 8209924 7768441 290898 3262263
7108656 1804980 2305714 7339660
Our specified RANDOM_SEED is:
12 12 12 12 12 12
12 12 12 12 12 12
12 12 12 12 12 12
12 12 12 12 12 12
12 12 12 12 12 12
12 12 12 12
2.8610232902792632E-006 2.8610232902792632E-006 2.8610232902792632E-006
2.8610232902792632E-006 2.8610232902792632E-006 4.2915349354188947E-006
4.2915349354188947E-006 4.2915349354188947E-006 4.2915349354188947E-006
4.2915349354188947E-006
Huh. That’s not that random…looking (as xkcd taught us, it still might be random). It turns out if you increase the size of rn (say, dimension 200), it does start “looking” random:
$ pgfortran randomtest.F90 && ./a.out
seed_size for RANDOM_SEED is:
34
System generated RANDOM_SEED is:
4193689 4181511 1633777 4958531 6587629 6505414
44629 3833330 1937304 2133152 2014503 6498151
3342430 7256452 3886428 1245379 6241516 6257177
887417 3292396 3900654 5533165 3650824 2752512
6321256 2223711 996759 5180850 230618 2486501
1362322 1675880 470661 4135704
Our specified RANDOM_SEED is:
12 12 12 12 12 12
12 12 12 12 12 12
12 12 12 12 12 12
12 12 12 12 12 12
12 12 12 12 12 12
12 12 12 12
2.8610232902792632E-006 2.8610232902792632E-006 2.8610232902792632E-006
2.8610232902792632E-006 2.8610232902792632E-006 4.2915349354188947E-006
4.2915349354188947E-006 4.2915349354188947E-006 4.2915349354188947E-006
4.2915349354188947E-006 5.7220465805585263E-006 5.7220465805585263E-006
5.7220465805585263E-006 5.7220465805585263E-006 5.7220465805585263E-006
7.1525582256981579E-006 7.1525582256981579E-006 8.5830698708377895E-006
8.5830698708377895E-006 8.5830698708377895E-006 1.0013581515977421E-005
1.0013581515977421E-005 1.2874604806256684E-005 1.2874604806256684E-005
1.2874604806256684E-005 1.4305116451396316E-005 1.4305116451396316E-005
1.8596651386815211E-005 1.8596651386815211E-005 1.8596651386815211E-005
2.0027163031954842E-005 2.0027163031954842E-005 2.5749209612513368E-005
2.5749209612513368E-005 2.7179721257653000E-005 2.8610232902792632E-005
2.8610232902792632E-005 3.5762791128490790E-005 3.5762791128490790E-005
4.0054326063909684E-005 4.1484837709049316E-005 4.1484837709049316E-005
5.0067907579887105E-005 5.0067907579887105E-005 5.8650977450724895E-005
6.0081489095864526E-005 6.0081489095864526E-005 7.0095070611841948E-005
7.0095070611841948E-005 8.4400187063238263E-005 8.5830698708377895E-005
8.7261210353517527E-005 9.8705303514634579E-005 9.8705303514634579E-005
1.2016297819172905E-004 1.2159348983686868E-004 1.2731553641742721E-004
1.4019014122368390E-004 1.4019014122368390E-004 1.7023088577161616E-004
1.7166139741675579E-004 1.8596651386815211E-004 2.0027163031954842E-004
2.0027163031954842E-004 2.4032595638345811E-004 2.4175646802859774E-004
2.7036670093139037E-004 2.8610232902792632E-004 2.8753284067306595E-004
3.3903125989809269E-004 3.4046177154323232E-004 3.9052967912311942E-004
4.0769581886479500E-004 4.1484837709049316E-004 4.7922140112177658E-004
4.8065191276691621E-004 5.6076056489473558E-004 5.7935721628155079E-004
6.0081489095864526E-004 6.7949303144132500E-004 6.8092354308646463E-004
8.0108652127819369E-004 8.2111368431014853E-004 8.7118159189003563E-004
9.6559536046925132E-004 9.6845638375953058E-004 1.1401177811762864E-003
1.1615754558533808E-003 1.2617112710131551E-003 1.3732911793340463E-003
1.3833047608500237E-003 1.6193391822980630E-003 1.6422273686202971E-003
1.8224718359078906E-003 1.9526483956155971E-003 1.9841196518086690E-003
2.2988322137393880E-003 2.3231509117067617E-003 2.6235583571860843E-003
2.7737620799257456E-003 2.8553012436987046E-003 3.2644275742086393E-003
3.2916072954662923E-003 3.7636761383623707E-003 3.9353375357791265E-003
4.1170125147118597E-003 4.6377187535426856E-003 4.6749120563163160E-003
5.3830153206604336E-003 5.5775649043994235E-003 5.9394843506197503E-003
6.5903671491582827E-003 6.6590317081249850E-003 7.6818475343998216E-003
7.9007158161061852E-003 8.5630427078058347E-003 9.3641292290840283E-003
9.5143329518236897E-003 1.0946275108608461E-002 1.1192323111572478E-002
1.2326718846168205E-002 1.3299466764863155E-002 1.3631345466535549E-002
1.5583993862151146E-002 1.5867235167888794E-002 1.7709734166828639E-002
1.8877031669262578E-002 1.9570829817155300E-002 2.2174361011309429E-002
2.2526266876013779E-002 2.5391581701228461E-002 2.6777747485368764E-002
2.8133872524961134E-002 3.1538490240393457E-002 3.2040599827837468E-002
3.6337856809836921E-002 3.7970070596941241E-002 4.0460591371129340E-002
4.4837957005256612E-002 4.5671945294373018E-002 5.1921850671988068E-002
5.3837305764830035E-002 5.8170325537957979E-002 6.3714988674519191E-002
6.5242775111528317E-002 7.4096211683297497E-002 7.6363572640843813E-002
8.3561907239186439E-002 9.0492736159887954E-002 9.3376647636489452E-002
0.1056347019236910 0.1084041724686813 0.1198997640490234
0.1284628067568292 0.1338372390076188 0.1504726589289476
0.1540761177630543 0.1718216147210114 0.1823001125216592
0.1920075645455768 0.2141876476034668 0.2193188928745826
0.2459178264043089 0.2586636851625030 0.2755694717847632
0.3046803837633547 0.3126955405110721 0.3515525283279999
0.3670678576311843 0.3954692358337866 0.4331431905201839
0.4465327795186909 0.5020251872569474 0.5211439753942386
0.5672908505547980 0.6154433030418431 0.6385403440642676
0.7162128348604142 0.7404628682688212 0.8132086769591069
0.8741069882043462 0.9141098158490308 2.0893218623768917E-002
5.3158408779893307E-002 0.1647612052871068 0.2411748458355305
0.3095790516828174 0.4540364091439528 0.4996911882985842
0.6667863925440543 0.7623188212297691 0.8768699022376154
6.9479712185795961E-002 0.1382315323628518 0.3829992274044685
0.5027816894985904 0.6900785791967223 0.9435867003901421
5.2341348211882632E-002 0.4038924460282374
After a while, the numbers aren’t the same anymore, which is good, but there seems to be a bit of a startup “bias” where the values monotonically increase until you reach 1.0 and then it “looks” random.
I’ll note that the generator in gfortran or ifort don’t seem to have the same issue:
$ gfortran randomtest.F90 && ./a.out
seed_size for RANDOM_SEED is:
8
System generated RANDOM_SEED is:
437395160 1404128605 572505362 -1187264075 454383258 525702629 973594203 1758310677
Our specified RANDOM_SEED is:
12 12 12 12 12 12 12 12
0.75896909164940396 0.48897481984547464 0.79513291277748022 0.98009539185978756 0.83467919508698862 0.85101238952040359 0.16483897501956624 0.14000430522166651 3.96003308392569586E-002 0.46935267743450682
$ ifort randomtest.F90 && a.out
randomtest.F90(23): (col. 1) remark: LOOP WAS VECTORIZED.
seed_size for RANDOM_SEED is:
2
System generated RANDOM_SEED is:
1351962852 156569157
Our specified RANDOM_SEED is:
12 12
0.999996210913953 0.694234916013649 0.769806140769980
0.997026373980212 0.443334298526595 0.940541781459957
0.975740233872978 1.607431767802543E-002 0.449636255958620
7.682652516768049E-003
I’m guessing the Standard doesn’t specify what the random number generator for an implementation is, just how it will be called, since each compiler gives different answers. Is there a “best practices” guide for calling the PGI RANDOM_NUMBER routine?
ETA: Oh, and I should have said this, but obviously, if your RANDOM_SEED is pretty random (like the system seed), pgfortran “looks” random:
$ pgfortran randomtest.F90 && ./a.out
seed_size for RANDOM_SEED is:
34
System generated RANDOM_SEED is:
8093102 2880621 3485957 6920646 7079926 8288461
7336175 2932069 7005763 4518683 2351199 4098661
2355125 3674405 7170886 2951714 306196 743387
188311 1938423 6027605 4458888 6263290 2424047
1007418 1494285 2611832 3318276 7764244 3167453
6653661 2773408 7100761 5707740
0.5215295654376177 0.2205747040406578 0.3656515800431919
0.6801460585258496 0.2190848209415321 1.0128064769659773E-002
0.6585979470808354 0.7175234108447057 0.7687646990084431
0.4501628522964864
That’s good, but usually when one specifies a seed, it’ll just be some simple number like 1. You’ll note from my code stub, I tried:
seed = [ (i,i=1,seed_size) ]
but even that still has the “spin-up” problem, just not as severe (rn(1) and rn(2) aren’t the same, say).
Thanks,
Matt