Efficient convergence check

I have a grid with thousands of double precision reals.

It's iterating through, and I need it to stop when it's reached convergence to 3 decimal places.

The target is to have it run as fast as possible, but needs to give the same result every (to 3 dp) every time.

At the minute I'm doing something like this

REAL(KIND=DP) :: TOL = 0.001_DP

DO WHILE(.NOT. CONVERGED)
    CONVERGED = .TRUE.
    DO I = 1, NUM_POINTS
        NEW POTENTIAL = !blah blah blah
        IF (CONVERGED) THEN
            IF (NEW_POTENTIAL < OLD_POTENTIAL - TOL .OR. NEW_POTENTIAL > OLD_POTENTIAL + TOL) THEN
                CONVERGED = .FALSE.
            END IF
        END IF
        OLD_POTENTIAL = NEW POTENTIAL
    END DO  
END DO

I'm thinking that many IF statements can't be too great for performance. I thought about checking for convergence at the end; finding the average value (summing the whole grid, divide by num_points), and checking if that has converged开发者_JAVA技巧 in the same way as above, but I'm not convinced this will always be accurate.

What is the best way of doing this?

If I understand correctly you've got some kind of time-stepping going on, where you create the values in new_potential by calculations on old_potential. Then make old equal to new and carry on.

You could replace your existing convergence tests with the single statement

converged = all(abs(new_potential - old_potential)<tol)

which might be faster. If the speed of the test is a major concern you could test only every other (or every third or fourth ...) iteration

A few comments:

1) If you used a potential array with 2 planes, instead of an old_ and new_potential, you could transfer new_ into old_ by swapping indices at the end of each iteration. As your code stands there's a lot of data movement going on.

2) While semantically you are right to have a while loop, I'd always use a do loop with a maximum number of iterations, just in case the convergence criterion is never met.

3) In your declaration REAL(KIND=DP) :: TOL = 0.001_DP the specification of DP on the numerical value of TOL is redundant, REAL(KIND=DP) :: TOL = 0.001 is adequate. I'd also make this a parameter, the compiler may be able to optimise its use if it knows that it is immutable.

4) You don't really need to execute CONVERGED = .TRUE. inside the outermost loop, set it before the first iteration -- this will save you a nanosecond or two.

Finally, if your convergence criterion is that every element in the potential array has converged to 3dp then that is what you should test for. It would be relatively easy to construct counterexamples for your suggested averages. However, my concern would be that your system will never converge on every element and that you should be using some matrix norm computation to determine convergence. SO is not the place for a lesson in that topic.

What are the calculations for the convergence criteria? Unless they are worse then the calculations to advance the potential it is probably better to have the IF statement to terminate the loop as soon as possible rather than guess a very large number of iterations to be sure to obtain a good solution.

Re High Performance Mark's suggestion #1, if the copying operation is a significant portion of the run time, you could also use pointers.

The only way to be sure about this stuff is to measure the run time ... Fortran provides intrinsic functions to measure both CPU and clock time. Otherwise you may modify your some portion of you code to make it faster, perhaps making it less easier to understand and possibly introducing a bug, possibly without much improvement in runtime ... if that portion was taking a small amount of the total runtime, no amount of cleverness will can make much difference.

As High Performance Mark says, though the current semantics are elegant, you probably want to guard against an infinite loop. One approach:

PotentialLoop: do i=1, MaxIter

  blah

  Converged = test...
  if (Converged) exit PotentialLoop

  blah

end do PotentialLoop

if (.NOT. Converged) write (*, *) "error, did not converge"

I = 1
DO
  NEWPOT = !bla bla bla
  IF (ABS(NEWPOT-OLDPOT).LT.TOL) EXIT
  OLDPOT = NEWPOT
  I = MOD(I,NUMPOINTS) + 1
END DO

Maybe better

I = 1
DO
  NEWPOT = !bla bla bla
  IF (ABS(NEWPOT-OLDPOT).LT.TOL) EXIT
  OLDPOT = NEWPOT
  IF (I.EQ.NUMPOINTS) THEN
    I = 1
  ELSE
    I = I + 1
  END IF
END DO