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