Line data    Source code

       1             : /* 
       2             :  * File:   OpComposition.cpp
       3             :  * Author: Pantelis Sopasakis
       4             :  * 
       5             :  * Created on September 14, 2015, 9:26 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 "OpComposition.h"
      22             : 
      23           7 : OpComposition::OpComposition(LinearOperator& A, LinearOperator& B) : LinearOperator(), m_A(A), m_B(B) {
      24             :     // check dimensions
      25           6 :     if (A.dimensionIn() != B.dimensionOut()) {
      26           1 :         throw std::invalid_argument("A and B have incompatible dimensions; AoB is not well defined.");
      27             :     }
      28           5 : }
      29             : 
      30           8 : OpComposition::~OpComposition() {
      31           8 : }
      32             : 
      33           3 : int OpComposition::call(Matrix& y, double alpha, Matrix& x, double gamma) {
      34           3 :     Matrix t(m_B.dimensionOut().first, m_B.dimensionOut().second);
      35           3 :     int status = m_B.call(t, 1.0, x, 0.0); // t = B(x)
      36           3 :     if (ForBESUtils::is_status_error(status)) {
      37           0 :         return status;
      38             :     }
      39           3 :     status = std::max(status, m_A.call(y, alpha, t, gamma));
      40           3 :     return status;
      41             : }
      42             : 
      43           1 : int OpComposition::callAdjoint(Matrix& y, double alpha, Matrix& x, double gamma) {
      44           1 :     Matrix t = m_A.callAdjoint(x);            
      45           1 :     return m_B.callAdjoint(y, alpha, t, gamma);
      46             : }
      47             : 
      48           5 : std::pair<size_t, size_t> OpComposition::dimensionIn() {
      49           5 :     return m_B.dimensionIn();
      50             : }
      51             : 
      52          23 : std::pair<size_t, size_t> OpComposition::dimensionOut() {
      53          23 :     return m_A.dimensionOut();
      54             : }
      55             : 
      56           1 : bool OpComposition::isSelfAdjoint() {
      57           1 :     return m_A.isSelfAdjoint() && m_B.isSelfAdjoint();
      58             : }
      59             : 
      60             :