GiNaC memory leak?

[Richard B. Kreckel]

Hi there,

You have just proven Richy's second law:

1) Heuristics can be made to work in 1-epsilon of all cases.
2) Real problems tend to concentrate in the epsilon corner.

:-)

Your problem is the heuristics in the routine matrix::solve().  I have
just run your problem with 16 as input value and found that you have a
sparse 456x456 matrix that you try to solve for 456 variables.  Most of
the matrix entries are small integer values with 83 exceptions where
integers multiplied with (-1+u)^(-1) appear.  The heuristics chooses
Gaussian elimination only if no symbolic entries appear in the matrix and
I am quite surprised to see Bareiss' elimination perform so poor for
your case.  Well, you hit the epsilon corner, it appears.  In any case,
Gaussian elimination can solve your problem, it took about 20 minutes here
(whereupon it returned an empty solution set).

You have three options:
1) Don't use lsolve().  Build a matrix from the equation system, test if
   the system turned out to be linear and fill the variable matrix just 
   like lsolve does.  Then call matrix::solve() with a hint for the
   algorithm as last argument.
2) Patch your personal library of GiNaC by adding an algorithmic hint
   as last argument in ginac/inifcns.cpp by changing the line
       solution = sys.solve(vars,rhs);
   to something like
       solution = sys.solve(vars,rhs,solve_algo::gauss);
   and use that library instead.  But then you have disabled Heuristics
   completely and you may screw up in other cases.
3) Let lsolve() accept an algorithmic hint.  Pass that along to
   matrix::solve().  Test it.  Send us a patch and we'll include it.

Regards
    -richy.

PS: BTW, for arguments 10 and 14 your program tells me something about
    inconsistencies.
PPS: You shouldn't write is_ex_of_type(c2, numeric).  Write
    is_a<numeric>(c2) instead.
-- 
Richard Kreckel
<Richard.Kreckel at Uni-Mainz.DE>
<http://wwwthep.physik.uni-mainz.de/~kreckel/>


On Wed, 8 Aug 2001, Marko Riedel wrote:
> 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();
> }