adding gmm

contrib/gmm/gmm.h

+// -*- c++ -*- (enables emacs c++ mode)
+//===========================================================================
+//
+// Copyright (C) 2002-2008 Yves Renard
+//
+// This file is a part of GETFEM++
+//
+// Getfem++  is  free software;  you  can  redistribute  it  and/or modify it
+// under  the  terms  of the  GNU  Lesser General Public License as published
+// by  the  Free Software Foundation;  either version 2.1 of the License,  or
+// (at your option) any later version.
+// This program  is  distributed  in  the  hope  that it will be useful,  but
+// WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+// License for more details.
+// You  should  have received a copy of the GNU Lesser General Public License
+// along  with  this program;  if not, write to the Free Software Foundation,
+// Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+// As a special exception, you may use this file as part of a free software
+// library without restriction.  Specifically, if other files instantiate
+// templates or use macros or inline functions from this file, or you compile
+// this file and link it with other files to produce an executable, this
+// file does not by itself cause the resulting executable to be covered by
+// the GNU General Public License.  This exception does not however
+// invalidate any other reasons why the executable file might be covered by
+// the GNU General Public License.
+//
+//===========================================================================
+
+/**@file gmm.h
+   @author  Yves Renard <Yves.Renard@insa-lyon.fr>
+   @date October 13, 2002.
+   @brief Include common gmm files.
+*/
+#ifndef GMM_H__
+#define GMM_H__
+
+#include "gmm_kernel.h"
+#include "gmm_dense_lu.h"
+#include "gmm_dense_qr.h"
+
+#include "gmm_iter_solvers.h"
+#include "gmm_condition_number.h"
+#include "gmm_inoutput.h"
+
+#include "gmm_lapack_interface.h"
+#include "gmm_superlu_interface.h"
+
+
+#include "gmm_domain_decomp.h"
+
+#endif //  GMM_H__

contrib/gmm/gmm_MUMPS_interface.h

+// -*- c++ -*- (enables emacs c++ mode)
+//===========================================================================
+//
+// Copyright (C) 2003-2008 Yves Renard
+//
+// This file is a part of GETFEM++
+//
+// Getfem++  is  free software;  you  can  redistribute  it  and/or modify it
+// under  the  terms  of the  GNU  Lesser General Public License as published
+// by  the  Free Software Foundation;  either version 2.1 of the License,  or
+// (at your option) any later version.
+// This program  is  distributed  in  the  hope  that it will be useful,  but
+// WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+// License for more details.
+// You  should  have received a copy of the GNU Lesser General Public License
+// along  with  this program;  if not, write to the Free Software Foundation,
+// Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+// As a special exception, you may use this file as part of a free software
+// library without restriction.  Specifically, if other files instantiate
+// templates or use macros or inline functions from this file, or you compile
+// this file and link it with other files to produce an executable, this
+// file does not by itself cause the resulting executable to be covered by
+// the GNU General Public License.  This exception does not however
+// invalidate any other reasons why the executable file might be covered by
+// the GNU General Public License.
+//
+//===========================================================================
+
+/**@file gmm_MUMPS_interface.h
+   @author Yves Renard <Yves.Renard@insa-lyon.fr>,
+   @author Julien Pommier <Julien.Pommier@insa-toulouse.fr>
+   @date December 8, 2005.
+   @brief Interface with MUMPS (LU direct solver for sparse matrices).
+*/
+#if defined(GMM_USES_MUMPS)
+
+#ifndef GMM_MUMPS_INTERFACE_H
+#define GMM_MUMPS_INTERFACE_H
+
+#include "gmm_kernel.h"
+
+
+extern "C" {
+
+#include <smumps_c.h>
+#undef F_INT
+#undef F_DOUBLE
+#undef F_DOUBLE2
+#include <dmumps_c.h>
+#undef F_INT
+#undef F_DOUBLE
+#undef F_DOUBLE2
+#include <cmumps_c.h>
+#undef F_INT
+#undef F_DOUBLE
+#undef F_DOUBLE2
+#include <zmumps_c.h>
+#undef F_INT
+#undef F_DOUBLE
+#undef F_DOUBLE2
+
+}
+
+namespace gmm {
+
+  template <typename T> struct ij_sparse_matrix {
+    std::vector<int> irn;
+    std::vector<int> jcn;
+    std::vector<T>        a;
+    
+    template <typename L> void store(const L& l, size_type i) {
+       typename linalg_traits<L>::const_iterator it = vect_const_begin(l),
+	 ite = vect_const_end(l);
+       for (; it != ite; ++it)
+	 { irn.push_back(i + 1); jcn.push_back(it.index() + 1); a.push_back(*it); }
+    }
+
+    template <typename L> void build_from(const L& l, row_major) {
+      for (size_type i = 0; i < mat_nrows(l); ++i)
+	store(mat_const_row(l, i), i);
+    }
+
+    template <typename L> void build_from(const L& l, col_major) {
+      for (size_type i = 0; i < mat_ncols(l); ++i)
+	store(mat_const_col(l, i), i);
+      irn.swap(jcn);
+    }
+
+    template <typename L> ij_sparse_matrix(const L& A) {
+      size_type nz = nnz(A);
+      irn.reserve(nz); jcn.reserve(nz); a.reserve(nz);
+      build_from(A,  typename principal_orientation_type<typename
+	       linalg_traits<L>::sub_orientation>::potype());
+    }
+  };
+
+  /* ********************************************************************* */
+  /*   MUMPS solve interface                                               */
+  /* ********************************************************************* */
+
+
+  template <typename T> struct mumps_interf {};
+
+  template <> struct mumps_interf<float> {
+    typedef SMUMPS_STRUC_C  MUMPS_STRUC_C;
+    typedef float value_type;
+
+    static void mumps_c(MUMPS_STRUC_C &id) { smumps_c(&id); }
+  };
+
+  template <> struct mumps_interf<double> {
+    typedef DMUMPS_STRUC_C  MUMPS_STRUC_C;
+    typedef double value_type;
+    static void mumps_c(MUMPS_STRUC_C &id) { dmumps_c(&id); }
+  };
+
+  template <> struct mumps_interf<std::complex<float> > {
+    typedef CMUMPS_STRUC_C  MUMPS_STRUC_C;
+    typedef mumps_complex value_type;
+    static void mumps_c(MUMPS_STRUC_C &id) { cmumps_c(&id); }
+  };
+
+  template <> struct mumps_interf<std::complex<double> > {
+    typedef ZMUMPS_STRUC_C  MUMPS_STRUC_C;
+    typedef mumps_double_complex value_type;
+    static void mumps_c(MUMPS_STRUC_C &id) { zmumps_c(&id); }
+  };
+
+
+  /** MUMPS solve interface  
+   *  Works only with sparse or skyline matrices
+   */
+  template <typename MAT, typename VECTX, typename VECTB>
+  void MUMPS_solve(const MAT &A, const VECTX &X_, const VECTB &B) {
+    VECTX &X = const_cast<VECTX &>(X_);
+
+    typedef typename linalg_traits<MAT>::value_type T;
+    typedef typename mumps_interf<T>::value_type MUMPS_T;
+    GMM_ASSERT2(gmm::mat_nrows(A) == gmm::mat_ncols(A), "Non square matrix");
+  
+    std::vector<T> rhs(gmm::vect_size(B)); gmm::copy(B, rhs);
+  
+    ij_sparse_matrix<T> AA(A);
+  
+    const int JOB_INIT = -1;
+    const int JOB_END = -2;
+    const int USE_COMM_WORLD = -987654;
+
+    typename mumps_interf<T>::MUMPS_STRUC_C id;
+
+#ifdef GMM_USES_MPI
+    int rank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+#endif
+    
+    id.job = JOB_INIT;
+    id.par = 1;
+    id.sym = 0;
+    id.comm_fortran = USE_COMM_WORLD;
+    mumps_interf<T>::mumps_c(id);
+    
+#ifdef GMM_USES_MPI
+    if (rank == 0) {
+#endif
+      id.n = gmm::mat_nrows(A);
+      id.nz = AA.irn.size();
+      id.irn = &(AA.irn[0]);
+      id.jcn = &(AA.jcn[0]);
+      id.a = (MUMPS_T*)(&(AA.a[0]));
+      id.rhs = (MUMPS_T*)(&(rhs[0]));
+#ifdef GMM_USES_MPI
+    }
+#endif
+
+#define ICNTL(I) icntl[(I)-1]
+#define INFO(I) info[(I)-1]
+    id.ICNTL(1) = -1; id.ICNTL(2) = -1; id.ICNTL(3) = -1; id.ICNTL(4) = 0;
+    id.job = 6;
+    
+    id.ICNTL(14) += 40; /* small boost to the workspace size as we have encountered some problem
+			   who did not fit in the default settings of mumps.. 
+			   by default, ICNTL(14) = 15 or 20
+		       */
+    //cout << "ICNTL(14): " << id.ICNTL(14) << "\n";
+
+    mumps_interf<T>::mumps_c(id);
+    if (id.INFO(1) < 0) {
+      switch (id.INFO(1)) {
+	case -6 : case -10 :
+	  GMM_ASSERT1(false, "Solve with MUMPS failed: matrix is singular");
+	case -13 :
+	  GMM_ASSERT1(false, "Solve with MUMPS failed: not enough memory");
+        case -9:
+	  GMM_ASSERT1(false, "Solve with MUMPS failed: error " << id.INFO(1)
+		      << ", increase ICNTL(14)");
+	default :
+	  GMM_ASSERT1(false, "Solve with MUMPS failed with error "
+		      << id.INFO(1));
+      }
+    }
+
+    id.job = JOB_END;
+    mumps_interf<T>::mumps_c(id);
+
+    gmm::copy(rhs, X);
+
+#undef ICNTL
+#undef INFO
+
+  }
+
+
+
+  /** MUMPS solve interface for distributed matrices 
+   *  Works only with sparse or skyline matrices
+   */
+  template <typename MAT, typename VECTX, typename VECTB>
+  void MUMPS_distributed_matrix_solve(const MAT &A, const VECTX &X_,
+				      const VECTB &B) {
+    VECTX &X = const_cast<VECTX &>(X_);
+
+    typedef typename linalg_traits<MAT>::value_type T;
+    typedef typename mumps_interf<T>::value_type MUMPS_T;
+    GMM_ASSERT2(gmm::mat_nrows(A) == gmm::mat_ncols(A), "Non-square matrix");
+  
+    std::vector<T> rhs(gmm::vect_size(B)); gmm::copy(B, rhs);
+  
+    ij_sparse_matrix<T> AA(A);
+  
+    const int JOB_INIT = -1;
+    const int JOB_END = -2;
+    const int USE_COMM_WORLD = -987654;
+
+    typename mumps_interf<T>::MUMPS_STRUC_C id;
+
+#ifdef GMM_USES_MPI
+    int rank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+#endif
+    
+    id.job = JOB_INIT;
+    id.par = 1;
+    id.sym = 0;
+    id.comm_fortran = USE_COMM_WORLD;
+    mumps_interf<T>::mumps_c(id);
+    
+    id.n = gmm::mat_nrows(A);
+    id.nz_loc = AA.irn.size();
+    id.irn_loc = &(AA.irn[0]);
+    id.jcn_loc = &(AA.jcn[0]);
+    id.a_loc = (MUMPS_T*)(&(AA.a[0]));
+
+#ifdef GMM_USES_MPI
+    if (rank == 0) {
+#endif
+      id.rhs = (MUMPS_T*)(&(rhs[0]));
+#ifdef GMM_USES_MPI
+    }
+#endif
+
+#define ICNTL(I) icntl[(I)-1]
+#define INFO(I) info[(I)-1]
+    id.ICNTL(1) = -1; id.ICNTL(2) = 6; // id.ICNTL(2) = -1;
+    id.ICNTL(3) = 6;
+    // id.ICNTL(3) = -1; 
+    id.ICNTL(4) = 2;
+    id.ICNTL(5)=0; id.ICNTL(18)=3;
+    id.job = 6;
+    mumps_interf<T>::mumps_c(id);
+    if (id.INFO(1) < 0) {
+      switch (id.INFO(1)) {
+      case -6 : case -10 :
+	GMM_ASSERT1(false, "Solve with MUMPS failed: matrix is singular");
+      case -13:
+	GMM_ASSERT1(false, "Solve with MUMPS failed: not enough memory");
+      case -9:
+	GMM_ASSERT1(false, "Solve with MUMPS failed: error " << id.INFO(1)
+		    << ", increase ICNTL(14)");
+      default :
+	GMM_ASSERT1(false, "Solve with MUMPS failed with error " <<id.INFO(1));
+      }
+    }
+
+    id.job = JOB_END;
+    mumps_interf<T>::mumps_c(id);
+#ifdef GMM_USES_MPI
+    MPI_Bcast(&(rhs[0]),id.n,gmm::mpi_type(T()),0,MPI_COMM_WORLD);
+#endif
+    gmm::copy(rhs, X);
+
+#undef ICNTL
+#undef INFO
+
+  }
+
+
+
+
+}
+
+  
+#endif // GMM_MUMPS_INTERFACE_H
+
+#endif // GMM_USES_MUMPS

