Dear All,
I am using Dgetri to invert a square matrix that I have. The matrix consists of the Radial Basis values of a grid that I have, i.e., RBF(Nnode,Nnode).
I have used two different compilers to invert the matrix. The first one is in the environment of Microsoft powerstaion 4 and the other is using inter fortran. The problem is that when I create the RBF matrix, all the components in two different codes are the same but after the inversion, they are not, though the values are close.

calculated after the inverting with intel:
89588.4169487697 -400462.914663885 993255.716366269
-1815020.42991486 2737745.90690058 -3621529.66069135
4359302.34907216 -4892909.07940900 5208790.59570593
-5324438.63545066

the old calculation, (I have trust in this results because it is verified) this one takes a much longer time for inversion, though.
92178.015799056420000 -414516.040633608200000 1035342.675828539000000
-1907313.242167841000000 2903705.975803736000000 -3881289.913598144000000
4726376.390789632000000 -5372746.814054917000000 5798863.224470832000000
-6015613.975911859000000 6052730.233410349000000 -5947828.868681677000000

I was thinking of the precision problem, but it cannot be the case because the calculation of the RBF function, which includes calculating the distance of each point in the grid with respect to the other points, is correct for both codes. They are equal to each other. Now, I think maybe the Intel compiler is somehow reducing the accuracy in the calculation in the inversion step. Any idea is highly welcome. Thank you in advance.

Cheers,

Farzad

Thank you for

On Sunday, September 3, 2023 at 1:43:40 PM UTC+2, Farzad Tatar wrote:
> Dear All,
> I am using Dgetri to invert a square matrix that I have. The matrix consists of the Radial Basis values of a grid that I have, i.e., RBF(Nnode,Nnode).
> I have used two different compilers to invert the matrix. The first one is in the environment of Microsoft powerstaion 4 and the other is using inter fortran. The problem is that when I create the RBF matrix, all the components in two different codes are the same but after the inversion, they are not, though the values are close.
>
> calculated after the inverting with intel:
> 89588.4169487697 -400462.914663885 993255.716366269
> -1815020.42991486 2737745.90690058 -3621529.66069135
> 4359302.34907216 -4892909.07940900 5208790.59570593
> -5324438.63545066
>
> the old calculation, (I have trust in this results because it is verified) this one takes a much longer time for inversion, though.
> 92178.015799056420000 -414516.040633608200000 1035342.675828539000000
> -1907313.242167841000000 2903705.975803736000000 -3881289.913598144000000
> 4726376.390789632000000 -5372746.814054917000000 5798863.224470832000000
> -6015613.975911859000000 6052730.233410349000000 -5947828.868681677000000
>
> I was thinking of the precision problem, but it cannot be the case because the calculation of the RBF function, which includes calculating the distance of each point in the grid with respect to the other points, is correct for both codes. They are equal to each other. Now, I think maybe the Intel compiler is somehow reducing the accuracy in the calculation in the inversion step. Any idea is highly welcome. Thank you in advance.
>
> Cheers,
>
> Farzad
>
> Thank you for

I do not quite understand the results you post: the result for Intel has 10 numbers, for MicroSoft powerstation there are 12 numbers - neither set of numbers fit in a square matrix ;). So, I have no idea if the two are really close as you say.

You should beware that the inverse of a matrix can be difficult to compute accurately, though these seem to be small. You can check the results by multiplying the original matrix with the inverse matrix (or even by multiplying the inverse matrix by the original). The answer should be the identity matrix or nearly so. Differences may occur because of rounding effects.

It would help if you could show your code.

One last remark: MicroSoft Powerstation has been out of support for several decades.

Regards,

Arjen

> I do not quite understand the results you post: the result for Intel has 10 numbers, for MicroSoft powerstation there are 12 numbers - > neither set of numbers fit in a square matrix ;). So, I have no idea if the two are really close as you say.
Sorry for the misleading information. The numbers are the first few elements of the first column just to be able to compare them. the matrix shape is 1400by 1400. for this case.

> You should beware that the inverse of a matrix can be difficult to compute accurately, though these seem to be small. You can check the results by multiplying the original matrix with the inverse matrix (or even by multiplying the inverse matrix by the original). The answer should be the identity matrix or nearly so. Differences may occur because of rounding effects.

I will try this though I already know that the result is not accurate. I actually tried REAL128 from the intrinsic kind function in the modern fortran and the results got closer to the one that I was receiving from the old code. It means the reason might be due to the precision. But aren't Real*8 and Real64 the same?

