CGAL users discussion list

[Zohar Levi <>]

Dear Laurent,

You misread what I wrote almost completely. I did use the TAUCS solver and on windows I have a precompiled libs. Again I'll stress my problems with the library:

1. Sometimes it crashes.

2. Using it on the attached code beneath which doesn't need solving a large matrix, it works really-really slow. Over 4 minutes for a polyhedron compared to less than a second with my lapack solver (attached). And yes, on large matrices it works fast enough. (If you are using CG for small matrices, you have a problem there.)

3. Linking with the debug libraries generates a corrupted executable. It's the first time I see such a thing, and I didn't look into it yet.

For your convenience I've attached some other linear solvers I wrote this morning, if you don't like my matlab solver:

ublas based - using normal equations, really sucks.
lapack - based on your CGAL::svd solver. Works great and the fastest.

Fragments of my code which uses the solvers. It implements a parameterization of every one ring neighborhood in a polyhedron:

// uv should already be resized for the source polygon vertices count
bool parameterize_poly(Polyhedron &poly, vector<Point_2> &uv)
{
    typedef CGAL::Parameterization_polyhedron_adaptor_3<Polyhedron>
        Parameterization_polyhedron_adaptor;
    Parameterization_polyhedron_adaptor mesh_adaptor(poly);

//    typedef CGAL::Square_border_arc_length_parameterizer_3<Parameterization_polyhedron_adaptor>
//        Border_parameterizer;
    // Circular border parameterizer (the default)
    typedef CGAL::Circular_border_arc_length_parameterizer_3<Parameterization_polyhedron_adaptor>
        Border_parameterizer;

    // TAUCS solver
//    typedef CGAL::Taucs_solver_traits<double> Solver;
//    typedef CGAL::Taucs_symmetric_solver_traits<double> Solver;
//    typedef cgal_matlab_solver<double> Solver;
    typedef cgal_lapack_solver Solver;

    // Floater Mean Value Coordinates parameterization is the default
    // (defaults are circular border and OpenNL solver)
    typedef CGAL::Mean_value_coordinates_parameterizer_3<Parameterization_polyhedron_adaptor,
        Border_parameterizer, Solver>
        Parameterizer;

    // crashes
//    typedef CGAL::Discrete_authalic_parameterizer_3<Parameterization_polyhedron_adaptor>
//        Parameterizer;

    Parameterizer::Error_code err = CGAL::parameterize(mesh_adaptor, Parameterizer());

//Border_parameterizer border;
//Parameterizer::Error_code err = border.parameterize_border(mesh_adaptor);

    if ( err != Parameterizer::OK ) {
        cout << "parameterize_poly() Parameterizer error: " << Parameterizer::get_error_message(err) << endl;
        return false;
    }

//return true;
    // (u,v) pairs
    for ( Polyhedron::Vertex_const_iterator pVertex = poly.vertices_begin() ;
        pVertex != poly.vertices_end() ; pVertex++ ) {
        // (u,v) pair is stored in any halfedge
        uv[pVertex->id] = Point_2(mesh_adaptor.info(pVertex->halfedge())->uv().x(),
            mesh_adaptor.info(pVertex->halfedge())->uv().y());
//cout << uv[pVertex->id] << endl;
    }

    return true;
}


void doVertex(Polyhedron::Vertex_handle &v, vector<Point_2> &param_uv)
{
    // Build umbrella patch
    vector<Polyhedron::Face_handle> faces;
    faces.resize(v->vertex_degree());
    int fCounter = 0;
    Polyhedron::Halfedge_around_vertex_circulator cir = v->vertex_begin();
    do {
        Polyhedron::Face_handle f = cir->opposite()->face();
        faces[fCounter++] = f;
        cir++;
    } while ( cir != v->vertex_begin() );

    Polyhedron umbrella;
    create_polyhedron_from_faces(faces, umbrella);
//cout << "Umb " << umbrella.size_of_facets() << endl;

    if ( !parameterize_poly(umbrella, param_uv) )
        return;
}

void test_parameterization(Polyhedron &poly, vector<Point_2> &param_uv)
{
    for ( Polyhedron::Vertex_iterator vit = poly.vertices_begin()
        ; vit != poly.vertices_end() ; vit++ ) {
        if ( !is_border_vertex(vit) ) {
            do_vertex(vit, vMesh, param_uv);
//return;
        }
    }
}

---

From: Laurent Saboret <>
To:
Sent: Thursday, August 27, 2009 11:59:18 AM
Subject: Re: [cgal-discuss] Parameterization solver

Hi Zohar,

Floater's mean coordinates (class Mean_value_coordinates_parameterizer_3) and Tutte's barycenter (class Barycentric_mapping_parameterizer_3) have a template parameter SparseLinearAlgebraTraits_d equal by default to OpenNL::DefaultLinearSolverTraits.
OpenNL::DefaultLinearSolverTraits is a general purpose BICGSTAB solver.