contrib/gmm/gmm_algobase.h

+// -*- c++ -*- (enables emacs c++ mode)
+//===========================================================================
+//
+// Copyright (C) 2000-2008 Yves Renard
+//
+// This file is a part of GETFEM++
+//
+// Getfem++  is  free software;  you  can  redistribute  it  and/or modify it
+// under  the  terms  of the  GNU  Lesser General Public License as published
+// by  the  Free Software Foundation;  either version 2.1 of the License,  or
+// (at your option) any later version.
+// This program  is  distributed  in  the  hope  that it will be useful,  but
+// WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+// License for more details.
+// You  should  have received a copy of the GNU Lesser General Public License
+// along  with  this program;  if not, write to the Free Software Foundation,
+// Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+// As a special exception, you may use this file as part of a free software
+// library without restriction.  Specifically, if other files instantiate
+// templates or use macros or inline functions from this file, or you compile
+// this file and link it with other files to produce an executable, this
+// file does not by itself cause the resulting executable to be covered by
+// the GNU General Public License.  This exception does not however
+// invalidate any other reasons why the executable file might be covered by
+// the GNU General Public License.
+//
+//===========================================================================
+
+/** @file gmm_algobase.h 
+   @author  Yves Renard <Yves.Renard@insa-lyon.fr>
+   @date September 28, 2000.
+    @brief Miscelleanous algorithms on containers.
+*/
+
+#ifndef GMM_ALGOBASE_H__
+#define GMM_ALGOBASE_H__
+#include "gmm_std.h"
+#include "gmm_except.h"
+#include <functional>
+
+namespace gmm {
+
+  /* ********************************************************************* */
+  /* Definitition de classes de comparaison.                               */
+  /* retournant un int.                                                    */
+  /* ********************************************************************* */
+  
+  template <class T>
+    struct less : public std::binary_function<T, T, int> {
+    inline int operator()(const T& x, const T& y) const
+    { return (x < y) ? -1 : ((y < x) ? 1 : 0); }
+  };
+
+  template<> struct less<int> : public std::binary_function<int, int, int>
+  { int operator()(int x, int y) const { return x-y; } };
+  template<> struct less<char> : public std::binary_function<char, char, int>
+  { int operator()(char x, char y) const { return int(x-y); } };
+  template<> struct less<short> : public std::binary_function<short,short,int>
+  { int operator()(short x, short y) const { return int(x-y); } };
+  template<> struct less<unsigned char>
+     : public std::binary_function<unsigned char, unsigned char, int> {
+    int operator()(unsigned char x, unsigned char y) const
+    { return int(x)-int(y); }
+  };
+  
+
+  template <class T>
+    struct greater : public std::binary_function<T, T, int> {
+    inline int operator()(const T& x, const T& y) const
+    { return (y < x) ? -1 : ((x < y) ? 1 : 0); }
+  };
+
+  template<> struct greater<int> : public std::binary_function<int, int, int>
+  { int operator()(int x, int y) const { return y-x; } };
+  template<> struct greater<char> : public std::binary_function<char,char,int>
+  { int operator()(char x, char y) const { return int(y-x); } };
+  template<> struct greater<short>
+      : public std::binary_function<short, short, int>
+  { int operator()(short x, short y) const { return int(y-x); } };
+  template<> struct greater<unsigned char>
+    : public std::binary_function<unsigned char, unsigned char, int> {
+    int operator()(unsigned char x, unsigned char y) const
+      { return int(y)-int(x); }
+  };
+
+  template <typename T> inline T my_abs(T a) { return (a < T(0)) ? T(-a) : a; }
+  
+  template <class T>
+    struct approx_less : public std::binary_function<T, T, int> { 
+    double eps;
+    inline int operator()(const T &x, const T &y) const
+    { if (my_abs(x - y) <= eps) return 0; if (x < y) return -1; return 1; }
+    approx_less(double e = 1E-13) { eps = e; }
+  };
+
+  template <class T>
+    struct approx_greater : public std::binary_function<T, T, int> { 
+    double eps;
+    inline int operator()(const T &x, const T &y) const
+    { if (my_abs(x - y) <= eps) return 0; if (x > y) return -1; return 1; }
+    approx_greater(double e = 1E-13) { eps = e; }
+  };
+
+  template<class ITER1, class ITER2, class COMP>
+    int lexicographical_compare(ITER1 b1, const ITER1 &e1,
+				ITER2 b2, const ITER2 &e2, const COMP &c)  {
+    int i;
+    for ( ; b1 != e1 && b2 != e2; ++b1, ++b2)
+      if ((i = c(*b1, *b2)) != 0) return i;
+    if (b1 != e1) return 1; if (b2 != e2) return -1; return 0; 
+  }
+
+  template<class CONT, class COMP = gmm::less<typename CONT::value_type> >
+    struct lexicographical_less : public std::binary_function<CONT, CONT, int>
+  { 
+    COMP c;
+    int operator()(const CONT &x, const CONT &y) const {
+      return gmm::lexicographical_compare(x.begin(), x.end(),
+					  y.begin(), y.end(), c);
+    }
+    lexicographical_less(const COMP &d = COMP()) { c = d; }
+  };
+
+  template<class CONT, class COMP = gmm::less<typename CONT::value_type> >
+  struct lexicographical_greater
+    : public std::binary_function<CONT, CONT, int> { 
+    COMP c;
+    int operator()(const CONT &x, const CONT &y) const {
+      return -gmm::lexicographical_compare(x.begin(), x.end(),
+					   y.begin(), y.end(), c);
+    }
+    lexicographical_greater(const COMP &d = COMP()) { c = d; }
+  };
+  
+
+  /* ********************************************************************* */
+  /* "Virtual" iterators on sequences.                                     */
+  /* The class T represent a class of sequence.                            */
+  /* ********************************************************************* */
+
+  template<class T> struct sequence_iterator {
+    
+    typedef T             value_type;
+    typedef value_type*   pointer;
+    typedef value_type&   reference;
+    typedef const value_type& const_reference;
+    typedef std::forward_iterator_tag iterator_category;
+
+    T Un;
+
+    sequence_iterator(T U0 = T(0)) { Un = U0; }
+    
+    sequence_iterator &operator ++()
+    { ++Un; return *this; }
+    sequence_iterator operator ++(int)
+    { sequence_iterator tmp = *this; (*this)++; return tmp; }
+	
+    const_reference operator *() const { return Un; }
+    reference operator *() { return Un; }
+    
+    bool operator ==(const sequence_iterator &i) const { return (i.Un==Un);}
+    bool operator !=(const sequence_iterator &i) const { return (i.Un!=Un);}
+  };
+
+  /* ********************************************************************* */
+  /* generic algorithms.                                                   */
+  /* ********************************************************************* */
+
+  template <class ITER1, class SIZE, class ITER2>
+  ITER2 copy_n(ITER1 first, SIZE count, ITER2 result) {
+    for ( ; count > 0; --count, ++first, ++result) *result = *first;
+    return result;
+  }
+
+  template<class ITER>
+    typename std::iterator_traits<ITER>::value_type
+      mean_value(ITER first, const ITER &last) {
+    GMM_ASSERT2(first != last, "mean value of empty container");
+    size_t n = 1;
+    typename std::iterator_traits<ITER>::value_type res = *first++;
+    while (first != last) { res += *first; ++first; ++n; }
+    res /= n;
+    return res;
+  }
+
+  template<class CONT>
+    typename CONT::value_type
+  mean_value(const CONT &c) { return mean_value(c.begin(), c.end()); }
+
+  template<class ITER> /* hum ... */
+    void minmax_box(typename std::iterator_traits<ITER>::value_type &pmin,
+		    typename std::iterator_traits<ITER>::value_type &pmax,
+		    ITER first, const ITER &last) {
+    typedef typename std::iterator_traits<ITER>::value_type PT;
+    if (first != last) { pmin = pmax = *first; ++first; }
+    while (first != last) {
+      typename PT::const_iterator b = (*first).begin(), e = (*first).end();
+      typename PT::iterator b1 = pmin.begin(), b2 = pmax.begin();
+      while (b != e)
+	{ *b1 = std::min(*b1, *b); *b2 = std::max(*b2, *b); ++b; ++b1; ++b2; }
+    }
+  }
+
+  template<typename VEC> struct sorted_indexes_aux {
+    const VEC &v;
+  public:
+    sorted_indexes_aux(const VEC& v_) : v(v_) {}
+    template <typename IDX>
+    bool operator()(const IDX &ia, const IDX &ib) const
+    { return v[ia] < v[ib]; }
+  };
+
+  template<typename VEC, typename IVEC> 
+  void sorted_indexes(const VEC &v, IVEC &iv) {
+    iv.clear(); iv.resize(v.size());
+    for (size_t i=0; i < v.size(); ++i) iv[i] = i;
+    std::sort(iv.begin(), iv.end(), sorted_indexes_aux<VEC>(v));
+  }
+
+}
+
+
+#endif /* GMM_ALGOBASE_H__ */

