/*
 * Copyright 2009-2019 The VOTCA Development Team (http://www.votca.org)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

#include <iostream>
#include <votca/xtp/matrixfreeoperator.h>

namespace votca {
namespace xtp {

Eigen::RowVectorXd MatrixFreeOperator::OperatorRow(Index) const {
  return Eigen::RowVectorXd::Zero(0);
}

Eigen::MatrixXd MatrixFreeOperator::OperatorBlock(Index, Index) const {
  return Eigen::MatrixXd::Zero(0, 0);
}

Eigen::VectorXd MatrixFreeOperator::diagonal() const {
  Eigen::VectorXd D = Eigen::VectorXd::Zero(_size);
  if (useRow()) {
#pragma omp parallel for schedule(guided)
    for (Index i = 0; i < _size; i++) {
      Eigen::RowVectorXd row_data = this->OperatorRow(i);
      D(i) = row_data(i);
    }
  }

  if (useBlock()) {
    Index blocksize = getBlocksize();
    if (size() % blocksize != 0) {
      throw std::runtime_error("blocksize is not a multiple of matrix size");
    }
    Index blocks = size() / blocksize;

// this is inefficient if blocks<num_ofthreads
#pragma omp parallel for schedule(guided)
    for (Index i = 0; i < blocks; i++) {
      Eigen::MatrixXd block = OperatorBlock(i, i);
      D.segment(i * blocksize, blocksize) += block.diagonal();
    }
  }

  return D;
}

// get the full matrix if we have to
Eigen::MatrixXd MatrixFreeOperator::get_full_matrix() const {
  Eigen::MatrixXd matrix = Eigen::MatrixXd::Zero(_size, _size);
  if (useRow()) {
#pragma omp parallel for schedule(guided)
    for (Index i = 0; i < _size; i++) {
      matrix.row(i) = this->OperatorRow(i);
    }
  }

  if (useBlock()) {
    Index blocksize = getBlocksize();
    if (size() % blocksize != 0) {
      throw std::runtime_error("blocksize is not a multiple of matrix size");
    }
    Index blocks = size() / blocksize;

// this is inefficient if blocks<num_ofthreads
#pragma omp parallel for schedule(guided)
    for (Index i_row = 0; i_row < blocks; i_row++) {
      for (Index i_col = 0; i_col < blocks; i_col++) {
        matrix.block(i_row * blocksize, i_col * blocksize, blocksize,
                     blocksize) += OperatorBlock(i_row, i_col);
      }
    }
  }

  return matrix;
}

// get the size
Index MatrixFreeOperator::size() const { return this->_size; }

// set the size
void MatrixFreeOperator::set_size(Index size) { this->_size = size; }

}  // namespace xtp
}  // namespace votca