Formatted Read Statements yet again.

Further to my last post, I want to add a code snippet for reading in large arrays with many columns. It has some nice implicit do-loops. "myfile.ans" contains 1 column of integers, followed by Ncols of data that I want to read into "Blam" and Ncols that I want to read into "Dlam".

integer :: i, Nrows, Ncols
real :: Blam(Nrows, Ncols), Dlam(Nrows, Ncols)
integer :: pixies(Nrows)
character(1) :: tmp

open(unit=2,file='myfile.ans',status='OLD',action='READ')
do i = 1,Nrows,
read(2,20) pixies(i), (tmp,Blam(i,j), j=1,Ncols), (tmp,Dlam(i,j), j=1,Ncols)
enddo
close(2)
20 format(I3, 36(A1,F6.4))

Annoyingly, I can't seem to replace 36 (=2*Ncols) with a variable name or my compiler shouts. A mystery for another time I think.

Formatted Read Statements again.

Grrr. I'm trying to read in a formatted file of integers such as:
1 1 1
10 10 10
100 100 100
(... three characters wide, one space between each column)
Fortran is giving me junk unless I format my read statement with
10 format(I3,A1,I3,A1,I3)
where I read the single space into a temporary variable declared as:
character(1) :: tmp
This is why I am in a bad mood with Fortran today. I put faith in it, and this is how it repays me. Plus, the fridge has frozen my tomatoes. I wonder if they are in league.

Compiling modules separately

When I have a module declared in mymodule.f90, I can compile it like so:

gfortran -m32 -c mymodule.f90

... I can also compile the main file this way:

gfortran -m32 -c mainfile.f90

Finally I must link it, possibly also with the NAG library:

gfortran -m32 mymodule.o mainfile.o /opt/NAG/fll3a21dfl/lib/libnag_nag.a -o mainfile

Fortran array initialization

This webpage here is very helpful for fortran arrays.

I'd forgotten that you can initialize vectors like this:

x = (/(i,i=0,maxx-1)/)

Instead of:

do i=1,maxx
x(i) = i-1
enddo

Addressing integer bits in Fortran

I want to find out about 8-bit integers (kind = 1) in gfortran.
Apparently they seem to be signed and stored as "two's complement". When you give the integer a value, you can view the bits like so:

integer(kind = 1) I
integer numbit(8) ! stores the bit values as 0 or 1's

I = 5
do j= 1,8
numbit(i) = ibits(I,j-1,1)
enddo
write(*,*) 'bits: ', numbit

...where here IBITS is extracting 1 bit from I at position j-1 and returning it as an integer with the bit value right-justified and the rest of the bits 0. So "1" at position j-1 is returned by ibits(I,j-1,1) as "00000001" which has integer value 1.

Alternatively, I could do this:

integer(kind = 1) I
integer numbit(8) ! stores the bit values as 0 or 1's

I = 5
do j= 1,8
numbit(i) = btest(I,j-1)
enddo
write(*,*) 'bits: ', numbit

...where BTEST returns .true. if the bit in position j-1 is "1", and .false. otherwise. This is implicitly converted to an integer by the gfortran compiler (but not without warnings at compile time), where .true. => 1 and .false. => 0.

Both of these return:
bits: 1 0 1 0 0 0 0 0

I guess the take-home message is that when you write out a bit-string, eg.
5 = [0101]
then the right-most bit is addressed as bit 0, and in general, bit number 1 is the second bit as you move left, etc.

So if you start with all bits set to zero, you can use IBSET to set certain bits to 1 and get whatever number you fancy...

I = 0 ! [00000000]
I = ibset(I,0) ! [00000001]
I = ibset(I,2) ! [00000101]
write(*,*) 'I = ', I

... and we get I = 5, as intended.

I guess this is a stupid post, but it wasn't obvious to me that the bit string would be numbered from the right starting at index 0. The default index for matrices in Fortran 95 starts at 1. Ho hum.

Incompatible ranks?!

Here's a fun problem that I have encountered. It might be a bug, or it might not (I've run out of patience with google-ing it now I have my program working).

I have a two dimensional matrix. I sum it down one dimension and then find the position of the maximum entry.

subroutine ionfinder(subframe,x_dim,y_dim, irow, icol)

! define imputs, and function type itself
integer, intent(in) :: x_dim, y_dim
integer, intent(in) :: subframe(x_dim,y_dim)
integer, intent(out) :: irow, icol(4) !ion positions

! the code follows...*********************
irow = maxloc(sum(subframe,1))
write(*,*) 'ions in row (y): ', irow

end subroutine ionfinder

...but this causes my compiler (gfortran) to shriek:

irow = maxloc(sum(subframe,1))
1
Error: Incompatible ranks 0 and 1 in assignment at (1)

Huuuurrrmmm. Seems a bit silly, but the following kludge fixes things:

subroutine ionfinder(subframe,x_dim,y_dim, irow, icol)

! define imputs, and function type itself
integer, intent(in) :: x_dim, y_dim
integer, intent(in) :: subframe(x_dim,y_dim)
integer, intent(out) :: irow, icol(4) !ion positions

! define other datum within function
integer tmp(1)

! the code follows...*********************
tmp = maxloc(sum(subframe,1))
irow = tmp(1)
write(*,*) 'ions in row (y): ', irow

end subroutine ionfinder

Well well well. I should probably investigate the cause further but I'm not going to because I have physics to do, not Fortran to fix.

Formatted Read Statements

This is not a very exciting problem or solution, but I shall post it anyway.
I was trying to read formatted data from a file like this:

0 9.99929445e-01
1 6.46409406e-05
2 3.80249651e-06
3 1.03999904e-06
4 4.24499610e-07
5 2.23999794e-07

I first tried the following:

integer :: myint(5)
real(kind = 10) :: myreal(5)

open(unit=2,file=filepath,status='OLD',action='READ')
do i = 1,5
read(2,20) myint(i),myreal(i)
end do
close(2)
20 format(I6,E15.8E3)

The I6 refers to the field width of the integer column. But alas, the output was this:

0 0.99992944E+001
1 0.64640941E+001
2 0.38024965E+001
3 0.10399990E+001
4 0.42449961E+001
5 0.22399979E+001

Oddly enough, everything was in order except that it would round powers from e-01 to e-09 up to E+001, and e-10 to e-19 became E+000 etc. Most odd. The solution was to change the format statement to skip the 2 blank spaces between the columns, which I had forgotten to do:

20 format(I6,2X,E15.8E3)

Simple really, but I don't know why it gives such a strange behaviour.

Oh Fortran, how I love thee.

Fortran's logical operators that operate on integers happen to be bitwise. I didn't care about this until now. Here's a summary for two integers, a and b, which can take values 0 or 1.

logical AND:

(ab)	iand(a,b)
00	0
01	0
10	0
11	1

logical OR:

(ab)	ior(a,b)
00	0
01	1
10	1
11	1

logical XOR:

(ab)	ieor(a,b)
00	0
01	1
10	1
11	0

logical NOT:

(a)	not(a)
0	-1
1	-2

WTF? I have been bitten by the bitwise bear! In order to get zeros and ones, I must use ieor(1,a) instead of not(a). Pah.

Also, the tables in this post look all shit in my browser but I can't be bothered to fix them when I have Fortran to do instead.

Today's post was brought to you by the numbers 0 and 1 and the letter F.

Another Fortran idiosyncracy bites me on the bum

I thought I'd have this program done by the end of the afternoon. So I start to write a function to calculate a nice simple Gaussian function. No NaNs and +Infs to contend with here, I thought. How wrong I was.
This time, I was finding that the intrinsic function exp(x) returned zero for x less than around -0.1E3. After some bashing, my compiler was sweet enough to tell me:

Result of EXP underflows its kind

Apparently, Fortran is not content with just "exp". It also has "dexp" and "cexp" for double-precision and complex floating point numbers respectively. So, dexp(x) it is then. But be careful how you define x, or you will get a shouty like this:

