Special Hints

1.

The following example shows how to use HDF image writing utilities in an HPF program:

PROGRAM jacobi

  IMPLICIT NONE

!HPF$ NOSEQUENCE

  ! All these are parallel variables with "field" and "mask" distributed
  ! blockwise.

  INTEGER, PARAMETER :: n = 200, iterations = 2000
  REAL, DIMENSION(n), PARAMETER :: north_boundary = 30.0, &
       east_boundary = 80.0, west_boundary = 80.0, south_boundary = 250.0
  REAL, DIMENSION(n, n) :: field = 30.0
  LOGICAL, DIMENSION(n, n) :: mask = .TRUE.

!HPF$ ALIGN mask(:,:) with field
!HPF$ DISTRIBUTE (*, BLOCK) :: field

  ! The following are HDF related variables.

  INTERFACE
     EXTRINSIC(hpf_serial) FUNCTION d8aimg(file_name, raster, nx, ny, mode)
       INTEGER :: d8aimg
       INTEGER :: nx, ny, mode
       INTEGER(kind=1) :: raster(nx,ny)
       CHARACTER(len=*) :: file_name
     END FUNCTION d8aimg
     EXTRINSIC(hpf_serial) FUNCTION d8pimg(file_name, raster, nx, ny, mode)
       INTEGER :: d8pimg
       INTEGER :: nx, ny, mode
       INTEGER(kind=1) :: raster(nx,ny)
       CHARACTER(len=*) :: file_name       
     END FUNCTION d8pimg
  END INTERFACE

  INTEGER(kind=1), DIMENSION(n, n) :: raster
  INTEGER :: i, j, status
  INTEGER, PARAMETER :: comp_rle = 11

  field(UBOUND(field, dim=1), :) = east_boundary
  mask(UBOUND(mask, dim=1), :) = .FALSE.
  field(LBOUND(field, dim=1), :) = west_boundary
  mask(LBOUND(mask, dim=1), :) = .FALSE.
  field(:, UBOUND(field, dim=2)) = north_boundary
  mask(:, UBOUND(mask, dim=2)) = .FALSE.
  field(:, LBOUND(field, dim=2)) = south_boundary
  mask(:, LBOUND(mask, dim=2)) = .FALSE.

  raster=INT(field, 1)
  status=d8pimg('movie.hdf', raster, n, n, COMP_RLE)

  DO i = 1, iterations
     WHERE (mask)
        field = (EOSHIFT(field, 1, dim=1) + EOSHIFT(field, -1, dim=1) &
             + EOSHIFT(field, 1, dim=2) + EOSHIFT(field, -1, dim=2)) * 0.25
     END WHERE
     IF (MOD(i, 10) .EQ. 0) THEN
        raster = INT(field, 1)
        status = d8aimg('movie.hdf', raster, n, n, comp_rle)
        WRITE(*,*) i, ' -> ', (INT(field(j, 10)), j = 1, 5), &
             (INT(field(j, 10)), j=96, 105), &
             (INT(field(j, 10)), j=196, 200)
     END IF
  END DO

END PROGRAM jacobi

    Here is how to compile, link, and run this program:

gustav@sp20:../jacobi 11:14:45 !523 $ make jacobi-raster
xlhpf90 -o jacobi-raster jacobi-raster.f \
    -L/afs/ovpit.indiana.edu/@sys/HDF/lib -lmfhdf -ldf -ljpeg -lz -lm
** jacobi   === End of Compilation 1 ===
1501-510  Compilation successful for file jacobi-raster.f.
gustav@sp20:../jacobi 11:14:51 !524 $ make run-raster
poe jacobi-raster -euidevice css0 -procs 8 -euilib ip -ilevel 3 \
    > jacobi-raster.log 2> jacobi-raster.err
gustav@sp20:../jacobi 11:25:33 !525 $

2.

There is a small but nasty bug in the IBM diffusion example,   which lives in

/usr/lpp/pessl/example/hpf

and which we've been working with. You must fix this bug, because otherwise you will not be able to write large images.

In main.f you will find the following code:

         do it = 1, dif_ntemps
          do i = 1, dif_npts
            iy = (i-1)/dif_nx+1
            ix = mod(i-1,dif_nx)+1
            do k = 1, dif_npts
              temp(k,it) = temp(k,it) + xsines(k,ix) * ysines(k,iy) &
                            * x(i,it)
            enddo
           enddo
          enddo

towards the end of the file. This is the part of the program where final temperatures are calculated for various times. The index i in the second DO loop runs from 1 through dif_npts, which is incorrect. It should run from 1 through n, because that index numbers harmonics, and we have $n = n_x \times n_y$ harmonics altogether. The variable that is used in the code, dif_npts, stores the numbers of points of the physical grid used, e.g., for writing images.

If your image is, say, $200\times200$, then dif_npts will be 40,000. On the other hand the default number of harmonics is 7 in each direction. So n = 49, which is a lot less than 40,000. Trying to run that second loop from 1 through 40,000 will invariably result in incorrect memory references.

In summary that part of the code should look as follows:

         do it = 1, dif_ntemps
          do i = 1, n
            iy = (i-1)/dif_nx+1
            ix = mod(i-1,dif_nx)+1
            do k = 1, dif_npts
              temp(k,it) = temp(k,it) + xsines(k,ix) * ysines(k,iy) &
                            * x(i,it)
            enddo
           enddo
          enddo

I'd like to thank Xin Hong,   for having drawn my attention to this problem.

Code fragments presented in this document and in

