Precision in matrix algebra applications

Hi,

I am new to Ada and I am trying to develop Finite Element programs for my own engineering use. In order to gain some insight into programming in Ada I developed and compared results with a Fortran program I made and trust (it compares very well with the well established software called Code-Aster, developed by Electrecite-de-France).

Although the results are very good I found big % discrepancies in the calculation of displacements in a handful of very small numbers (4 bad results in 774 with magnitude of about 1e-9 to 1e-11, which are 3 orders of magnitude smaller than the maximum results), which I suspect is related to the precision of the numerical calculations. I searched the web and books but could not found any solution therefore here I am.

Anyway, my question is: How do I declare Real Vectors and Real Matrix types with 18 digits of precision? Any help is greatly appreciated!

Regards,
Gilberto

2 Likes

If you need 18 digits, then say so:

type Real is digits 18;

package Real_Matrices is new Ada.Numerics.Generic_Real_Arrays (Real => Real);

Your compiler may reject this if it cannot provide 18 digits. If you want the maximum precision for your compiler, you can say that, too:

type Real is digits System.Max_Digits;

This will give you more than 18 digits if the compiler supports it.

5 Likes

Thank you JC001 for helping me, I appreciate.

That said, I gave it a try and there was no change at all in the results. I am getting to the conclusion that the problem might be something else. Perhaps it is the solver I am using…

Thank you again,
Gilberto

Hi Gilberto,

it is great to see another FEM/Fortran/Code_Aster user in the Ada forums! I am also Irvise in the CA forums and the guy that did the tutorials about OpenTURNS-Persalys in YouTube.

Here are my two cents, as I am no Ada expert. Afaik, GNAT uses IEEE 754 for floats, so the behaviour should be very similar (same) as Fortran’s/CA’s. However, I would check whether the assembly generated by the routines is the same. I believe the best tool to check for that is https://godbolt.org/ It has Ada and GFortran support. I would also make sure the optimisation level -Ox (with x in 0, 1, 2, 3, fast, but please, do not use fast for FEM :slight_smile: ) is the same for both compilers. I would not be surprised if Ada’s treatment of Floats has a few checks/less optimisations in them with lower optimisation levels.

And how are you comparing the results? Ada’s printing of floats is not as standard as it is in other more common languages. Though it is quite close to what Fortran does.

Also, in case you are not aware, Ada/GNAT, sports an official binding support for Fortran code. You can read more about it in the Ada Reference Manual, Annex B.5 Interfacing with Fortran.

Cheers,
Fer

Cool that you start doing FE in Ada!
I did the same a few decades ago, which reminds me I should push more of my FE & Ada resources on the Web. A few bits of them are in an open-source library called “mathpaqs”.

As an update I managed to improve the results. It turns out the problem was a “human interface” kind of problem (me!). Anyways, I was working with Floats and I was ignorant of the fact that using floats I only get 6 digits… (or so it appears; I might be wrong again since I am a neophyte)…

I changed my little amateur program to work with Reals and thanks to help given by Jeff I managed to get results very close to my Fortran code (and therefore to Code-Aster). The problem then became “raised STORAGE_ERROR : stack overflow or erroneous memory access”…

Now I have to work on the memory problem…Well, but at least I got to solve the precision problem…

Thank you again Jeff.

1 Like

Hi Irvise,

It is great to see you here and thank you for the help provided! As mentioned in my updated message above I managed to fix the problem I had created myself (I got laugh at myself for that…).

Anyways, your help is also greatly appreciated!

Thank you,
Gilberto

Hi Gautier,

My FE programs are little amateur hobby things as I am not a professional in numerical methods. It is more a way of exercising the brain with linear algebra, programming, continuum mechanics, etc. As it turns out I got to find your Ada codes so thank you for that as FE in Ada seems rare in the web…

Regards,
Gilberto

1 Like

ARM 3.5.7 (14) defines the required precision of Float: “In an implementation that supports floating point types with 6 or more digits of precision, the requested decimal precision for Float shall be at least 6.” In plain English, this means that Float may have any precision at all. Again, the message is clear: if you need a specific precision, don’t use Float.

I’m not sure whether my GNAT Math Extensions project (crate gnat_math_extensions) is at all relevant – just the equivalent of ARM Annex G, without the restrictions. I wrote it after a question on StackOverflow caught my interest; not an actual user.

1 Like

It usually comes from attempt to allocate huge matrix of the stack. So, avoid return matrix as function return value and avoid to declare auxiliary matricies in declarative part of the subprogram/code block.

If you are using Ada.Numerics… It is time to rewrite it, unless API was changed last years.

PS. I’ve own package for huge matrix operations because of this. :disappointed:

2 Likes

My two cents here, but I would recommend against that unless there is a good reason. I would personally use/consume an optimised matrix library such as *BLAS or an specialised solver, such as MUMPS or PETSc :slight_smile: