Line data Source code
1 : /*
2 : * File: OpDCT2.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on September 15, 2015, 3:36 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 "OpDCT2.h"
22 : #include <cmath>
23 : #include <iostream>
24 :
25 :
26 : void update_y_helper_n_even(Matrix& y, double alpha, Matrix& x, double gamma, size_t n);
27 : void update_y_helper_n_odd(Matrix& y, double alpha, Matrix& x, double gamma, size_t n);
28 : double power_of_minus_one(size_t k);
29 :
30 0 : OpDCT2::OpDCT2() : LinearOperator(), m_dimension(_EMPTY_OP_DIM) {
31 :
32 0 : }
33 :
34 51 : OpDCT2::OpDCT2(size_t n) : m_dimension(_VECTOR_OP_DIM(n)) {
35 51 : }
36 :
37 100 : OpDCT2::~OpDCT2() {
38 100 : }
39 :
40 : static const double FOO[4] = {1.0, 0.0, -1.0, 0.0};
41 :
42 18831 : inline double power_of_minus_one(size_t k) {
43 18831 : return (k % 2 == 0) ? 1.0 : -1.0;
44 : }
45 :
46 577 : void update_y_helper_n_even(Matrix& y, double alpha, Matrix& x, double gamma, size_t n) {
47 577 : size_t nu = n / 2;
48 577 : bool apply_trick = (n >= 8) && (n % 4 == 0);
49 6488 : for (size_t k = 1; k < n; k++) {
50 5911 : double yk = 0.0;
51 5911 : if (apply_trick && k % 2 == 0) {
52 2090 : for (size_t i = 0; i < nu / 2; i++) {
53 : double aik;
54 1903 : size_t mu = k / 2;
55 1903 : aik = std::cos(M_PI * (static_cast<double> (i) + 0.5) * static_cast<double> (k) / static_cast<double> (n));
56 : yk += (
57 1903 : x[i]
58 1903 : + power_of_minus_one(mu) * (x[nu - i - 1] + x[nu + i])
59 1903 : + x[n - i - 1]
60 1903 : ) * aik;
61 187 : }
62 : } else {
63 46902 : for (size_t i = 0; i < nu; i++) {
64 : double aik;
65 41178 : aik = std::cos(M_PI * (static_cast<double> (i) + 0.5) * static_cast<double> (k) / static_cast<double> (n));
66 41178 : if (k % 2 == 1) {
67 24114 : yk += (x[i] - x[n - i - 1]) * aik;
68 : } else {
69 17064 : yk += (x[i] + x[n - i - 1]) * aik;
70 : }
71 : }
72 : }
73 5911 : y[k] = gamma * y[k] + alpha * yk;
74 : }
75 577 : }
76 :
77 23 : void update_y_helper_n_odd(Matrix& y, double alpha, Matrix& x, double gamma, size_t n) {
78 851 : for (size_t k = 1; k < n; k++) {
79 828 : double yk = 0.0;
80 17756 : for (size_t i = 0; i < n / 2; i++) {
81 : double aik;
82 16928 : aik = std::cos(M_PI * (static_cast<double> (i) + 0.5) * static_cast<double> (k) / static_cast<double> (n));
83 16928 : yk += (x[i] + power_of_minus_one(k) * x[n - 1 - i]) * aik;
84 : }
85 828 : yk += FOO[k % 4] * x[n / 2];
86 828 : y[k] = gamma * y[k] + alpha * yk;
87 : }
88 23 : }
89 :
90 600 : int OpDCT2::call(Matrix& y, double alpha, Matrix& x, double gamma) {
91 600 : size_t n = x.getNrows();
92 600 : double y0 = gamma * y[0];
93 7939 : for (size_t i = 0; i < n; i++) {
94 7339 : y0 += alpha * x[i];
95 : }
96 600 : y[0] = y0;
97 600 : if (n % 2 == 0) { // if n is even
98 577 : update_y_helper_n_even(y, alpha, x, gamma, n);
99 : } else {
100 23 : update_y_helper_n_odd(y, alpha, x, gamma, n);
101 : }
102 600 : return ForBESUtils::STATUS_OK;
103 : }
104 :
105 554 : int OpDCT2::callAdjoint(Matrix& y, double alpha, Matrix& x, double gamma) {
106 554 : size_t n = x.getNrows();
107 6844 : for (size_t k = 0; k < n; k++) {
108 6290 : double v = 0.0;
109 217390 : for (size_t i = 0; i < n; i++) {
110 : double aki;
111 211100 : aki = std::cos(M_PI * (static_cast<double> (k) + 0.5) * static_cast<double> (i) / static_cast<double> (n));
112 211100 : v += x[i] * aki;
113 : }
114 6290 : y[k] = gamma * y[k] + alpha * v;
115 : }
116 554 : return ForBESUtils::STATUS_OK;
117 : }
118 :
119 1115 : std::pair<size_t, size_t> OpDCT2::dimensionIn() {
120 1115 : return m_dimension;
121 : }
122 :
123 1207 : std::pair<size_t, size_t> OpDCT2::dimensionOut() {
124 1207 : return m_dimension;
125 : }
126 :
127 1 : bool OpDCT2::isSelfAdjoint() {
128 1 : return false;
129 3 : }
|
