Line data    Source code

       1             : #include "PowerIterations.h"
       2             : #include "DataContainer.h"
       3             : #include "Error.h"
       4             : #include "LinearOperator.h"
       5             : #include <cmath>
       6             : 
       7             : #include <iostream>
       8             : 
       9             : namespace elsa
      10             : {
      11             :     template <class data_t>
      12             :     data_t powerIterations(const LinearOperator<data_t>& op, index_t maxiters, real_t rtol,
      13             :                            real_t atol)
      14          12 :     {
      15          12 :         if (op.getDomainDescriptor().getNumberOfCoefficients()
      16          12 :             != op.getRangeDescriptor().getNumberOfCoefficients()) {
      17           0 :             throw LogicError("powerIterations: Operator for power iterations must be symmetric");
      18           0 :         }
      19             : 
      20          12 :         DataContainer<data_t> u(
      21          12 :             op.getDomainDescriptor(),
      22          12 :             Vector_t<data_t>::Random(op.getDomainDescriptor().getNumberOfCoefficients()));
      23             : 
      24          12 :         data_t norm = u.l2Norm();
      25          12 :         data_t normOld = norm;
      26             : 
      27          12 :         u /= norm;
      28             : 
      29          24 :         for (index_t i = 0; i < maxiters; i++) {
      30          24 :             u = op.apply(u);
      31             : 
      32          24 :             norm = u.l2Norm();
      33          24 :             if (norm == 0.0) {
      34           0 :                 throw Error("powerIterations: Reached 0 vector after {} iterations", i + 1);
      35           0 :             }
      36             : 
      37          24 :             if (std::isinf(norm)) {
      38           0 :                 throw Error("powerIterations: Reached non-finite vector after {} iterations",
      39           0 :                             i + 1);
      40           0 :             }
      41             : 
      42          24 :             if (std::abs(norm - normOld) < atol + rtol * std::abs(norm)) {
      43          12 :                 break;
      44          12 :             } else {
      45          12 :                 u /= norm;
      46          12 :                 normOld = norm;
      47          12 :             }
      48          24 :         }
      49             : 
      50          12 :         return norm;
      51          12 :     }
      52             : 
      53             :     template float powerIterations(const LinearOperator<float>&, index_t, real_t, real_t);
      54             :     template double powerIterations(const LinearOperator<double>&, index_t, real_t, real_t);
      55             : } // namespace elsa