fortran and columns 
Author Message
 fortran and columns

Hi, i'm trying to get the following program to work properly. I believe
it's not working right, because i don't know how fortran works with
respect to columns. i know the program has to start off in column 7, but
are there specific rules for do statements and if statement with regards
to columns or any other column rules other than the programming starting
in column 7? The program is from a book by the way.

Thanks

      DIMENSION V(0:20,0:20)
      DATA PIE/3.14159/
      DATA A,B/1.0,1.0/
      DATA V1,V2,V3,V4/0.0,10.0,20.0,-10.0/
      NX=4
      NY=4
      H=A/FLOAT(NX)
      DO 10 I=1,NX-1
      DO 10 J=1,NY-1
      V(I,J)=(V1+V2+V3+V4)/4.0
   10 CONTINUE
      DO 20 I=1,NX-1
      V(I,0)=V1
      V(I,NY)=V3
   20 CONTINUE
      DO 30 J=1,NY-1
      V(0,J)=V4
      V(NX,J)=V2
   30 CONTINUE
      V(0,0)=(V1+V4)/2.0
      V(NX,0)=(V1+V2)/2.0
      V(0,NY)=(V3+V4)/2.0
      V(NX,NY)=(V2+V3)/2.0
      T=COS(PIE/NX)+COS(PIE/NY)
      W=(8.0-SQRT(64.0-16.0*T*T))/(T*T)
      WRITE(6,40) W
   40 FORMAT(2X,'SOR FACTOR OMEGA=',F10.6)
      W4=4/4.0
      NCOUNT=0
   50 RMIN=0.0
      DO 70 I=1,NX-1
      X=H*FLOAT(I)
      DO 70 J=1,NY-1
      Y=H*FLOAT(J)
      G=-36.0*PIE*X*(Y-1.0)
      R=W4*(V(I+1,J)+V(I-1,J)+V(I,J+1)+V(I,J-1)-4.0*V(I,J)-G*H*H)
      RMIN=RMIN+ABS(R)
      V(I,J)=V(I,J)+R
   70 CONTINUE
      RMIN=RMIN/FLOAT(NX*NY)
         IF(RMIN.GE.0.0001) THEN
         NCOUNT=NCOUNT+1
         IF(NCOUNT.LT.100) THEN
         GO TO 50
         ELSE
         WRITE(6,80)
   80 FORMAT(2X,'SOLUTION DOES NOT CONVERGE IN 100 ITERATIONS')
         GO TO 100
         ENDIF
      ENDIF
      WRITE(6,90) NCOUNT
   90 FORMAT(2X,'SOLUTION CONVERGES IN',2X,I3,2X,'ITERATIONS',/)
  100 CONTINUE
      DO 120 I=1,NX-1
      DO 120 J=1,NY-1
      WRITE(6,110) I,J,V(I,J)
  110 FORMAT(2X,'I=',I3,2X,'J=',I3,2X,'V=',E12.6,/)
  120 CONTINUE
      DO 150 I=1,NX-1
      X=H*FLOAT(I)
      DO 150 J=1,NY-1
      Y=H*FLOAT(J)
      SUM=0.0
      DO 130 M=1,10
      FM=FLOAT(M)
      DO 130 N=1,10
      FN=FLOAT(N)
      FACTOR1=(FM*PIE/A)**2+(FN*PIE/B)**2
      FACTOR2=((-1.0)**(M+N))*144.0*A*B/(PIE*FM*FN)
      FACTOR3=1.0-(1.0-(-1.0)**N)/B
      FACTOR=FACTOR2*FACTOR3/FACTOR1
      SUM=SUM+FACTOR*SIN(FM*PIE*X/A)*SIN(FN*PIE*Y/B)
 130  CONTINUE
      VH=SUM
      C1=4.0*V1/PIE
      C2=4.0*V2/PIE
      C3=4.0*V3/PIE
      C4=4.0*V4/PIE
      SUM=0.0
      DO 140 K=1,10
      N=2*K-1
      AN=FLOAT(N)
      A1=SIN(AN*PIE*X/B)
      A2=SINH(AN*PIE*(A-Y)/B)
      A3=AN*SINH(AN*PIE*A/B)
      TERM1=C1*A1*A2/A3
      B1=SINH(AN*PIE*X/A)
      B2=SIN(AN*PIE*Y/A)
      B3=AN*SINH(AN*PIE*B/A)
      TERM2=C2*B1*B2/B3
      D1=SIN(AN*PIE*X/B)
      D2=SINH(AN*PIE*Y/B)
      D3=AN*SINH(AN*PIE*A/B)
      TERM3=C3*D1*D2/D3
      E1=SINH(AN*PIE*(B-X)/A)
      E2=SIN(AN*PIE*Y/A)
      E3=AN*SINH(AN*PIE*B/A)
      TERM4=C4*E1*E2/E3
      TERM=TERM1+TERM2+TERM3+TERM4
      SUM=SUM+TERM
  140 CONTINUE
      VI=SUM
      V(I,J)=VH+VI
  150 CONTINUE
      DO 160 I=1,NX-1
      DO 160 J=1,NY-1
      WRITE(6,110) I,J,V(I,J)
  160 CONTINUE
      STOP
      END    



Tue, 18 Jul 2006 05:24:19 GMT  
 fortran and columns

"jakeh"  wrote

