Line data    Source code

       1             : /* 
       2             :  * File:   OpReverseVector.cpp
       3             :  * Author: Pantelis Sopasakis
       4             :  * 
       5             :  * Created on September 15, 2015, 12:57 PM
       6             :  * 
       7             :  * ForBES is free software: you can redistribute it and/or modify
       8             :  * it under the terms of the GNU Lesser General Public License as published by
       9             :  * the Free Software Foundation, either version 3 of the License, or
      10             :  * (at your option) any later version.
      11             :  *  
      12             :  * ForBES is distributed in the hope that it will be useful,
      13             :  * but WITHOUT ANY WARRANTY; without even the implied warranty of
      14             :  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
      15             :  * GNU Lesser General Public License for more details.
      16             :  * 
      17             :  * You should have received a copy of the GNU Lesser General Public License
      18             :  * along with ForBES. If not, see <http://www.gnu.org/licenses/>.
      19             :  */
      20             : 
      21             : #include "OpReverseVector.h"
      22             : #include <algorithm>
      23             : #include <iterator>
      24             : #include <limits>
      25             : #include <cmath>
      26             : 
      27           1 : OpReverseVector::OpReverseVector() {
      28           1 :     m_vectorDim = 0;
      29           1 : }
      30             : 
      31          16 : OpReverseVector::~OpReverseVector() {
      32          16 : }
      33             : 
      34          12 : OpReverseVector::OpReverseVector(size_t n) {
      35          12 :     this->m_vectorDim = n;
      36          12 : }
      37             : 
      38          16 : int OpReverseVector::call(Matrix& y, double alpha, Matrix& x, double gamma) {
      39          16 :     if (y.getNrows() == 0) {
      40           2 :         y = Matrix(x.getNrows(), 1);
      41             :     }
      42          16 :     bool is_gamma_zero = (std::abs(gamma) < std::numeric_limits<double>::epsilon());
      43          16 :     size_t n = x.length();
      44         201 :     for (size_t i = 0; i < n / 2; ++i) {
      45             :         double temp;
      46         185 :         temp = x.get(n - i - 1, 0);
      47         185 :         if (is_gamma_zero) {
      48         185 :             y[n - i - 1] = alpha * x[i];
      49         185 :             y[i] = alpha * temp;
      50             :         } else {
      51           0 :             y[n - i - 1] = gamma * y[n - i - 1] + alpha * x[i];
      52           0 :             y[i] = gamma * y[i] + alpha * temp;
      53             :         }
      54             :     }
      55          16 :     if (n % 2 == 1) {
      56           5 :         size_t middle_idx = n / 2;
      57           5 :         if (is_gamma_zero) {
      58           5 :             y[middle_idx] = alpha * x[middle_idx];
      59             :         } else {
      60           0 :             y[middle_idx] = gamma * y[middle_idx] + alpha * x[middle_idx];
      61             :         }
      62             :     }
      63          16 :     return ForBESUtils::STATUS_OK;
      64             : }
      65             : 
      66           1 : int OpReverseVector::callAdjoint(Matrix& y, double alpha, Matrix& x, double gamma) {
      67           1 :     return call(y, alpha, x, gamma);
      68             : }
      69             : 
      70          14 : std::pair<size_t, size_t> OpReverseVector::dimensionIn() {
      71          14 :     return _VECTOR_OP_DIM(m_vectorDim);
      72             : }
      73             : 
      74          43 : std::pair<size_t, size_t> OpReverseVector::dimensionOut() {
      75          43 :     return _VECTOR_OP_DIM(m_vectorDim);
      76             : }
      77             : 
      78           1 : bool OpReverseVector::isSelfAdjoint() {
      79           1 :     return true;
      80             : }
      81             : 
      82             : 
      83             : 
      84             : 
      85             : 
      86             : 
      87             :