IMSL library: Crash in function surface_fitting

[Daylis__Samuil]

Have a big issue with surface_fitting routine: It doesn't accept fine grid of spline breakpoints. I attached simplest possible project to demonstrate this issue. IF 11.0 Project contains file surface_fitting_ex1.f90 that ships with IMSL. In this code Integer parameter ngrid specify number of knots in a grid of breakpoints. Code ships with ngrid=9 on 2X2 square grid. 

Code works if ngrid<15. At ngrid>=15 surface_fitting routine throws an unhandled exception. Sure this value is too low to cause any computational problems but still it does. I tried IF 13.3 with updated version of IMSL6 but the only benefit was that surface_fitting started crashing from ngrid=16 and up. If there is any solution to this problem please let me know.

[Steven_L_Intel1]

We'll take a look - thanks.

[mecej4]

I found that your program runs with no errors until ngrid=40 when compiled with the 32-bit Ifort 13.1 compiler. At ngrid=50 the program failed with stack overflow, rather than an unhandled exception as you stated.

I suggest that you consider different ways of allocating the two-dimensional array coeff and other large arrays in the subroutine rather than have them allocated on the stack, or increase the stack limit when linking. One has to look into such issues when adopting an example code from the IMSL manual for use with a large scale computation.

[John_Campbell]

I was puzzled by your code example for:

[fortran]
      real(kind(1d0)), target :: spline_data (4, ndata),   & 
...
! example exponential function at these values.
          spline_data(1:2,:) =  two*rand(spline_data(1:2,:))
[/fortran]

Does rand support this undefined argument ?

[Daylis__Samuil]

To make things more clear added very simple test. All data (25 points) explicitely  entered , all memory is allocated explicitely. Visual Studio 2008, IF 11.0.066. Same story, but stack overflow at ngrid>=16.

BTW, tried VS 2012 + IF XE 13.1 update 3. Again,  stack overflow at ngrid>=16.

[Steven_L_Intel1]

Code inside the IMSL routine is using a stack temporary. You can get around this by increasing the linker stack limit (Linker > System > Stack Reserve Size). I started with 100000000 and it worked fine with ngrid=16 - that's 100 times the default.

[Daylis__Samuil]

Thanks Steve. This stack size fits even ngrid=50. How could it be that documentation keeps mum about this great "feature"? Or maybe not a feature but a development mistake?

[Steven_L_Intel1]

It's not a development mistake, though it can be aggravating. Stack overflow is a generic problem. You might be interested in an article I wrote a while back on the topic: Doctor Fortran - Don't Blow Your Stack!

That said, I'm going to recommend to Rogue Wave that they build IMSL with the /heap-arrays option so that this sort of problem will be lessened when using IMSL.

[Bernard]

As Steve pointed it out try to increase default values for /STACK reserve and commit options.The default option is 1MB and could be the reason of stack overflow in your case.

P.s

Had such a issue few days ago when allocated large structures of arrays for my multithreaded programme.Resolved it by increasing default /STACK  VALUES.

[John_Campbell]

Would the following change improve things ?

[fortran]
      integer, parameter :: ngrid   = 15,               &
                            nord    = 4,                &
                            ndegree = nord-1,           & ! 3
                            ncoef   = ngrid+ndegree-1,  & ! 17
                            nbkpt   = ngrid+2*ndegree,  & ! 21
                            ndata   = 25
      real(kind(1d0)), parameter   :: zero=0d0, one=1d0, two=2d0
      real(kind(1d0)), target      :: bkpt(nbkpt), delta
      real(kind(1d0)) ,allocatable :: inputx(:),inputy(:),inputf(:)
      real(kind(1d0)) ,allocatable :: coeff(:,:), spline_data(:,:)
      real(kind(1d0)) :: val,xa,ya      
      type (d_spline_knots) knotsx, knotsy
!
      allocate ( spline_data(4,ndata) )
      allocate ( coeff(ncoef,ncoef) )
      allocate ( inputx(ndata), inputy(ndata), inputf(ndata) )
[/fortran]

I always find thet the solution to stack overflow is to not increase the stack size, but fix the problem..

John

[Daylis__Samuil]

Evidently, this is not user data problem, this is this routine design/build problem. See, all input data is allocated explicitely except those two structures knotsx, knotsy. Sizes of arrays passing via those structures are very small, 1 MB stack should suffice... but it is not.

[Bernard]

>>>I always find thet the solution to stack overflow is to not increase the stack size, but fix the problem..>>>

In my case increasing the stack size worked.I did not have stack overflow caused by deep recursive function calls.

[Steven_L_Intel1]

In a case where the error is happening in code you don't control, the only available option is to increase the stack size in the linker properties (or EDITBIN the EXE). If you don't control the EXE either, such as your code is a DLL called from VB.NET, for example, you're out of luck,

[Bernard]

Completely agree Steve.In may case the overflow was caused by large allocation of AoS on stack.So increasing stack size helped.