Type of argument 'x' in call to 'dexp' at (1) should be REAL(8), not REAL(4)

So one should try dexp(-0.1D3) (note the letter D!) rather than exp(-0.1E3), which will return zero, possibly grumbling in the process.

NaN, +Inf, -Inf, horrid things...

This is the answer.
Based on such, I have a sample function to detect a +Inf:

function ispinf(x)
! function to detect if x is +Inf
integer :: ispinf
double precision :: x

!local variables
integer :: IPInf
real :: PInf
data IPInf/B'01111111100000000000000000000000'/ ! +Infinity
PInf = transfer(IPinf,Pinf)

if (x.eq.PInf) then
ispinf = 1
else
ispinf = 0
endif

end function ispinf

...in which I pacify the fussy NAG routine

As usual, this is of interest to probably no-one but me, but I have just found the bug in my program.
I am trying to call NAG routine G05MZF which generates random integers from a probability distribution supplied by you. I have found routines give unexpected answers if you don't pass them variables of the correct type, so I was careful about this. The error I was getting at run-time was this:

At line 106 of file /scratch/zohair/FLL3A21DF/fll3a21df/source/g/g05mzfn.f
Fortran runtime error: Missing initial left parenthesis in format

I have had errors like this before which have been due to the reference vector (which contains the CDF) being too small, but this wasn't the case here. Instead, it turned out to be angry about being passed what it considered to be a non-normalized PDF. For example, passing a PDF which was 307 entries long with each entry set to 1.0D0/307.0D0 (ie. a flat distribution) worked fine. However, setting each entry to the actual value, 0.32573290D-02 was not accurately normalised enough, and it threw a hissy fit.

So the moral of the story is:
If you think you've normalized your PDF, the NAG library knows better. Expressed in code, this is:
P = P/sum(P)

In other news, I find it incredible how easy it is to destroy my code. Simply compiling program.f90 to the output file program.f90 will replace my code with executable junk. Surely they should have idiot-proofed this?

*********************************

Edit: 20th March 2008

Actually there's more to it than this. Firstly, the bizarre and unhelpful error is due to using gcc version 4.1.2 rather than the recommended version 4.2 which I don't yet have. If you use version 4.3, the error is more helpful, and gives you the value of IFAIL and results in a NAG hard failure telling you the deviation of the probability sum from unity. There are other issues too, but using the right compiler is helpful!

Keyboard shortcuts for kwrite

As usual, this will be of very little interest to anyone but me, but it's actually made my day. I have accidentally discovered the keyboard shortcuts useful for editing fortran code in kwrite.

Ctrl+I - indents
Ctrl+Shift+I - un-indents
Ctrl + D - comments out code (by ! marks)
Ctrl + Shift + D - un-comments code

How exciting my day is today. I actually could have looked this up in the kwrite menu, so I am stupid too!

File I/O in Fortran 95

Can a "unit" be re-used after disconnecting from a previous file?
I believe the answer is yes! Two different units are only needed when trying to access two files at the same time. Once a file has been disconnected, one can re-use the number. Hoorah!

I can't even remember how to compile!

My previous programs which don't require linking to a NAG library seem to compile thus:

gfortran -m32 program.f90 -o program

If they do require linkage, I compile like this:

gfortran -m32 program.f90 /opt/NAG/fll3a21dfl/lib/libnag_nag.a -o program

Silly me. Won't forget that now, will I?

Converting integers to strings in Fortran

Assisted by this page, I present the result:

character(11), parameter :: fstart = "Xfiles/Xbin"
character(4), parameter :: fend =".txt"
character(3) :: fnum
character(18) :: fname
integer :: blob

blob = 2
write (fnum,1000) blob
1000 format(I3)
fname = fstart//trim(adjustl(fnum))//fend

The double slash "//" concatenates strings. The strange form of the write command (rather than just char(blob) gives you the value of blob as a character, rather than the ascii character for that number. The adjustl(fnum) function moves the contents of string "fnum" to the left, to leave all trailing blanks to the right. The call to trim(string) then removes the trailing blanks. This is necessary to avoid spaces appearing in "fname" if the value of blob is only 2 digits or less, since I defined "fnum" to be three characters long.

More NAG problemos

So, to re-cap, the NAG library wants me to have gfortran 4.2.0 or better. I installed the only ready-built new version of gfortran I could find for linux (since Mandriva doesn't offer a package newer than version 4.1.2), which was 4.3.0. I have the NAG library working fine with this new version of gfortran, and it cries if I try to use the old fortran compiler. To re-cap, I now have two compilers on my computer:

"New"
gcc version 4.3.0 20071121 (experimental)
"Old"
gcc version 4.1.2 20070302 (prerelease)

So all is splendid (when I use the routine S14BAF to find the incomplete gamma function) after a few angry words about data-types. You must use double-precision variables (not real(kind = 10) or anything) or else it just returns junk. This will only work if I compile with the new gfortran version.

However , (and there's always a "however" with me and linux!), now I try to compile a program which uses the random number routines G05MKF and similar. The new compiler now gets angry:

/opt/NAG/fll3a21dfl/lib/libnag_nag.a(g05mkxn.o): In function `g05mkxn_':
g05mkxn.f:(.text+0xde): undefined reference to `_gfortran_pow_r8_i4'
/opt/NAG/fll3a21dfl/lib/libnag_nag.a(g05mkzn.o): In function `g05mkzn_':
g05mkzn.f:(.text+0xb1): undefined reference to `_gfortran_pow_r8_i4'
collect2: ld returned 1 exit status

This subroutine "_gfortran_pow_r8_i4" is apparently supposed to be in libgfortran.a and many people have solved their problems by explicitly telling their compiler where to look for this when linking, eg.
/etc/irun/bin/gfortran -m32 simtest.f90 /opt/NAG/fll3a21dfl/lib/libnag_nag.a -o simtest -L/etc/irun/lib -lgfortran
But this still doesn't work. Anyway, I added this directory to the path before hand, so this probably is not the problem. I look at the contents of the libgfortran library using the command:
nm libgfortran.a |more
where I've piped it through the "more" doo-dah so it doesn't all flash by at once. There are lots of similar routines to "_gfortran_pow_r8_i4", but no actual version of the one I want. So my compiler is not lying! I have not so far found a new version of the library which contains this, but my old compiler library does! So for all my random number programs, I must still use the old compiler, which seems to work just fine. Perhaps I can add the subroutine from the old library into the new one? Or will there be other issues?
This page is quite helpful about libraries, I found...

NAG library- getting the thing to work!

So, I receive my license key, and I've installed the library that I needed (FLL3A21DFL, x86-32, Linux, gfortran compiler, Double Precision) and followed the installers instructions on this page. It installed in:

/opt/NAG/fll3a21dfl/

It's now time to run a test! The User's Note suggests to run an example, like thus:
nag_example e04ucf
...for which I moved to my own directory where I had write permission. However, then I need to tell it where the nag_example script is, so I actually had to run:
/opt/NAG/fll3a21dfl/scripts/nag_example e04ucf
The computer said:

Copying e04ucfe.f to current directory
cp /opt/NAG/fll3a21dfl/examples/source/e04ucfe.f .
Compiling and linking e04ucfe.f to produce executable e04ucfe.exe
gfortran -m32 e04ucfe.f /opt/NAG/fll3a21dfl/lib/libnag_nag.a -o e04ucfe.exe
Copying e04ucfe.d to current directory
cp /opt/NAG/fll3a21dfl/examples/data/e04ucfe.d .
Running e04ucfe.exe with data from e04ucfe.d
./e04ucfe.exe < e04ucfe.d > e04ucfe.r
/opt/NAG/fll3a21dfl/scripts/nag_example: line 58: 12681 Segmentation fault ./${example}.exe <${example}.d >${example}.r

...as it happily seg-faulted. I guessed this is because it called "gfortran" which still refers to the old version of the compiler on my machine (I still haven't worked out how to change that!) But now I guess the correct commands to re-compile it:
/etc/irun/bin/gfortran -m32 e04ucfe.f /opt/NAG/fll3a21dfl/lib/libnag_nag.a -o e04ucfe.exe
Then to run it:
./e04ucfe.exe < e04ucfe.d > e04ucfe.r
And it cried:

./e04ucfe.exe: error while loading shared libraries: libgfortran.so.3: cannot open shared object file: No such file or directory

No matter, I had forgotten to type
export LD_LIBRARY_FLAG=/etc/irun/lib
export LD_LIBRARY_PATH=/etc/irun/lib
...because I still haven't managed to sort this problem out a better way. Perhaps making a symbolic link to the missing libraries as suggested in the NAG Installer's notes might work:

To simplify access to the libraries you may wish to create symbolic links in a system location such as /usr/lib pointing at the installed libraries. They would then be in the default search path of the linker during the link phase, and be available for execution at run time (in the case of shareable libraries).

Anyway, re-compiling and running now works fine. Time to call a routine from my program. Since this is the NAG library (and not the NAG f90 library), I must call functions with strange names like S14AAF. After adding the compilers own libraries to LD_LIBRARY_PATH (see above!), I additionally added the NAG libraries:
LD_LIBRARY_PATH=/opt/NAG/fll3a21dfl/lib:/opt/NAG/fll3a21dfl/acml:
${LD_LIBRARY_PATH}
but really all on one line, and with no spaces. Then...
export LD_LIBRARY_PATH
Such that typing:
echo $LD_LIBRARY_PATH
Gives me:

/opt/NAG/fll3a21dfl/lib:/opt/NAG/fll3a21dfl/acml:/etc/irun/lib

I've no idea if that was strictly necessary. (Edit- it doesn't seem to be!) In my fortran code, I must declare the types of the NAG functions I call, as well as the fact that they are external:
real(kind = 10) :: factorial, S14AAF, S14BAF
external S14AAF, S14BAF
Where factorial is an internal function I wrote myself. Functions are then called as you would expect...
gam = S14AAF(n+1.0,ifail)
or somesuch line. The program is compiled using:
/etc/irun/bin/gfortran -m32 gammatest.f90 /opt/NAG/fll3a21dfl/lib/libnag_nag.a -o gammatest
...since I had to tell the compiler where to find the external functions.

Fun and games with gfortran and the NAG library

I am trying to install the NAG library at work to run with gfortran. The problem is that it won't work with the old version of gfortran (4.1.2), which is the only package available for my version of the Mandriva linux kernel (2007). I had to update gfortran to 4.2.0 or better. First I needed to find out which kernel version I had (32 bit or 64 bit on a 64-bit processor?). I typed:
cat /etc/issue
which seemed to tell me that I have a 32 bit kernel. An alternative is to try:
uname =a
I then downloaded the 32 bit gfortran binaries for Linux, and saved them in a temporary file. I then changed to the directory in which I wished to install them:
cd /etc
Next I un-tarred them:
tar -zxvf $temp_folder_path/gfortran-linux.tar.gz
Typing:
/etc/irun/bin/gfortran -v
told me all about the new version I had installed. When trying to compile my test code, I had to explicitly give the path to the new version, otherwise it used the old version (haven't worked out how to change this yet!)
/etc/irun/bin/gfortran -o test test.f90
So far, so good, but when I tried to run it:
./test
The computer screamed:

./test: error while loading shared libraries: libgfortran.so.3: cannot open shared object file: No such file or directory

Typing:
ldd test
showed me which libraries the program calls, and which (might) be missing. To solve this:
export LD_LIBRARY_FLAG=/etc/irun/lib
and similarly...
export LD_LIBRARY_PATH=/etc/irun/lib
To see what the variables contain:
echo $LD_LIBRARY_PATH
More fun to come, as I try and get the NAG library to work...

Fortran 95

Here I am, installing Fortran 95 compilers at work on Mandriva One. It didn't seem to be installed, so I just installed the two packages pertaining to Fortran from rpmdrake, and now my file test.f90 compiles when I type:

f95 -o test test.f90

and runs with:

./test

Hoorah for that.