forums.silverfrost.com

Ronan · Joined: 27 Jun 2007 Posts: 19

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( Smile

INTEGER,INTENT(IN)::kdiag( Smile

INTEGER::n,i,ki,l,kj,j,ll,m,k
REAL(iwp): Mad

! 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.

JohnCampbell · Joined: 16 Feb 2006 Posts: 2554 Location: Sydney

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

Ronan · Joined: 27 Jun 2007 Posts: 19

sorry
In main program kv and kdiag are declared as:

IanLambley · Joined: 17 Dec 2006 Posts: 490 Location: Sunderland

Have you actually allocated them before use?

Ronan · Joined: 27 Jun 2007 Posts: 19

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.

Ronan · Joined: 27 Jun 2007 Posts: 19

Any sugestion please ?

JohnCampbell · Joined: 16 Feb 2006 Posts: 2554 Location: Sydney

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

Ronan · Joined: 27 Jun 2007 Posts: 19

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

JohnHorspool · Joined: 26 Sep 2005 Posts: 270 Location: Gloucestershire UK

JohnCampbell · Joined: 16 Feb 2006 Posts: 2554 Location: Sydney

Ronan · Joined: 27 Jun 2007 Posts: 19

The codes are FEM for structural analyse, I've made mistake with unit conversion Smile

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

JohnHorspool · Joined: 26 Sep 2005 Posts: 270 Location: Gloucestershire UK

Ronan · Joined: 27 Jun 2007 Posts: 19

I suppose it's default single precision

declaration is :

JohnHorspool · Joined: 26 Sep 2005 Posts: 270 Location: Gloucestershire UK

I suggest you try double precision for all your reals.

JohnCampbell · Joined: 16 Feb 2006 Posts: 2554 Location: Sydney

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