Invalid floating point operation !

Dear community members,

I'm experiencing the difficulties when trying to assign real type variable to another real type variable, traversed through the previous post in the forums and found that its conversant and some of them have come acrros with such issue before (sligthly different case when assigning integer to real) and reported as underflow bug . In my case it shouldn't be a problem as long as i've both vars with same types. To test it I've downloaded the official FTN95 v5.5 personal version where developers says that this problem has been eliminated but result for me is same.

Any suggestions how to overcome this impediment will be appreciated!

sub-unit is as follows:

IMPLICIT NONE INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15) REAL(iwp),INTENT(IN OUT)::kv(:) INTEGER,INTENT(IN)kdiag(:) INTEGERn,i,ki,l,kj,j,ll,m,k REAL(iwp)::x

! code ommited here

kv(ki+i)=SQRT(x) !! HERE COMPLAINS THE COMPILER IF X VALUE IS AROUND or smaller than 1..E-16 END DO RETURN

Kind Regards,

PS. I'm a bit dubious that I've downloaded the latest version (don't know whether it is for free or commercial affiliates) , cause Plato IDE about screen shows the version 4.3.

I tried to identify your problem as best I could, but I could not reproduce the problem. Try using the Code button when copying in sample. My best estimate of your code (ronan.f95) was IMPLICIT NONE INTEGER,PARAMETERiwp=SELECTED_REAL_KIND(15) REAL(iwp) :: kv(8) INTEGER :: kdiag(8) ! write (,) 'iwp = ',iwp call ronan (kv, kdiag) end subroutine ronan (kv, kdiag) ! IMPLICIT NONE INTEGER,PARAMETERiwp=SELECTED_REAL_KIND(15) REAL(iwp),INTENT(IN OUT)kv(8) INTEGER,INTENT(IN)kdiag(8) ! INTEGERn,i,ki,l,kj,j,ll,m,k REAL(iwp) x

! code ommited here 
i = 2
x = 1.e-18
!
kv = 0
do ki = 1,6
   kv(ki+i)=SQRT(x) !! HERE COMPLAINS THE COMPILER IF X VALUE IS AROUND or smaller than 1..E-16 
END DO 
RETURN 
end

How was x defined ? I compiled and ran with ftn95 ronan.f95 /check /link sdbg ronan.exe with no identified problem ??

sorry In main program kv and kdiag are declared as:

program Ronan
IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 INTEGER,ALLOCATABLE::kdiag(:)
 REAL(iwp),ALLOCATABLE::kv(:)
 .........
  call sparin(kv, kdiag)
 .......
end program Ronan

!and here is the routine

SUBROUTINE sparin(kv,kdiag)
IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::kv(:)
 INTEGER,INTENT(IN)::kdiag(:)
 INTEGER::n,i,ki,l,kj,j,ll,m,k
 REAL(iwp)::x
 n=UBOUND(kdiag,1)  
 kv(1)=SQRT(kv(1))
 DO i=2,n
   ki=kdiag(i)-i
   l=kdiag(i-1)-ki+1
   DO j=l,i
     x=kv(ki+j)
     kj=kdiag(j)-j
     IF(j/=1)THEN
       ll=kdiag(j-1)-kj+1
       ll=max(l,ll)
       IF(ll/=j)THEN
         m=j-1
         DO k=ll,m 
           x=x-kv(ki+k)*kv(kj+k) 
         END DO
       END IF
     END IF	   
     !IF (kv(kj+j)/=0) THEN 
         kv(ki+j)=x/kv(kj+j) 
     !END IF
   END DO
   kv(ki+i)=SQRT(x)
 END DO
RETURN
END SUBROUTINE sparin

Compiles normally w/o errors

Have you actually allocated them before use?

Yes I'm sure

because that problem doesn't arrise always, for example if I change the input values then program shows the correct outputs flawlessly with no complaint from compiler.

Any sugestion please ?

Ronan,

Your code looks like a Cholesky-Crout decomposition. If you are having problems with SQRT, why not change to not taking the square root of the diagonal. This solution, which I use, is well documented, and more stable. Alternatively, if the diagonal is becomiong zero or -ve, you could set the value to 1, effectively fixing that equation, as it is not well defined by the data. You should be testing the value of the diagonal at each step. As a start test if x < 1.e-16 ( or a small value for your problem), note the occurrence and set x = 1 and continue. You are effectively adding a further equation 'F = 1 x X' to the problem definition.

John

Thanks in advance John,

I think I misleaded you, compiler doesn't rise the exception when values are smaller than 1.e-16 or something, the problem is with negative value in sqrt function, as you guessed, I've added

    if (x > 0) then  
      kv(ki+i)=SQRT(x) 
    else 
      kv(ki+i)=1  
    end if

but outputs are wrong! Does that altering in algorithm change the whole solution of equations? If it's so, how could possibly the original writer thought to deal with negative numbers?

Regards,

how could possibly the original writer thought to deal with negative numbers?

Depends on the original intention for use of the code by the writer. For instance with a stiffness matirx you should never encounter this problem. If you do then the input data is wrong, not the code.

but outputs are wrong!

Not necessarily, do a back substitution and check the error. You should find that values associated with the equation that had a negative sqrt may not have a significant effect on the correct answer. If the system of equations should be symmetric positive definite, then large negatives say that the equations do not look good and you should check the values you are using. Are all calculations in real*8 ?

The codes are FEM for structural analyse, I've made mistake with unit conversion 😃 corrected the values (EA, EIx, EIy and GJ) results are still near to the correct but not the exact and have some precision problem (e.g. my result = -0.7952; correct =-0.08425 ! precision problem here my result = -6.093; correct =-0.6100 ! precision problem here my result = -110.4; correct =-111.06 ! but this one is near

Are all calculations in real*8 ?

For stifness matrix and its solutions they are real, but for node nums, element nums, connectivity are integers.

Thanks for patience,

For stifness matrix and its solutions they are real

real4 or real8 ?

or in other words have you specified double precision or left it at the default single precision?

I suppose it's default single precision

declaration is :

IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 INTEGER::fixed_freedoms,i,iel,k,loaded_nodes,ndim,ndof,nels,neq,nod=2,   &
   nodof,nn,nprops,np_types,nr,nlen
 REAL(iwp)::penalty=1.0e20_iwp,zero=0.0_iwp
!-----------------------dynamic arrays------------------------------------
 INTEGER,ALLOCATABLE::etype(:),g(:),g_g(:,:),g_num(:,:),kdiag(:),nf(:,:), &
   no(:),node(:),num(:),sense(:)
 REAL(iwp),ALLOCATABLE::action(:),coord(:,:),eld(:),gamma(:),g_coord(:,:),&
   km(:,:),kv(:),loads(:),prop(:,:),value(:)
 CHARACTER(LEN=15)::argv, load_type

I suggest you try double precision for all your reals.

You said you were writing a FEM program. If the problem is that your stiffness matrix is not positive definite, it probably is that the accumulated values in your stiffness matrix are wrong. Sources of this error may be not using real*8, as John is suggesting, or your error may be more basic, that the element stiffness matrices are wrong. You should check their formulation. My suggestion of setting the diagonal value to 1 is adding a 'spring' stiffness to that freedom, the stiffness being = '1 - diagonal value' As the diagonal value should not be -ve, this is only masking the problem, so that a general solution can be found. This is typical of modelling of shell bending elements, where only 5 degrees of freedom are provided by the shell element, while the node has 6 general freedoms, 3 displacments and 3 rotations. If the elements are beams, they are defined for all 6 freedoms, so should not have zero diagonal, unless your boundary conditions are wrong (check torsion), so back to checking the FE problem. If the problem persists, check this as a possibility.

John