Line data    Source code

       1             : /* 
       2             :  * File:   OpDCT3.cpp
       3             :  * Author: Pantelis Sopasakis
       4             :  * 
       5             :  * Created on September 15, 2015, 6:40 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 "OpDCT3.h"
      22             : #include <cmath>
      23             : #include <limits>
      24             : 
      25           3 : OpDCT3::OpDCT3(size_t dimension) :
      26             : LinearOperator(),
      27           3 : m_dimension(_VECTOR_OP_DIM(dimension)) {
      28           3 : }
      29             : 
      30             : //LCOV_EXCL_START
      31             : OpDCT3::OpDCT3() : m_dimension(_EMPTY_OP_DIM) {
      32             : }
      33             : //LCOV_EXCL_STOP
      34             : 
      35           6 : OpDCT3::~OpDCT3() {
      36           6 : }
      37             : 
      38          53 : int OpDCT3::call(Matrix& y, double alpha, Matrix& x, double gamma) {
      39          53 :     size_t n = x.length();
      40          53 :     if (m_dimension.first != 0 && n != m_dimension.first) {
      41           0 :         throw std::invalid_argument("x-dimension is invalid");
      42             :     }
      43          53 :     double x0_2 = x[0] / 2.0;
      44         953 :     for (size_t k = 0; k < n; k++) {
      45         900 :         double yk = x0_2;
      46       18750 :         for (size_t i = 1; i < n; i++) {
      47       17850 :             yk += (x[i] * std::cos(i * M_PI * (k + 0.5) / n));
      48             :         }
      49         900 :         y[k] = gamma * y[k] + alpha * yk;
      50             :     }
      51          53 :     return ForBESUtils::STATUS_OK;
      52             : }
      53             : 
      54          53 : int OpDCT3::callAdjoint(Matrix& y, double alpha, Matrix& x, double gamma) {
      55          53 :     size_t n = x.length();
      56          53 :     if (m_dimension.first != 0 && n != m_dimension.first) {
      57           0 :         throw std::invalid_argument("x-dimension is invalid");
      58             :     }
      59          53 :     double tk = 0.0;
      60         983 :     for (size_t i = 0; i < n; i++) {
      61         930 :         tk += x[i] / 2.0;
      62             :     }
      63          53 :     y[0] = gamma * y[0] + alpha * tk;
      64         930 :     for (size_t k = 1; k < n; k++) {
      65         877 :         tk = 0.0;
      66       21997 :         for (size_t i = 0; i < n; i++) {
      67       21120 :             tk += (x[i] * std::cos(k * M_PI * (i + 0.5) / n));
      68             :         }
      69         877 :         y[k] = gamma * y[k] + alpha * tk;
      70             :     }
      71          53 :     return ForBESUtils::STATUS_OK;
      72             : }
      73             : 
      74         114 : std::pair<size_t, size_t> OpDCT3::dimensionIn() {
      75         114 :     return m_dimension;
      76             : }
      77             : 
      78         114 : std::pair<size_t, size_t> OpDCT3::dimensionOut() {
      79         114 :     return m_dimension;
      80             : }
      81             : 
      82           1 : bool OpDCT3::isSelfAdjoint() {
      83           1 :     return false;
      84             : }
      85             : 
      86             : 
      87             : 
      88             : 
      89             : 
      90             :