Quote:
> Hi, i'm trying to get the following program to work properly. I believe
> it's not working right, because i don't know how fortran works with
> respect to columns. i know the program has to start off in column 7, but
> are there specific rules for do statements and if statement with regards
> to columns or any other column rules other than the programming starting
> in column 7?

For input and output it is the format statements in the code that determine
columns.  You will probably need a book on fortran. What exactly is your
problem?


Tue, 18 Jul 2006 06:20:53 GMT  
 fortran and columns

Quote:

> Hi, i'm trying to get the following program to work properly. I believe
> it's not working right, because i don't know how fortran works with
> respect to columns. i know the program has to start off in column 7, but
> are there specific rules for do statements and if statement with regards
> to columns or any other column rules other than the programming starting
> in column 7? The program is from a book by the way.

> Thanks

|-|-----|-|-------- .........--|
 1 2     6 7                 72

1-5 : statement numbers
6   : continuation mark
7-72: fortran code

You don't have a newer version of Fortran to avoid this old fixed form?
Fortran 90 and on is free form.

/Sakis



Tue, 18 Jul 2006 05:49:47 GMT  
 fortran and columns

Quote:

> i know the program has to start off in column 7, but
> are there specific rules for do statements and if statement with regards
> to columns or any other column rules other than the programming starting
> in column 7?

As of f90, there are 2 different sets of rules.  The old fixed-source
form ones and the new free-form source ones.  There are several
differences, but the code you showed can work ok in either one
(simplified by the fact that it has no comments or continuation lines,
which are some of the areas of difference).

The free-form source rules are the easiest - the only limitation
relating to columns is that you can't go past column 132 (which
this code is in no danger of doing).  So if you are using an f90
compiler (I don't recall that you said what compiler you were
using), easiest thing might be to tell the compiler that the code
is free source form.  For most compilers, you can do that just by
naming the file with a .f90 extension.

If you need to stick with the old fixed source form (possibly because
you are using an older compiler that only does f77), then the
relevant rules are.

1. The statement must start in column 7 or later.  You can indent after
   column 7 or not, completely to your taste - doesn't matter.  But
   don't start before the 7th column.

2. The statement numbers do *NOT* go after column 7.  They must go
   in columns 1-5.  Doesn't matter where in 1-5, as long as they are
   in that range.

   The statement numbers are those numbers before a few of the
   statements, like the 10 in

      10 CONTINUE

   If you misunderstood the column 7 business and put the statement
   numbers after columnn 7, then that won't work at all.

3. Don't let anything go past column 72.  That can be a confusing
   one that is hard to spot.  Doesn't look like any of your lines
   are dangerously close yet, but be careful of it anyway.

4. You might have noticed that column 6 got skipped over.  Yep.
   Don't let anything at all get in column 6 or you'll likely get
   very confusing error messages.  Column 6 is only for continuations,
   which you don't have any of.

5. The code doesn't have any comments, so you don't have to worry
   about that (though ahem, the lack of comments is not a good sign.)

Oh, and no, there are no rules about DO indentation or anything like
that.  There are things that people consider to be good style and to
help readability, but no actual rules.

--
Richard Maine                       |  Good judgment comes from experience;
email: my first.last at org.domain  |  experience comes from bad judgment.
org: nasa, domain: gov              |        -- Mark Twain



Tue, 18 Jul 2006 06:07:14 GMT  
 fortran and columns

Quote:


>> Hi, i'm trying to get the following program to work properly. I believe
>> it's not working right, because i don't know how fortran works with
>> respect to columns. i know the program has to start off in column 7, but
>> are there specific rules for do statements and if statement with regards
>> to columns or any other column rules other than the programming starting
>> in column 7? The program is from a book by the way.

>> Thanks

> |-|-----|-|-------- .........--|
>  1 2     6 7                 72

> 1-5 : statement numbers
> 6   : continuation mark
> 7-72: fortran code

> You don't have a newer version of Fortran to avoid this old fixed form?
> Fortran 90 and on is free form.

> /Sakis

By the way, I almost forgot it. Comments require a C (or a *) in column 1.

/Sakis



Tue, 18 Jul 2006 05:54:43 GMT  
 fortran and columns

Quote:

> Hi, i'm trying to get the following program to work properly. I believe
> it's not working right, because i don't know how fortran works with

The compiler liked your typing skills and even produced a result
claiming convergence in 77 tries. Attached is a lower case version
that might be easier on your eyes.

