codac2_IntvFullPivLU.h

Go to the documentation of this file.

 
#pragma once
 
#include <ostream>
#include "codac2_Matrix.h"
#include "codac2_Row.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;

        IntervalMatrix cokernel() const;

        template<typename Derived>
        Derived image(const Eigen::MatrixBase<Derived> &M) const;

        template<typename Derived>
        Derived coimage(const Eigen::MatrixBase<Derived> &M) const;

        IntervalMatrix solve(const IntervalMatrix& rhs) const;

        void solve(const IntervalMatrix& rhs, IntervalMatrix& B) 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 Row& transformation() const;

        const IntervalMatrix& matrix_LU() const;
 
 
    private:
 
        Eigen::FullPivLU<Matrix> _LU;
        Row transform; 
        IntervalMatrix matrixLU_; 
     
        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;
}
inline const Row &IntvFullPivLU::transformation() const {
     return this->transform;
}
 
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;
}
 
template<typename Derived>
inline  Derived IntvFullPivLU::coimage
  (const Eigen::MatrixBase<Derived> &M) const
{
   int rk = this->rank().lb();
   if (rk==0) {
      return Derived::Zero(1,M.cols()); 
   }
   Derived ret = 
    Derived::Zero(rk,M.cols());
   Index p = 0;
   Index dim = std::min(matrixLU_.rows(),matrixLU_.cols());
   auto P = this->_LU.permutationP();
   for (Index i = 0; i<dim; i++) {
      if (!matrixLU_(i,i).contains(0.0)) {
         ret.row(p) = M.row(P.indices().coeff(i));
         p++;
      }
   }
   return ret;
}
 
}