contrib/gmm/gmm_condition_number.h

+// -*- c++ -*- (enables emacs c++ mode)
+//===========================================================================
+//
+// Copyright (C) 2003-2008 Yves Renard
+//
+// This file is a part of GETFEM++
+//
+// Getfem++  is  free software;  you  can  redistribute  it  and/or modify it
+// under  the  terms  of the  GNU  Lesser General Public License as published
+// by  the  Free Software Foundation;  either version 2.1 of the License,  or
+// (at your option) any later version.
+// This program  is  distributed  in  the  hope  that it will be useful,  but
+// WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+// License for more details.
+// You  should  have received a copy of the GNU Lesser General Public License
+// along  with  this program;  if not, write to the Free Software Foundation,
+// Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+// As a special exception, you may use this file as part of a free software
+// library without restriction.  Specifically, if other files instantiate
+// templates or use macros or inline functions from this file, or you compile
+// this file and link it with other files to produce an executable, this
+// file does not by itself cause the resulting executable to be covered by
+// the GNU General Public License.  This exception does not however
+// invalidate any other reasons why the executable file might be covered by
+// the GNU General Public License.
+//
+//===========================================================================
+
+/**@file gmm_condition_number.h
+   @author  Yves Renard <Yves.Renard@insa-lyon.fr>, Julien Pommier <Julien.Pommier@insa-toulouse.fr>
+   @date August 27, 2003.
+   @brief computation of the condition number of dense matrices.
+*/
+#ifndef GMM_CONDITION_NUMBER_H__
+#define GMM_CONDITION_NUMBER_H__
+
+#include "gmm_dense_qr.h"
+
+namespace gmm {
+
+  /** computation of the condition number of dense matrices using SVD.
+
+      Uses symmetric_qr_algorithm => dense matrices only.
+
+      @param M a matrix.
+      @param emin smallest (in magnitude) eigenvalue
+      @param emax largest eigenvalue.
+   */
+  template <typename MAT> 
+  typename number_traits<typename 
+  linalg_traits<MAT>::value_type>::magnitude_type
+  condition_number(const MAT& M, 
+	  typename number_traits<typename
+	  linalg_traits<MAT>::value_type>::magnitude_type& emin,
+	  typename number_traits<typename
+	  linalg_traits<MAT>::value_type>::magnitude_type& emax) {
+    typedef typename linalg_traits<MAT>::value_type T;
+    typedef typename number_traits<T>::magnitude_type R;
+
+    size_type m = mat_nrows(M), n = mat_ncols(M);
+    emax = emin = R(0);
+    std::vector<R> eig(m+n);
+
+    if (m+n == 0) return R(0);
+    if (is_hermitian(M)) {
+      eig.resize(m);
+      gmm::symmetric_qr_algorithm(M, eig);
+    }
+    else {
+      dense_matrix<T> B(m+n, m+n); // not very efficient ??
+      gmm::copy(conjugated(M), sub_matrix(B, sub_interval(m, n), sub_interval(0, m)));
+      gmm::copy(M, sub_matrix(B, sub_interval(0, m),
+					  sub_interval(m, n)));
+      gmm::symmetric_qr_algorithm(B, eig);
+    }
+    emin = emax = gmm::abs(eig[0]);
+    for (size_type i = 1; i < eig.size(); ++i) {
+      R e = gmm::abs(eig[i]); 
+      emin = std::min(emin, e);
+      emax = std::max(emax, e);
+    }
+    // cout << "emin = " << emin << " emax = " << emax << endl;
+    if (emin == R(0)) return gmm::default_max(R());
+    return emax / emin;
+  }
+
+  template <typename MAT> 
+  typename number_traits<typename 
+  linalg_traits<MAT>::value_type>::magnitude_type
+  condition_number(const MAT& M) { 
+    typename number_traits<typename
+      linalg_traits<MAT>::value_type>::magnitude_type emax, emin;
+    return condition_number(M, emin, emax);
+  }
+
+  template <typename MAT> 
+  typename number_traits<typename 
+  linalg_traits<MAT>::value_type>::magnitude_type
+  Frobenius_condition_number_sqr(const MAT& M) { 
+    typedef typename linalg_traits<MAT>::value_type T;
+    typedef typename number_traits<T>::magnitude_type R;
+    size_type m = mat_nrows(M), n = mat_ncols(M);
+    dense_matrix<T> B(std::min(m,n), std::min(m,n));
+    if (m < n) mult(M,gmm::conjugated(M),B);
+    else       mult(gmm::conjugated(M),M,B);
+    R trB = abs(mat_trace(B));
+    lu_inverse(B);
+    return trB*abs(mat_trace(B));
+  }
+
+  template <typename MAT> 
+  typename number_traits<typename 
+  linalg_traits<MAT>::value_type>::magnitude_type
+  Frobenius_condition_number(const MAT& M)
+  { return sqrt(Frobenius_condition_number_sqr(M)); }
+
+  /** estimation of the condition number (TO BE DONE...)
+   */
+  template <typename MAT> 
+  typename number_traits<typename 
+  linalg_traits<MAT>::value_type>::magnitude_type
+  condest(const MAT& M, 
+	  typename number_traits<typename
+	  linalg_traits<MAT>::value_type>::magnitude_type& emin,
+	  typename number_traits<typename
+	  linalg_traits<MAT>::value_type>::magnitude_type& emax) {
+    return condition_number(M, emin, emax);
+  }
+  
+  template <typename MAT> 
+  typename number_traits<typename 
+  linalg_traits<MAT>::value_type>::magnitude_type
+  condest(const MAT& M) { 
+    typename number_traits<typename
+      linalg_traits<MAT>::value_type>::magnitude_type emax, emin;
+    return condest(M, emin, emax);
+  }
+}
+
+#endif

contrib/gmm/gmm_conjugated.h