--
T.Silva
------------------------------------------------------
Applied Algorithms    http://sdynamix.com

      dimension v(0:20,0:20)
      data pie/3.14159/
      data a,b/1.0,1.0/
      data v1,v2,v3,v4/0.0,10.0,20.0,-10.0/
      nx=4
      ny=4
      h=a/float(nx)
      do 10 i=1,nx-1
      do 10 j=1,ny-1
      v(i,j)=(v1+v2+v3+v4)/4.0
   10 continue
      do 20 i=1,nx-1
      v(i,0)=v1
      v(i,ny)=v3
   20 continue
      do 30 j=1,ny-1
      v(0,j)=v4
      v(nx,j)=v2
   30 continue
      v(0,0)=(v1+v4)/2.0
      v(nx,0)=(v1+v2)/2.0
      v(0,ny)=(v3+v4)/2.0
      v(nx,ny)=(v2+v3)/2.0
      t=cos(pie/nx)+cos(pie/ny)
      w=(8.0-sqrt(64.0-16.0*t*t))/(t*t)
      write(6,40) w
   40 format(2x,'sor factor omega=',f10.6)
      w4=4/4.0
      ncount=0
   50 rmin=0.0
      do 70 i=1,nx-1
      x=h*float(i)
      do 70 j=1,ny-1
      y=h*float(j)
      g=-36.0*pie*x*(y-1.0)
      r=w4*(v(i+1,j)+v(i-1,j)+v(i,j+1)+v(i,j-1)-4.0*v(i,j)-g*h*h)
      rmin=rmin+abs(r)
      v(i,j)=v(i,j)+r
   70 continue
      rmin=rmin/float(nx*ny)
         if(rmin.ge.0.0001) then
         ncount=ncount+1
         if(ncount.lt.100) then
         go to 50
         else
         write(6,80)
   80 format(2x,'solution does not converge in 100 iterations')
         go to 100
         endif
      endif
      write(6,90) ncount
   90 format(2x,'solution converges in',2x,i3,2x,'iterations',/)
  100 continue
      do 120 i=1,nx-1
      do 120 j=1,ny-1
      write(6,110) i,j,v(i,j)
  110 format(2x,'i=',i3,2x,'j=',i3,2x,'v=',e12.6)
  120 continue
      do 150 i=1,nx-1
      x=h*float(i)
      do 150 j=1,ny-1
      y=h*float(j)
      sum=0.0
      do 130 m=1,10
      fm=float(m)
      do 130 n=1,10
      fn=float(n)
      factor1=(fm*pie/a)**2+(fn*pie/b)**2
      factor2=((-1.0)**(m+n))*144.0*a*b/(pie*fm*fn)
      factor3=1.0-(1.0-(-1.0)**n)/b
      factor=factor2*factor3/factor1
      sum=sum+factor*sin(fm*pie*x/a)*sin(fn*pie*y/b)
 130  continue
      vh=sum
      c1=4.0*v1/pie
      c2=4.0*v2/pie
      c3=4.0*v3/pie
      c4=4.0*v4/pie
      sum=0.0
      do 140 k=1,10
      n=2*k-1
      an=float(n)
      a1=sin(an*pie*x/b)
      a2=sinh(an*pie*(a-y)/b)
      a3=an*sinh(an*pie*a/b)
      term1=c1*a1*a2/a3
      b1=sinh(an*pie*x/a)
      b2=sin(an*pie*y/a)
      b3=an*sinh(an*pie*b/a)
      term2=c2*b1*b2/b3
      d1=sin(an*pie*x/b)
      d2=sinh(an*pie*y/b)
      d3=an*sinh(an*pie*a/b)
      term3=c3*d1*d2/d3
      e1=sinh(an*pie*(b-x)/a)
      e2=sin(an*pie*y/a)
      e3=an*sinh(an*pie*b/a)
      term4=c4*e1*e2/e3
      term=term1+term2+term3+term4
      sum=sum+term
  140 continue
      vi=sum
      v(i,j)=vh+vi
  150 continue
      do 160 i=1,nx-1
      do 160 j=1,ny-1
      write(6,110) i,j,v(i,j)
  160 continue
      stop
      end



Tue, 18 Jul 2006 10:05:30 GMT  
 fortran and columns
Jake,
Your code is in very, very old Fortran.
Since the 1990 standard, it has not been necessary for commands to start in
or beyond column 7.
Attached is a more modern version of your code - together with the output.

Cheers
--
Alan Miller
http://users.bigpond.net.au/amiller
Retired Statistician


Quote:
> Hi, i'm trying to get the following program to work properly. I believe
> it's not working right, because i don't know how fortran works with
> respect to columns. i know the program has to start off in column 7, but
> are there specific rules for do statements and if statement with regards
> to columns or any other column rules other than the programming starting
> in column 7? The program is from a book by the way.

> Thanks

>       DIMENSION V(0:20,0:20)
>       DATA PIE/3.14159/
>       DATA A,B/1.0,1.0/
>       DATA V1,V2,V3,V4/0.0,10.0,20.0,-10.0/
>       NX=4
>       NY=4
>       H=A/FLOAT(NX)
>       DO 10 I=1,NX-1
>       DO 10 J=1,NY-1
>       V(I,J)=(V1+V2+V3+V4)/4.0
>    10 CONTINUE
>       DO 20 I=1,NX-1
>       V(I,0)=V1
>       V(I,NY)=V3
>    20 CONTINUE
>       DO 30 J=1,NY-1
>       V(0,J)=V4
>       V(NX,J)=V2
>    30 CONTINUE
>       V(0,0)=(V1+V4)/2.0
>       V(NX,0)=(V1+V2)/2.0
>       V(0,NY)=(V3+V4)/2.0
>       V(NX,NY)=(V2+V3)/2.0
>       T=COS(PIE/NX)+COS(PIE/NY)
>       W=(8.0-SQRT(64.0-16.0*T*T))/(T*T)
>       WRITE(6,40) W
>    40 FORMAT(2X,'SOR FACTOR OMEGA=',F10.6)
>       W4=4/4.0
>       NCOUNT=0
>    50 RMIN=0.0
>       DO 70 I=1,NX-1
>       X=H*FLOAT(I)
>       DO 70 J=1,NY-1
>       Y=H*FLOAT(J)
>       G=-36.0*PIE*X*(Y-1.0)
>       R=W4*(V(I+1,J)+V(I-1,J)+V(I,J+1)+V(I,J-1)-4.0*V(I,J)-G*H*H)
>       RMIN=RMIN+ABS(R)
>       V(I,J)=V(I,J)+R
>    70 CONTINUE
>       RMIN=RMIN/FLOAT(NX*NY)
>          IF(RMIN.GE.0.0001) THEN
>          NCOUNT=NCOUNT+1
>          IF(NCOUNT.LT.100) THEN
>          GO TO 50
>          ELSE
>          WRITE(6,80)
>    80 FORMAT(2X,'SOLUTION DOES NOT CONVERGE IN 100 ITERATIONS')
>          GO TO 100
>          ENDIF
>       ENDIF
>       WRITE(6,90) NCOUNT
>    90 FORMAT(2X,'SOLUTION CONVERGES IN',2X,I3,2X,'ITERATIONS',/)
>   100 CONTINUE
>       DO 120 I=1,NX-1
>       DO 120 J=1,NY-1
>       WRITE(6,110) I,J,V(I,J)
>   110 FORMAT(2X,'I=',I3,2X,'J=',I3,2X,'V=',E12.6,/)
>   120 CONTINUE
>       DO 150 I=1,NX-1
>       X=H*FLOAT(I)
>       DO 150 J=1,NY-1
>       Y=H*FLOAT(J)
>       SUM=0.0
>       DO 130 M=1,10
>       FM=FLOAT(M)
>       DO 130 N=1,10
>       FN=FLOAT(N)
>       FACTOR1=(FM*PIE/A)**2+(FN*PIE/B)**2
>       FACTOR2=((-1.0)**(M+N))*144.0*A*B/(PIE*FM*FN)
>       FACTOR3=1.0-(1.0-(-1.0)**N)/B
>       FACTOR=FACTOR2*FACTOR3/FACTOR1
>       SUM=SUM+FACTOR*SIN(FM*PIE*X/A)*SIN(FN*PIE*Y/B)
>  130  CONTINUE
>       VH=SUM
>       C1=4.0*V1/PIE
>       C2=4.0*V2/PIE
>       C3=4.0*V3/PIE
>       C4=4.0*V4/PIE
>       SUM=0.0
>       DO 140 K=1,10
>       N=2*K-1
>       AN=FLOAT(N)
>       A1=SIN(AN*PIE*X/B)
>       A2=SINH(AN*PIE*(A-Y)/B)
>       A3=AN*SINH(AN*PIE*A/B)
>       TERM1=C1*A1*A2/A3
>       B1=SINH(AN*PIE*X/A)
>       B2=SIN(AN*PIE*Y/A)
>       B3=AN*SINH(AN*PIE*B/A)
>       TERM2=C2*B1*B2/B3
>       D1=SIN(AN*PIE*X/B)
>       D2=SINH(AN*PIE*Y/B)
>       D3=AN*SINH(AN*PIE*A/B)
>       TERM3=C3*D1*D2/D3
>       E1=SINH(AN*PIE*(B-X)/A)
>       E2=SIN(AN*PIE*Y/A)
>       E3=AN*SINH(AN*PIE*B/A)
>       TERM4=C4*E1*E2/E3
>       TERM=TERM1+TERM2+TERM3+TERM4
>       SUM=SUM+TERM
>   140 CONTINUE
>       VI=SUM
>       V(I,J)=VH+VI
>   150 CONTINUE
>       DO 160 I=1,NX-1
>       DO 160 J=1,NY-1
>       WRITE(6,110) I,J,V(I,J)
>   160 CONTINUE
>       STOP
>       END


Subject: fortran and columns
Date: Friday, 30 January 2004 8:24 AM

Hi, i'm trying to get the following program to work properly. I believe
it's not working right, because i don't know how fortran works with
respect to columns. i know the program has to start off in column 7, but
are there specific rules for do statements and if statement with regards
to columns or any other column rules other than the programming starting
in column 7? The program is from a book by the way.

