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.

Installing an Epson scanner under Mandriva

Another how-to, brought to you by Mr B.:

To install and use the Epson Perfection 2480 scanner under Mandriva:

1) install "sane" (in package manager, install "libsane", "sane-backends" and perhaps "sane-frontends").

2) login as root, go to /etc/sane.d/ and edit the file "dll.conf" and ensure the line "snapscan" is uncommented (without a # at the beginning of the line).

3) from the CD supplied with the scanner, find the file "ESCAN/ModUsd.cab" and copy to your desktop (or elsewhere!).

4) go to the directory on your hard drive containing "ModUsd.cab" and run

$ cabextract ModUsd.cab

cabextract is a program to extract files from a micro$oft .cab format, which I have described in a previous post.

5) login as root and copy the firmware file "Esfw41.bin" to /etc/sane.d/esfw41.bin

6) run the following (as root):

$ scanimage -L

...and the scanner should be detected!

7) use "Kooka" or "Xscanimage" to use scanner (Xscanimage is from the package "sane-frontends")

Extracting Micro$oft .cab files

Here it is, courtesy of Mr B.:

To extract Microshite .cab files, install "cabextract" (from package manager).
Copy blah.cab into a directory and:

$ cabextract blah.cab

OR

$ cabextract -l blah.cab

to list the files in the .cab.

Fixing dvips

We have a lovely shiney (the case really is!) new computer named Lewis. (This is after Lewis Hamilton, because he is so fast).

Anyway, it was super fast to build and amazingly the whole kit works, even the processor fan started so our new core2duo didn't melt. Mandriva 2008.1 installed quicker than it takes our old laptops to boot up. Now I have to iron out all the little problems getting things to work. (Like why Firefox is thinking every word in the english language that I'm typing is spelt wrong!!!).

First up, I can't seem to convert dvi to ps in Kile. Actually, it doesn't work on the command line either- it says:

dvips: ! Couldn't find header file cm-super-t1.enc.
Note that an absolute path or a relative path with .. are denied in -R2 mode.

R2 is some sort of secure mode. I am using texlive. Part of the problem is that dvips is correct- I really didn't have this file. I installed the following package:

texlive-texmf-cmsuper

But dvips still couldn't see it where it was hiding (which was in /usr/share/texmf-texlive/fonts/enc/dvips/cm-super). It seems that dvips uses the same search-y doo-dah as TeX, so all I needed to do was update the lists by running

texhash

as a superuser. It's happy now.

Povray

Povray is amazing. This post is really just to remind myself of how to render an image with anti-aliasing from the command line:

povray myimage.pov -geometry 640x480 +A0.1

Hoorah!

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!

Installing a .sty file!

At last, I succeed in one thing!
I needed bbding.sty, so I downloaded the whole zip from CTAN. It included the .ins and .dtx files. I stuck it in a local temporary place and followed the instructions on this handy page.
I ran latex on the .ins file and then on the .dtx file. I then moved the whole folder and everything in it to my local texmf tree:

/usr/local/share/texmf/tex/latex/

into a folder of the same name, ie. bbding. Next, as a superuser I ran:
texhash
and we all lived happily ever after.

Merging PDF files

Assisted by this page, I managed to create a handy pdf file from the six .jpg files of sheet music I'd been sent. First I opened them in showFoto and printed them to six .pdf files, and then merged them on the linux command-line thus:
gs -dNOPAUSE -sDEVICE=pdfwrite -sOUTPUTFILE=Merged.pdf -dBATCH 1.pdf 2.pdf 3.pdf 3.pdf
etcetera, etcetera, etcetera...

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?