+// -*- c++ -*- (enables emacs c++ mode)
+//===========================================================================
+//
+// Copyright (C) 2003-2008 Yves Renard
+//
+// This file is a part of GETFEM++
+//
+// Getfem++  is  free software;  you  can  redistribute  it  and/or modify it
+// under  the  terms  of the  GNU  Lesser General Public License as published
+// by  the  Free Software Foundation;  either version 2.1 of the License,  or
+// (at your option) any later version.
+// This program  is  distributed  in  the  hope  that it will be useful,  but
+// WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+// License for more details.
+// You  should  have received a copy of the GNU Lesser General Public License
+// along  with  this program;  if not, write to the Free Software Foundation,
+// Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+// As a special exception, you may use this file as part of a free software
+// library without restriction.  Specifically, if other files instantiate
+// templates or use macros or inline functions from this file, or you compile
+// this file and link it with other files to produce an executable, this
+// file does not by itself cause the resulting executable to be covered by
+// the GNU General Public License.  This exception does not however
+// invalidate any other reasons why the executable file might be covered by
+// the GNU General Public License.
+//
+//===========================================================================
+
+/**@file gmm_conjugated.h
+   @author  Yves Renard <Yves.Renard@insa-lyon.fr>
+   @date September 18, 2003.
+   @brief handle conjugation of complex matrices/vectors.
+*/
+#ifndef GMM_CONJUGATED_H__
+#define GMM_CONJUGATED_H__
+
+#include "gmm_def.h"
+
+namespace gmm {
+  ///@cond DOXY_SHOW_ALL_FUNCTIONS
+
+  /* ********************************************************************* */
+  /*		Conjugated references on vectors            		   */
+  /* ********************************************************************* */
+
+  template <typename IT> struct conjugated_const_iterator {
+    typedef typename std::iterator_traits<IT>::value_type      value_type;
+    typedef typename std::iterator_traits<IT>::pointer         pointer;
+    typedef typename std::iterator_traits<IT>::reference       reference;
+    typedef typename std::iterator_traits<IT>::difference_type difference_type;
+    typedef typename std::iterator_traits<IT>::iterator_category
+    iterator_category;
+
+    IT it;
+    
+    conjugated_const_iterator(void) {}
+    conjugated_const_iterator(const IT &i) : it(i) {}
+    
+    inline size_type index(void) const { return it.index(); }
+    conjugated_const_iterator operator ++(int)
+    { conjugated_const_iterator tmp = *this; ++it; return tmp; }
+    conjugated_const_iterator operator --(int) 
+    { conjugated_const_iterator tmp = *this; --it; return tmp; }
+    conjugated_const_iterator &operator ++() { ++it; return *this; }
+    conjugated_const_iterator &operator --() { --it; return *this; }
+    conjugated_const_iterator &operator +=(difference_type i)
+      { it += i; return *this; }
+    conjugated_const_iterator &operator -=(difference_type i)
+      { it -= i; return *this; }
+    conjugated_const_iterator operator +(difference_type i) const
+      { conjugated_const_iterator itb = *this; return (itb += i); }
+    conjugated_const_iterator operator -(difference_type i) const
+      { conjugated_const_iterator itb = *this; return (itb -= i); }
+    difference_type operator -(const conjugated_const_iterator &i) const
+      { return difference_type(it - i.it); }
+    
+    value_type operator  *() const { return gmm::conj(*it); }
+    value_type operator [](size_type ii) const { return gmm::conj(it[ii]); }
+    
+    bool operator ==(const conjugated_const_iterator &i) const
+      { return (i.it == it); }
+    bool operator !=(const conjugated_const_iterator &i) const
+      { return (i.it != it); }
+    bool operator < (const conjugated_const_iterator &i) const
+      { return (it < i.it); }
+  };
+
+  template <typename V> struct conjugated_vector_const_ref {
+    typedef conjugated_vector_const_ref<V> this_type;
+    typedef typename linalg_traits<V>::value_type value_type;
+    typedef typename linalg_traits<V>::const_iterator iterator;
+    typedef typename linalg_traits<this_type>::reference reference;
+    typedef typename linalg_traits<this_type>::origin_type origin_type;
+
+    iterator begin_, end_;
+    const origin_type *origin;
+    size_type size_;
+
+    conjugated_vector_const_ref(const V &v)
+      : begin_(vect_const_begin(v)), end_(vect_const_end(v)),
+	origin(linalg_origin(v)),
+	size_(vect_size(v)) {}
+
+    reference operator[](size_type i) const
+    { return gmm::conj(linalg_traits<V>::access(origin, begin_, end_, i)); }
+  };
+
+  template <typename V> struct linalg_traits<conjugated_vector_const_ref<V> > {
+    typedef conjugated_vector_const_ref<V> this_type;
+    typedef typename linalg_traits<V>::origin_type origin_type;
+    typedef linalg_const is_reference;
+    typedef abstract_vector linalg_type;
+    typedef typename linalg_traits<V>::value_type value_type;
+    typedef value_type reference;
+    typedef abstract_null_type iterator;
+    typedef conjugated_const_iterator<typename
+                   linalg_traits<V>::const_iterator> const_iterator;
+    typedef typename linalg_traits<V>::storage_type storage_type;
+    typedef typename linalg_traits<V>::index_sorted index_sorted;
+    static size_type size(const this_type &v) { return v.size_; }
+    static iterator begin(this_type &v) { return iterator(v.begin_); }
+    static const_iterator begin(const this_type &v)
+    { return const_iterator(v.begin_); }
+    static iterator end(this_type &v)
+    { return iterator(v.end_); }
+    static const_iterator end(const this_type &v)
+    { return const_iterator(v.end_); }
+    static value_type access(const origin_type *o, const const_iterator &it,
+			     const const_iterator &ite, size_type i)
+    { return gmm::conj(linalg_traits<V>::access(o, it.it, ite.it, i)); }
+    static const origin_type* origin(const this_type &v) { return v.origin; }
+  };
+
+  template<typename V> std::ostream &operator <<
+    (std::ostream &o, const conjugated_vector_const_ref<V>& m)
+  { gmm::write(o,m); return o; }
+
+  /* ********************************************************************* */
+  /*		Conjugated references on matrices            		   */
+  /* ********************************************************************* */
+
+  template <typename M> struct conjugated_row_const_iterator {
+    typedef conjugated_row_const_iterator<M> iterator;
+    typedef typename linalg_traits<M>::const_row_iterator ITER;
+    typedef typename linalg_traits<M>::value_type value_type;
+    typedef ptrdiff_t difference_type;
+    typedef size_t size_type;
+
+    ITER it;
+
+    iterator operator ++(int) { iterator tmp = *this; it++; return tmp; }
+    iterator operator --(int) { iterator tmp = *this; it--; return tmp; }
+    iterator &operator ++()   { it++; return *this; }
+    iterator &operator --()   { it--; return *this; }
+    iterator &operator +=(difference_type i) { it += i; return *this; }
+    iterator &operator -=(difference_type i) { it -= i; return *this; }
+    iterator operator +(difference_type i) const 
+    { iterator itt = *this; return (itt += i); }
+    iterator operator -(difference_type i) const
+    { iterator itt = *this; return (itt -= i); }
+    difference_type operator -(const iterator &i) const
+    { return it - i.it; }
+
+    ITER operator *() const { return it; }
+    ITER operator [](int i) { return it + i; }
+
+    bool operator ==(const iterator &i) const { return (it == i.it); }
+    bool operator !=(const iterator &i) const { return !(i == *this); }
+    bool operator < (const iterator &i) const { return (it < i.it); }
+
+    conjugated_row_const_iterator(void) {}
+    conjugated_row_const_iterator(const ITER &i) : it(i) { }
+
+  };
+
+  template <typename M> struct  conjugated_row_matrix_const_ref {
+    
+    typedef conjugated_row_matrix_const_ref<M> this_type;
+    typedef typename linalg_traits<M>::const_row_iterator iterator;
+    typedef typename linalg_traits<M>::value_type value_type;
+    typedef typename linalg_traits<this_type>::origin_type origin_type;
+
+    iterator begin_, end_;
+    const origin_type *origin;
+    size_type nr, nc;
+
+    conjugated_row_matrix_const_ref(const M &m)
+      : begin_(mat_row_begin(m)), end_(mat_row_end(m)),
+	origin(linalg_origin(m)), nr(mat_ncols(m)), nc(mat_nrows(m)) {}
+
+    value_type operator()(size_type i, size_type j) const
+    { return gmm::conj(linalg_traits<M>::access(begin_+j, i)); }
+  };
+
+  template <typename M>
+  struct linalg_traits<conjugated_row_matrix_const_ref<M> > {
+    typedef conjugated_row_matrix_const_ref<M> this_type;
+    typedef typename linalg_traits<M>::origin_type origin_type;
+    typedef linalg_const is_reference;
+    typedef abstract_matrix linalg_type;
+    typedef typename linalg_traits<M>::value_type value_type;
+    typedef value_type reference;
+    typedef typename linalg_traits<M>::storage_type storage_type;
+    typedef typename linalg_traits<M>::const_sub_row_type vector_type;
+    typedef conjugated_vector_const_ref<vector_type> sub_col_type;
+    typedef conjugated_vector_const_ref<vector_type> const_sub_col_type;
+    typedef conjugated_row_const_iterator<M> col_iterator;
+    typedef conjugated_row_const_iterator<M> const_col_iterator;
+    typedef abstract_null_type const_sub_row_type;
+    typedef abstract_null_type sub_row_type;
+    typedef abstract_null_type const_row_iterator;
+    typedef abstract_null_type row_iterator;
+    typedef col_major sub_orientation;
+    typedef typename linalg_traits<M>::index_sorted index_sorted;
+    static inline size_type ncols(const this_type &m) { return m.nc; }
+    static inline size_type nrows(const this_type &m) { return m.nr; }
+    static inline const_sub_col_type col(const const_col_iterator &it)
+    { return conjugated(linalg_traits<M>::row(it.it)); }
+    static inline const_col_iterator col_begin(const this_type &m)
+    { return const_col_iterator(m.begin_); }
+    static inline const_col_iterator col_end(const this_type &m)
+    { return const_col_iterator(m.end_); }
+    static inline const origin_type* origin(const this_type &m)
+    { return m.origin; }
+    static value_type access(const const_col_iterator &it, size_type i)
+    { return gmm::conj(linalg_traits<M>::access(it.it, i)); }
+  };
+
+  template<typename M> std::ostream &operator <<
+  (std::ostream &o, const conjugated_row_matrix_const_ref<M>& m)
+  { gmm::write(o,m); return o; }
+
+
+  template <typename M> struct conjugated_col_const_iterator {
+    typedef conjugated_col_const_iterator<M> iterator;
+    typedef typename linalg_traits<M>::const_col_iterator ITER;
+    typedef typename linalg_traits<M>::value_type value_type;
+    typedef ptrdiff_t difference_type;
+    typedef size_t size_type;
+
+    ITER it;
+
+    iterator operator ++(int) { iterator tmp = *this; it++; return tmp; }
+    iterator operator --(int) { iterator tmp = *this; it--; return tmp; }
+    iterator &operator ++()   { it++; return *this; }
+    iterator &operator --()   { it--; return *this; }
+    iterator &operator +=(difference_type i) { it += i; return *this; }
+    iterator &operator -=(difference_type i) { it -= i; return *this; }
+    iterator operator +(difference_type i) const 
+    { iterator itt = *this; return (itt += i); }
+    iterator operator -(difference_type i) const
+    { iterator itt = *this; return (itt -= i); }
+    difference_type operator -(const iterator &i) const
+    { return it - i.it; }
+
+    ITER operator *() const { return it; }
+    ITER operator [](int i) { return it + i; }
+
+    bool operator ==(const iterator &i) const { return (it == i.it); }
+    bool operator !=(const iterator &i) const { return !(i == *this); }
+    bool operator < (const iterator &i) const { return (it < i.it); }
+
+    conjugated_col_const_iterator(void) {}
+    conjugated_col_const_iterator(const ITER &i) : it(i) { }
+
+  };
+
+  template <typename M> struct  conjugated_col_matrix_const_ref {
+    
+    typedef conjugated_col_matrix_const_ref<M> this_type;
+    typedef typename linalg_traits<M>::const_col_iterator iterator;
+    typedef typename linalg_traits<M>::value_type value_type;
+    typedef typename linalg_traits<this_type>::origin_type origin_type;
+
+    iterator begin_, end_;
+    const origin_type *origin;
+    size_type nr, nc;
+
+    conjugated_col_matrix_const_ref(const M &m)
+      : begin_(mat_col_begin(m)), end_(mat_col_end(m)),
+	origin(linalg_origin(m)), nr(mat_ncols(m)), nc(mat_nrows(m)) {}
+
+    value_type operator()(size_type i, size_type j) const
+    { return gmm::conj(linalg_traits<M>::access(begin_+i, j)); }
+  };
+
+  template <typename M>
+  struct linalg_traits<conjugated_col_matrix_const_ref<M> > {
+    typedef conjugated_col_matrix_const_ref<M> this_type;
+    typedef typename linalg_traits<M>::origin_type origin_type;
+    typedef linalg_const is_reference;
+    typedef abstract_matrix linalg_type;
+    typedef typename linalg_traits<M>::value_type value_type;
+    typedef value_type reference;
+    typedef typename linalg_traits<M>::storage_type storage_type;
+    typedef typename linalg_traits<M>::const_sub_col_type vector_type;
+    typedef conjugated_vector_const_ref<vector_type> sub_row_type;
+    typedef conjugated_vector_const_ref<vector_type> const_sub_row_type;
+    typedef conjugated_col_const_iterator<M> row_iterator;
+    typedef conjugated_col_const_iterator<M> const_row_iterator;
+    typedef abstract_null_type const_sub_col_type;
+    typedef abstract_null_type sub_col_type;
+    typedef abstract_null_type const_col_iterator;
+    typedef abstract_null_type col_iterator;
+    typedef row_major sub_orientation;
+    typedef typename linalg_traits<M>::index_sorted index_sorted;
+    static inline size_type nrows(const this_type &m) { return m.nr; }
+    static inline size_type ncols(const this_type &m) { return m.nc; }
+    static inline const_sub_row_type row(const const_row_iterator &it)
+    { return conjugated(linalg_traits<M>::col(it.it)); }
+    static inline const_row_iterator row_begin(const this_type &m)
+    { return const_row_iterator(m.begin_); }
+    static inline const_row_iterator row_end(const this_type &m)
+    { return const_row_iterator(m.end_); }
+    static inline const origin_type* origin(const this_type &m)
+    { return m.origin; }
+    static value_type access(const const_row_iterator &it, size_type i)
+    { return gmm::conj(linalg_traits<M>::access(it.it, i)); }
+  };
+
+  template<typename M> std::ostream &operator <<
+  (std::ostream &o, const conjugated_col_matrix_const_ref<M>& m)
+  { gmm::write(o,m); return o; }
+
+
+  template <typename L, typename SO> struct conjugated_return__ {
+    typedef conjugated_row_matrix_const_ref<L> return_type;
+  };
+  template <typename L> struct conjugated_return__<L, col_major> {
+    typedef conjugated_col_matrix_const_ref<L> return_type;
+  };
+  template <typename L, typename T, typename LT> struct conjugated_return_ {
+    typedef const L & return_type;
+  };
+  template <typename L, typename T>
+  struct conjugated_return_<L, std::complex<T>, abstract_vector> {
+    typedef conjugated_vector_const_ref<L> return_type;
+  };
+  template <typename L, typename T>
+  struct conjugated_return_<L, T, abstract_matrix> {
+    typedef typename conjugated_return__<L,
+    typename principal_orientation_type<typename
+    linalg_traits<L>::sub_orientation>::potype
+    >::return_type return_type;
+  };
+  template <typename L> struct conjugated_return {
+    typedef typename
+    conjugated_return_<L, typename linalg_traits<L>::value_type,
+		       typename linalg_traits<L>::linalg_type		       
+		       >::return_type return_type;
+  };
+
+  ///@endcond
+  /** return a conjugated view of the input matrix or vector. */
+  template <typename L> inline
+  typename conjugated_return<L>::return_type
+  conjugated(const L &v) {
+    return conjugated(v, typename linalg_traits<L>::value_type(),
+		      typename linalg_traits<L>::linalg_type());
+  }
+  ///@cond DOXY_SHOW_ALL_FUNCTIONS
+
+  template <typename L, typename T, typename LT> inline
+  const L & conjugated(const L &v, T, LT) { return v; }
+
+  template <typename L, typename T> inline
+  conjugated_vector_const_ref<L> conjugated(const L &v, std::complex<T>,
+					    abstract_vector)
+  { return conjugated_vector_const_ref<L>(v); }
+
+  template <typename L, typename T> inline
+  typename conjugated_return__<L,
+    typename principal_orientation_type<typename
+    linalg_traits<L>::sub_orientation>::potype>::return_type
+  conjugated(const L &v, T, abstract_matrix) {
+    return conjugated(v, typename principal_orientation_type<typename
+		      linalg_traits<L>::sub_orientation>::potype());
+  }
+
+  template <typename L> inline
+  conjugated_row_matrix_const_ref<L> conjugated(const L &v, row_major)
+  { return conjugated_row_matrix_const_ref<L>(v); }
+
+  template <typename L> inline
+  conjugated_col_matrix_const_ref<L> conjugated(const L &v, col_major)
+  { return conjugated_col_matrix_const_ref<L>(v); }
+  
+  ///@endcond
+  
+
+}
+
+#endif //  GMM_CONJUGATED_H__

contrib/gmm/gmm_dense_Householder.h

+// -*- c++ -*- (enables emacs c++ mode)
+//===========================================================================
+//
+// Copyright (C) 2003-2008 Yves Renard
+//
+// This file is a part of GETFEM++
+//
+// Getfem++  is  free software;  you  can  redistribute  it  and/or modify it
+// under  the  terms  of the  GNU  Lesser General Public License as published
+// by  the  Free Software Foundation;  either version 2.1 of the License,  or
+// (at your option) any later version.
+// This program  is  distributed  in  the  hope  that it will be useful,  but
+// WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+// License for more details.
+// You  should  have received a copy of the GNU Lesser General Public License
+// along  with  this program;  if not, write to the Free Software Foundation,
+// Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+// As a special exception, you may use this file as part of a free software
+// library without restriction.  Specifically, if other files instantiate
+// templates or use macros or inline functions from this file, or you compile
+// this file and link it with other files to produce an executable, this
+// file does not by itself cause the resulting executable to be covered by
+// the GNU General Public License.  This exception does not however
+// invalidate any other reasons why the executable file might be covered by
+// the GNU General Public License.
+//
+//===========================================================================
+
+/**@file gmm_dense_Householder.h
+   @author Caroline Lecalvez <Caroline.Lecalvez@gmm.insa-toulouse.fr>
+   @author Yves Renard <Yves.Renard@insa-lyon.fr>
+   @date June 5, 2003.
+   @brief Householder for dense matrices.
+*/
+
+#ifndef GMM_DENSE_HOUSEHOLDER_H
+#define GMM_DENSE_HOUSEHOLDER_H
+
+#include "gmm_kernel.h"
+
+namespace gmm {
+  ///@cond DOXY_SHOW_ALL_FUNCTIONS
+
+  /* ********************************************************************* */
+  /*    Rank one update  (complex and real version)                        */
+  /* ********************************************************************* */
+
+  template <typename Matrix, typename VecX, typename VecY>
+  inline void rank_one_update(Matrix &A, const VecX& x,
+			      const VecY& y, row_major) {
+    typedef typename linalg_traits<Matrix>::value_type T;
+    size_type N = mat_nrows(A);
+    GMM_ASSERT2(N <= vect_size(x) && mat_ncols(A) <= vect_size(y),
+	       "dimensions mismatch");
+    typename linalg_traits<VecX>::const_iterator itx = vect_const_begin(x);
+    for (size_type i = 0; i < N; ++i, ++itx) {
+      typedef typename linalg_traits<Matrix>::sub_row_type row_type;
+      row_type row = mat_row(A, i);
+      typename linalg_traits<row_type>::iterator
+	it = vect_begin(row), ite = vect_end(row);
+      typename linalg_traits<VecY>::const_iterator ity = vect_const_begin(y);
+      T tx = *itx;
+      for (; it != ite; ++it, ++ity) *it += conj_product(*ity, tx);
+    }
+  }
+
+  template <typename Matrix, typename VecX, typename VecY>
+  inline void rank_one_update(Matrix &A, const VecX& x,
+			      const VecY& y, col_major) {
+    typedef typename linalg_traits<Matrix>::value_type T;
+    size_type M = mat_ncols(A);
+    GMM_ASSERT2(mat_nrows(A) <= vect_size(x) && M <= vect_size(y),
+		"dimensions mismatch");
+    typename linalg_traits<VecY>::const_iterator ity = vect_const_begin(y);
+    for (size_type i = 0; i < M; ++i, ++ity) {
+      typedef typename linalg_traits<Matrix>::sub_col_type col_type;
+      col_type col = mat_col(A, i);
+      typename linalg_traits<col_type>::iterator
+	it = vect_begin(col), ite = vect_end(col);
+      typename linalg_traits<VecX>::const_iterator itx = vect_const_begin(x);
+      T ty = *ity;
+      for (; it != ite; ++it, ++itx) *it += conj_product(ty, *itx);
+    }
+  }
+  
+  ///@endcond
+  template <typename Matrix, typename VecX, typename VecY>
+  inline void rank_one_update(const Matrix &AA, const VecX& x,
+			      const VecY& y) {
+    Matrix& A = const_cast<Matrix&>(AA);
+    rank_one_update(A, x, y, typename principal_orientation_type<typename
+		    linalg_traits<Matrix>::sub_orientation>::potype());
+  }
+  ///@cond DOXY_SHOW_ALL_FUNCTIONS
+
+  /* ********************************************************************* */
+  /*    Rank two update  (complex and real version)                        */
+  /* ********************************************************************* */
+
+  template <typename Matrix, typename VecX, typename VecY>
+  inline void rank_two_update(Matrix &A, const VecX& x,
+			      const VecY& y, row_major) {
+    typedef typename linalg_traits<Matrix>::value_type value_type;
+    size_type N = mat_nrows(A);
+    GMM_ASSERT2(N <= vect_size(x) && mat_ncols(A) <= vect_size(y),
+		"dimensions mismatch");
+    typename linalg_traits<VecX>::const_iterator itx1 = vect_const_begin(x);
+    typename linalg_traits<VecY>::const_iterator ity2 = vect_const_begin(y);
+    for (size_type i = 0; i < N; ++i, ++itx1, ++ity2) {
+      typedef typename linalg_traits<Matrix>::sub_row_type row_type;
+      row_type row = mat_row(A, i);
+      typename linalg_traits<row_type>::iterator
+	it = vect_begin(row), ite = vect_end(row);
+      typename linalg_traits<VecX>::const_iterator itx2 = vect_const_begin(x);
+      typename linalg_traits<VecY>::const_iterator ity1 = vect_const_begin(y);
+      value_type tx = *itx1, ty = *ity2;
+      for (; it != ite; ++it, ++ity1, ++itx2)
+	*it += conj_product(*ity1, tx) + conj_product(*itx2, ty);
+    }
+  }
+
+  template <typename Matrix, typename VecX, typename VecY>
+  inline void rank_two_update(Matrix &A, const VecX& x,
+			      const VecY& y, col_major) {
+    typedef typename linalg_traits<Matrix>::value_type value_type;
+    size_type M = mat_ncols(A);
+    GMM_ASSERT2(mat_nrows(A) <= vect_size(x) && M <= vect_size(y),
+		"dimensions mismatch");
+    typename linalg_traits<VecX>::const_iterator itx2 = vect_const_begin(x);
+    typename linalg_traits<VecY>::const_iterator ity1 = vect_const_begin(y);
+    for (size_type i = 0; i < M; ++i, ++ity1, ++itx2) {
+      typedef typename linalg_traits<Matrix>::sub_col_type col_type;
+      col_type col = mat_col(A, i);
+      typename linalg_traits<col_type>::iterator
+	it = vect_begin(col), ite = vect_end(col);
+      typename linalg_traits<VecX>::const_iterator itx1 = vect_const_begin(x);
+      typename linalg_traits<VecY>::const_iterator ity2 = vect_const_begin(y);
+      value_type ty = *ity1, tx = *itx2;
+      for (; it != ite; ++it, ++itx1, ++ity2)
+	*it += conj_product(ty, *itx1) + conj_product(tx, *ity2); 
+    }
+  }
+  
+  ///@endcond
+  template <typename Matrix, typename VecX, typename VecY>
+  inline void rank_two_update(const Matrix &AA, const VecX& x,
+			      const VecY& y) {
+    Matrix& A = const_cast<Matrix&>(AA);
+    rank_two_update(A, x, y, typename principal_orientation_type<typename
+		    linalg_traits<Matrix>::sub_orientation>::potype());
+  }
+  ///@cond DOXY_SHOW_ALL_FUNCTIONS
+
+  /* ********************************************************************* */
+  /*    Householder vector computation (complex and real version)          */
+  /* ********************************************************************* */
+
+  template <typename VECT> void house_vector(const VECT &VV) {
+    VECT &V = const_cast<VECT &>(VV);
+    typedef typename linalg_traits<VECT>::value_type T;
+    typedef typename number_traits<T>::magnitude_type R;
+    
+    R mu = vect_norm2(V), abs_v0 = gmm::abs(V[0]);
+    if (mu != R(0))
+      gmm::scale(V, (abs_v0 == R(0)) ? T(R(1) / mu)
+		 : (safe_divide(T(abs_v0), V[0]) / (abs_v0 + mu)));
+    if (gmm::real(V[vect_size(V)-1]) * R(0) != R(0)) gmm::clear(V);
+    V[0] = T(1);
+  }
+
+  template <typename VECT> void house_vector_last(const VECT &VV) {
+    VECT &V = const_cast<VECT &>(VV);
+    typedef typename linalg_traits<VECT>::value_type T;
+    typedef typename number_traits<T>::magnitude_type R;
+
+    size_type m = vect_size(V);
+    R mu = vect_norm2(V), abs_v0 = gmm::abs(V[m-1]);
+    if (mu != R(0))
+      gmm::scale(V, (abs_v0 == R(0)) ? T(R(1) / mu)
+		 : ((abs_v0 / V[m-1]) / (abs_v0 + mu)));
+    if (gmm::real(V[0]) * R(0) != R(0)) gmm::clear(V);
+    V[m-1] = T(1);
+  }
+  
+  /* ********************************************************************* */
+  /*    Householder updates  (complex and real version)                    */
+  /* ********************************************************************* */
+
+  // multiply A to the left by the reflector stored in V. W is a temporary.
+  template <typename MAT, typename VECT1, typename VECT2> inline
+    void row_house_update(const MAT &AA, const VECT1 &V, const VECT2 &WW) {
+    VECT2 &W = const_cast<VECT2 &>(WW); MAT &A = const_cast<MAT &>(AA);
+    typedef typename linalg_traits<MAT>::value_type value_type;
+    typedef typename number_traits<value_type>::magnitude_type magnitude_type;
+
+    gmm::mult(conjugated(A),
+	      scaled(V, value_type(magnitude_type(-2)/vect_norm2_sqr(V))), W);
+    rank_one_update(A, V, W);
+  }
+
+  // multiply A to the right by the reflector stored in V. W is a temporary.
+  template <typename MAT, typename VECT1, typename VECT2> inline
+    void col_house_update(const MAT &AA, const VECT1 &V, const VECT2 &WW) {
+    VECT2 &W = const_cast<VECT2 &>(WW); MAT &A = const_cast<MAT &>(AA);
+    typedef typename linalg_traits<MAT>::value_type value_type;
+    typedef typename number_traits<value_type>::magnitude_type magnitude_type;
+    
+    gmm::mult(A,
+	      scaled(V, value_type(magnitude_type(-2)/vect_norm2_sqr(V))), W);
+    rank_one_update(A, W, V);
+  }
+
+  ///@endcond
+
+  /* ********************************************************************* */
+  /*    Hessemberg reduction with Householder.                             */
+  /* ********************************************************************* */
+
+  template <typename MAT1, typename MAT2>
+    void Hessenberg_reduction(const MAT1& AA, const MAT2 &QQ, bool compute_Q){
+    MAT1& A = const_cast<MAT1&>(AA); MAT2& Q = const_cast<MAT2&>(QQ);
+    typedef typename linalg_traits<MAT1>::value_type value_type;
+    if (compute_Q) gmm::copy(identity_matrix(), Q);
+    size_type n = mat_nrows(A); if (n < 2) return;
+    std::vector<value_type> v(n), w(n);
+    sub_interval SUBK(0,n);
+    for (size_type k = 1; k+1 < n; ++k) {
+      sub_interval SUBI(k, n-k), SUBJ(k-1,n-k+1);
+      v.resize(n-k);
+      for (size_type j = k; j < n; ++j) v[j-k] = A(j, k-1);
+      house_vector(v);
+      row_house_update(sub_matrix(A, SUBI, SUBJ), v, sub_vector(w, SUBJ));
+      col_house_update(sub_matrix(A, SUBK, SUBI), v, w);
+      // is it possible to "unify" the two on the common part of the matrix?
+      if (compute_Q) col_house_update(sub_matrix(Q, SUBK, SUBI), v, w);
+    }
+  }
+
+  /* ********************************************************************* */
+  /*    Householder tridiagonalization for symmetric matrices              */
+  /* ********************************************************************* */
+
+  template <typename MAT1, typename MAT2> 
+  void Householder_tridiagonalization(const MAT1 &AA, const MAT2 &QQ,
+				     bool compute_q) {
+    MAT1 &A = const_cast<MAT1 &>(AA); MAT2 &Q = const_cast<MAT2 &>(QQ);
+    typedef typename linalg_traits<MAT1>::value_type T;
+    typedef typename number_traits<T>::magnitude_type R;
+
+    size_type n = mat_nrows(A); if (n < 2) return;
+    std::vector<T> v(n), p(n), w(n), ww(n);
+    sub_interval SUBK(0,n);
+
+    for (size_type k = 1; k+1 < n; ++k) { // not optimized ...
+      sub_interval SUBI(k, n-k);
+      v.resize(n-k); p.resize(n-k); w.resize(n-k); 
+      for (size_type l = k; l < n; ++l) 
+	{ v[l-k] = w[l-k] = A(l, k-1); A(l, k-1) = A(k-1, l) = T(0); }
+      house_vector(v);
+      R norm = vect_norm2_sqr(v);
+      A(k-1, k) = gmm::conj(A(k, k-1) = w[0] - T(2)*v[0]*vect_hp(w, v)/norm);
+
+      gmm::mult(sub_matrix(A, SUBI), gmm::scaled(v, T(-2) / norm), p);
+      gmm::add(p, gmm::scaled(v, -vect_hp(v, p) / norm), w);
+      rank_two_update(sub_matrix(A, SUBI), v, w);
+      // it should be possible to compute only the upper or lower part
+
+      if (compute_q) col_house_update(sub_matrix(Q, SUBK, SUBI), v, ww);
+    }
+  }
+
+  /* ********************************************************************* */
+  /*    Real and complex Givens rotations                                  */
+  /* ********************************************************************* */
+
+  template <typename T> void Givens_rotation(T a, T b, T &c, T &s) {
+    typedef typename number_traits<T>::magnitude_type R;
+    R aa = gmm::abs(a), bb = gmm::abs(b);
+    if (bb == R(0)) { c = T(1); s = T(0);   return; }
+    if (aa == R(0)) { c = T(0); s = b / bb; return; }
+    if (bb > aa)
+      { T t = -safe_divide(a,b); s = T(R(1) / (sqrt(R(1)+gmm::abs_sqr(t)))); c = s * t; }
+    else
+      { T t = -safe_divide(b,a); c = T(R(1) / (sqrt(R(1)+gmm::abs_sqr(t)))); s = c * t; }
+  }
+
+  // Apply Q* v
+  template <typename T> inline
+  void Apply_Givens_rotation_left(T &x, T &y, T c, T s)
+  { T t1=x, t2=y; x = gmm::conj(c)*t1 - gmm::conj(s)*t2; y = c*t2 + s*t1; }
+
+  // Apply v^T Q
+  template <typename T> inline
+  void Apply_Givens_rotation_right(T &x, T &y, T c, T s)
+  { T t1=x, t2=y; x = c*t1 - s*t2; y = gmm::conj(c)*t2 + gmm::conj(s)*t1; }
+
+  template <typename MAT, typename T>
+  void row_rot(const MAT &AA, T c, T s, size_type i, size_type k) {
+    MAT &A = const_cast<MAT &>(AA); // can be specialized for row matrices
+    for (size_type j = 0; j < mat_ncols(A); ++j)
+      Apply_Givens_rotation_left(A(i,j), A(k,j), c, s);
+  }
+
+  template <typename MAT, typename T>
+  void col_rot(const MAT &AA, T c, T s, size_type i, size_type k) {
+    MAT &A = const_cast<MAT &>(AA); // can be specialized for column matrices
+    for (size_type j = 0; j < mat_nrows(A); ++j)
+      Apply_Givens_rotation_right(A(j,i), A(j,k), c, s);
+  }
+  
+}
+
+#endif
+

contrib/gmm/gmm_dense_lu.h

+// -*- c++ -*- (enables emacs c++ mode)
+//===========================================================================
+//
+// Copyright (C) 2003-2008 Yves Renard
+//
+// This file is a part of GETFEM++
+//
+// Getfem++  is  free software;  you  can  redistribute  it  and/or modify it
+// under  the  terms  of the  GNU  Lesser General Public License as published
+// by  the  Free Software Foundation;  either version 2.1 of the License,  or
+// (at your option) any later version.
+// This program  is  distributed  in  the  hope  that it will be useful,  but
+// WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+// License for more details.
+// You  should  have received a copy of the GNU Lesser General Public License
+// along  with  this program;  if not, write to the Free Software Foundation,
+// Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+// As a special exception, you may use this file as part of a free software
+// library without restriction.  Specifically, if other files instantiate
+// templates or use macros or inline functions from this file, or you compile
+// this file and link it with other files to produce an executable, this
+// file does not by itself cause the resulting executable to be covered by
+// the GNU General Public License.  This exception does not however
+// invalidate any other reasons why the executable file might be covered by
+// the GNU General Public License.
+//
+//===========================================================================
+
+// This file is a modified version of lu.h from MTL.
+// See http://osl.iu.edu/research/mtl/
+// Following the corresponding Copyright notice.
+//===========================================================================
+// Copyright (c) 2001-2003 The Trustees of Indiana University.
+// All rights reserved.
+// Copyright (c) 1998-2001 University of Notre Dame. All rights reserved.
+// Authors: Andrew Lumsdaine, Jeremy G. Siek, Lie-Quan Lee
+//
+// This file is part of the Matrix Template Library
+//
+// Indiana University has the exclusive rights to license this product under
+// the following license.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//   1.  All redistributions of source code must retain the above  copyright
+//       notice, the list of authors in the original source code,  this list
+//       of conditions and the disclaimer listed in this license;
+//   2.  All redistributions in binary form must reproduce the above
+//       copyright notice, this list of conditions and the disclaimer listed
+//       in this license in the documentation and/or other materials provided
+//       with the distribution;
+//   3.  Any documentation included with all redistributions must include the
+//       following acknowledgement:
+//	"This product includes software developed at the University of
+//       Notre Dame and the Pervasive Technology Labs at Indiana University.
+//       For technical information contact Andrew Lumsdaine at the Pervasive
+//       Technology Labs at Indiana University. 
+//	 For administrative and license questions contact the Advanced
+//       Research and Technology Institute at 1100 Waterway Blvd.
+//       Indianapolis, Indiana 46202, phone 317-274-5905, fax  317-274-5902."
+//	 Alternatively, this acknowledgement may appear in the software
+//       itself, and wherever such third-party acknowledgments normally appear.
+//   4.  The name "MTL" shall not be used to endorse or promote  products
+//       derived from this software without prior written  permission from
+//       Indiana University. For written permission, please  contact Indiana
+//       University Advanced Research & Technology  Institute.
+//   5.  Products derived from this software may not be called "MTL", nor
+//       may "MTL" appear in their name, without  prior written permission
+//       of Indiana University Advanced Research & Technology Institute.
+//
+// Indiana University provides no reassurances that the source code provided
+// does not infringe the patent or any other intellectual property rights of
+// any other entity. Indiana University disclaims any liability to any
+// recipient for claims brought by any other entity based on infringement of 
+// intellectual property rights or otherwise.
+//
+// LICENSEE UNDERSTANDS THAT SOFTWARE IS PROVIDED "AS IS" FOR WHICH NO
+// WARRANTIES AS TO CAPABILITIES OR ACCURACY ARE MADE. INDIANA UNIVERSITY
+// GIVES NO WARRANTIES AND MAKES NO REPRESENTATION THAT SOFTWARE IS FREE OF
+// INFRINGEMENT OF THIRD PARTY PATENT, COPYRIGHT, OR OTHER PROPRIETARY RIGHTS. 
+// INDIANA UNIVERSITY MAKES NO WARRANTIES THAT SOFTWARE IS FREE FROM "BUGS",
+// "VIRUSES", "TROJAN HORSES", "TRAP DOORS", "WORMS", OR OTHER HARMFUL CODE.
+// LICENSEE ASSUMES THE ENTIRE RISK AS TO THE PERFORMANCE OF SOFTWARE AND/OR
+// ASSOCIATED MATERIALS, AND TO THE PERFORMANCE AND VALIDITY OF 
+// INFORMATION GENERATED USING SOFTWARE.
+//===========================================================================
+
+/**@file gmm_dense_lu.h
+   @author  Andrew Lumsdaine, Jeremy G. Siek, Lie-Quan Lee, Y. Renard
+   @date June 5, 2003.
+   @brief LU factorizations and determinant computation for dense matrices.
+*/
+#ifndef GMM_DENSE_LU_H
+#define GMM_DENSE_LU_H
+
+#include "gmm_dense_Householder.h"
+#include "gmm_opt.h"
+
+namespace gmm {
+
+
+  /** LU Factorization of a general (dense) matrix (real or complex).
+  
+  This is the outer product (a level-2 operation) form of the LU
+  Factorization with pivoting algorithm . This is equivalent to
+  LAPACK's dgetf2. Also see "Matrix Computations" 3rd Ed.  by Golub
+  and Van Loan section 3.2.5 and especially page 115.
+  
+  The pivot indices in ipvt are indexed starting from 1
+  so that this is compatible with LAPACK (Fortran).
+  */
+  template <typename DenseMatrix, typename Pvector>
+  size_type lu_factor(DenseMatrix& A, Pvector& ipvt) {
+    typedef typename linalg_traits<DenseMatrix>::value_type T;
+    typedef typename number_traits<T>::magnitude_type R;
+    size_type info(0), i, j, jp, M(mat_nrows(A)), N(mat_ncols(A));
+    size_type NN = std::min(M, N);
+    std::vector<T> c(M), r(N);
+    
+    GMM_ASSERT2(ipvt.size()+1 >= NN, "IPVT too small");
+    for (i = 0; i+1 < NN; ++i) ipvt[i] = i;
+      
+    if (M || N) {
+      for (j = 0; j+1 < NN; ++j) {
+	R max = gmm::abs(A(j,j)); jp = j;
+	for (i = j+1; i < M; ++i)		   /* find pivot.          */
+	  if (gmm::abs(A(i,j)) > max) { jp = i; max = gmm::abs(A(i,j)); }
+	ipvt[j] = jp + 1;
+	
+	if (max == R(0)) { info = j + 1; break; }
+        if (jp != j) for (i = 0; i < N; ++i) std::swap(A(jp, i), A(j, i));
+	
+        for (i = j+1; i < M; ++i) { A(i, j) /= A(j,j); c[i-j-1] = -A(i, j); }
+        for (i = j+1; i < N; ++i) r[i-j-1] = A(j, i);  // avoid the copy ?
+	rank_one_update(sub_matrix(A, sub_interval(j+1, M-j-1),
+				 sub_interval(j+1, N-j-1)), c, conjugated(r));
+      }
+      ipvt[j] = j + 1;
+    }
+    return info;
+  }
+  
+  /** LU Solve : Solve equation Ax=b, given an LU factored matrix.*/
+  //  Thanks to Valient Gough for this routine!
+  template <typename DenseMatrix, typename VectorB, typename VectorX,
+	    typename Pvector>
+  void lu_solve(const DenseMatrix &LU, const Pvector& pvector, 
+		VectorX &x, const VectorB &b) {
+    typedef typename linalg_traits<DenseMatrix>::value_type T;
+    copy(b, x);
+    for(size_type i = 0; i < pvector.size(); ++i) {
+      size_type perm = pvector[i]-1;     // permutations stored in 1's offset
+      if(i != perm) { T aux = x[i]; x[i] = x[perm]; x[perm] = aux; }
+    }
+    /* solve  Ax = b  ->  LUx = b  ->  Ux = L^-1 b.                        */
+    lower_tri_solve(LU, x, true);
+    upper_tri_solve(LU, x, false);
+  }
+
+  template <typename DenseMatrix, typename VectorB, typename VectorX>
+  void lu_solve(const DenseMatrix &A, VectorX &x, const VectorB &b) {
+    typedef typename linalg_traits<DenseMatrix>::value_type T;
+    dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
+    std::vector<int> ipvt(mat_nrows(A));
+    gmm::copy(A, B);
+    size_type info = lu_factor(B, ipvt);
+    GMM_ASSERT1(!info, "Singular system, pivot = " << info);
+    lu_solve(B, ipvt, x, b);
+  }
+  
+  template <typename DenseMatrix, typename VectorB, typename VectorX,
+	    typename Pvector>
+  void lu_solve_transposed(const DenseMatrix &LU, const Pvector& pvector, 
+			   VectorX &x, const VectorB &b) {
+    typedef typename linalg_traits<DenseMatrix>::value_type T;
+    copy(b, x);
+    lower_tri_solve(transposed(LU), x, false);
+    upper_tri_solve(transposed(LU), x, true);
+    for(size_type i = pvector.size(); i > 0; --i) {
+      size_type perm = pvector[i-1]-1;    // permutations stored in 1's offset
+      if(i-1 != perm) { T aux = x[i-1]; x[i-1] = x[perm]; x[perm] = aux; }
+    }
+
+  }
+
+
+  ///@cond DOXY_SHOW_ALL_FUNCTIONS
+  template <typename DenseMatrixLU, typename DenseMatrix, typename Pvector>
+  void lu_inverse(const DenseMatrixLU& LU, const Pvector& pvector,
+		  DenseMatrix& AInv, col_major) {
+    typedef typename linalg_traits<DenseMatrixLU>::value_type T;
+    std::vector<T> tmp(pvector.size(), T(0));
+    std::vector<T> result(pvector.size());
+    for(size_type i = 0; i < pvector.size(); ++i) {
+      tmp[i] = T(1);
+      lu_solve(LU, pvector, result, tmp);
+      copy(result, mat_col(AInv, i));
+      tmp[i] = T(0);
+    }
+  }
+
+  template <typename DenseMatrixLU, typename DenseMatrix, typename Pvector>
+  void lu_inverse(const DenseMatrixLU& LU, const Pvector& pvector,
+		  DenseMatrix& AInv, row_major) {
+    typedef typename linalg_traits<DenseMatrixLU>::value_type T;
+    std::vector<T> tmp(pvector.size(), T(0));
+    std::vector<T> result(pvector.size());
+    for(size_type i = 0; i < pvector.size(); ++i) {
+      tmp[i] = T(1); // to be optimized !!
+      // on peut sur le premier tri solve reduire le systeme
+      // et peut etre faire un solve sur une serie de vecteurs au lieu
+      // de vecteur a vecteur (accumulation directe de l'inverse dans la
+      // matrice au fur et a mesure du calcul ... -> evite la copie finale
+      lu_solve_transposed(LU, pvector, result, tmp);
+      copy(result, mat_row(AInv, i));
+      tmp[i] = T(0);
+    }
+  }
+  ///@endcond  
+
+  /** Given an LU factored matrix, build the inverse of the matrix. */
+  template <typename DenseMatrixLU, typename DenseMatrix, typename Pvector>
+  void lu_inverse(const DenseMatrixLU& LU, const Pvector& pvector,
+		  const DenseMatrix& AInv_) {
+    DenseMatrix& AInv = const_cast<DenseMatrix&>(AInv_);
+    lu_inverse(LU, pvector, AInv, typename principal_orientation_type<typename
+	       linalg_traits<DenseMatrix>::sub_orientation>::potype());
+  }
+
+  /** Given a dense matrix, build the inverse of the matrix, and
+      return the determinant */
+  template <typename DenseMatrix>
+  typename linalg_traits<DenseMatrix>::value_type
+  lu_inverse(const DenseMatrix& A_) {
+    typedef typename linalg_traits<DenseMatrix>::value_type T;
+    DenseMatrix& A = const_cast<DenseMatrix&>(A_);
+    dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
+    std::vector<int> ipvt(mat_nrows(A));
+    gmm::copy(A, B);
+    size_type info = lu_factor(B, ipvt);
+    GMM_ASSERT1(!info, "Non invertible matrix, pivot = " << info);
+    lu_inverse(B, ipvt, A);
+    return lu_det(B, ipvt);
+  }
+
+  /** Compute the matrix determinant (via a LU factorization) */
+  template <typename DenseMatrixLU, typename Pvector>
+  typename linalg_traits<DenseMatrixLU>::value_type
+  lu_det(const DenseMatrixLU& LU, const Pvector &pvector) {
+    typedef typename linalg_traits<DenseMatrixLU>::value_type T;
+    T det(1);
+    for (size_type j = 0; j < std::min(mat_nrows(LU), mat_ncols(LU)); ++j)
+      det *= LU(j,j);
+    for(size_type i = 0; i < pvector.size(); ++i)
+      if (i != size_type(pvector[i]-1)) { det = -det; }
+    return det;
+  }
+
+  template <typename DenseMatrix>
+  typename linalg_traits<DenseMatrix>::value_type
+  lu_det(const DenseMatrix& A) {
+    typedef typename linalg_traits<DenseMatrix>::value_type T;
+    dense_matrix<T> B(mat_nrows(A), mat_ncols(A));
+    std::vector<int> ipvt(mat_nrows(A));
+    gmm::copy(A, B);
+    lu_factor(B, ipvt);
+    return lu_det(B, ipvt);
+  }
+
+}
+
+#endif
+