Thanks

      DIMENSION V(0:20,0:20)
      DATA PIE/3.14159/
      DATA A,B/1.0,1.0/
      DATA V1,V2,V3,V4/0.0,10.0,20.0,-10.0/
      NX=4
      NY=4
      H=A/FLOAT(NX)
      DO 10 I=1,NX-1
      DO 10 J=1,NY-1
      V(I,J)=(V1+V2+V3+V4)/4.0
   10 CONTINUE
      DO 20 I=1,NX-1
      V(I,0)=V1
      V(I,NY)=V3
   20 CONTINUE
      DO 30 J=1,NY-1
      V(0,J)=V4
      V(NX,J)=V2
   30 CONTINUE
      V(0,0)=(V1+V4)/2.0
      V(NX,0)=(V1+V2)/2.0
      V(0,NY)=(V3+V4)/2.0
      V(NX,NY)=(V2+V3)/2.0
      T=COS(PIE/NX)+COS(PIE/NY)
      W=(8.0-SQRT(64.0-16.0*T*T))/(T*T)
      WRITE(6,40) W
   40 FORMAT(2X,'SOR FACTOR OMEGA=',F10.6)
      W4=4/4.0
      NCOUNT=0
   50 RMIN=0.0
      DO 70 I=1,NX-1
      X=H*FLOAT(I)
      DO 70 J=1,NY-1
      Y=H*FLOAT(J)
      G=-36.0*PIE*X*(Y-1.0)
      R=W4*(V(I+1,J)+V(I-1,J)+V(I,J+1)+V(I,J-1)-4.0*V(I,J)-G*H*H)
      RMIN=RMIN+ABS(R)
      V(I,J)=V(I,J)+R
   70 CONTINUE
      RMIN=RMIN/FLOAT(NX*NY)
         IF(RMIN.GE.0.0001) THEN
         NCOUNT=NCOUNT+1
         IF(NCOUNT.LT.100) THEN
         GO TO 50
         ELSE
         WRITE(6,80)
   80 FORMAT(2X,'SOLUTION DOES NOT CONVERGE IN 100 ITERATIONS')
         GO TO 100
         ENDIF
      ENDIF
      WRITE(6,90) NCOUNT
   90 FORMAT(2X,'SOLUTION CONVERGES IN',2X,I3,2X,'ITERATIONS',/)
  100 CONTINUE
      DO 120 I=1,NX-1
      DO 120 J=1,NY-1
      WRITE(6,110) I,J,V(I,J)
  110 FORMAT(2X,'I=',I3,2X,'J=',I3,2X,'V=',E12.6,/)
  120 CONTINUE
      DO 150 I=1,NX-1
      X=H*FLOAT(I)
      DO 150 J=1,NY-1
      Y=H*FLOAT(J)
      SUM=0.0
      DO 130 M=1,10
      FM=FLOAT(M)
      DO 130 N=1,10
      FN=FLOAT(N)
      FACTOR1=(FM*PIE/A)**2+(FN*PIE/B)**2
      FACTOR2=((-1.0)**(M+N))*144.0*A*B/(PIE*FM*FN)
      FACTOR3=1.0-(1.0-(-1.0)**N)/B
      FACTOR=FACTOR2*FACTOR3/FACTOR1
      SUM=SUM+FACTOR*SIN(FM*PIE*X/A)*SIN(FN*PIE*Y/B)
 130  CONTINUE
      VH=SUM
      C1=4.0*V1/PIE
      C2=4.0*V2/PIE
      C3=4.0*V3/PIE
      C4=4.0*V4/PIE
      SUM=0.0
      DO 140 K=1,10
      N=2*K-1
      AN=FLOAT(N)
      A1=SIN(AN*PIE*X/B)
      A2=SINH(AN*PIE*(A-Y)/B)
      A3=AN*SINH(AN*PIE*A/B)
      TERM1=C1*A1*A2/A3
      B1=SINH(AN*PIE*X/A)
      B2=SIN(AN*PIE*Y/A)
      B3=AN*SINH(AN*PIE*B/A)
      TERM2=C2*B1*B2/B3
      D1=SIN(AN*PIE*X/B)
      D2=SINH(AN*PIE*Y/B)
      D3=AN*SINH(AN*PIE*A/B)
      TERM3=C3*D1*D2/D3
      E1=SINH(AN*PIE*(B-X)/A)
      E2=SIN(AN*PIE*Y/A)
      E3=AN*SINH(AN*PIE*B/A)
      TERM4=C4*E1*E2/E3
      TERM=TERM1+TERM2+TERM3+TERM4
      SUM=SUM+TERM
  140 CONTINUE
      VI=SUM
      V(I,J)=VH+VI
  150 CONTINUE
      DO 160 I=1,NX-1
      DO 160 J=1,NY-1
      WRITE(6,110) I,J,V(I,J)
  160 CONTINUE
      STOP
      END

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

PROGRAM columns

! Code converted using TO_F90 by Alan Miller
! Date: 2004-01-30  Time: 09:26:00

IMPLICIT NONE
REAL     :: an, a1, a2, a3, b1, b2, b3, c1, c2, c3, c4, d1, d2, d3,   &
    e1, e2, e3, factor, factor1, factor2, factor3, fm, fn, g, h,   &
    r, rmin, sum, t, term, term1, term2, term3, term4, vh, vi, w, w4, x, y
INTEGER  :: i, j, k, ncount, m, n, nx, ny
REAL     :: v(0:20,0:20)
REAL, PARAMETER  :: pie = 3.141593, a = 1.0, b = 1.0
REAL, PARAMETER  :: v1 = 0.0, v2 = 10.0, v3 = 20.0, v4 = -10.0

nx=4
ny=4
h=a/nx
DO  i=1,nx-1
  DO  j=1,ny-1
    v(i,j)=(v1+v2+v3+v4)/4.0
  END DO
END DO
DO  i=1,nx-1
  v(i,0)=v1
  v(i,ny)=v3
END DO
DO  j=1,ny-1
  v(0,j)=v4
  v(nx,j)=v2
END DO
v(0,0)=(v1+v4)/2.0
v(nx,0)=(v1+v2)/2.0
v(0,ny)=(v3+v4)/2.0
v(nx,ny)=(v2+v3)/2.0
t=COS(pie/nx) + COS(pie/ny)
w=(8.0-SQRT(64.0-16.0*t*t))/(t*t)
WRITE(6,40) w
40 FORMAT('  SOR FACTOR OMEGA=', f10.6)
w4=4/4.0
ncount=0
50 rmin=0.0
DO  i=1,nx-1
  x=h*i
  DO  j=1,ny-1
    y=h*j
    g=-36.0*pie*x*(y-1.0)
    r=w4*(v(i+1,j) + v(i-1,j) + v(i,j+1) + v(i,j-1) - 4.0*v(i,j) - g*h*h)
    rmin=rmin+ABS(r)
    v(i,j)=v(i,j) + r
  END DO
