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

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

Bild des Benutzers GoldenRetriever

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!

  

13 Beiträge / 0 neu
Letzter Beitrag
Nähere Informationen zur Compiler-Optimierung finden Sie in unserem Optimierungshinweis.
Bild des Benutzers Tim Prince

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.

Bild des Benutzers FortranFan

Zitat:

GoldenRetriever schrieb:

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,

Bild des Benutzers John Campbell

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

Bild des Benutzers Steve Lionel (Intel)

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.

Steve
Bild des Benutzers GoldenRetriever

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

Bild des Benutzers GoldenRetriever

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

Anlagen: 

AnhangGröße
Herunterladen GESV_0.ZIP2.34 KB
Bild des Benutzers Steve Lionel (Intel)

I am moving this thread to the MKL forum.

Steve
Bild des Benutzers John Campbell

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

Bild des Benutzers GoldenRetriever

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.

 

 

Anlagen: 

AnhangGröße
Herunterladen GESV_140304.zip2.55 KB
Bild des Benutzers John Campbell

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

!    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

Bild des Benutzers GoldenRetriever

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

Bild des Benutzers John Campbell

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

Anlagen: 

AnhangGröße
Herunterladen gesv03.f906.34 KB
Herunterladen gesv03_0.log3.77 KB

Melden Sie sich an, um einen Kommentar zu hinterlassen.