Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
28456 Discussions

Memory

JohnNichols
Valued Contributor III
‎03-01-2016 11:01 AM
194 Views

Back in the old days - Conte and deBoore was a nice textbook on Numerical Analysis. I remember using this textbook in my Mechanical Engineering course on Numerical Analysis - of course one could only use Fortran then.

Skip forward to 1987 and I needed a quick inverter to solve a minor problem, so I used the old Conte and de Boore method. 

I still have the code and the examples. 

The structural analysis package I lifted from 1973 is giving me one or two minor and interesting answers that I wanted to check - nothing major just annoyances as I work out the program.  Combined with FEAST it is a treat to use compared to a commerical package, one can get at all the numbers quickly and can play with them.

 

So I thought I would run out the old program,  found the floppy drive, load it into VS 2013 - compiles first time - run out a 30 by 30 test case from the 73 code and get the C and deB to invert the matrix, to see if the results are the same. 

Answer was a big NO.

Load up the C and deBoore test case, one can still find a PDF of the textbook online -- must be illegal - but what the heck. Still getting the wrong answer on reading in the data.  After a bit of exploration with Watch on VS 2013, I learnt something I had completely forgotten, one can pass a vector through a subroutine call and then declare it as an array in the subroutine. Here was the sticking point - it is in column and not row order.   Fixing the input routines solved a simple problem and now I can check the 73 code.  Of course one should load up the MKL inverter, but what the heck life is short play hard.  I need to put the MKL inverter into the 73 code, that is next.

One of the problems with eigenvectors is visuallying the result.  I have DISLIN - but I was darned if I could get it to run at the weekend - just ornery I guess and a 9 year old yellling in my ear, I had a simple bitmap drawer from Japan, so I thought - what the heck it is only simple line sketches.  I thought I need to draw lines, circles and some text.  I was fun -- text looks like MS DOS prompts after heated to 1200 C, but who cares.  It is possible to create bitmaps from a standard FORTRAN program -- of course it is not pretty - but then again neither am I.

Here is the code - very old Fortran - but it still works. and is often said - if it is not broken do not fix it

C *********************************************************************
      PROGRAM MATRIX
C *********************************************************************
C      This is a matrix inversion progamme taken from Conte and deBoor
C      Elementary Numerical Analysis.
C      It has been amended to run on an IBM PC under Microsoft Fortran
C      by John Nichols.
C *********************************************************************
      implicit character*60(z)

      DIMENSION Aorg(10000)
      Dimension A(10000),ainv(10000),b(100),ipivot(100),c(100),d(100)
      write(*,1)
1     format(1x,'This is an IBM MICROSOFT FORTRAN programme',/
     11x,'from CONTE and de BOOR,amended by John Nichols',/
     21x,'in July 1987.',/
     31x,' Ver. 1.0 --Ph (049 25464).'//)
C
C *********************************************************************
C
C      Read input and output devices for files.
C
C *********************************************************************
C

      write(*,20)
20    format(1x,'File name for input data :   '\)

      read(*,30)zFILE1
30    format(a12)

      write(*,40)
40    format(1x,'File name for output data :   '\)

      read(*,50)zFILE2
50    format(a12)

C
C *********************************************************************
C
C      Open files for input ouput
C
C *********************************************************************
C

      write(*,60)
60    format(1x,'Opening the input file to access data!'/)

      open(1,file=zFILE1,status='OLD')

      write(*,70)
70    format(1x,'Opening the output file for data!'/)

      open(2,file=zFILE2,status='UNKNOWN')

C
C *********************************************************************
C
C      Read input data from file 1
C
C *********************************************************************
C

100   continue
      c(1:100)=0.0

      read(1,80)N
80    format(I3)

      write(*,90)N
90    format(1x,'The matrix size is ',i2/)

      read(1,*)(d(j),j=1,n)
      read(1,*)N1

      nsq=N*N

      read(1,*)(a(j),j=1,nsq)
      do 102 i = 1,nsq
          aorg(i) = a(i)
102   end do

      write(*,101)(a(j),j=1,nsq)
  101 format(3(2x,F10.6))
      write(*,101)(d(j),j=1,n)
C
C *********************************************************************
C
C      Call factor subroutine
C
C *********************************************************************
C

      write(*,150)
150   format(1x,'At the factor subroutine!',/)

      call factor(a,a,ipivot,b,n,iflag)
      go to (120,110),iflag

110   continue

C
C *********************************************************************
C
C      The matrix is singular
C
C *********************************************************************
C

      write(*,130)
130   format(1x,'This matrix is singular have another shot!',/)
      go to 210

120   continue

C
C *********************************************************************
C
C      Call subst subroutine
C
C *********************************************************************
C

      do 140 I=1,N
          b(i)=0.0
140   continue

      ibeg=1

      do  160 j=1,N
          b(j)=1.0

          call subst(a,b,ainv(ibeg),ipivot,n)

          b(j)=0.0
          ibeg=ibeg+N
160   continue

C
C *********************************************************************
C
C      Inverted matrix results
C
C *********************************************************************
C

      write(2,170)
      write(*,170)
170   format(1x,'The computed inverse is ',//)

      write(2,101)(ainv(j),j=1,nsq)
      write(*,101)(ainv(j),j=1,nsq)
180   format((3F10.5))
210   continue

      call MatrixResult(aorg, ainv,d,c,n)
C
C *********************************************************************
C
C      End of programme
C
C *********************************************************************
C

      write(*,200)
200   format(1x,'The end of the programme',/)


      end

C
C *********************************************************************
C
      Subroutine MatrixResult(a,ainv,d,c, n)
C 
C *********************************************************************       
C
      dimension a(n,n),d(n),c(n),ainv(n,n), f(n,n)
      integer n,i,j

      do 100 i = 1,n
          do 200 j =1,n

              c(i) = a(i,j)*d(i)+c(i)
200       Continue
100   Continue
      f(1:n,1:n)=0.0

      do 400 i=1, n
          do 300 j=1, n
              do 250 k=1, n
                  f(i,j)=f(i,j)+a(j,k)*ainv(k,i)
250   enddo
300   end do
400   end do
      do 500 i = 1,n
          write(*,180)(F(i,j),j=1,n)
500   end do
180   format((30F4.1))
      write(*,70)
70    format(1x,'Opening the C  data!'/)
      write(*,*)(c(i),i=1,n)

      return

      end

C
C
C *********************************************************************
C
      Subroutine Factor(a,w,ipivot,d,n,iflag)
C 
C *********************************************************************       
C
      dimension a(n,n),w(n,n),ipivot(n),d(n)

C
C *********************************************************************
C
C      initialize vectors
C
C *********************************************************************
C

      iflag=1

      do 10 i=1,n
          ipivot(i)=i
          rowmax=0.0
          do 20 j=1,n
              w(i,j)=a(i,j)
              rowmax=amax1(rowmax,abs(w(i,j)))
20        continue
          if(rowmax .eq. 0.0) go to 999
          d(i)=rowmax
10    continue

C
C *********************************************************************
C
C      Gaussian elimination with scaled partial pivotting
C
C *********************************************************************
C

      nm1=n-1
      if(nm1 .eq. 0) return
      do 30 k=1,nm1
          j=k
          kp1=k+1
          ip=ipivot(k)
          colmax=abs(w(ip,k))/d(ip)
          do 40 i=kp1,n
              ip=ipivot(i)
              awikov=abs(w(ip,k))/d(ip)
              if (awikov .le. colmax ) go to 40
              colmax=awikov
              j=i
40        continue
          if (colmax .eq. 0.0) go to 999
          ipk=ipivot(j)
          ipivot(j)=ipivot(k)
          ipivot(k)=ipk
          do 30 i=kp1,n
              ip=ipivot(i)
              w(ip,k)=w(ip,k)/w(ipk,k)
              ratio=-w(ip,k)
              do 30 j=kp1,n
                  w(ip,j)=ratio*w(ipk,j)+w(ip,j)
30    continue
      if (w(ip,n) .eq. 0.0) go to 999
      return
999   iflag=2
      return
      end

C
C *********************************************************************
C
      subroutine subst (w,b,x,ipivot,n)
C
C *********************************************************************
C
      dimension w(n,n),b(n),x(n),ipivot(n)
      if (n .gt. 1) go to 10
      x(1)=b(1)/w(1,1)
      return
10    continue
      ip=ipivot(1)
      x(1)=b(ip)
      do 15 k=2,n
          ip=ipivot(k)
          km1=k-1
          sum=0.0
          do 14 j=1,km1
              sum=w(ip,j)*x(j)+sum
14        continue
          x(k)=b(ip)-sum
15    continue
      x(n)=x(n)/w(ip,n)
      k=n
      do 20 np1mk=2,n
          kp1=k
          k=k-1
          ip=ipivot(k)
          sum=0.0
          do 19 j=kp1,n
              sum=w(ip,j)*x(j)+sum
19        continue
          x(k)=(x(k)-sum)/w(ip,k)
20    continue
      return
      end

 

 

 

0 Kudos

Link Copied

0 Replies
Reply

Community support is provided during standard business hours (Monday to Friday 7AM - 5PM PST). Other contact methods are available here.

Intel does not verify all solutions, including but not limited to any file transfers that may appear in this community. Accordingly, Intel disclaims all express and implied warranties, including without limitation, the implied warranties of merchantability, fitness for a particular purpose, and non-infringement, as well as any warranty arising from course of performance, course of dealing, or usage in trade.

For more complete information about compiler optimizations, see our Optimization Notice.