Hello,
In 4.8 Multidimensional integrals, there is a method to calculate 3D integrals which has been further modified to eliminate the use of global variables as in http://www.nr.com/forum/showthread.php?t=1663 which works perfectly.
I have to calculate 2D integral multiple times (at least 256*256*9*200) in a loop.
As the process continues in the loop,its seems that the program consumes memory in the increasing order(say 3 or 5MB per sec).
So, i get either adrees voilation error or out of memory error even if i am using 3GB ram.
I think that since this is a function and structures the memeory is free up after each call.But it doesnot happens so.

What is the real problem and how can i solve this issue?
If i do it for large number of loops may be i have to do parrallel implementation of code.How can i do that?



#include "qromb.h"
#include <math.h>

//--------------------------------------------------------------------
// After debugging/testing, you can put the template stuff in
// its own header

//--------------------------------------------------------------------
// Start of template stuff
// Derived from the quad3d_t in post number 3 of this thread:
//http://www.nr.com/forum/showthread.php?t=1663
//
template < class T >
struct NRf2
{
Double xsav;
T *func2d;
Double operator() (const Double y) {
return (*func2d) (xsav, y);
}
};

template < class T, class Y1, class Y2 >
struct NRf1 {
Y1 & y1;
Y2 & y2;
NRf2 < T > f2;
NRf1(Y1 & yy1, Y2 & yy2):y1(yy1), y2(yy2) { }
Double operator() (const Double x) {
f2.xsav = x;
return qromb(f2, y1(x), y2(x));
}
};

template < class T, class Y1, class Y2 >
Double quad2d_t(T & func, const Double x1, const Double x2, Y1 & y1, Y2 & y2)
{
NRf1 < T, Y1, Y2 > f1(y1, y2);
f1.f2.func2d = &func;
return qromb(f1, x1, x2);
}
// End of template stuff
//--------------------------------------------------------------------

//Structures for calculating the volume of intersection of sphere and cylinder


struct Ftor
{ double row;
Ftor(const Double & sphr):row(sphr) {}
Double operator() (const Double & x, const Double & y) const
{
double b;
if ((row*row - x*x - y*y)<0)
{
// ShowMessage(row*row - x*x - y*y);
b= sqrt(fabs(row*row - x*x - y*y));
}
else
b= sqrt(row*row - x*x - y*y);


return b;
}
};

// For a given value of x, this calculates
// the lower limit for integration over y.
//
struct Fty1
{
// double row,r,d;

// Fty1(const Double & r):radius(r) {}

Double operator() (const Double & x) const
{
return 0;
}
};


//
//For a given value of x, this calculates
// the upper limit for integration over y.
//

struct Fty2
{
double row,r,d;

Fty2(const Double & cylr,const Double & sphr,const Double & dcsphacyl):row(sphr),r(cylr),d(dcsphacyl) {}

Double operator() (const Double &x) const
{
return GetMin(sqrt(fabs(r*r - (x-d)*(x-d))),sqrt(fabs(row*row - x*x)));
}
};

//Functions for calculating the volume of intersection of sphere and cylinder which is called in a loop
// sphr : the radius of sphere
// cylr : the radius of cylinder
// dcsphacyl : the distance from the center of the sphere to the axis of the cylinder
//reference: F. Lamarche and C. Leroy, Evaluation of the volume of intersection
//of a sphere with a cylinder by elliptic integrals, Computer Phys. Comm. 59 (1990) 359-369.
double Cyl_Sphere_Intersection(double cylr, double sphr,double dcsphacyl){

double Intersection;
sphr=1;cylr=1;dcsphacyl=1;
Ftor func(sphr);
Fty1 lymin;
Fty2 lymax(cylr,sphr,dcsphacyl);
Intersection=4*quad2d_t(func, GetMax(dcsphacyl-cylr,-sphr), GetMin(dcsphacyl+cylr,sphr),lymin, lymax);

return Intersection;
}

...
What is the real problem and how can i solve this issue?...The template class in your code snippet seems to be pretty much the same as in this thread: http://www.nr.com/forum/showthread.php?t=1920

The example that I showed in my initial response on that thread does not have any memory leaks. (Tested and verified with valgrind on my Linux workstations.) Since we can't possibly know what else is in your code in addition to what you posted (and I won't guess), I can't really help.


Regards,

Dave