topic GoldenRetriever, in Intel® oneAPI Math Kernel Library

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

GoldenRetriever — Sun, 16 Feb 2014 10:22:16 GMT

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!

If you're looking for a

TimP — Sun, 16 Feb 2014 11:37:16 GMT

If you're looking for a solution within MKL, the MKL forum would be a more likely choice to get you an answer. Otherwise, you could look into trying the iterative refinement by building the iterative refinement method using a combination of double real(8) and quad real(16) data types.

Quote:GoldenRetriever wrote:

FortranFan — Mon, 17 Feb 2014 22:12:00 GMT

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,

Potentially you can itterate

John_Campbell — Mon, 17 Feb 2014 23:34:35 GMT

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

16 digits would be especially

Steven_L_Intel1 — Tue, 18 Feb 2014 02:34:23 GMT

16 digits would be especially difficult if you're using double precision, since the type itself is good to only about 15 digits. 11 digits for a complex calculation is about right.

Thank you all for the reply.

GoldenRetriever — Wed, 19 Feb 2014 10:39:59 GMT

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.

https://www.dropbox.com/s/m9snqr3cddkx1g2/GESV.ZIP

Thank you all for the reply.

GoldenRetriever — Wed, 19 Feb 2014 10:43:31 GMT

Thank you all for the reply. Appreciated that.

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

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.

https://www.dropbox.com/s/m9snqr3cddkx1g2/GESV.ZIP

I am moving this thread to

Steven_L_Intel1 — Wed, 19 Feb 2014 17:08:01 GMT

I am moving this thread to the MKL forum.

Did you try what I suggested

John_Campbell — Mon, 03 Mar 2014 04:23:44 GMT

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

Thanks John for your follow

GoldenRetriever — Tue, 04 Mar 2014 12:34:54 GMT

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.

GoldenRetriever,

John_Campbell — Wed, 05 Mar 2014 04:40:00 GMT

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]

Thanks John for that,

GoldenRetriever — Wed, 05 Mar 2014 10:51:49 GMT

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. ;)

Attached is the hybrid

John_Campbell — Thu, 06 Mar 2014 05:57:33 GMT

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