/afs/ovpit.indiana.edu/common/www/htdocs/gustav/diffusion

are correct in this respect.

3.

When   you generate animations based on rasters it is important that colours on all generated images conform to the same colour map. What this means is that a pixel colour 255 must correspond to the highest temperature ever seen in the system (this is likely to be the highest temperature recorded on the first snapshot, because this whole system cools down), and some very low pixel colour, e.g., 10, should correspond to the lowest temperature seen in the system, i.e., 0. Once you find a correct scaling factor, you must then scale all data using the same scaling factor before or while writing the images.

  You should avoid pixel colours 0 and 1, because these often correspond to black and white, and you want a neatly graded display from red through dark blue, without any white and black spots or lines. Other low value colours may be bad too, e.g., they may or may not correspond to pure yellow, green, red, etc. You can find the lowest "safe" number simply by experimenting (or reading HDF manuals). In our jacobi.f example, I've used 30, which is safe enough.

The movie should look as follows. Initially there should be a reddish region somewhere inside the square dissolving gradually towards blue around the edges. Then as time goes on, the red should become gradually bluer until towards the end you'll have blue everywhere, i.e., all the heat should flow away from the beam.

Observe that in program jacobi.f I didn't have to do any scaling, because I have forced my temperatures to stay between 30 and 250 by the virtue of imposed boundary conditions. I could therefore map from temperature onto the raster image without any scaling at all. Of course I did it on purpose, to avoid cluttering the program with that kind of stuff.

But you will have to scale your data in the Spectral Method example.

4.

You will find PESSL manual on:  

/usr/lpp/pessl/doc/postscript/pessl.00.ps

on the SP. You can view it on-line with

/afs/ovpit.indiana.edu/@sys/X11R6.4/bin/ghostview

You can also copy it to a more convenient location, e.g., to your AFS directory.

5.

You will find HPF manuals in:  

/usr/lpp/xlhpf/postscript

There are three files there, which are all compressed. Use the command uncompress to uncompress them. You will have to copy them to your AFS directory to do that. There may not be enough space in your SP directory. The files are:

xlhpflr.ps.Z

HPF Reference Manual

xlhpfug.ps.Z

HPF User Guide

xlhpftune.ps.Z

HPF Tuning Guide

6.

You will find Fortran-90 manuals in  

/usr/lpp/xlf/postscript

Most importantly:

xlflr.ps.Z

Fortran Reference

xlfug.ps.Z

Fortran User Guide

xlftune.ps.Z

Fortran Tuning Guide

7.

You will find HDF manuals on  

/afs/ovpit.indiana.edu/common/HDF/doc/Users_Guide/ug41r2.pdf
/afs/ovpit.indiana.edu/common/HDF/doc/Ref_Manual/RM_4.1r2b.pdf

These have to be read with Acrobat Reader, acroread.

8.

It turns   out that you can run POE jobs from your AFS directory. This way you can avoid problems with insufficient disk space on the NFS and make it easier for yourselves to view the data generated on the SP with tools installed in the IRIX fileset of the AFS.

Sometimes AFS doesn't come up properly on some of the nodes that you have declared in your host.list file. If POE hits any of those nodes during execution, and the run has been started up from AFS, it'll bump out and complain. As it happened, AFS didn't run on nodes 40 and 42 the last couple of days and that has caused some problems.

Anyhow, here is how it works for me:

In my AFS home directory

/afs/ovpit.indiana.edu/home/gustav

I have created a directory called "jacobi" and made it AFS-open to all:

 $ cd /afs/ovpit.indiana.edu/home/gustav
 $ fs setacl jacobi -acl system:anyuser all
 $ fs listacl jacobi
    Access list for jacobi is
    Normal rights:
      system:anyuser rlidwka
      gustav rlidwka
 $ cd jacobi
 $ make clean
 $ make jacobi-raster
 $ make run-raster
    poe jacobi-raster -euidevice css0 -procs 8 -euilib ip -ilevel 3 \
        > jacobi-raster.log 2> jacobi-raster.err
 $ cd ..
 $ fs setacl jacobi -acl system:anyuser none
 $ fs listacl jacobi
    Access list for jacobi is
    Normal rights:
      gustav rlidwka

Perhaps a better way would be to create a directory called work in your AFS home and make it -acl system:anyuser none, and then create a directory jacobi inside work and make it -acl system:anyuser all. I've tested that too and it worked just fine. That way you don't have to be so paranoid about running fs setacl on your jacobi directory after you have finished your POE run.

Your AFS quota is currently 100MB per student. There is plenty of space on the AFS and that quota can be easily increased if you need more.

Please feel free to see me if you would like me to help you set up your AFS and SP environments.

9.

 Here is a yet another hint that relates to the second part of the exercise, that is, explaining the two code fragments   in subroutines get_diffusion_matrix and expand_temp_profile.

You need to study the discussion of the FFT subroutine in the PESSL manual

/afs/ovpit.indiana.edu/common/www/htdocs/gustav/SP-docs/pessl/pessl.00.ps

View that manual with ghostview and go to pages 623-632. Observe that the Fourier transform of the diffusion coefficient function is going to be complex in general (in the lecture I have demonstrated that the matrix F is real and symmetric, but that is not the same as Fourier transform of D(x,y) alone).

On the other hand, a Fourier transform of the initial condition T₀(x,y)is going to be real (can you explain why?). The sizing condition in this case is somewhat different than in a real-to-complex case. Compare with notes on page 626 of the PESSL manual.