> It would help if you could show your code.
This subroutine is a very small part of my project. Here is the part I am dealing with.

SUBROUTINE RBF_inv(iflag)
USE COMMON_1

USE, INTRINSIC :: ISO_FORTRAN_ENV, DP => REAL64

INTEGER, INTENT(IN) :: iflag
INTEGER :: ii, jj, i1, j1, INFO, NRHS, M, N, ising, lwork
INTEGER, DIMENSION(NNODE_LS) :: IPIV
REAL (DP):: xi_, yi_, xj_, yj_, dist_, temp, determ
REAL (DP):: F_intp
REAL (DP):: xni, yni, xnj, ynj,xw_ls_norm
REAL(DP), ALLOCATABLE:: work(:)

! Determine the coeff of the Radial Basis Functions
! This approach is used to interpolate the function everywhere in the domain

! 1. Invert RBF_LS only at the 1st iteration (iflag==0) , then I will have the RBF coeffs

IF (iflag == 0) THEN
WRITE(*,*) ' Solve global RBF iter', iter, Nnode_LS !nnode_ls is the number of node in each row.

! Matrix of the RBF function values
DO ii = 1, Nnode_LS
xi_ = xy_lsgrid(1, ii)
yi_ = xy_lsgrid(2, ii)
DO jj = 1, Nnode_LS
xj_ = xy_lsgrid(1, jj)
yj_ = xy_lsgrid(2, jj)
dist_ = SQRT((xj_ - xi_)**2.0_dp + (yj_ - yi_)**2.0_dp)
RBF_LS(ii, jj) = EXP(-((dist_ * epsw)**2.0_dp))
!IF (RBF_LS(II,JJ) <= 1.0E-58_DP) RBF_LS(II,JJ)=0.0_DP
END DO
END DO

! Solve the linear system
NRHS=1
print*, "1:", rbf_ls(1:10,1)
! Perform LU decomposition of RBF_LS
CALL DGETRF(NNODE_LS, NNODE_LS, RBF_LS, NNODE_LS, ipiv, info)
print*, "2:composed LU", rbf_ls(1:10,1)
IF (info == 0) THEN
! RBF_LS was successfully factorized, now compute its inverse
ALLOCATE(work(lwork)) !this is work vector
CALL DGETRI(NNODE_LS, RBF_LS, NNODE_LS, ipiv, work, nnode_ls, info)
print*, "3:", rbf_ls(1:10,1)
IF (info == 0) THEN
! A_inv contains the inverse of A
! Use A_inv as needed
ELSE
PRINT *, "Error computing the inverse of A."
END IF
DEALLOCATE(work)
ELSE
PRINT *, "Error in LU decomposition of A."
END IF

> One last remark: MicroSoft Powerstation has been out of support for several decades.
The point is that this old method is giving a more accurate result in comparison with the new one.

Thank you,

Farzad