contrib/gmm/gmm_interface_bgeot.h

+// -*- c++ -*- (enables emacs c++ mode)
+//===========================================================================
+//
+// Copyright (C) 2002-2008 Yves Renard
+//
+// This file is a part of GETFEM++
+//
+// Getfem++  is  free software;  you  can  redistribute  it  and/or modify it
+// under  the  terms  of the  GNU  Lesser General Public License as published
+// by  the  Free Software Foundation;  either version 2.1 of the License,  or
+// (at your option) any later version.
+// This program  is  distributed  in  the  hope  that it will be useful,  but
+// WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+// License for more details.
+// You  should  have received a copy of the GNU Lesser General Public License
+// along  with  this program;  if not, write to the Free Software Foundation,
+// Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+// As a special exception, you may use this file as part of a free software
+// library without restriction.  Specifically, if other files instantiate
+// templates or use macros or inline functions from this file, or you compile
+// this file and link it with other files to produce an executable, this
+// file does not by itself cause the resulting executable to be covered by
+// the GNU General Public License.  This exception does not however
+// invalidate any other reasons why the executable file might be covered by
+// the GNU General Public License.
+//
+//===========================================================================
+
+/**@file gmm_interface_bgeot.h
+   @author  Yves Renard <Yves.Renard@insa-lyon.fr>
+   @date October 13, 2002.
+   @brief interface for bgeot::small_vector
+*/
+#ifndef GMM_INTERFACE_BGEOT_H__
+#define GMM_INTERFACE_BGEOT_H__
+
+
+namespace gmm {
+
+  /* ********************************************************************* */
+  /*		                                         	 	   */
+  /*		Traits for bgeot objects                     		   */
+  /*		                                         		   */
+  /* ********************************************************************* */
+
+  template <typename T> struct linalg_traits<bgeot::small_vector<T> > {
+    typedef bgeot::small_vector<T> this_type;
+    typedef this_type origin_type;
+    typedef linalg_false is_reference;
+    typedef abstract_vector linalg_type;
+    typedef T value_type;
+    typedef T& reference;
+    typedef typename this_type::iterator iterator;
+    typedef typename this_type::const_iterator const_iterator;
+    typedef abstract_dense storage_type;
+    typedef linalg_true index_sorted;
+    static size_type size(const this_type &v) { return v.size(); }
+    static iterator begin(this_type &v) { return v.begin(); }
+    static const_iterator begin(const this_type &v) { return v.begin(); }
+    static iterator end(this_type &v) { return v.end(); }
+    static const_iterator end(const this_type &v) { return v.end(); }
+    static origin_type* origin(this_type &v) { return &v; }
+    static const origin_type* origin(const this_type &v) { return &v; }
+    static void clear(origin_type* o, const iterator &it, const iterator &ite)
+    { std::fill(it, ite, value_type(0)); }
+    static void do_clear(this_type &v)
+    { std::fill(v.begin(), v.end(), value_type(0)); }
+    static value_type access(const origin_type *, const const_iterator &it,
+			     const const_iterator &, size_type i)
+    { return it[i]; }
+    static reference access(origin_type *, const iterator &it,
+			    const iterator &, size_type i)
+    { return it[i]; }
+    static void resize(this_type &v, size_type n) { v.resize(n); }
+  };
+
+}
+
+
+#endif //  GMM_INTERFACE_BGEOT_H__

