Matrix sparsitiy pattern in eps format

Hi John,

thanks for your excellent comments and ideas!

I have taken a look at our code and your comment below applies. In many cases the code and the documentation is the other way arround. This is really very confusing and something ons should take care of. The stiffness matrix for our air gap element is such an example. Since this is a dense matrix your suggestion to make use of array syntax on the inner loop might improve the performance.

The other main issue to consider is documentation and readibility. If the algorithm on which the calculation is based, is row based, it might make sense to write the code so that it looks like the algorithm. ( or change the algorithm so that it looks column based). Documenting and using more readable names for variables reduces bugs and makes the code easier to support. Coming back to code that is confusing wastes a lot of time.

By the way the Zienkiewicz code for the fourth edition,pcfeap, must be the upgrade to the miniFEM in the third edition. This is also a good example of improvements due to new language options and possibilities. Unfortunately one needs to change the code a little bit to be able to use with FTN95.

Quoted from JohnCampbell

You are correct to identify that column based coding is more suitable for Fortran.

Using arrays was the first thing I got totally wrong with Fortran. After writing about 300 lines of code for my Tetris game, I then realized my game blocks looked a little bit funny (my data was in row-major order). 😉

I think array initialization is quite complex on Fortran.

In PL/I, I can just write: dcl A (3,4) float init (1,2,3,4,5,6,7,8,9,10,11,12);

and let compiler do the math.

Does Fortran have something similar to PL/I's iSUB defining?

dcl A (m,n) float;
dcl B (n,m) float defined (A(2sub, 1sub));

Above, matrix B would be defined to be transpose of the matrix A.

Hi John,

how is your visualization with %gr coming on? Unfortunately I am not that fit in Clearwin. Since I use eps in my documents I prefer to do all figures in pure eps rather than converting it. However, I think it should be easy to do it in Clearwin.

Coming back to the numbering: For the same example I only changed the node numbering manually. In this small example one can clearly see that the solution is independent of the node numbering. The profile for the second numbering scheme is larger than the first one. For problems with many nodes this is, as you know, critical for a direct solver - compact storage scheme by Alan Jennings.

[URL=http://imageshack.us/photo/my-images/696/spy1u2.png/][/URL]

      subroutine spy(a,m,n)
      implicit none
C
      include 'split_file.ins'
C
      integer,dimension(m,n) :: a
      integer,dimension(m)   :: numbering

      real,parameter :: cell_size=10.0
      real,parameter :: font_size=0.8*cell_size
      real,parameter :: r= 0.35*font_size
      integer,parameter :: io=33,tex=35

      real :: x_lo,x_hi,y_lo,y_hi,xdiag
      real :: scale,x,y
      integer :: m,n,i,j
    
      x_lo = 0;
      y_lo = 0;
      x_hi = m*cell_size;
      y_hi = m*cell_size;
C
C Open the eps file and define the bounding box 
C
      FileToWrite = trim(path)//trim(name)//'-spy.eps'
      open(io,file=FileToWrite);
      write(io,'(A)') '%!PS-Adobe-1.0'
      if (cell_size*real(n) > 595.0) then
        scale = 595.0/(cell_size*real(n));
        write(io,'(A,4I5)') '%%BoundingBox:',0,0,595,595
        write(io,'(2F13.6,1X,A)') 0.0,0.0,'translate'
        write(io,'(2F13.6,1X,A)') scale,scale,'scale'
      else
        write(io,'(A,4I5)') '%%BoundingBox:',
     *  int(x_lo+0.5),int(y_lo+0.5),int(x_hi+0.5),int(y_hi+0.5)
      endif
      write(io,'(2F13.6,2x,A)') x_lo,y_lo,'newpath moveto'
      write(io,'(2F13.6,2x,A)') x_hi,y_lo,'lineto'
      write(io,'(2F13.6,2x,A)') x_hi,y_hi,'lineto'
      write(io,'(2F13.6,2x,A)') x_lo,y_hi,'lineto'
      write(io,'(2F13.6,2x,A)') x_lo,y_lo,'lineto'
      write(io,'(A)') 'stroke'
C
C Plot the nonzero matrix entries
C
      do i=1,m
C First the diagonal element is plotted.
        write(io,'(A)') '0 setgray'
        x = x_lo + (real(i) - 0.5)*cell_size;
        y = y_hi - (real(i) - 0.5)*cell_size;
        write(io,'(3F13.6,2x,A)') x,y,r,'0 360 arc closepath'
        write(io,'(A)') 'gsave 1 0 0 setrgbcolor fill grestore'
        write(io,'(A)') 'stroke'
        xdiag = x
        do j=1,n
          if (a(i,j).NE. 0) then
            x = x_lo + (real(a(i,j)) - 0.5)*cell_size;
            y = y_hi - (real(i)      - 0.5)*cell_size;
            write(io,'(3F13.6,2I5,2x,A)') x,y,r,0,360,'arc closepath'
            if (x.LT.xdiag) then
              write(io,'(A,1X,3I3,1X,A)') 'gsave',0,0,1,
     *          ' setrgbcolor fill grestore'
            else
              write(io,'(A,1X,3I3,1X,A)') 'gsave',0,1,0,
     *          ' setrgbcolor fill grestore'
            endif
            write(io,'(A)') 'stroke'
          endif
        enddo
      enddo
      write(io,'(A)') 'showpage'
      write(io,'(A)') '%%EOF'
      close(io);    
      return
      end subroutine spy

Jalih, You said:

I think array initialization is quite complex on Fortran.

Although I am not familiar with PL/I, I am not sure why you find initialisation of a fortran array to be complex? There are a number of ways to initialise a 2D array. You could consider using: DATA statement, or array syntax, or DO loops, or FORALL structure, or RESHAPE intrinsic syntax. The following example is taken from Lahey documentation, which shows how a single vector can be re-shaped into a 2D array.

Example
x = reshape((/1,2,3,4/), (/3,2/), pad=(/0/))
! x is assigned |1 4|
!               |2 0|
!               |3 0|

Depending on either a column or row based storage, I don't see how this contributes to the complexity. Just use what is required.

If you need to think in columns, rather than rows, it is not that hard to copy a vector to a 2D array, as in k = 0 do i = 1,n do j = 1,m k = k+1 if (k > size(vector) then array(i,j) = default else array(i,j) = vector(k) end if end do end do I hope this gives you some options to investigate. John