Bajib
02-02-2011, 02:20 AM
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;
}
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;
}