vanand@molotov.eng.clemson.edu (Anand Lakshmikumaran) (06/05/90)
Hi Folks !
I encountered a very odd error in a couple of my FORTRAN programs.
I have debugged the programs to the point that i know where the error occurs
but can't understand why that error should occur. I have simplified one of
the programs and am attaching it below so that if anyone have faced a similar
problem before or know how to avoid this error could kindly reply to me by email.
This program is just a part of a bigger program.
The problem is:
I am trying to invert a matrix 'b' (which in this example is of order 4
but could be much bigger actually) using a subroutine. 'write statement' no. 1
gives the correct value but 'write statement' no. 2 which SHOULD give the same value
doesn't. I have attached the output also.
The odd part of it is that if i change the dimension of 'b' in my first statement
to 'b(4,4)' (i.e. which is the actual array size for this particular example) then
the output is accurate. The correct values of 'b inverse' is also attached.
P.S. I have tried to run it on both SUN Fortran and VAX Fortran. Both of them
give the same error.
My email address is vanand@eng.clemson.edu
Thanks for reading through this.
-ANAND V. L.
vanand@eng.clemson.edu
dimension b(50,50)
n = 4
do i = 1,n,1
b(i,i) = 4.0
if(i.ne.(n))then
b(i,i+1) = 1.
b(i+1,i) = 1.
endif
enddo
c*** WRITE STATEMENT NO.1 *******************************************
do i = 1,n,1
write(*,*)(b(i,j),j=1,n,1)
enddo
C**************************************************************************
call minv(n,b)
C*** WRITE STATEMENT NO. 4 *******************************************
do i = 1,n,1
write(*,*)(b(i,j),j=1,n,1)
enddo
C**************************************************************************
end
C ********************************************************
C MINV IS MATRIX INVERSE.
C ********************************************************
SUBROUTINE MINV(N,A)
dimension a(n,n)
C**********************************************************************
C modified from p. 197 of James,Smith,& Wolford
C matrix inversion using Gauss-Jordan reduction
C inverted matrix overlays original matrix in memory
C
C input: matrix a to be inverted, N=order of matrix [A]
C
C output: a is replaced by its inverse (original destroyed)
C
C calculate elements of new matrix
C**********************************************************************
C*** WRITE STATEMENT NO. 2 ****************************************
do 100 i = 1,n,1
write(*,*)(a(i,j),j=1,n,1)
100 continue
C************************************************************************
DO 600 K=1,N
C calculate new elements of pivot row
do 20 J=1,N
IF (J.ne.K) THEN
A(K,J)=A(K,J)/A(K,K)
endif
20 continue
C calculate element replacing pivot element
A(K,K)=1.d0/A(K,K)
C calculate new elements not in pivot row or pivot column
do 40 I=1,N
do 30 J=1,N
if ((I.ne.K).and.(J.ne.K)) then
A(i,j)=A(i,j)-A(k,j)*A(i,k)
endif
30 continue
40 CONTINUE
C calculate replacement elements for pivot column
C (except pivot element)
do 50 i=1,N
if(i.ne.K) a(i,K)=-a(i,K)*a(K,K)
50 CONTINUE
600 continue
C*** WRITE STATEMENT NO. 3 ********************************************
do 200 i = 1,n,1
write(*,*)(a(i,j),j=1,n,1)
200 continue
C****************************************************************************
return
end
THE OUTPUT ON SUN FORTRAN WAS AS FOLLOWS
WRITE STATEMENT NO. 1
4.00000 1.00000 0. 0.
1.00000 4.00000 1.00000 0.
0. 1.00000 4.00000 1.00000
0. 0. 1.00000 4.00000
WRITE STATEMENT NO. 2
4.00000 0. 0. 0.
1.00000 0. 0. 0.
0. 0. 0. 0.
0. 0. 0. 0.
WRITE STATEMENT NO. 3
-\NA\N \NA\N -\NA\N -\NA\N
\NA\N \NA\N \NA\N \NA\N
-\NA\N -\NA\N -\NA\N -\NA\N
-\NA\N -\NA\N -\NA\N -\NA\N
WRITE STATEMENT NO. 4
-\NA\N 1.00000 0. 0.
\NA\N 4.00000 1.00000 0.
-\NA\N 1.00000 4.00000 1.00000
-\NA\N 0. 1.00000 4.00000
THE CORRECT VALUES OF THE MATRIX B AND ITS INVERSE ARE
MATRIX B
4.0 1.0 0.0 0.0
1.0 4.0 1.0 0.0
0.0 1.0 4.0 1.0
0.0 0.0 1.0 4.0
B INVERSE
0.267943 -7.17703E-02 1.91388E-02 -4.78469E-03
-7.17703E-02 0.287081 -7.65550E-02 1.91388E-02
1.91388E-02 -7.65550E-02 0.287081 -7.17703E-02
-4.78469E-03 1.91388E-02 -7.17703E-02 0.267943sun@me.utoronto.ca (Andy Sun Anu-guest) (06/05/90)
In article <9206@hubcap.clemson.edu> vanand@molotov.eng.clemson.edu (Anand Lakshmikumaran) writes: >Hi Folks ! > I encountered a very odd error in a couple of my FORTRAN programs. >I have debugged the programs to the point that i know where the error occurs >but can't understand why that error should occur. I have simplified one of >the programs and am attaching it below so that if anyone have faced a similar >problem before or know how to avoid this error could kindly reply to me by email. >This program is just a part of a bigger program. > The problem is: > I am trying to invert a matrix 'b' (which in this example is of order 4 >but could be much bigger actually) using a subroutine. 'write statement' no. 1 >gives the correct value but 'write statement' no. 2 which SHOULD give the same value >doesn't. I have attached the output also. > The odd part of it is that if i change the dimension of 'b' in my first statement >to 'b(4,4)' (i.e. which is the actual array size for this particular example) then >the output is accurate. The correct values of 'b inverse' is also attached. > >P.S. I have tried to run it on both SUN Fortran and VAX Fortran. Both of them > give the same error. > >My email address is vanand@eng.clemson.edu > >Thanks for reading through this. >-ANAND V. L. >vanand@eng.clemson.edu > > > dimension b(50,50) > n = 4 > do i = 1,n,1 > b(i,i) = 4.0 > if(i.ne.(n))then > b(i,i+1) = 1. > b(i+1,i) = 1. > endif > enddo > >c*** WRITE STATEMENT NO.1 ******************************************* > > do i = 1,n,1 > write(*,*)(b(i,j),j=1,n,1) > enddo > >C************************************************************************** > > call minv(n,b) > >C*** WRITE STATEMENT NO. 4 ******************************************* > > do i = 1,n,1 > write(*,*)(b(i,j),j=1,n,1) > enddo > >C************************************************************************** > > end > The rest of the stuff deleted... A lot of th first year students I tutored made this mistake. Nothing wrong with the subroutine, what's wrong is the way you pass the parameters. You are partly correct in saying that changing the dimension to (4,4) works. In fact, if you declared an (m,n) array in the main program, no matter how small a portion that you are passing to a subroutine (e.g. you've declared dimension (100,100) but just use (2,2)), you HAVE TO pass the maximum ROW element over. In my example, you have to pass (100,2) as parameters, not (2,2) or else you array got all messed up. In your case, you b is ALWAYS 50 if (50,50) is what you've declared. Fortran stores your 2D array (say, a 3x3) sequentially like this: 1 (1,1) 2 (1,2) 3 (1,3) (1,1) (1,2) (1,3) 4 (2,1) (2,1) (2,2) (2,3) 5 (2,2) (3,1) (3,2) (3,3) 6 (2,3) 7 (3,1) 8 (3,2) 9 (3,3) if you only pass a 2x2 over to a subroutine, this is what you are actually passing: 1 (1,1) 2 (1,2) (1,1) (1,2) 3 (1,3) (1,3) (2,1) 4 (2,1) You see why it get all messed up. If you are a good programmer and initialized the array to zeros first once you delcared it, you will find zeros all over the place if you pass the array the way you did. If you don't initialize it, you can get ANYTHING: random numbers depending on what those memory locations got in them. Andy _______________________________________________________________________________ Andy Sun | Internet: sun@me.utoronto.ca University of Toronto, Canada | UUCP : ...!utai!me!sun Dept. of Mechanical Engineering | BITNET : sun@me.utoronto.BITNET