codac2_IntvFullPivLU.h

Go to the documentation of this file.

 
#pragma once
 
#include <ostream>
#include "codac2_Matrix.h"
#include "codac2_IntervalMatrix.h"
#include "codac2_BoolInterval.h"
 
namespace codac2
{
  class IntvFullPivLU
  {
     public:
     explicit IntvFullPivLU(const Matrix &M);

     explicit IntvFullPivLU(const IntervalMatrix &M);

     BoolInterval is_injective() const;
     BoolInterval is_invertible() const;
     BoolInterval is_surjective() const;
     Interval determinant() const;
     Interval rank() const;
     Interval dimension_of_kernel() const;
     IntervalMatrix kernel() const;
     template <typename Derived> Derived
    image(const Eigen::MatrixBase<Derived> &M) const;
     IntervalMatrix solve(const IntervalMatrix &rhs) const;
 

     IntervalMatrix reconstructed_matrix() const;
     double max_pivot() const;
     
     const Eigen::FullPivLU<Matrix>::PermutationPType 
            &permutation_P () const;
     const Eigen::FullPivLU<Matrix>::PermutationQType 
            &permutation_Q () const;
     const Eigen::FullPivLU<Matrix> &eigen_LU() const;

     const IntervalMatrix &matrix_LU() const;
 
 
     private:
        Eigen::FullPivLU<Matrix> _LU;
 
        IntervalMatrix matrixLU_;   /* LU matrix */
     
        void compute_matrix_LU(const IntervalMatrix &M, double nonzero);
        static IntervalMatrix build_LU_bounds(const IntervalMatrix &E);
  };
 
inline const IntervalMatrix &IntvFullPivLU::matrix_LU() const {
     return this->matrixLU_;
}
inline const Eigen::FullPivLU<Matrix>::PermutationPType 
            &IntvFullPivLU::permutation_P() const {
    return this->_LU.permutationP();
}
inline const Eigen::FullPivLU<Matrix>::PermutationQType 
            &IntvFullPivLU::permutation_Q() const {
    return this->_LU.permutationQ();
}
inline const Eigen::FullPivLU<Matrix> &IntvFullPivLU::eigen_LU() const {
    return this->_LU;
}
 
template<typename Derived>
inline  Derived IntvFullPivLU::image
    (const Eigen::MatrixBase<Derived> &M) const
{
   int rk = this->rank().lb();
   if (rk==0) {
      return Derived::Zero(M.rows(),1); 
        /* NdDamien : où est-ce qu'on dit que Derived est une
           matrice ??? je crois que je hais le C++ ;)  */
   }
   Derived ret = 
        Derived::Zero(M.rows(),rk);
   Index p = 0;
   Index dim = std::min(matrixLU_.rows(),matrixLU_.cols());
   auto Q = this->_LU.permutationQ();
   for (Index i = 0; i<dim; i++) {
      if (!matrixLU_(i,i).contains(0.0)) {
         ret.col(p) = M.col(Q.indices().coeff(i));
         p++;
      }
   }
   return ret;
}
 
 
}