Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.

Maximum number of steps exceeded DIVPAG

Divya_K_
Beginner
981 Views

Hello,

I am running a source code using the Intel® Parallel Studio XE Composer Edition for Fortran with Rogue Wave IMSL 7.0 for Windows*

I am using a solver called the IVPAG and it returns an error-"Maximum number of iterations exceeded"

I have investigated the IMSL manual and documentation to obtain the parameter that will allow me to increase the number of steps, but there is no such parameter available to edit. Several other solvers allow to modify this parameter , but I do not see any option to change this parameter in IVPAG.

How can I overcome this problem? 

Thanks!

0 Kudos
14 Replies
mecej4
Honored Contributor III
981 Views

You may find it worthwhile to use DAESL instead of IVPAG; see the comments in #2 of https://forums.roguewave.com/showthread.php?1542-Ivpag .

A pertinent question: are the ODEs that you are integrating stiff?

0 Kudos
Divya_K_
Beginner
981 Views

Hello. Thank you for your response. Yes, they are stiff. However, I have had good success with using IVPAG in terms of the quality of the solution with the other input files. So, ideally I would have liked to continue using IVPAG if that is possible. 

Typically, when I have come across this error, I was able to solve it using a trial and error combination of step size and error criteria. However, it doesn't work this time and I had to see if there is any other parameter that I can change toward fixing this error? For instance, I've seen other solvers that were able to set the parameter for the 

I could try using the new solver you suggested and see how it compares to the previous results. I will keep things posted in the forum. Thank you for your suggestions.

0 Kudos
jimdempseyatthecove
Honored Contributor III
981 Views

Would it be worth showing us your user supplied function?

Jim Dempsey

0 Kudos
Divya_K_
Beginner
981 Views

Hi Jim,

What do you hope to find out from the user supplied function? The only reason I am asking this is because as you may have already guessed, there are several lines of code and some of the parameters have been passed on from another subroutine where the IVPAG is called. If I know why you may want to look at the user supplied function, I can figure out which sections of the code to show here? 

Thanks!

0 Kudos
mecej4
Honored Contributor III
981 Views

It is characteristic of many algorithms that they work well on most problems, but fail on a small set of pathological problems. In such cases, a different solver may succeed with no special effort required. All that we are interested in is to obtain from you the exact specification (ODE, initial conditions, range of independent variable) of one pathological problem instance.

When the ODE contains some parameters, we need just the values of the parameters, even if some  long calculation had to be done in order to obtain those values.

0 Kudos
Divya_K_
Beginner
981 Views

The first call is made to DIVPAG as follows inside a DO and CONTINUE loop. 

    DO 1000 K = 1,NSTEPS

        TOUT=TOUT+SSIZE         ------SSIZE is the inputted step size, in this case 0.4
        TIME(K)=TOUT   -----------TOUT=0 before it enters the DO loop
        TIMEMIN(K)=TOUT/60.0

      CALL DIVPAG (IDO,N,DIFFUN,FCNJ,A,T0,TOUT,EPS,PARAM,YO)

These parameters are:

IDO=1

N=18

A=A(1,1)  //Not very important I believe. It uses default value. 

T0=0

TOUT=0

EPS=0.0000000000000001     

PARAM is defined as 

PARAM=0.0D0
      PARAM(4) = 500000
      PARAM(6) = 5
      PARAM(10)=MABSE  //This value is 1
      PARAM(12)=MBDF   //This value is 2

      PARAM(13)=MSOLVE  //This value is 2
      PARAM(19) = 0 
      IDO=1    

User supplied function=DIFFUN. DIFFUN is a separate subroutine and has parameters as

SUBROUTINE DIFFUN(N,T,YO,YDOT)

YO is the matrix that contains the initial values for the ODEs and the new values at t+delta t are stored in YDOT. Then, it continues until convergence.

I do want to mention that in some of the cases, the program never seems to come back out of the solver. Normally, I would expect that it either gives the solution or "maximum number of iterations exceeded" error. However, there are times when neither of these show up. And it seems stuck in the solver. The reason I know this is because I have a print statement after the program comes out of the solver. I wonder why this happens.

 

             

0 Kudos
mecej4
Honored Contributor III
981 Views

There is something seriously wrong in the way that you use IVPAG, or I am misreading what you wrote in #7.

  • You are integrating the ODE from 0 to SSIZE, then from 0 to 2*SSIZE, ..., 0 to n.SSIZE. This is very wasteful and full of needless repetition of calculations. You should, instead, integrate from 0 to SSIZE, then from SSIZE to 2.SSIZE, ..., (n-1)SSIZE to n.SSIZE. In other words, in each iteration of the DO loop, do not restart the integration from t = 0, but continue from the value of t reached in the previous iteration of the DO loop.
  • Your description of the function DIFFUN is wrong. According to the description in the IMSL manual for IVPAG, DIFFUN is supposed to return the derivatives y_dot, given t and y. As far as this subroutine is concerned, it is irrelevant to talk about "initial" and "final" values. Given the differential equation set dy/dt = f(t,y), this subroutine is expected to return the vector of values of f(t,y) corresponding to the input values t (scalar) and y (vector).
  • Unless your ODE system is implicit, you do not need to worry about the subroutine FCNJ for the Jacobian.

Please study the example ODE problems and codes in the IMSL documentation and construct your program in a similar manner. Many IMSL subroutines have rather complex interfaces, and care is needed to use them correctly.

