|
forums.silverfrost.com Welcome to the Silverfrost forums
|
View previous topic :: View next topic |
Author |
Message |
Ronan
Joined: 27 Jun 2007 Posts: 19
|
Posted: Tue May 25, 2010 8:18 am Post subject: 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(
INTEGER::n,i,ki,l,kj,j,ll,m,k
REAL(iwp):
! 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. |
|
Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006 Posts: 2554 Location: Sydney
|
Posted: Tue May 25, 2010 8:51 am Post subject: |
|
|
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 Code: | IMPLICIT NONE
INTEGER,PARAMETER::iwp=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,PARAMETER::iwp=SELECTED_REAL_KIND(15)
REAL(iwp),INTENT(IN OUT)::kv(8)
INTEGER,INTENT(IN)::kdiag(8)
!
INTEGER::n,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 ?? |
|
Back to top |
|
|
Ronan
Joined: 27 Jun 2007 Posts: 19
|
Posted: Tue May 25, 2010 9:07 am Post subject: |
|
|
sorry
In main program kv and kdiag are declared as:
Code: |
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
Code: | 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 |
|
Back to top |
|
|
IanLambley
Joined: 17 Dec 2006 Posts: 490 Location: Sunderland
|
Posted: Tue May 25, 2010 12:16 pm Post subject: |
|
|
Have you actually allocated them before use? |
|
Back to top |
|
|
Ronan
Joined: 27 Jun 2007 Posts: 19
|
Posted: Tue May 25, 2010 12:25 pm Post subject: Invalid floating point operation |
|
|
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. |
|
Back to top |
|
|
Ronan
Joined: 27 Jun 2007 Posts: 19
|
Posted: Tue May 25, 2010 3:14 pm Post subject: Invalid floating point operation |
|
|
Any sugestion please ? |
|
Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006 Posts: 2554 Location: Sydney
|
Posted: Wed May 26, 2010 2:34 am Post subject: |
|
|
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 |
|
Back to top |
|
|
Ronan
Joined: 27 Jun 2007 Posts: 19
|
Posted: Wed May 26, 2010 8:18 am Post subject: |
|
|
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
Code: | 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, |
|
Back to top |
|
|
JohnHorspool
Joined: 26 Sep 2005 Posts: 270 Location: Gloucestershire UK
|
Posted: Wed May 26, 2010 10:04 am Post subject: |
|
|
Quote: | 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. |
|
Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006 Posts: 2554 Location: Sydney
|
Posted: Wed May 26, 2010 12:49 pm Post subject: |
|
|
Quote: | 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 ? |
|
Back to top |
|
|
Ronan
Joined: 27 Jun 2007 Posts: 19
|
Posted: Thu May 27, 2010 10:59 am Post subject: |
|
|
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
Quote: | 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, |
|
Back to top |
|
|
JohnHorspool
Joined: 26 Sep 2005 Posts: 270 Location: Gloucestershire UK
|
Posted: Thu May 27, 2010 12:18 pm Post subject: |
|
|
Quote: | For stifness matrix and its solutions they are real |
real*4 or real*8 ?
or in other words have you specified double precision or left it at the default single precision? |
|
Back to top |
|
|
Ronan
Joined: 27 Jun 2007 Posts: 19
|
Posted: Thu May 27, 2010 12:35 pm Post subject: |
|
|
I suppose it's default single precision
declaration is :
Code: | 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 |
|
|
Back to top |
|
|
JohnHorspool
Joined: 26 Sep 2005 Posts: 270 Location: Gloucestershire UK
|
Posted: Thu May 27, 2010 1:18 pm Post subject: |
|
|
I suggest you try double precision for all your reals. |
|
Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006 Posts: 2554 Location: Sydney
|
Posted: Sun May 30, 2010 1:35 pm Post subject: |
|
|
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 |
|
Back to top |
|
|
|
|
You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|