On Monday, September 4, 2023 at 8:34:45 AM UTC+12, Farzad Tatar wrote:
> > I do not quite understand the results you post: the result for Intel has 10 numbers, for MicroSoft powerstation there are 12 numbers - > neither set of numbers fit in a square matrix ;). So, I have no idea if the two are really close as you say.
> Sorry for the misleading information. The numbers are the first few elements of the first column just to be able to compare them. the matrix shape is 1400by 1400. for this case.
> > You should beware that the inverse of a matrix can be difficult to compute accurately, though these seem to be small. You can check the results by multiplying the original matrix with the inverse matrix (or even by multiplying the inverse matrix by the original). The answer should be the identity matrix or nearly so. Differences may occur because of rounding effects.
> I will try this though I already know that the result is not accurate. I actually tried REAL128 from the intrinsic kind function in the modern fortran and the results got closer to the one that I was receiving from the old code. It means the reason might be due to the precision. But aren't Real*8 and Real64 the same?
> > It would help if you could show your code.
> This subroutine is a very small part of my project. Here is the part I am dealing with.
>
> SUBROUTINE RBF_inv(iflag)
> USE COMMON_1
>
> USE, INTRINSIC :: ISO_FORTRAN_ENV, DP => REAL64
>
> INTEGER, INTENT(IN) :: iflag
> INTEGER :: ii, jj, i1, j1, INFO, NRHS, M, N, ising, lwork
> INTEGER, DIMENSION(NNODE_LS) :: IPIV
> REAL (DP):: xi_, yi_, xj_, yj_, dist_, temp, determ
> REAL (DP):: F_intp
> REAL (DP):: xni, yni, xnj, ynj,xw_ls_norm
> REAL(DP), ALLOCATABLE:: work(:)
>
> ! Determine the coeff of the Radial Basis Functions
> ! This approach is used to interpolate the function everywhere in the domain
>
>
> ! 1. Invert RBF_LS only at the 1st iteration (iflag==0) , then I will have the RBF coeffs
>
> IF (iflag == 0) THEN
> WRITE(*,*) ' Solve global RBF iter', iter, Nnode_LS !nnode_ls is the number of node in each row.
>
> ! Matrix of the RBF function values
> DO ii = 1, Nnode_LS
> xi_ = xy_lsgrid(1, ii)
> yi_ = xy_lsgrid(2, ii)
> DO jj = 1, Nnode_LS
> xj_ = xy_lsgrid(1, jj)
> yj_ = xy_lsgrid(2, jj)
> dist_ = SQRT((xj_ - xi_)**2.0_dp + (yj_ - yi_)**2.0_dp)
> RBF_LS(ii, jj) = EXP(-((dist_ * epsw)**2.0_dp))
> !IF (RBF_LS(II,JJ) <= 1.0E-58_DP) RBF_LS(II,JJ)=0.0_DP
> END DO
> END DO
>
> ! Solve the linear system
> NRHS=1
> print*, "1:", rbf_ls(1:10,1)
> ! Perform LU decomposition of RBF_LS
> CALL DGETRF(NNODE_LS, NNODE_LS, RBF_LS, NNODE_LS, ipiv, info)
> print*, "2:composed LU", rbf_ls(1:10,1)
> IF (info == 0) THEN
> ! RBF_LS was successfully factorized, now compute its inverse
> ALLOCATE(work(lwork)) !this is work vector
> CALL DGETRI(NNODE_LS, RBF_LS, NNODE_LS, ipiv, work, nnode_ls, info)
> print*, "3:", rbf_ls(1:10,1)
> IF (info == 0) THEN
> ! A_inv contains the inverse of A
> ! Use A_inv as needed
> ELSE
> PRINT *, "Error computing the inverse of A."

> END IF
> DEALLOCATE(work)
> ELSE
> PRINT *, "Error in LU decomposition of A."
> END IF
> > One last remark: MicroSoft Powerstation has been out of support for several decades.
> The point is that this old method is giving a more accurate result in comparison with the new one.
>
How do you know the old method is more accurate? I have heard of cases where people thought a result was right because it agreed with a previous one, but both turned out to be wrong. We have not been shown what either matmul(A_inv ,A) or matmul(A,A_inv) is. Nor have we been told how ill-conditioned A is. I recall years ago that one of my students thought of a way to do a problem in quantum mechanics that involved solving a linear system of the form A*x=y. His pure mathematics was right, but he had never learnt about ill-conditioned matrices. I said that the elements of y were obtained from physical experiments and were known only to within about 1 part in 10**6 and so I asked him to try changing some of them by that amount and see what happened to x. That showed him why his method was not in general use for his problem.

> How do you know the old method is more accurate?

Thank you for sharing your valuable experience with me. As a student, I am always dubious about the results I obtain. So I do the double check all the time. The reason why I think the old method is correct is that I calculate the area of my grid based on the matrix that I obtained later. Then with a simple comparison, I realized that the new code was wrong. Probably I am doing something wrong with Lapack subroutines. I think of a precision problem because when I change the intrinsic function from real64 to real128, the results get improved, not enough to be compared with the old code though. If the problem was rooted in the ill-conditioning of the RBF__LS matrix, it should not have worked properly with the old code as well. Please correct me if I am wrong.
Thank you for taking the time.

Cheers,

Farzad

> I will try this though I already know that the result is not accurate. I actually tried REAL128 from the intrinsic kind function in the modern fortran and the results got closer to the one that I was receiving from the old code. It means the reason might be due to the precision. But aren't Real*8 and Real64 the same?
> How do you know the old method is more accurate?

Thanks for the hint of inverting the inverted matrix. I did it and realized that the result of the mentioned calculation was not accurate. Indeed it has the accuracy until the first four decimals. It means it is using a single precision instead of a double. Nevertheless, I am sure the LAPack subroutines and files I am using start with "D" letter instead of "S". Can anyone give me a clue?

Regards,

Farzad Tatar

