0

In a simple program I want to build a matrix by defining a function. The problem is that a variable declared and initialized in the program has the exact assinged value of zero (zero_test)and some components of the matrix(D(4,1), D(1,4) etc.) which are assigned to 0., are not exactly zero. The latter have a value of order E-308 or E-291etc. I wonder why is there a difference.

Due to some articles I read, like this, the other components(D(1,1), D(1,2) etc.) are imprecise because of the transformation to the binary format.

System info: IA-32, Intel Visual Fortran 11.1.048 integrated with Microsoft Visual Studio 2008

The code:

program test

implicit none

real(8),parameter :: LAMBDA=75.e9,MU=50.e9
integer,parameter :: ndi=3,ntens=4
real(8) :: D(ntens,ntens),zero_test=0.


D = clcElasticStiffnessMatrix(LAMBDA,MU,ndi,ntens)

contains

        function clcElasticStiffnessMatrix(LAMBDA,MU,ndi,ntens)

        implicit none

        integer :: ndi,ntens,i,j
        real(8) :: clcElasticStiffnessMatrix(ntens,ntens),LAMBDA,MU


        do i=1,ndi    
        do j=i,ndi
            if(i .eq. j) then
                clcElasticStiffnessMatrix(i,j) = 2.*MU + LAMBDA
            else 
                clcElasticStiffnessMatrix(i,j) = LAMBDA
                clcElasticStiffnessMatrix(j,i) = clcElasticStiffnessMatrix(i,j)
            end if
        end do
        end do

        do i=ndi+1,ntens
        do j=i,ntens
            if(i .eq. j) then
                clcElasticStiffnessMatrix(i,j) = MU
            else 
                clcElasticStiffnessMatrix(i,j) = 0.
                clcElasticStiffnessMatrix(j,i) = clcElasticStiffnessMatrix(i,j)
            end if

        end do
        end do      

        end function

end program

Matrix D in break mode:

 D:
   174999994368.000        74999996416.0000        74999996416.0000       2.641443384627243E-308
   74999996416.0000        174999994368.000        74999996416.0000       2.640433316727162E-308
   74999996416.0000        74999996416.0000        174999994368.000      -1.051992669438322E-291
  2.640110775815455E-308  0.000000000000000E+000  6.151018477594351E-318   49999998976.0000
Community
  • 1
  • 1
user3410012
  • 53
  • 2
  • 6
  • Is there any specific question you need an answer about that behavior? – innoSPG Mar 22 '14 at 17:58
  • Yes I want to know why the assigned variables in the matrix have a different value from zero. – user3410012 Mar 22 '14 at 22:32
  • You should not look for exact equality in floatting point calculations. It is strongly discouraged. Modern compilers will warn you of such comparison, instead, you should use the epsilon function. – innoSPG Mar 23 '14 at 02:14
  • if i understand the question its not about calculations. When a value is *assigned* to 0. it should in fact pass a `.eq.0.` test. See what you get if you set lambda an mu to zero as well. – agentp Mar 23 '14 at 13:21
  • the error in this case is in your loop construct. you are not assigning *anything* to those matrix components that are showing up with the `E-308` junk. My advice here initialise the entire matrix to zero, then assign the nonzero components in the loop. – agentp Mar 23 '14 at 13:27
  • You are right. Something was wrong with my loop index range. It should be `j=1,i`. Thanks. – user3410012 Mar 23 '14 at 21:50
  • Is there any debugging option (flag) that could throw an error when a variable, which is not assigned to a value, is used in a statement? – user3410012 Mar 23 '14 at 22:05
  • You can use ` -finit-real=nan` in gfortran or the valgrind debugger. – Vladimir F Героям слава Mar 24 '14 at 09:10

1 Answers1

0

The problem is that you are assigning real(4) values to real(8) variables when you write

LAMBDA=75.e9,MU=50.e9

or

clcElasticStiffnessMatrix(i,j) = 0.

Try with d instead e and specify always:

LAMBDA=75.d9,MU=50.d9
clcElasticStiffnessMatrix(i,j) = 0.d0
Alfonso Sancho
  • 136
  • 6