In hindsight, had a means been provided in IVPAG to increase the number of iterations, you could have found yourself in a worse situation.

It would be useful if you gave a statement of the mathematical problem at this juncture. That would set the stage for how to proceed with obtaining a numerical solution.

0 Kudos
Divya_K_
Beginner
981 Views

I will provide two separate comments. In this comment, I will provide responses to the issues you have mentioned in your comment. In the second comment, I will describe the problem we are trying to solve.

  • Perhaps, I should have mentioned that t0 is re-initialized at every time step. So, that the integration is happening in the way you have suggested it should. The integration begins from 0 to ssize, ssize to 2*ssize, 2*ssize to 3*ssize....and so on. I have confirmed this in the program by providing a print statement every time the DIVPAG is called

 

  • I believe I have provided DIFFUN as is suggested in the manual. See here. http://docs.roguewave.com/imsl/fortran/7.1/html/fnlmath/index.html#page/FNLMath/mch5.08.05.html  DIFFUN is FCN which is to be provided in the form of FCN(N,T,Y,YPRIME). Those are the parameters provided by DIFFUN in the same order. In the same link that I have provided, if you scroll down you will see an Euler equation based problem in which YPRIME is generated i.e the equations are essentially provided using the values in Y. In my code, DIFFUN calls another subroutine ODEQUATN that performs this action of generating equations. So, I am not sure what it is you were pointing.

 

  • We do not have FCNJ in our code.

 

0 Kudos
Divya_K_
Beginner
981 Views

The problem we are trying to solve is as follows-

We have a chemical system. This chemical system has different species namely C1,C2,C3, and C4. For now, we will use 4 species but we end up having around 150 species in our system.  They interact with one another. What I mean by this is simply that,

C1+C2 -> C3    k1=rate constant

C2+C4 -> C1    k2=rate constant  and so on. Therefore we take into account their interaction which changes the nature of the system. So, if we started with 10 of C1, after a certain time, t, we would have 6 of C1.

The goal is to quantify the values of C1,C2,C3 and C4 after a "time" desired by the user. We are able to describe their interaction by a parameter called "rate of reaction" given by dC1/dt = -k1*C1*C2 + k2*C2*C4. So that, here C1,C2,C3 and C4 are independent variables and t is the dependent variable.

Similarly, we write such simultaneous ODEs and solve for each of these variables. They are stiff likely because their values are in a wide range (I believe this is the reason for their stiffness but it could be anything else). For instance, C1 could be 0.001 and C3 could be 10^-11.

It is to solve these simultaneous ODEs that we use DIVPAG. 

0 Kudos
mecej4
Honored Contributor III
981 Views

Ah, I see.

In #7 you wrote:

     YO is the matrix that contains the initial values for the ODEs and the new values at t+delta t are stored in YDOT.

That statement is what caused me alarm. Likewise, I saw no updating of T0 in the code fragment in #7.

Simultaneous chemical reactions are a well-known source of stiff ODEs. Robertson's example is famous. For comprehensive coverage of algorithms, software and comparison of results, see Francesca Mazzia's Web page at https://archimede.dm.uniba.it/~testset/testsetivpsolvers/ and Hairer's Web page at https://www.unige.ch/~hairer/testset/testset.html .

Since the characteristic response times of the simultaneous reactions can span many orders of magnitude, using a uniform time step is not a good choice. See, for example, the time axes in the figures at https://archimede.dm.uniba.it/~testset/problems/rober.php .

0 Kudos
Divya_K_
Beginner
981 Views

Thankyou for the useful links.

I am not sure how to choose the right step size at every iteration, especially when solving computational problems. While solving it by hand, one can make a guess to achieve the solution because we know the values at each iteration. Of course, these were really small problems and therefore, less iterations, so it's feasible.

However, I am not sure how I can do this in a computer program, especially, if it has 1000+ iterations, so it is obviously not practical to look at the value after every iteration and provide a new step size. A uniform step size mostly works or has worked so far for all our systems. Any thoughts about this?

And again, are we able to go back and re-think on modifying the parameter that can change the maximum number of iterations?

 

0 Kudos
mecej4
Honored Contributor III
981 Views

Divya K. wrote:

And again, are we able to go back and re-think on modifying the parameter that can change the maximum number of iterations?

Only Roguewave can answer that.

0 Kudos
jimdempseyatthecove
Honored Contributor III
981 Views

>>However, I am not sure how I can do this in a computer program, especially, if it has 1000+ iterations, so it is obviously not practical to look at the value after every iteration and provide a new step size.

Why not?

Assume you have a candidate guess at useable step size. You can partition the guess into 2 or 4 sizes that bound the guess and run those in parallel. Upon return of all partitions determine a new guess.

Please note, if the number of cores/HT's are underutilized when running without partitioning the step size guess, then the extra work run with partitioning may come mostly free. And, by improving the step size estimations, you  can reduce the total number of steps, thus potentially reducing runtime.

Jim Dempsey

0 Kudos
Castor1
Beginner
981 Views

Hi Divya K.,

You can change the number of steps at the array "param(4)".  The standard value is 500. I usually put something very high, like 2d9. I believe this is the limit value.   Also, another thing you should check it is the tolerance. If I were you, I would print the time variable in the FCN (DIFFUN in your case) subroutine in order to have an idea about the step size.

Hi Jim Dempsey,

Related to you last suggestion, it is very interesting. But, I'm not following you. :-(     How can I do the partition in the guess step?

Thanks

 

0 Kudos
Reply