Il 04/09/23 15:38, Farzad Tatar ha scritto:
>... Indeed it has the accuracy until the first four decimals. It means it is using a single precision instead of a double. Nevertheless, I am sure the LAPack subroutines and files I am using start with "D" letter instead of "S". Can anyone give me a clue?

Before making uncontrolled assumptions, you should check the condition
number of your matrix (see the definition of condition numebr for
matices in https://en.wikipedia.org/wiki/Condition_number ).

On Monday, September 4, 2023 at 6:38:22 AM UTC-7, Farzad Tatar wrote:

(snip)

> Thanks for the hint of inverting the inverted matrix. I did it and realized that
> the result of the mentioned calculation was not accurate.
> Indeed it has the accuracy until the first four decimals. It means it is using a
> single precision instead of a double.

> Nevertheless, I am sure the LAPack subroutines and files I am using start with "D" letter instead of "S". Can anyone give me a clue?

It is extremely easy to lose precision with matrix inversion. LU decomposition is much better.

Many calculations only need single precision, but many matrix operations need to be done with
double precision intermediates to get single precision.

And if you are a little less lucky, only four digits.

In 1968, IBM added 128 bit floating point hardware to the 360/85, and continued
it with S/370 and later models. (Except the divide operation, DXR, was done in software
emulation until recently, as studies showed that it was rarely enough used.)

And yet most others didn't follow along. VAX has a defined H-float, but for most it
is either software emulation, or extra cost option.

RISC-V has it. I don't know which implementations include it.

> Before making uncontrolled assumptions, you should check the condition
> number of your matrix (see the definition of condition numebr for
> matices in https://en.wikipedia.org/wiki/Condition_number ).

Thank you for your insightful advice. I calculated the condition number of my matrix using the Forbenus norm and here are the norms and the condition number of my RBF matrix.

Norm , Norm of the inverse matrix and Condition number are: 142.049233380836
1.829619706901719E+016 2.598960767438591E+018

Considering the condition number, I may lose the accuracy up to 18 decimals.. Am I right? Please let me know if you have any textbook or reference on how to improve the condition number. I have read the RBF interpolation is well-conditioned, or maybe I am wrong and that has been compactly supported RBF interpolation. I need to delve more into this topic. Thank you anyway for your valuable advice.

Cheers,

F Tatar

> Many calculations only need single precision, but many matrix operations need to be done with
> double precision intermediates to get single precision.
Thank you, I will remember this to avoid unnecessary memory allocation to the operations where matrix inversion is not included.

> In 1968, IBM added 128 bit floating point hardware to the 360/85, and continued
> it with S/370 and later models. (Except the divide operation, DXR, was done in software
> emulation until recently, as studies showed that it was rarely enough used.)

Since I realized that the matrix I was trying to invert was ill-conditioned, I will fix this problem first, then I will tackle the accuracy if it was needed. Thank you anyway for the recourse.
Regards,

Farzad Tatar

Farzad Tatar wrote:

>
> > Before making uncontrolled assumptions, you should check the
> > condition number of your matrix (see the definition of condition
> > numebr for matices in
> > https://en.wikipedia.org/wiki/Condition_number ).
>
> Thank you for your insightful advice. I calculated the condition
> number of my matrix using the Forbenus norm and here are the norms
> and the condition number of my RBF matrix.
>
> Norm , Norm of the inverse matrix and Condition number are:
> 142.049233380836 1.829619706901719E+016 2.598960767438591E+018
>
> Considering the condition number, I may lose the accuracy up to 18
> decimals. Am I right? Please let me know if you have any textbook or
> reference on how to improve the condition number. I have read the RBF
> interpolation is well-conditioned, or maybe I am wrong and that has
> been compactly supported RBF interpolation. I need to delve more into
> this topic. Thank you anyway for your valuable advice.
>
> Cheers,
>
> F Tatar

There are two direct solutions:
(1) reduce the effective correlation that you are assuming, which will
move the matrix in the direction of an identity matrix. Ensure that the
largest off-diagonal element isn't too close to one.
(2) reduce the size of the field over which you are doing
interpolation. That is, for any given target point, limit the number of
neighbours being used for fitting: maybe only 20 rather than 1400.

You might also want to reconsider you choice of basis-function. You
could also look at the problem from the point of view of statistical
optimal interpolation rather than basis-functions: this gives a very
similar algorithm but with a different interpretation, which may be
helpful.

On Tuesday, September 5, 2023 at 12:52:32 PM UTC+2, David Jones wrote:
> Farzad Tatar wrote:
>
> >
> > > Before making uncontrolled assumptions, you should check the
> > > condition number of your matrix (see the definition of condition
> > > numebr for matices in
> > > https://en.wikipedia.org/wiki/Condition_number ).
> >
> > Thank you for your insightful advice. I calculated the condition
> > number of my matrix using the Forbenus norm and here are the norms
> > and the condition number of my RBF matrix.
> >
> > Norm , Norm of the inverse matrix and Condition number are:
> > 142.049233380836 1.829619706901719E+016 2.598960767438591E+018
> >
> > Considering the condition number, I may lose the accuracy up to 18
> > decimals. Am I right? Please let me know if you have any textbook or
> > reference on how to improve the condition number. I have read the RBF
> > interpolation is well-conditioned, or maybe I am wrong and that has
> > been compactly supported RBF interpolation. I need to delve more into
> > this topic. Thank you anyway for your valuable advice.
> >
> > Cheers,
> >
> > F Tatar
> There are two direct solutions:
> (1) reduce the effective correlation that you are assuming, which will
> move the matrix in the direction of an identity matrix. Ensure that the
> largest off-diagonal element isn't too close to one.
> (2) reduce the size of the field over which you are doing
> interpolation. That is, for any given target point, limit the number of
> neighbours being used for fitting: maybe only 20 rather than 1400.
>
> You might also want to reconsider you choice of basis-function. You
> could also look at the problem from the point of view of statistical
> optimal interpolation rather than basis-functions: this gives a very
> similar algorithm but with a different interpretation, which may be
> helpful.

Thank you. I did not totally understand your first point. Regarding your second advice, I have used a similar job, though I cannot change the size of RBF matrix, I have set far values to be zero. However, this did not help. I was expecting a meaningful improvement in reducing the condition number of my problem, but it is still in the order of 10E+13.

I tried a different basis function but I did not improve the condition number.

A solution that I discovered on the internet is to use a Singular Value Decomposition algorithm for the matrix inversion. Do you think it could help me in my case?

Cheers,

Farzad

Farzad Tatar wrote:

> On Tuesday, September 5, 2023 at 12:52:32 PM UTC+2, David Jones wrote:
> > Farzad Tatar wrote:
> >
> > >
> > > > Before making uncontrolled assumptions, you should check the
> > > > condition number of your matrix (see the definition of
> > > > condition numebr for matices in
> > > > https://en.wikipedia.org/wiki/Condition_number ).
> > >
> > > Thank you for your insightful advice. I calculated the condition
> > > number of my matrix using the Forbenus norm and here are the
> > > norms and the condition number of my RBF matrix.
> > >
> > > Norm , Norm of the inverse matrix and Condition number are:
> > > 142.049233380836 1.829619706901719E+016 2.598960767438591E+018
> > >
> > > Considering the condition number, I may lose the accuracy up to
> > > 18 decimals. Am I right? Please let me know if you have any
> > > textbook or reference on how to improve the condition number. I
> > > have read the RBF interpolation is well-conditioned, or maybe I
> > > am wrong and that has been compactly supported RBF interpolation.
> > > I need to delve more into this topic. Thank you anyway for your
> > > valuable advice.
> > >
> > > Cheers,
> > >
> > > F Tatar
> > There are two direct solutions:
> > (1) reduce the effective correlation that you are assuming, which
> > will move the matrix in the direction of an identity matrix. Ensure
> > that the largest off-diagonal element isn't too close to one.
> > (2) reduce the size of the field over which you are doing
> > interpolation. That is, for any given target point, limit the
> > number of neighbours being used for fitting: maybe only 20 rather
> > than 1400.
> >
> > You might also want to reconsider you choice of basis-function. You
> > could also look at the problem from the point of view of
> > statistical optimal interpolation rather than basis-functions: this
> > gives a very similar algorithm but with a different interpretation,
> > which may be helpful.
>
> Thank you. I did not totally understand your first point. Regarding
> your second advice, I have used a similar job, though I cannot change
> the size of RBF matrix, I have set far values to be zero. However,
> this did not help. I was expecting a meaningful improvement in
> reducing the condition number of my problem, but it is still in the
> order of 10E+13.
>
> I tried a different basis function but I did not improve the
> condition number.
>
> A solution that I discovered on the internet is to use a Singular
> Value Decomposition algorithm for the matrix inversion. Do you think
> it could help me in my case?
>
> Cheers,
>
> Farzad

Rather than playing with the "far values", I think you should be
looking at any close pairs to make sure that none of the underlying
locations are not effectively identical. Otherwise you could try
replacing the values on the diagonal of your matrix which should be
being evaluated as one, by one plus a small number. The effect of this
is to change the exact interpolation to an approximate interpolation.
You may find this referred to as a "nugget effect".

As to the singular value decomposition, this may just identify the fact
that you have a difficult matrix. I don't think it would help with
matrix inversion, but there may be some other way to use the theory to
do the interpolation.

You say that you cannot change the size of RBF matrix, but this should
be possible by splitting the interpolation problem into a set of
overlapping sub-regions and doing interpolation for each. This is a
generalised nearest-neighbour approach. Alternatively you could treat
the idea as developing your understanding of what is going-on, by
looking at the results for a 1 by 1 matrix, a 2 by 2 matrix,... etc.

On 9/3/23 3:34 PM, Farzad Tatar wrote:

I think you are on the right track with examination of the condition
number of the matrix. As you already note, the solution to this problem
may involve reformulations using different basis functions, or maybe the
same type of basis functions but different orders of the function, or
different grid spacing, and so on.

Here are a couple of comments about your posted code.

> dist_ = SQRT((xj_ - xi_)**2.0_dp + (yj_ - yi_)**2.0_dp)

You should almost never use a floating point exponent with an integer
value in fortran. You should use an integer exponent here. The floating
point exponent expression is much more expensive (it requires a LOG()
and and EXP() evaluation internally), and it is potentially less accurate.

> ALLOCATE(work(lwork)) !this is work vector

I do not see where lwork is set in your code. Normally with LAPACK
procedures, you make a workspace query call first, and then use the
returned value to allocate the workspace, and then call the procedure to
do the work.

$.02 -Ron Shepard

Ron Shepard wrote:

> On 9/3/23 3:34 PM, Farzad Tatar wrote:
>
> I think you are on the right track with examination of the condition
> number of the matrix. As you already note, the solution to this
> problem may involve reformulations using different basis functions,
> or maybe the same type of basis functions but different orders of the
> function, or different grid spacing, and so on.
>
> Here are a couple of comments about your posted code.
>
> > dist_ = SQRT((xj_ - xi_)**2.0_dp + (yj_ - yi_)**2.0_dp)
>
>
> You should almost never use a floating point exponent with an integer
> value in fortran. You should use an integer exponent here. The
> floating point exponent expression is much more expensive (it
> requires a LOG() and and EXP() evaluation internally), and it is
> potentially less accurate.
>
.. and in the next statement the result is effectively squared, so the
SQRT calculation is unnecessary.

I suggest to have a look at: "Multivariate interpolation of large sets of scattered data", Robert J. Renka

https://dl.acm.org/doi/abs/10.1145/45054.45055

And the implementations: QSHEP2D, QSHEP3D (TOMS 660, TOMS 661)

On Wednesday, September 6, 2023 at 10:52:12 AM UTC+2, Edmondo Giovannozzi wrote:
> I suggest to have a look at: "Multivariate interpolation of large sets of scattered data", Robert J. Renka
Thank you Edmondo,
I will read it for a better understanding as I am not expert in this field.

> > You should almost never use a floating point exponent with an integer
> > value in fortran. You should use an integer exponent here. The
> > floating point exponent expression is much more expensive (it
> > requires a LOG() and and EXP() evaluation internally), and it is
> > potentially less accurate.
Thank you for your advice. Yes, I did a short search and now I see that even though it solves the problem, it is better to not use decimal for that. In fact, I am using Fortran to improve the numerical efficiency of an extensive problem. Therefore, following such useful hints is mandatory for me.

Cheers,

Farzad

On Wednesday, September 6, 2023 at 6:00:47 AM UTC-7, Farzad Tatar wrote:
> > > You should almost never use a floating point exponent with an integer
> > > value in fortran. You should use an integer exponent here. The
> > > floating point exponent expression is much more expensive (it
> > > requires a LOG() and and EXP() evaluation internally), and it is
> > > potentially less accurate.

> Thank you for your advice. Yes, I did a short search and now I see that
> even though it solves the problem, it is better to not use decimal for that.

Note that besides efficiency, the 2.0 exponent won't work with negative values.

> Note that besides efficiency, the 2.0 exponent won't work with negative values.
Thank you for that, I did not know it. Apart from getting the response to your specific question, you will learn many other things, when you ask questions in this group. I like it. Thank you all for the information you are providing.