CGAL users discussion list

[Zohar Levi <>]

Okay, I read the code more carefully. I'm not sure who put the comment of solving with CG. I looked in taucs_solver_traits.h and it seems to use the old interface (and I'm not sure it uses CG by default). There is also a class which implements the new taucs interface but it is not in use: taucs_symmetric_solver_traits. I tried playing with the options like "taucs.factor.lu=true", but to no avail.

Just so you won't say I'm just complaining all the time, I implemented a new solver traits using Matlab api (file attached). Although it took me most of the night to write it, it is stable (why taucs sometimes crashes while solving a simple linear system??), and parameterizing every face in a 15,000 faces mesh took less than a second comparing to the taucs which took more than 4 minutes!!

I wouldn't like to sully taucs reputation and all, especially since I like Sivan very much. Maybe we are doing something wrong with it, or it is good for very large systems. In this case, at least my own, I think one should add also an LAPACK simple linear solver to cgal.

Cheers.
 

BTW.

The Voronoi works great.

---

From: Zohar Levi <>
To:
Sent: Wednesday, August 26, 2009 9:17:18 PM
Subject: Re: [cgal-discuss] Parameterization solver

Ok, I found the problem with the parameterization. For Floatter's mean coordinates or Tutte's barycenter you need to calculate the lambdas and solve a simple linear system. You solved the linear system using conjugate gradients instead of a simple LU and forward substitution. Could anyone explain this please??

Could anyone hint on the simplest and fastest direction to solve this?

---

From: Laurent Saboret <>
To:
Sent: Tuesday, August 25, 2009 10:59:55 AM
Subject: Re: [cgal-discuss] Parameterization solver

>
> Hi,
>
> I used the parameterizer to parameterize a vertex with its first ring of neighbors (using default values). We are talking about 7 vertices in average. On a larger neighborhood the parameterizer worked fine. On this small neighborhood the default solver almost always complained it couldn't solve the system. Using the Tau solver it worked almost always - there was a time when it crashed during LU factorization - but there's a bigger problem: It took it almost a second to solve the system. This renders it to be impractical.
>
> I think you should also add precompiled debug libraries of the tau solver. I understand it crashes much more than the release version, which makes it unusable, but at least I could compile my application in debug, and operate other parts of it. Adding the libraries would help me with unnecessarily clattering my code with #ifdefs.
>
> Thanks.
>

Hi (again) Zohar,

CGAL 3.4 contains precompiled debug libraries of TAUCS and a few bug fixes in the parameterization code. Please try again with CGAL 3.4 in debug mode.

If the problem is still present, please post your code and input file to this mailing list.

Best regards,
-- Laurent Saboret
INRIA Sophia-Antipolis
-- You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://lists-sop.inria.fr/wws/info/cgal-discuss

#pragma once

#include <CGAL/Taucs_matrix.h>
#include <CGAL/Taucs_vector.h>

#include "utils.h"

template<class T>
class my_solver {
public:
        typedef CGAL::Taucs_matrix<T> Matrix;
        typedef CGAL::Taucs_vector<T> Vector;

        my_solver()     
        {
        }

        /// Solve the sparse linear system "A*x = b".
        /// Return true on success. The solution is then (1/D) * x.
        ///
        /// Preconditions:
        /// - A.row_dimension()    == b.dimension().
        /// - A.column_dimension() == x.dimension().
        bool linear_solver (const Matrix& A, const Vector& b, Vector& x, T& D)
        {
                D = 1;          // does not support homogeneous coordinates

                // A
                mwSize dim[2];
                dim[0] = A.row_dimension(); dim[1] = A.column_dimension();
                mxArray *arr = mxCreateNumericArray(2, dim, mxDOUBLE_CLASS, 
mxREAL);
                double *dat = (double*)mxGetData(arr);
                for ( unsigned int i = 0 ; i < dim[0] ; i++ )
                        for ( unsigned int j = 0 ; j < dim[1] ; j++ ) {
                                mwSize ind[2];
                                ind[0] = (mwSize)i; ind[1] = (mwSize)j;
                                int imat = mxCalcSingleSubscript(arr, 2, ind);
                                dat[imat] = A.get_coef(i, j);
                        }
                engPutVariable(getMatlab(), "cgal_A", arr);
                mxDestroyArray(arr);

                // b
                dim[0] = b.dimension();
                arr = mxCreateNumericArray(1, dim, mxDOUBLE_CLASS, mxREAL);
                dat = (double*)mxGetData(arr);
                for ( unsigned int i = 0 ; i < dim[0] ; i++ )
                        dat[i] = b[i];
                engPutVariable(getMatlab(), "cgal_b", arr);
                mxDestroyArray(arr);

                // x
                engEvalString(getMatlab(), "cgal_x = cgal_A \\ cgal_b;");
                const char *res = matlab_output_buffer();
                if ( res ) {
                        cout << res << endl;
                        return false;
                }

                arr = engGetVariable(getMatlab(), "cgal_x");
                dat = (double*)mxGetData(arr);
                for ( unsigned int i = 0 ; i < x.dimension() ; i++ )
                        x[i] = dat[i];
                mxDestroyArray(arr);

                return true;
        }
};