Page 1 of 2 12 LastLast
Results 1 to 10 of 17

Thread: Ivpag

  1. #1

    Ivpag

    test tetst
    Last edited by prashantsondkar; 10-14-2011 at 01:25 PM.

  2. #2
    Senior Member
    Join Date
    Sep 2005
    Location
    Albuquerque, NM USA
    Posts
    140
    There should be no issue with a problem of this size. Write out all the values in PARAM(1:20) before the call.
    Dickie T. Bird

  3. #3
    Senior Member
    Join Date
    Sep 2005
    Location
    Albuquerque, NM USA
    Posts
    140
    Yes, the routine has the internal work arrays allocated with the SAVE attribute. So call with IDO=3 to release the work store. This is important if you continue on with another model for example. You could also bring the defns of PARAM outside the loop and cut the size of PARAM(4) back to a smaller value. I would leave it at the default unless you get an error message that the integrator took too many steps. You have repeated the comment in defns of PARAM(12) and PARAM(13). Look at the model example 3 in the documentation for the integrator, FNL version 7. Follow the same steps except use your own model.
    Dickie T. Bird

  4. #4
    Senior Member
    Join Date
    Sep 2005
    Location
    Albuquerque, NM USA
    Posts
    140
    Attach your driver, perhaps leaving out the routines f1, jac1. The short story is that storage is released after making your integration steps and before starting again with IGO=1. There is no need to release the storage during the intervening steps. Perhaps it will clear up when your code is made available. A full running driver would be best, of course. You can send this with a private message, attaching the file.
    Dickie T. Bird

  5. #5
    Senior Member
    Join Date
    Sep 2005
    Location
    Albuquerque, NM USA
    Posts
    140
    Do 199 K = 1, NPTS
    Ttm(K) = tout
    tout = tout + DT
    IF(K .eq. NPTS) IDO=3 !<= New line


    This new line is important here. You are calling the integrator several times in your outer loop. Each time you finish an inner set of calls, IDO=3 releases storage. Your code was not executed so there may still be issues.
    Dickie T. Bird

  6. #6
    Senior Member
    Join Date
    Sep 2005
    Location
    Albuquerque, NM USA
    Posts
    140
    That new line should come just before the call to DIVPAG. The way you wrote this it is hard to tell. You had best get in touch with post-sales support.
    Dickie T. Bird

  7. #7
    Senior Member mecej4's Avatar
    Join Date
    Dec 2009
    Posts
    139
    The answer to your memory problems is quite simple.

    In the subroutiness f1() and jac1(), you allocate lots of arrays, which are never deallocated, unless you are using features of Fortran later than Fortran 95 that perform automatic deallocation of locals and your compiler supports such usage.

    Since these subroutines are called thousands of times by the integrator, there are thousands of allocations, which keep piling up. Thus, your program is tying up memory, most of which has no further use (a.k.a. memory leak). Eventually, your computer will run out of memory.

    I question the need for these arrays to be allocatable. If they are needed to be such, and the compiler does not automatically deallocate them, it is your responsibility to deallocate them after the need to retain them is gone.

    You also have a number of bugs in your code: uninitialized variables used in expressions, unused subroutines and lots of unused variables and statements. It is a mess. Until it is cleaned up, there is little point in discussing the results, which at this point are probably quite incorrect because of the bugs.

    Given that your code is buggy, why do you jump to the conclusion that any run time errors must be attributable to inadequacies in IMSL? Debug your code first.
    Last edited by mecej4; 09-30-2011 at 09:45 AM.

  8. #8
    Senior Member mecej4's Avatar
    Join Date
    Dec 2009
    Posts
    139
    But I have one more question, is there way to speed up computation with this routine. Like parallel processing. MATLAB seem to be salving same problem in really less time (order or 10 times fast).
    That is an invalid comparison. Your Fortran code solves ODEs using time-step control based on accuracy. Unfortunately, the same code uses uninitialized variables in several places. Therefore, the time-steps needed to match your excessively tight accuracy requirement are unpredictable, as is the CPU time for the computation.

    Before considering parallelization, remove the bugs in your code, and compute quantities only once if they do not change (e.g., mass matrix, perhaps?).

  9. #9
    Senior Member mecej4's Avatar
    Join Date
    Dec 2009
    Posts
    139
    I need to compute mass matrix couple of times
    You compute not only the mass matrix but, worse, its inverse, every time your subroutine F1() is called by the integrator. That is several thousands of times, which does not qualify as "a couple".

    Make the inverse mass matrix EMIN a SAVE, ALLOCATABLE array in Subroutine F1. Allocate and compute it only during the first run through F1. This change by itself can give you an order of magnitude increase in speed. If you can do likewise for the stiffness matrix, you could pick up another order of magnitude in speed. I estimate that your code, after making these changes, should run in about 10 seconds on a current X86/X64 CPU.

    Secondly, your ODEs are probably not stiff. If so, it would be more efficient to use a non-stiff integrator such as DIVPRK, which has the additional advantage of not needing to compute the jacobian.

    Get rid of unnecessary COMMON blocks. You do not need them and they make complex code difficult to optimize and even more difficult to debug.

    An entire matrix (or any entire array of ranks 1,2,3,...) passed as an argument consumes just as much stack as a scalar: one address. It is local arrays that can consume stack space; to alleviate this problem, some compilers allow you to specify that local arrays beyond a certain size be allocated on the heap.
    Last edited by mecej4; 10-01-2011 at 08:08 PM.

  10. #10
    Senior Member mecej4's Avatar
    Join Date
    Dec 2009
    Posts
    139
    Code:
    double precision, allocatable, save :: EMin(:,:)
    ...
    LOGICAL,SAVE :: INITIAL = .true.
    ..
    IF (INITIAL) THEN
       --- compute inverse or L-U factors of mass matrix into array EMIN ---
       INITIAL=.false.
    ENDIF

Tags for this Thread

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •