GiNaC memory leak?

[Marko Riedel]

Greetings.

I am sending some code that I wrote in order to solve a differential
equation. You run the compiled program with a command-line argument
telling it what terms to use. Interesting values are e.g. 10, 14, 16
and 19. (The higher values take some time to compute.)

Warning! When you run the program on the higher values, it starts
consuming huge amounts of memory. On a SuSE Linux machine the OS
eventually kills the process. What is going on here? The program
solves a system of equations with a 1000 x 1000 approx. coefficient
matrix. What is all the memory being used for? Take, say 32 bytes per
entry; the system should then allocate about 32M of memory. In fact it 
tries to allocate more than 512M of memory!

Best regards,

Marko Riedel

----------------------------------------------------------------------

#include <stdio.h>
#include <fstream>
#include <ginac/ginac.h>

using namespace std;
using namespace GiNaC;

static symbol u("u");
static symbol z("z");
static symbol l1("l1");
static symbol l2("l2");
static symbol phi("phi");

static void extracteqs(ex eq, lst *l)
{
  int i, j;
  ex c1, c2;

  for(i=eq.ldegree(z); i<=eq.degree(z); i++){
    c1=eq.coeff(z, i);
    for(j=eq.ldegree(u); j<=eq.degree(u); j++){
      c2=c1.coeff(u, j);
      if(is_ex_of_type(c2, numeric)){
	if(c2!=0){
	  cerr << "inconsistency: " 
	       << i << " " << j << endl;
	  exit(-1);
	}
      }
      else{
	(*l).append(c2==0);
	cerr << (*l).nops() << endl;
      }
    }
  }
}

int main(int argc, char **argv)
{
  int maxcoeff, i, j;

  if(argc!=2){
    cout << argv[0] << " <maxcoeff>" << endl;
    exit(-1);
  }

  sscanf(argv[1], "%d", &maxcoeff);
  if(maxcoeff<1){
    cout << "argument " << maxcoeff 
	 << " out of range" << endl;
    exit(-2);
  }

  ex rt, lf1, lf2;
  lst vs;

  rt=0; lf1=0; lf2=0;
  for(i=0; i<maxcoeff; i++){
    for(j=0; j<maxcoeff; j++){
      symbol a;

      rt+=a*pow(z, i)*pow(u, j);
      vs.append(a);

      if(i<maxcoeff-6 && j<maxcoeff-6){
	symbol b, c;

	lf1+=b*pow(z, i)*pow(u, j);
	lf2+=c*pow(z, i)*pow(u, j);
	vs.append(b);
	vs.append(c);
      }
    }
  }

  cerr << rt << endl;
  cerr << lf1 << endl;
  cerr << lf2 << endl;

  cerr << vs << endl;

  ex s, sf, sf1, sf2;

  sf=pow(1-u, 11)*pow(1-z, 7)*pow(1-z*u, 7);
  sf1=pow(1-u, 9)*pow(1-z, 5)*pow(1-z*u, 5);
  sf2=pow(1-u, 10)*pow(1-z, 2)*pow(1-z*u, 2);
  s=rt/sf1+
    lf1*u*log(1/(1-z))/sf2+
    lf2*log(1/(1-z*u))/sf2;

  ex f, rz, rzv, rzu, rzvu;

  f=(u-u*u*z*z)/pow(1-z, 2)/pow(1-z*u, 2);
  rz=1/pow(1-z, 2);
  rzv=1/pow(1-z, 2);
  rzu=rz.subs(z==z*u);
  rzvu=rzv.subs(z==z*u);

  ex t1, t2, t3, t4, t5, t6, t7, t8, t9;

  t1=6*u*u*rzu *f;
  t2=6*u*u*rzvu*f;
  t3=6*u*u*rzu *phi;

  t4=6*u*u*rzu *rz;
  t5=6*u*u*rzvu*rz;
  t6=6*u*u*rzu *rzv;

  t7=6*    rz  *f;
  t8=6*    rzv *f;
  t9=6*    rz  *phi;
  
  ex rhs1, rhs2;
  
  rhs1=(t1+t2+t4+t5+t6+t7+t8).normal();
  rhs2=(t3+t9).normal();

  ex eq, eq1;

  eq=(s.diff(z, 2)-rhs1-rhs2.subs(phi==s))
    .subs(lst(log(1/(1-z))==l1, log(1/(1-z*u))==l2))
    .expand();
  eq1=0;
  
  for(i=0; i<eq.nops(); i++){
    eq1+=(eq.op(i)*sf).normal().expand();
    cerr << "expand: " << i << " " << eq.nops() << endl;
  }

  ex leq1=eq1.coeff(l1, 1), leq2=eq1.coeff(l2, 1);
  eq1=eq1.subs(lst(l1==0, l2==0));

  lst ceqs;

  ex s0=(s.subs(z==0)*pow(1-u, 9)).normal();
  ex sz0=(s.diff(z).subs(z==0)*pow(1-u, 9)).normal();

  ex i0=(u*pow(1-u, 9)).expand();
  ex iz0=((4*u+4*u*u)*pow(1-u, 9)).expand();

  cerr << "extract initial: 1" << endl;
  extracteqs(s0-i0, &ceqs);
  cerr << "extract initial: 2" << endl;
  extracteqs(sz0-iz0, &ceqs);

  cerr << "extract: 1" << endl;
  extracteqs(eq1, &ceqs);
  cerr << "eqs: " << ceqs.nops() << endl;
  cerr << "extract: 2" << endl;
  extracteqs(leq1, &ceqs);
  cerr << "eqs: " << ceqs.nops() << endl;
  cerr << "extract: 3" << endl;
  extracteqs(leq2, &ceqs);
  cerr << "eqs: " << ceqs.nops() << endl;

  ex sl;
  
  cerr << "solving: " 
       << ceqs.nops() << " " 
       << vs.nops() << " ..." << endl;
  try{
    sl=lsolve(ceqs, vs);
  } 
  catch(exception &e){
    cerr << e.what() << endl;
    exit(-3);
  }
  cout << sl << endl;
  
  lst zlist;
  for(i=0; i<vs.nops(); i++){
    zlist.append(vs.op(i)==0);
  }

  ex result;
  result=s.subs(sl).subs(zlist);
  cout << result << endl;
  
  archive a;
  a.archive_ex(result, "result");

  char fname[80];
  sprintf(fname, "fixedpoints%d.gar", maxcoeff);

  ofstream out(fname);
  out << a;
  out.close();
}