contrib/gmm/gmm_kernel.h

+// -*- c++ -*- (enables emacs c++ mode)
+//===========================================================================
+//
+// Copyright (C) 2002-2008 Yves Renard
+//
+// This file is a part of GETFEM++
+//
+// Getfem++  is  free software;  you  can  redistribute  it  and/or modify it
+// under  the  terms  of the  GNU  Lesser General Public License as published
+// by  the  Free Software Foundation;  either version 2.1 of the License,  or
+// (at your option) any later version.
+// This program  is  distributed  in  the  hope  that it will be useful,  but
+// WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+// License for more details.
+// You  should  have received a copy of the GNU Lesser General Public License
+// along  with  this program;  if not, write to the Free Software Foundation,
+// Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+// As a special exception, you may use this file as part of a free software
+// library without restriction.  Specifically, if other files instantiate
+// templates or use macros or inline functions from this file, or you compile
+// this file and link it with other files to produce an executable, this
+// file does not by itself cause the resulting executable to be covered by
+// the GNU General Public License.  This exception does not however
+// invalidate any other reasons why the executable file might be covered by
+// the GNU General Public License.
+//
+//===========================================================================
+
+/**@file gmm_kernel.h 
+   @author  Yves Renard <Yves.Renard@insa-lyon.fr>
+   @date November 15, 2003.
+   @brief Include the base gmm files.
+ */
+
+#ifndef GMM_KERNEL_H__
+#define GMM_KERNEL_H__
+
+#include "gmm_def.h"
+#include "gmm_blas.h"
+#include "gmm_real_part.h"
+#include "gmm_interface.h"
+#include "gmm_sub_vector.h"
+#include "gmm_sub_matrix.h"
+#include "gmm_vector_to_matrix.h"
+#include "gmm_vector.h"
+#include "gmm_matrix.h"
+#include "gmm_tri_solve.h"
+#include "gmm_blas_interface.h"
+
+
+#endif //  GMM_KERNEL_H__