END DO
rmin=rmin/(nx*ny)
IF(rmin >= 0.0001) THEN
  ncount=ncount+1
  IF(ncount < 100) THEN
    GO TO 50
  ELSE
    WRITE(6,80)
    80 FORMAT('  SOLUTION DOES NOT CONVERGE IN 100 ITERATIONS')
    GO TO 100
  END IF
END IF
WRITE(6,90) ncount
90 FORMAT('  SOLUTION CONVERGES IN  ', i3, '  ITERATIONS'/)

100 DO  i=1,nx-1
  DO  j=1,ny-1
    WRITE(6,110) i,j,v(i,j)
    110 FORMAT('  I=', i3, '  J=', i3, '  V=', e12.6/)
  END DO
END DO
DO  i=1,nx-1
  x=h*i
  DO  j=1,ny-1
    y=h*j
    sum=0.0
    DO  m=1,10
      fm=m
      DO  n=1,10
        fn=n
        factor1=(fm*pie/a)**2 + (fn*pie/b)**2
        factor2=((-1.0)**(m+n))*144.0*a*b/(pie*fm*fn)
        factor3=1.0-(1.0-(-1.0)**n)/b
        factor=factor2*factor3/factor1
        sum=sum + factor*SIN(fm*pie*x/a)*SIN(fn*pie*y/b)
      END DO
    END DO
    vh=sum
    c1=4.0*v1/pie
    c2=4.0*v2/pie
    c3=4.0*v3/pie
    c4=4.0*v4/pie
    sum=0.0
    DO  k=1,10
      n=2*k-1
      an=n
      a1=SIN(an*pie*x/b)
      a2=SINH(an*pie*(a-y)/b)
      a3=an*SINH(an*pie*a/b)
      term1=c1*a1*a2/a3
      b1=SINH(an*pie*x/a)
      b2=SIN(an*pie*y/a)
      b3=an*SINH(an*pie*b/a)
      term2=c2*b1*b2/b3
      d1=SIN(an*pie*x/b)
      d2=SINH(an*pie*y/b)
      d3=an*SINH(an*pie*a/b)
      term3=c3*d1*d2/d3
      e1=SINH(an*pie*(b-x)/a)
      e2=SIN(an*pie*y/a)
      e3=an*SINH(an*pie*b/a)
      term4=c4*e1*e2/e3
      term=term1+term2+term3+term4
      sum=sum+term
    END DO
    vi=sum
    v(i,j)=vh+vi
  END DO
END DO
DO  i=1,nx-1
  DO  j=1,ny-1
    WRITE(6,110) i,j,v(i,j)
  END DO
END DO
STOP
END PROGRAM columns

Output:

 SOR FACTOR OMEGA=  1.171573
 SOLUTION CONVERGES IN   77  ITERATIONS

 I=  1  J=  1  V=-.530195E+38

 I=  1  J=  2  V=-.128261E+39

 I=  1  J=  3  V=-.949419E+38

 I=  2  J=  1  V=-.125229E+39

 I=  2  J=  2  V=-.242635E+39

 I=  2  J=  3  V=************

 I=  3  J=  1  V=-.136810E+39

 I=  3  J=  2  V=************

 I=  3  J=  3  V=************

 I=
...

read more »



Tue, 18 Jul 2006 12:43:46 GMT  
 fortran and columns


Quote:
> Jake,
> Your code is in very, very old Fortran.

I'm not sure whether to call it Fortran.  Does it converge only with data
limited to single precision range?

      real*8 v(0:20,0:20)
      real *8 pie,sum,vh,rmin,x,h,y,t,vi,g,r
      real*8 factor1,factor2,factor3,factor
      real a,b
      parameter( pie=3.1415927,a=1,b=1)
      real v1,v2,v3,v4
      parameter( v1=0,v2=10,v3=20,v4=-10)
      real*8 c1,c2,c3,c4
      parameter(c1=4*v1/pie,c2=4*v2/pie,c3=4*v3/pie,c4=4*v4/pie)
      real*8 term1,term2,term3,term4,term
      parameter(nx= 4, ny= 4, h= a/nx)
      real*8 a1,a2,a3,b1,b2,b3,d1,d2,d3,e1,e2,e3
      do i= 1,nx-1
        do j= 1,ny-1
          v(i,j)= (v1+v2+v3+v4)/4.0
          enddo
          v(i,0)= v1
          v(i,ny)= v3
        enddo
      do j= 1,ny-1
          v(0,j)= v4
          v(nx,j)= v2
        enddo
      v(0,0)= (v1+v4)/2.0
      v(nx,0)= (v1+v2)/2.0
      v(0,ny)= (v3+v4)/2.0
      v(nx,ny)= (v2+v3)/2.0
      t= cos(pie/nx)+cos(pie/ny)
      w= (8.0-sqrt(64.0-16.0*t*t))/(t*t)
      write(6,10)w
      w4= 1
      do ncount= 0,99
          rmin= 0.0
          do i= 1,nx-1
              x= h*(i)
              do j= 1,ny-1
                  y= h*(j)
                  g= -36.0*pie*x*(y-1.0)
                  r= w4*(v(i+1,j)+v(i-1,j)+v(i,j+1)+v(i,j-1)-4.0*v(i,j)-
     &g*h*h)
                  rmin= rmin+abs(r)
                  v(i,j)= v(i,j)+r
                enddo
            enddo
          rmin= rmin/(nx*ny)
          if(rmin <  0.0001)then
              write(6,30)ncount
             GO TO 100
            endif
