Solving sparse NxM linear eqs

sunaj · **Joined:** Wed Mar 16, 2011 1:11 pm **Posts:** 5

Hi,

I have a non square system of linear equations and I'm currently using rmatrixsolvels to solve it. This works fine.
But my system is very large and computationally expensive to solve. It is also a sparse system (it has relatively few values outside the diagonal).
Can I solve the system more efficiently using sparse matrix functions (http://www.alglib.net/matrixops/sparse.php) and if so - how?

Any help will be much appreciated:)

sunaj · **Joined:** Wed Mar 16, 2011 1:11 pm **Posts:** 5

I found a solution using http://www.stanford.edu/group/SOL/software/lsqr/ - much faster and much more memory efficient.

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

Sorry for replying too late, but ALGLIB has LSQR solver :)

sunaj · **Joined:** Wed Mar 16, 2011 1:11 pm **Posts:** 5

Hi again,

I've had a second look at the alglib LSQR solver (because the other leaks memory and I can't figure out why). But when I use linlsqrsolvesparse, it does not give the same results as when I use rmatrixsolvels ( the other LSQR sovler gives exactly the same as rmatrixsolvels). Am I doing something wrong (sorry for my messy test code)?

Thanks in advance!

Code:

   string tmpstr = JFile::readFile("..\\testLinAlg\\vb.csv");
   alglib::real_2d_array mac;
   int n,m;
   double *b = JStr::toNumArr<double>(tmpstr,n,newline);
   double **a = new double*[n];
   vector<string> al = JFile::readFileLines("..\\testLinAlg\\ma.csv",true);
   if (al.size()!=n) throw JException("Files do not match.");
   for (int r=0;r<n;r++) {
      al[r] = al[r].substr(0,al[r].length()-2);
      a[r] = JStr::toNumArr<double>(al[r],m,",");
   }

   //4 check
   double *d = new double[n*m];
   memset(d,0,n*m*sizeof(double));
   mac.setcontent(n,m,d);

   //pop data
   sparsematrix ma;
    sparsecreate(n, m, ma);
   for (int r=0;r<n;r++){
      for (int c=0;c<m;c++){
         if (a[r][c]!=0) {
            mac[r][c] = a[r][c];
            sparseset(ma, r, c, a[r][c]);
         }
      }
   }
   sparseconverttocrs(ma);
   real_1d_array rhs;
   //rhs.setlength(n);
   //for (int r=0;r<n;r++)rhs[r] = b[r];
   rhs.setcontent(n,b);

   //check with full solver
   alglib::real_1d_array check;
   alglib::ae_int_t      info;
   alglib::densesolverlsreport repc;
   alglib::rmatrixsolvels(mac,n,m,rhs,0.0,info,repc,check);
   cout << "velAdj:" << endl; showMatrix(check,10); //last value: 0.000263801


   //solve
    linlsqrstate s;
    linlsqrreport rep;
    real_1d_array x;
    linlsqrcreate(n, m, s);
   linlsqrsetcond(s, 1.0e-10,1.0e-10,n+m+50);
    linlsqrsolvesparse(s, ma, rhs);
    linlsqrresults(s, x, rep);

   printf("%d\n", int(rep.terminationtype)); // EXPECTED: 4
   //printf("%s\n", x.tostring(6).c_str());
   double *res = x.getcontent();
    cout << JArray::join<double>(res,min(m,10),newline) << endl; //last value: 0.000123301

forum.alglib.net

Forum rules

Solving sparse NxM linear eqs

Who is online