symmetric_matrix.cc

Go to the documentation of this file.

/*
 *            Copyright 2009-2020 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.
 *
 */
 
// Standard includes
#include <cassert>
#include <iostream>
 
// Local VOTCA includes
#include "votca/xtp/symmetric_matrix.h"
 
namespace votca {
namespace xtp {
 
Symmetric_Matrix::Symmetric_Matrix(const Eigen::MatrixXd& full)
    : dimension(full.rows()) {
  assert(full.rows() == full.cols() && "Input matrix not quadratic");
  data.resize((dimension + 1) * dimension / 2);
  for (Index i = 0; i < full.rows(); ++i) {
    for (Index j = 0; j <= i; ++j) {
      this->operator()(i, j) = full(i, j);
    }
  }
}
 
std::ostream& operator<<(std::ostream& out, const Symmetric_Matrix& a) {
 
  out << "[" << a.dimension << "," << a.dimension << "]\n";
  for (Index i = 0; i < a.dimension; ++i) {
    for (Index j = 0; j <= i; ++j) {
      out << a(i, j);
      if (i == j) {
        out << "\n";
      } else {
        out << " ";
      }
    }
  }
  return out;
}
 
double Symmetric_Matrix::TraceofProd(const Symmetric_Matrix& a) const {
  assert(data.size() == a.data.size() && "Matrices do not have same size");
  double result = 0.0;
 
  for (Index i = 0; i < dimension; ++i) {
    const Index index = (i * (i + 1)) / 2 + i;
    result += +data[index] * a.data[index];
  }
  for (Index i = 0; i < dimension; ++i) {
    const Index start = (i * (i + 1)) / 2;
    for (Index j = 0; j < i; ++j) {
      result += 2 * data[start + j] * a.data[start + j];
    }
  }
  return result;
}
 
void Symmetric_Matrix::AddtoEigenMatrix(Eigen::MatrixXd& full,
                                        double factor) const {
  for (Index j = 0; j < full.cols(); ++j) {
    const Index start = (j * (j + 1)) / 2;
    for (Index i = 0; i <= j; ++i) {
      full(i, j) += factor * data[start + i];
    }
 
    for (Index i = j + 1; i < full.rows(); ++i) {
      const Index index = (i * (i + 1)) / 2 + j;
      full(i, j) += factor * data[index];
    }
  }
  return;
}
 
void Symmetric_Matrix::AddtoEigenUpperMatrix(
    Eigen::SelfAdjointView<Eigen::MatrixXd, Eigen::Upper>& upper,
    double factor) const {
  for (Index j = 0; j < upper.cols(); ++j) {
    const Index start = (j * (j + 1)) / 2;
    for (Index i = 0; i <= j; ++i) {
      upper(i, j) += factor * data[start + i];
    }
  }
  return;
}
 
Eigen::MatrixXd Symmetric_Matrix::FullMatrix() const {
  Eigen::MatrixXd result = Eigen::MatrixXd(dimension, dimension);
  for (Index j = 0; j < result.cols(); ++j) {
    const Index start = (j * (j + 1)) / 2;
    for (Index i = 0; i <= j; ++i) {
      result(i, j) = data[start + i];
    }
 
    for (Index i = j + 1; i < result.rows(); ++i) {
      const Index index = (i * (i + 1)) / 2 + j;
      result(i, j) = data[index];
    }
  }
  return result;
}
 
void Symmetric_Matrix::FillFullMatrix(Eigen::MatrixXd& full) const {
  assert(full.rows() == dimension && full.cols() == dimension &&
         "FillFullMatrix needs a dim x dim matrix");
  // column j of the upper triangle is stored contiguously
  for (Index j = 0; j < dimension; ++j) {
    const Index start = (j * (j + 1)) / 2;
    full.col(j).head(j + 1) =
        Eigen::Map<const Eigen::VectorXd>(data.data() + start, j + 1);
  }
  full.triangularView<Eigen::StrictlyLower>() = full.transpose();
}
 
Eigen::MatrixXd Symmetric_Matrix::UpperMatrix() const {
  Eigen::MatrixXd result = Eigen::MatrixXd::Zero(dimension, dimension);
  for (Index j = 0; j < result.cols(); ++j) {
    const Index start = (j * (j + 1)) / 2;
    for (Index i = 0; i <= j; ++i) {
      result(i, j) = data[start + i];
    }
  }
  return result;
}
 
Index Symmetric_Matrix::index(Index i, Index j) const {
  Index index;
  if (i >= j) {
    index = (i * (i + 1)) / 2 + j;
  } else {
    index = (j * (j + 1)) / 2 + i;
  }
  return index;
}
 
}  // namespace xtp
}  // namespace votca