23009   enddo
        write(6,20)
  100 CONTINUE
      do i= 1,nx-1
        do j= 1,ny-1
          write(6,40)i,j,v(i,j)
          enddo
        enddo
      do i= 1,nx-1
          x= h*(i)
          do j= 1,ny-1
              y= h*(j)
              sum= 0.0
              m1mpn=1
              do m= 1,10
                  fm= (m)
                  m1n=-1
                  do n= 1,10
                      fn= (n)
                      factor1= (fm*pie/a)**2+(fn*pie/b)**2
                      factor2= m1mpn*144.0*a*b/(pie*fm*fn)
                      m1mpn=-m1mpn
                      factor3= 1.0-(1.0-m1n)/b
                      m1n=-m1n
                      factor= factor2*factor3/factor1
                      sum= sum+factor*sin(fm*pie/a*x)*sin(fn*pie/b*y)
                    enddo
                  m1mpn=-m1mpn
                enddo
              vh= sum
              sum= 0.0
              do k= 1,10
                  n= 2*k-1
                  an= (n)
                  a1= sin(an*pie/b*x)
                  a2= sinh(an*pie/b*(a-y))
                  a3= an*sinh(an*pie/b*a)
                  term1= c1*a1*a2/a3
                  b1= sinh(an*pie/a*x)
                  b2= sin(an*pie/a*y)
                  b3= an*sinh(an*pie/a*b)
                  term2= c2*b1*b2/b3
                  d1= sin(an*pie/b*x)
                  d2= sinh(an*pie/b*y)
                  d3= an*sinh(an*pie/b*a)
                  term3= c3*d1*d2/d3
                  e1= sinh(an*pie/a*(b-x))
                  e2= sin(an*pie/a*y)
                  e3= an*sinh(an*pie/a*b)
                  term4= c4*e1*e2/e3
                  term= term1+term2+term3+term4
                  sum= sum+term
                enddo
              vi= sum
              v(i,j)= vh+vi
            enddo
        enddo
      do i= 1,nx-1
        do j= 1,ny-1
          write(6,40)i,j,v(i,j)
          enddo
        enddo

10      format(2x,'SOR FACTOR OMEGA=',f10.6)
20      format(2x,'SOLUTION DOES NOT CONVERGE IN 100 ITERATIONS')
30      format(2x,'SOLUTION CONVERGES IN',2x,i3,2x,'ITERATIONS',/)
40      format(2x,'I=',i3,2x,'J=',i3,2x,'V=',e12.6,/)
      end



Tue, 18 Jul 2006 15:51:18 GMT  
 fortran and columns
Others have commented on the significance of columns and indentation.
The program as given produces a floating point overflow in iteration 76
(ish) (Salford ftn95, open Watcom). RMIN has become very large. Maybe
you should go back to the underlying theory and check the algorithm. Or
check that it was transcribed correctly from the book.
Quote:
>      W=(8.0-SQRT(64.0-16.0*T*T))/(T*T)
>      WRITE(6,40) W
>   40 FORMAT(2X,'SOR FACTOR OMEGA=',F10.6)
>      W4=4/4.0

       ^^^^^^^^
This is another way of saying w4 = 1.0. Was w4 = w/4.0 intended?
--
to reply, interpret the following:
John Mansell:  john at wcompsys dot co dot uk


Tue, 18 Jul 2006 17:59:03 GMT  
 fortran and columns

Quote:

> Hi, i'm trying to get the following program to work properly.  I
> believe it's not working right, because i don't know how fortran
> works with respect to columns ...
...
> The program is from a book by the way.

The columns look okay, but ...

Quote:
>       W4=4/4.0

That statement looks suspicious.  I smell a typo.  Maybe others
lurking in there as well.

A problem with transcribing mysterious code from books is that
sometimes even the books have typos.

Maybe that should be W4=W/4.0, as John Mansell suggested, or maybe
it should be W4=A/4.0 (A looks kinda similar to 4).

Quote:
>       FACTOR1=(FM*PIE/A)**2+(FN*PIE/B)**2
>       FACTOR2=((-1.0)**(M+N))*144.0*A*B/(PIE*FM*FN)
>       FACTOR3=1.0-(1.0-(-1.0)**N)/B

Another potential "gotcha" is that old FORTRAN limited symbol
names to 6 characters, and FACTOR1, FACTOR2, and FACTOR3 are
7 characters each.  Many old compilers allowed longer names as
an extension.  I seem to remember one that allowed 16 characters,
but ignored everything after the first 6, essentially truncating
all of these symbols to FACTOR.  Modern compilers allow much
longer names without truncating, but I'm surprised to see long
names in FORTRAN of this vintage.  It violates the '77 standard.
If you compile in strict ANSI FORTRAN 77 mode, I'm not sure what
will happen.  You should at least get a warning from the compiler.
You could try shortening the names:

      FACTR1=(FM*PIE/A)**2+(FN*PIE/B)**2
      FACTR2=((-1.0)**(M+N))*144.0*A*B/(PIE*FM*FN)
      FACTR3=1.0-(1.0-(-1.0)**N)/B
      FACTOR=FACTR2*FACTR3/FACTR1

--

Ted Hall



Tue, 18 Jul 2006 21:55:36 GMT  
 fortran and columns


Quote:
> Others have commented on the significance of columns and indentation.
> The program as given produces a floating point overflow in iteration 76
> (ish) (Salford ftn95, open Watcom). RMIN has become very large. Maybe
> you should go back to the underlying theory and check the algorithm. Or
> check that it was transcribed correctly from the book.
> >      W=(8.0-SQRT(64.0-16.0*T*T))/(T*T)
> >      WRITE(6,40) W
> >   40 FORMAT(2X,'SOR FACTOR OMEGA=',F10.6)
> >      W4=4/4.0
>        ^^^^^^^^
> This is another way of saying w4 = 1.0. Was w4 = w/4.0 intended?

Apparently, a good catch.  Output looks reasonably stable now:

  SOR FACTOR OMEGA=  1.171573
  SOLUTION CONVERGES IN    8  ITERATIONS

  I=  1  J=  1  V=-.324739E+01

  I=  1  J=  2  V=-.170334E+01

  I=  1  J=  3  V=0.430575E+01

  I=  2  J=  1  V=0.393097E-01

  I=  2  J=  2  V=0.301193E+01

  I=  2  J=  3  V=0.936813E+01

  I=  3  J=  1  V=0.304350E+01

  I=  3  J=  2  V=0.611078E+01

  I=  3  J=  3  V=0.110384E+02

  I=  1  J=  1  V=-.342926E+01

  I=  1  J=  2  V=-.202854E+01

  I=  1  J=  3  V=0.427733E+01

  I=  2  J=  1  V=-.118218E+00

  I=  2  J=  2  V=0.291348E+01

  I=  2  J=  3  V=0.959263E+01

  I=  3  J=  1  V=0.290189E+01

  I=  3  J=  2  V=0.606522E+01

  I=  3  J=  3  V=0.111330E+02



Tue, 18 Jul 2006 22:13:54 GMT  
 fortran and columns

Quote:
> I seem to remember one that allowed 16 characters,
> but ignored everything after the first 6,

Yes! The people who implemented that abomination without telling
their customers should have been drawn and quartered publicly!

Nothing like having a several-million-line program, put together
from contributions of more than hundred people from several dozen
organisations all over the world, crash at some strange place because
a COMMON variable nobody expects to touch is corrupted, only to find
out that an update to a subroutine used a seven-letter variable that
"aliased" at six letters with that COMMON variable, and the author of
routine including the definition of the COMMON "just to be on the safe
side". Thank God I didn't need to debug that one.

        Jan



Tue, 18 Jul 2006 22:31:34 GMT  
 fortran and columns

(snip)

Quote:
>>This is another way of saying w4 = 1.0. Was w4 = w/4.0 intended?

> Apparently, a good catch.  Output looks reasonably stable now:

>   SOR FACTOR OMEGA=  1.171573
>   SOLUTION CONVERGES IN    8  ITERATIONS

Numerical Recipes has a good explanation of Successive OverRelaxation,
or SOR, including some programs.   You might look at those.

-- glen



Wed, 19 Jul 2006 02:05:48 GMT  
 fortran and columns


Quote:

> (snip)

> >>This is another way of saying w4 = 1.0. Was w4 = w/4.0 intended?

> > Apparently, a good catch.  Output looks reasonably stable now:

> >   SOR FACTOR OMEGA=  1.171573
> >   SOLUTION CONVERGES IN    8  ITERATIONS

> Numerical Recipes has a good explanation of Successive OverRelaxation,
> or SOR, including some programs.   You might look at those.

> -- glen

Understanding SOR is only a small part of understanding this code snippet
and correcting the bugs in it.  Most of us probably don't care to know what
this code is intended to demonstrate. If the author cared to help in that
respect, it certainly doesn't show.  It might be intended (if expanded to a
finer mesh) to compare SOR solution with an analytical separation of
variables solution, judging by all those sinh() and sin() functions.


Wed, 19 Jul 2006 05:33:03 GMT  
 fortran and columns

Quote:
>> I seem to remember one that allowed 16 characters,
>> but ignored everything after the first 6,

> Yes! The people who implemented that abomination without telling
> their customers should have been drawn and quartered publicly!

> .... Thank God I didn't need to debug that one.

I did.  :-(

Well, not the several million lines part, but I've debugged code
hit by such abominations.

The bodies are well hidden.  :-)

--
Richard Maine
email: my last name at domain
domain: sumertriangle dot net



Wed, 19 Jul 2006 09:35:12 GMT  
 
 [ 16 post ]  Go to page: [1] [2]

 Relevant Pages 

1. New to FORTRAN - Reformatting One Column to Multiple Columns

2. how do i convert a multi column into a single column in the output in fortran

3. sum all columns based on column 1

4. Help: using hidden columns with columns resizable in a list box

5. Column selection in multi-column report.

6. CW5: Column selection in multi-column report.

7. Plotting 1st column vs 2nd column in Excel from LabVIEW

8. Row/Column in Fortran

9. reformatting Fortran code with characters beyond column 72

10. About the column-dominant storage of Fortran

11. Tab delimited columns in FORTRAN

12. why column arrays in Fortran (was opinions on computer languages)

 

 
Powered by phpBB® Forum Software