OpenNL is the default as it is shipped with CGAL. Unfortunately, it cannot solve large linear systems.
TAUCS is very fast although it is 6 years old. As we use TAUCS out-of-core solvers, it can also solve very large systems.
If you managed to compile TAUCS, I recommend you to use it, ie set SparseLinearAlgebraTraits = Taucs_solver_traits. See Taucs_parameterization.cpp for an example.
You can interactively compare both solvers using polyhedron_ex_parameterization.cpp example.

Anyway, both solvers should parameterize a 15000 facets mesh in a few seconds.

Here are few numbers on a 3GHz/3Gb of RAM PC:
Model            #facets      Method        Solver      Time (sec)
one_ring          8            Floater      OpenNL      5e-20
oni              2845          Floater      OpenNL      0.1
casting          10224        Floater      OpenNL      0.4
camel            19536        Floater      OpenNL      2.5
david            47753        Floater      OpenNL      9
one_ring          8            Floater      TAUCS        0.26
oni              2845          Floater      TAUCS        0.27
casting          10224        Floater      TAUCS        0.6
camel            19536        Floater      TAUCS        1.1
david            47753        Floater      TAUCS        13.6

Please post your code and your input example if you wish a more specific answer.

Best regards,
Laurent Saboret
INRIA


Zohar Levi a écrit :
>
> 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.
> 
> ------------------------------------------------------------------------
> *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??


-- 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/Lapack/Linear_algebra_lapack.h>

typedef CGAL::Lapack_svd svd;

class Matrix_cgal_svd : public svd::Matrix {
public:
        Matrix_cgal_svd(int size, int size2) : svd::Matrix(size, size2) {}

        void set_coef(int row, int col, svd::FT val, bool flag) {
                set(row, col, val);
        }

        void print() {
                cout << "[";
                for ( unsigned i = 0 ; i < number_of_rows() ; i++ ) {
                        for ( unsigned j = 0 ; j < number_of_columns() ; j++ 
) {
                                cout << operator()(i, j);
                                if ( j+1 < number_of_columns() )
                                        cout << ", ";
                        }
                        if ( i+1 < number_of_rows() )
                                cout << ";\n";
                }
                cout << "]" << endl;
        }
};

class Vector_cgal_svd : public svd::Vector {
public:
        Vector_cgal_svd(int size) : svd::Vector(size) {}

        svd::FT& operator[](int ind) {
                return m_vector[ind];
        }

        svd::FT operator[](int ind) const {
                return m_vector[ind];
        }

        void print() {
                cout << "[";
                for ( unsigned i = 0 ; i < size() ; i++ ) {
                        cout << m_vector[i];
                        if ( i+1 < size() )
                                cout << "; ";
                }
                cout << "]" << endl;
        }
};

class cgal_lapack_solver {
public:
        typedef Matrix_cgal_svd Matrix;
        typedef Vector_cgal_svd Vector;

        cgal_lapack_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().
        // Mustn't change A,b
        bool linear_solver(Matrix& A, Vector& b, Vector& x, svd::FT& D)
        {
                D = 1;          // does not support homogeneous coordinates

                for ( unsigned i = 0 ; i < x.size() ; i++ )
                        x[i] = b[i];

                svd::Matrix A_(A.number_of_rows(), A.number_of_columns());
                for ( unsigned i = 0 ; i < A.number_of_rows() ; i++ )
                        for ( unsigned j = 0 ; j < A.number_of_columns() ; 
j++ )
                                A_.set(i, j, A(i, j));
                svd::solve(A_, x);

                return true;
        }

};




#pragma once

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp> 

namespace ublas = boost::numeric::ublas; 

template<class T>
class Matrix_cgal_ublas : public ublas::matrix<T> {
public:
        Matrix_cgal_ublas(int size, int size2) : ublas::matrix<T>(size, 
size2) {}

        void set_coef(int row, int col, T val, bool flag) {
                at_element(row, col) = val;
        }
};

template<class T>
class cgal_ublas_solver {
public:
        typedef Matrix_cgal_ublas<T> Matrix;
        typedef ublas::vector<T> Vector;

        cgal_ublas_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

                ublas::matrix<T> AtA = prod(trans(A),A);
//cout << A << endl;
//cout << AtA << endl;
                ublas::permutation_matrix<std::size_t> pm(AtA.size1());

                lu_factorize(AtA, pm);

//cout << pm << endl;

                ublas::matrix<T> invAtA = 
ublas::identity_matrix<T>(AtA.size1());

                lu_substitute(AtA, pm, invAtA); 
//cout << invAtA << endl;
                x = prod(invAtA, b);

                // is too big
                for ( int i = 0 ; i < x.size() ; i++ ) {
                        if ( fabs(x[i])-1.0 > 1.0e-8 || isnanf(x[i]) ) {
                                cout << x[i] << " is too big" << endl;
                                return false;
                        }
                }


                return true;
        }

};