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