Prev: call c lib from fortran in ms visual studio 2005
Next: What exactly does the SAVE statement do ?
From: Gordon Sande on 23 May 2010 18:10
On 2010-05-23 18:28:19 -0300, "e p chandler" <epc8(a)juno.com> said:

> "eric948470" wrote:
>> (Richard Maine) wrote:
>> You should normally avoid unit numbers in the single digits.
>> They are often used for special purposes.
>
> On additional suggestion. Your original program contained
>
> a(i,j)=((mu(j)**2.0D0)/(delta**2.0D0))
> d(i,j)=1.0D0+(2.0D0*((mu(j)**2.0D0)/(delta**2.0D0)))
> c(i,j)=(mu(j)**2.0D0)/(delta**2.0D0)
>
> While one change probably will not save you much run time, it's far
> better when raising a real number to an integer power to use an integer
> as an exponent.
>
> so you would have **2 instead of **2.0D0.

As well as working to square negative numbers.

> There are a number of other suggestions I would make in terms of style,
> etc. but I'll save those in case the OP returns.
>
> --- e


From: jfh on 23 May 2010 18:41
On May 23, 6:14 pm, nos...(a)see.signature (Richard Maine) wrote:
> eric948470 <eric948...(a)gmail.com> wrote:
>
> ...
>
> > But the above are displayed on screen for only 6 iterations. Then the
> > output to the screen stops. But the calculation seems to go on until
> > 'f' converges to the required tolerance (I know because there are a
> > total of 206 txt files written, not 6). It does what I want it to do.
> > But I don't what's happening to the screen output. Can anybody help?
>
> > And btw, there's also an unwanted file named 'fort.6' that is made by
> > the program.
> ...
>
> > PS: I checked the contents of 'fort.6' just now, and the rest of the
> > output I wanted to be written to screen has been written to that file.
>
> You should normally avoid unit numbers in the single digits. They are
> often used for special purposes. In particular, unit 6 is very commonly
> used for the screen output (aka unit *). By opening unit 6 with a file
> name, you have messed with the screen output. The exact details of what
> will result from opening and closing unit 6 are likely to vary from
> compiler to compiler. What happened to your code sounds plausible. After
> you closed unit 6, subsequent output to unit * reopened it with the
> fort.6 file name (but no longer to the screen). But you don't really
> care about those details. Just know to avoid such unit numbers.
>
> You really don't need multiple unit numbers at all. Although you are
> writing multiple files, you are doing so one at a time, opening one
> file, writing to it, and then closing it before opening the next one.
> After you have closed the file, it is perfectly fine and normal to reuse
> the unit number for the next file. I would recommend you use just a
> single unit number, independent of the iteration number..... and make
> that unit number something with 2 digits.

One may even need to be fussy about which 2-digit number to use:-(.
Some years ago we had a graphics package that gave trouble with WRITE
statements in a program that also drew a diagram.
INQUIRE(UNIT=N,OPENED=OP) in a suitable loop showed that the graphics
subroutines and the operating system were already using these unit
numbers: 0 5 6 19 21 23 44. And the Graphics Device Guide suggested
that 31 32 33 51 52 might also be in use.

-- John Harper

From: robin on 23 May 2010 23:04
"e p chandler" <epc8(a)juno.com> wrote in message news:htc6ne$d26$1(a)news.eternal-september.org...
| "eric948470" wrote:
| On additional suggestion. Your original program contained
|
| a(i,j)=((mu(j)**2.0D0)/(delta**2.0D0))
| d(i,j)=1.0D0+(2.0D0*((mu(j)**2.0D0)/(delta**2.0D0)))
| c(i,j)=(mu(j)**2.0D0)/(delta**2.0D0)
|
| While one change probably will not save you much run time, it's far better
| when raising a real number to an integer power to use an integer as an
| exponent.
|
| so you would have **2 instead of **2.0D0.

Even better is to do one squaring instead of two, e.g.,

c(i,j)=(mu(j)/delta)**2


From: e p chandler on 23 May 2010 23:52

"robin" <robin51(a)dodo.com.au> wrote in message
news:4bf9ecbe$0$11951$c30e37c6(a)exi-reader.telstra.net...
> "e p chandler" <epc8(a)juno.com> wrote in message
> news:htc6ne$d26$1(a)news.eternal-september.org...
> | "eric948470" wrote:
> | On additional suggestion. Your original program contained
> |
> | a(i,j)=((mu(j)**2.0D0)/(delta**2.0D0))
> | d(i,j)=1.0D0+(2.0D0*((mu(j)**2.0D0)/(delta**2.0D0)))
> | c(i,j)=(mu(j)**2.0D0)/(delta**2.0D0)
> |
> | While one change probably will not save you much run time, it's far
> better
> | when raising a real number to an integer power to use an integer as an
> | exponent.
> |
> | so you would have **2 instead of **2.0D0.
>
> Even better is to do one squaring instead of two, e.g.,
>
> c(i,j)=(mu(j)/delta)**2
>

Yes, along with other possible improvements. The OP apparently is not
familiar with operator precedence in Fortran.

Of course it would make things more pleasant to have free form source, F95+
do/end do, replacing the goto with do..end do plus if () exit, whole array
operations and so on.













From: e p chandler on 28 May 2010 16:47

"eric948470" wrote

[snip section related to conflict of unit numbers with pre-connected units]

> PROGRAM hmwrk6
[snip]
> do 2 j=1,10
> a(1,j)=0.0D0
> d(1,j)=1.0D0+(mu(j)/delta)+(delta/(2.0D0*mu(j)))
> c(1,j)=mu(j)/delta
> beta(2,j)=c(1,j)/d(1,j)
> do 3 i=2,NSTEPS
> a(i,j)=((mu(j)**2.0D0)/(delta**2.0D0))
> d(i,j)=1.0D0+(2.0D0*((mu(j)**2.0D0)/(delta**2.0D0)))
> c(i,j)=(mu(j)**2.0D0)/(delta**2.0D0)
> beta(i+1,j)=c(i,j)/(d(i,j)-(beta(i,j)*a(i,j)))
> 3 continue
> a(NSTEPS+1,j)=(mu(j)/delta)-0.5D0
> d(NSTEPS+1,j)=0.5D0+(mu(j)/delta)
> c(NSTEPS+1,j)=0.0D0
> beta(NSTEPS+2,j)=c(NSTEPS+1,j)/(d(NSTEPS+1,j)-(beta(NSTEPS
> +1,j)
> $*a(NSTEPS+1,j)))
> 2 continue
>
> adfmax=1.0D0
> n=0
> 4 if(adfmax.ge.1.0D-6)then
> n=n+1
> write(*,*)'iteration',n
> do 5 j=1,10
> e(1,j)=(b(1)*delta)/(2.0D0*mu(j))
> alpha(2,j)=e(1,j)/d(1,j)
> do 6 i=2,NSTEPS
> e(i,j)=b(i)
> alpha(i+1,j)=(e(i,j)+(alpha(i,j)*a(i,j)))/(d(i,j)-
> (beta(i
> $,j)*a(i,j)))
> 6 continue
> e(NSTEPS+1,j)=(b(NSTEPS+1)*(0.5D0+(mu(j)/delta)))+
> (b(NSTEPS)
> $*(0.5D0-(mu(j)/delta)))
> alpha(NSTEPS+2,j)=(e(NSTEPS+1,j)+(alpha(NSTEPS
> +1,j)*a(NSTEPS
> $+1,j)))/(d(NSTEPS+1,j)-(beta(NSTEPS+1,j)*a(NSTEPS+1,j)))
>
> u(NSTEPS+1,j)=alpha(NSTEPS+2,j)
> do 7 i=NSTEPS,1,-1
> u(i,j)=alpha(i+1,j)+(beta(i+1,j)*u(i+1,j))
> 7 continue
>
> uprime(1,j)=((4.0D0*u(2,j))-(3.0D0*u(1,j))-u(3,j))/
> (2.0D0*de
> $lta)
> do 8 i=2,NSTEPS
> uprime(i,j)=(u(i+1,j)-u(i-1,j))/(2.0D0*delta)
> 8 continue
> uprime(NSTEPS+1,j)=((3.0D0*u(NSTEPS+1,j))+u(NSTEPS-1,j)-
> $(4.0D0*u(NSTEPS,j)))/(2.0D0*delta)
> 5 continue
[snip]

Others have already suggested to use a single two digit unit number for all
file output.

In addition to replacing **2.0D0 with **2, there are many other ways to
simplify this program. This program did compile and run at a reasonable
speed, once the I/O problem was fixed. I still suggest you do this for
pedagogical reasons.

In particular, I would convert this program to Fortran 95 (automated
conversion programs are available) and use its freatures (implicit none,
free form source, whole array operations, do .. end do loops, if () exit,
etc.

In particular the arrays A,C,D and E can be eliminated. A,C and D are almost
constant for the same value of J.
Once you combine expressions on successive lines, ALPHA and BETA only run
from (2 to nsteps+1,:).

If you set up two new variables for B and F, you can calculate contributions
from U and UPRIME for each value of J instead of later. Thus is possible to
make your all your arrays one dimensional.

Of course, a reasonable value of NSTEPS (lower than 10,000) might be in
order too...

If someone is interested I could work this through and post a revised
program.

--- Elliot





First  |  Prev  |  Next  |  Last
Pages: 1 2 3
Prev: call c lib from fortran in ms visual studio 2005
Next: What exactly does the SAVE statement do ?