Line data    Source code

       1             : #include "JacobiPreconditioner.h"
       2             : #include "VolumeDescriptor.h"
       3             : #include "Scaling.h"
       4             : 
       5             : namespace elsa
       6             : {
       7             :     template <typename data_t>
       8           0 :     JacobiPreconditioner<data_t>::JacobiPreconditioner(const LinearOperator<data_t>& op,
       9             :                                                        bool inverse)
      10             :         : LinearOperator<data_t>{op.getDomainDescriptor(), op.getRangeDescriptor()},
      11           0 :           _inverseDiagonal(op.getDomainDescriptor(),
      12           0 :                            inverse ? 1 / diagonalFromOperator(op) : diagonalFromOperator(op))
      13             :     {
      14           0 :     }
      15             : 
      16             :     template <typename data_t>
      17           0 :     JacobiPreconditioner<data_t>::JacobiPreconditioner(const JacobiPreconditioner<data_t>& other)
      18           0 :         : LinearOperator<data_t>{*other._domainDescriptor, *other._rangeDescriptor},
      19           0 :           _inverseDiagonal{other._inverseDiagonal.getDomainDescriptor(),
      20           0 :                            other._inverseDiagonal.getScaleFactors()}
      21             :     {
      22           0 :     }
      23             : 
      24             :     template <typename data_t>
      25           0 :     auto JacobiPreconditioner<data_t>::cloneImpl() const -> JacobiPreconditioner<data_t>*
      26             :     {
      27           0 :         return new JacobiPreconditioner<data_t>(*this);
      28             :     }
      29             : 
      30             :     template <typename data_t>
      31           0 :     bool JacobiPreconditioner<data_t>::isEqual(const LinearOperator<data_t>& other) const
      32             :     {
      33           0 :         if (!LinearOperator<data_t>::isEqual(other))
      34           0 :             return false;
      35             : 
      36             :         // static_cast as type checked in base comparison
      37           0 :         auto otherJacobiPrecond = static_cast<const JacobiPreconditioner<data_t>*>(&other);
      38             : 
      39           0 :         if (_inverseDiagonal != otherJacobiPrecond->_inverseDiagonal)
      40           0 :             return false;
      41             : 
      42           0 :         return true;
      43             :     }
      44             : 
      45             :     template <typename data_t>
      46           0 :     void JacobiPreconditioner<data_t>::applyImpl(const DataContainer<data_t>& x,
      47             :                                                  DataContainer<data_t>& Ax) const
      48             :     {
      49           0 :         _inverseDiagonal.apply(x, Ax);
      50           0 :     }
      51             : 
      52             :     template <typename data_t>
      53           0 :     void JacobiPreconditioner<data_t>::applyAdjointImpl(const DataContainer<data_t>& y,
      54             :                                                         DataContainer<data_t>& Aty) const
      55             :     {
      56           0 :         _inverseDiagonal.applyAdjoint(y, Aty);
      57           0 :     }
      58             : 
      59             :     template <typename data_t>
      60             :     DataContainer<data_t>
      61           0 :         JacobiPreconditioner<data_t>::diagonalFromOperator(const LinearOperator<data_t>& op)
      62             :     {
      63           0 :         DataContainer<data_t> e(op.getDomainDescriptor());
      64           0 :         e = 0;
      65           0 :         DataContainer<data_t> diag(op.getDomainDescriptor());
      66           0 :         for (index_t i = 0; i < e.getSize(); i++) {
      67           0 :             e[i] = 1;
      68           0 :             diag[i] = op.apply(e)[i];
      69           0 :             e[i] = 0;
      70             :         }
      71             : 
      72           0 :         return diag;
      73           0 :     }
      74             : 
      75             :     // ------------------------------------------
      76             :     // explicit template instantiation
      77             :     template class JacobiPreconditioner<float>;
      78             :     template class JacobiPreconditioner<double>;
      79             : } // namespace elsa