Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- Subroutine to solve system of equation with the accuracy better than double precision

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

GoldenRetriever

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

02-16-2014
02:22 AM

76 Views

Subroutine to solve system of equation with the accuracy better than double precision

I'm currently using the GESV subroutine from MKL library to solve the system of equations (A.x=B). The number I used is Double Precision.

I see the results are accurate up to 11 digits after the decimal point. Now I would like to have the accuracy is up to 16 digits but I can't find the solution yet.

Anyone suggests a solution for that would be highly appreciated.

Thanks a heap!

Link Copied

12 Replies

TimP

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

02-16-2014
03:37 AM

76 Views

FortranFan

Honored Contributor I

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

02-17-2014
02:12 PM

76 Views

GoldenRetriever wrote:

I'm currently using the GESV subroutine from MKL library to solve the system of equations (A.x=B). The number I used is Double Precision.

I see the results are accurate up to 11 digits after the decimal point. Now I would like to have the accuracy is up to 16 digits but I can't find the solution yet.

Anyone suggests a solution for that would be highly appreciated.

Thanks a heap!

It seems to me what you're trying to achieve, "have the accuracy .. up to 16 digits ", would only be possible under a very specific set of conditions for a system of equations, non-linear I assume. You would need to study your A matrix thoroughly, analyze its determinant, its eigenvalues, etc. and look at transformations to minimize any loss of precision during decompostion. And then you can make use of QUAD precision intelligently to improve accuracy. To me, this would be akin to developing a "custom solver" for a given problem set, a departure from using a general-purpose solver from MKL library.

But should you figure out how to do this with a library solver like in MKL, please post back on this forum if you can do so. I, for one, will be very keen to learn from your experience and would greatly appreciate your contribution!

Cheers,

John_Campbell

New Contributor II

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

02-17-2014
03:34 PM

76 Views

Potentially you can itterate on the error, by:

calculate x for A.x = B

calculate the error e = B - A.x

calculate x' for A.x' = e

Potentially then x+x' is a better solution.

You should investigate why A.x-B is only accurate to 10 significant figures, when up to 16 could be possible with real(8)

Often due to rounding during the calculation, the above approach does not improve the result.

You can try using real(16) (especially for the calculation of e using a real(16) accumulator), but this will only work if you can calculate the coefficients of A and B to real(16) precision. If you populate A and B as real(16) arrays with values calculated as real(8) you can achieve only limited improvement in the "accuracy" of the result, after consideration of the accuracy of the values used in A and B. This is due to the sensitivity of the error e to the inaccuracy of the coeffcicents of A and B.

Poor precision is probably due to an "ill-conditioned" A, which can only be corrected by changing the problem definition.

John

Steven_L_Intel1

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

02-17-2014
06:34 PM

76 Views

GoldenRetriever

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

02-19-2014
02:39 AM

76 Views

Thank you all for the reply. Appreciated that.

FortranFan: at this stage I prefer to use any available general-purpose solver such as from MKL library. To develop a kind of "custom solver" would be a bit out of my ability at the moment. Certainly will let you know if I find a way to get around with this.

John Campbell: what you suggested would be interesting to try out. I will have a go with it.

Tim Prince & Steve Lionel (Intel): Using double precision, my arrays A and B are accurate up to 15 digits; however, the results give out from "GESV" are only 11 digits of accuracy. I include a simple example of my case in the attached file for reference. In this example, the results are the displacements of a structure, they are supposed to be symmetric but some values are not accurately symmetric as expected.

I notice MKL has an Expert Driver named gesvx which provides error bounds on the solution. I am not sure if it is something can help in this case. Anyone has experience using that please advise.

GoldenRetriever

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

02-19-2014
02:43 AM

76 Views

Thank you all for the reply. Appreciated that.

FortranFan: at this stage I prefer to use any available general-purpose solver such as from MKL library. To develop a kind of "custom solver" would be a bit out of my ability at the moment. Certainly will let you know if I find a way to get around with this.

John Campbell: what you suggested would be interesting to try out. I will have a go with it.

Tim Prince & Steve Lionel (Intel): Using double precision, my arrays A and B are accurate up to 15 digits; however, the results give out from "GESV" are only 11 digits of accuracy. I include a simple example of my case in the attached file for reference. In this example, the results are the displacements of a structure, they are supposed to be symmetric but some values are not accurately symmetric as expected.

I notice MKL has an Expert Driver named gesvx which provides error bounds on the solution. I am not sure if it is something can help in this case. Anyone has experience using that please advise.

Steven_L_Intel1

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

02-19-2014
09:08 AM

76 Views

I am moving this thread to the MKL forum.

John_Campbell

New Contributor II

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-02-2014
08:23 PM

76 Views

Did you try what I suggested on solving the error ?

Often this does not work, but in some circumstances it can help.

The basic approach is:

given B = A . x

a first solution estimate is x1

this gives e = B - A . x1

solving e = A . xe

an improved solution is B = e + A . x1 = A . ( x1 + xe)

and e2 = B - A . (x1 + xe)

You can check the accuracy by either the maximum absolure error or the sum of the square error terms. ( e'.e )

Calculating e and e2 to higher precision can help.

I'd be interested if it proved effective or ineffective, if you had a chance to test it out.

John

GoldenRetriever

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-04-2014
04:34 AM

76 Views

Thanks John for your follow up.

I could not get it work even your method does sound make sense to me.

I include a simple example with array of 16x16 in attached file if you interest to have a try with it. Please let me know if you can spot any problem with what I have done.

Thanks very much mate.

John_Campbell

New Contributor II

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-04-2014
08:40 PM

76 Views

GoldenRetriever,

I removed all the libraries and replaced the solver with a real*16 solver.

I also placed code to check the size of the coefficients against the size of the error terms.

After cycle 1, the max coefficient term when calculating the error matrix is 6.45052935008452D+9

The max error is 7.666822057217360E-007, which implies an accuracy of about 16 significant figures.

The determinant of the AA matrix is 2.085791801338600E+123.

This AA matrix appears to be a well conditioned matrix and there is minimal improvement in error after 2 cycles. There is little opportunity to improve or demonstrate improvement of the error.

I think that the basic problem you have is the error you are looking for an error magnitude of 1.e-10, where you should be looking for an error in terms of significant figures of error. For the changes I made, I am getting 16 significant figures, while I expect your original code provided something similar.

John

For some reason the attachment did not work. I will paste the code I changed below. Excuse the expected layout. Will this IDZ ever get fixed !!

[fortran]! USE MKL95_LAPACK

! USE mkl95_blas

IMPLICIT NONE

INTEGER*4, parameter :: N = 16 ! size of problem

integer*4, parameter :: LIMIT = 11 ! number of itterations

!

INTEGER I, J, k, ERROR_CODE

real*8 A(N,N), B(N,1), AA(N,N), X(N,1) , XX(N,LIMIT), E(N,1)

real*8 ErrorMax, s, determinant, v

!

real*16 acum

!

INTERFACE

SUBROUTINE MATRIX_SOLVE (MATRIX_IN, RHS, determinant, ERROR_CODE)

!

real*8, dimension(:,:), intent (in) :: matrix_in

real*8, dimension(:,:), intent (inout) :: rhs

real*8, intent (out) :: determinant

integer*4, intent (out) :: error_code

!

END SUBROUTINE MATRIX_SOLVE

END INTERFACE

!

! BEGIN Read Input File

OPEN(111,FILE='TDReview_FEM-K_REDUCE.csv',STATUS='UNKNOWN', ACTION='READ')

REWIND(111)

READ(111,*) ((A(i,j), j=1,N), i=1,N)

CLOSE(111)

OPEN(112,FILE='TDReview_FEM-LOAD_REDUCE.csv',STATUS='UNKNOWN', ACTION='READ')

REWIND(112)

READ(112,*) ((B(i,j), j=1,1), i=1,N)

CLOSE(112)

! END Read Input File

! check A for symmetry

s = 0

v = 0

do i = 1,n

do j = i+1,n

s = max (s, abs(a(i,j)-a(j,i)) )

v = max (v, abs(a(i,j)) )

end do

end do

write (*,*) 'Max symmetry error =',s

write (*,*) 'Max coefficient =',v

!

XX=0D0

DO i=1,LIMIT

!

AA = A ! Reasign AA which is modified when calling GESV

!

! X(:,1) = SUM(XX,DIM=2)

do j = 1,n

acum = 0

do k = 1,i

acum = acum + xx(j,k)

end do

X(j,1) = acum

end do

!

! E = B - MATMUL(A,X)

v = 0

do j = 1,n

acum = B(j,1)

do k = 1,n

acum = acum - A(j,k)*X(k,1)

v = max ( v, abs(A(j,k)*X(k,1)) )

end do

E(j,1) = acum

end do

ErrorMax = MAXVAL(ABS(E))

s = 0

do k = 1,N

s = s + ABS (E(k,1))

end do

WRITE (*,*) 'LOOP',I

write (*,*) ' Max error term :', ErrorMax

write (*,*) ' Avg error term :', s/n

write (*,*) ' Max coefficient :', v

IF (ErrorMax<1.D-10) EXIT

!

! CALL GESV (AA, E)

call MATRIX_SOLVE (AA, E, determinant, ERROR_CODE)

write (*,*) ' Solve error =',error_code

write (*,*) ' determinant =',determinant

XX(:,i) = E(:,1)

END DO

! OPEN(111,FILE='TDReview_FEM-DISP_VEC.csv',STATUS='UNKNOWN', ACTION='WRITE')

DO i=1,N

WRITE (*,"(1000(1x,F55.19,','))") X(i,:)

END DO

END[/fortran]

GoldenRetriever

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-05-2014
02:51 AM

76 Views

Thanks John for that,

I could not see the mentioned "real*16 solver" in the code. I could not reproduced your code to make it works as well. Did I miss anything?

PS: the procedure to attach file on this forum is a bit different from the normal way of attaching file as most email webpage does. It is a bit confusing; after choosing to add the file we actually need one more step to upload them to get it TRULY added to the post. 75% of the time I failed to attach the file; sadly that's true. Should it be something Intel Forum to fix? Yes. ;)

John_Campbell

New Contributor II

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-05-2014
09:57 PM

76 Views

Attached is the hybrid version which mixes real*8 and real*16 for the equation solution. If you want a fully real*16 solution, change most real*8 to real*16. For a fully real*8 solution, change real*16 to real*8 in the solver, while keeping acum as real*16.

You should see that the error does not change significantly and given the maximum value of the coefficients and load case, the errors reported are quite acceptable.

I adapted this code to inspect the errors during the calculation. MKL does provide better routines for refining the error, but I think they would also find that as the matrix is well conditioned, there is not much scope for improvement.

John

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

For more complete information about compiler optimizations, see our Optimization Notice.