MatrixRepresentationEvaluationContainer.hpp

Go to the documentation of this file.

// This file is part of GQCG-GQCP.
//
// Copyright (C) 2017-2020  the GQCG developers
//
// GQCG-GQCP is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// GQCG-GQCP is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with GQCG-GQCP.  If not, see <http://www.gnu.org/licenses/>.
 
#pragma once
 
 
#include "Mathematical/Representation/SquareMatrix.hpp"
 
#include <Eigen/Sparse>
 
 
namespace GQCP {
 
 
template <typename _Matrix>
class MatrixRepresentationEvaluationContainer {
public:
    // The type that is used as a matrix representation, i.e. a container that stores the matrix representation.
    using Matrix = _Matrix;
 
    // The scalar type of one of the matrix elements: real or complex.
    using Scalar = typename Matrix::Scalar;
 
 
public:
    // The current index of the iterator, which is an address in a particular ONV basis.
    size_t index = 0;
 
    // The expected last position of the iterator, i.e. the dimension of the particular ONV basis.
    size_t end;
 
    // The matrix representation.
    Matrix matrix;
 
 
public:
    /*
     *  MARK: Constructors
     */
 
    MatrixRepresentationEvaluationContainer(const size_t dimension) :
        matrix {Matrix::Zero(dimension)},
        end {dimension} {}
 
 
    /*
     *  MARK: Additions and evaluations
     */
 
    void addDiagonally(const Scalar value) {
        this->matrix(this->index, this->index) += value;
    }
 
    void addColumnwise(const size_t column, const Scalar value) {
        this->matrix(this->index, column) += value;
    }
 
 
    void addRowwise(const size_t row, const Scalar value) {
        this->matrix(row, this->index) += value;
    }
 
 
    const Matrix& evaluation() const { return this->matrix; }
 
 
    /*
     *  MARK: Iterator
     */
 
    void increment() { this->index++; }
 
    bool isFinished() {
 
        if (this->index == this->end) {
            this->index = 0;
            return true;
        } else {
            return false;
        }
    }
};
 
 
template <typename _Scalar>
class MatrixRepresentationEvaluationContainer<Eigen::SparseMatrix<_Scalar>> {
public:
    // The scalar type of one of the matrix elements: real or complex.
    using Scalar = _Scalar;
 
 
public:
    size_t index = 0;  // Current position of the iterator in the dimension of the ONV basis.
    size_t end;        // The total dimension.
 
    Eigen::SparseMatrix<Scalar> matrix;                  // Matrix containing the evaluations.
    std::vector<Eigen::Triplet<Scalar>> triplet_vector;  // Vector which temporarily contains the added values.
 
 
public:
    /*
     *  CONSTRUCTORS
     */
 
    MatrixRepresentationEvaluationContainer(const size_t dimension) :
        matrix {Eigen::SparseMatrix<Scalar>(dimension, dimension)},
        end(dimension) {}
 
 
    /*
     *  PUBLIC METHODS
     */
 
    void addColumnwise(const size_t column, const Scalar value) { this->triplet_vector.emplace_back(this->index, column, value); }
 
    void addRowwise(const size_t row, const Scalar value) { this->triplet_vector.emplace_back(row, this->index, value); }
 
    void addToMatrix() {
        this->matrix.setFromTriplets(this->triplet_vector.begin(), this->triplet_vector.end());
        this->triplet_vector = {};
    }
 
    const Eigen::SparseMatrix<Scalar>& evaluation() const { return this->matrix; }
 
    void increment() { this->index++; }
 
    bool isFinished() {
        if (this->index == this->end) {
            this->index = 0;
            return true;
        } else {
            return false;
        }
    }
 
    void reserve(const size_t n) {
        this->triplet_vector.reserve(n);
        this->matrix.reserve(n);
    }
 
    const std::vector<Eigen::Triplet<Scalar>>& triplets() const { return triplet_vector; }
};
 
 
template <typename _Scalar>
class MatrixRepresentationEvaluationContainer<VectorX<_Scalar>> {
public:
    // The scalar type of one of the matrix elements: real or complex.
    using Scalar = _Scalar;
 
public:
    size_t index = 0;  // Current position of the iterator in the dimension of the ONV basis.
    size_t end;        // The total dimension.
 
    VectorX<Scalar> matvec;                     // Matvec containing the evaluations.
    const VectorX<Scalar>& coefficient_vector;  // Vector with which is multiplied.
    double sequential_double = 0;               // Double which temporarily contains the sum of added values, it corresponds to the value in the matvec of the current index and allows a for a single write access each iteration instead of multiple.
    double nonsequential_double = 0;            // Double gathered from the coefficient for non-sequential matvec additions, this corresponds to the value of the coefficient vector of the current index, it allows for a single read operation each iteration.
 
 
public:
    /*
     *  CONSTRUCTORS
     */
 
    MatrixRepresentationEvaluationContainer(const VectorX<Scalar>& coefficient_vector, const VectorX<Scalar>& diagonal) :
        end {static_cast<size_t>(coefficient_vector.rows())},
        coefficient_vector {coefficient_vector},
        matvec {diagonal.cwiseProduct(coefficient_vector)} {}
 
    MatrixRepresentationEvaluationContainer(const VectorX<Scalar>& coefficient_vector) :
        end {static_cast<size_t>(coefficient_vector.rows())},
        coefficient_vector {coefficient_vector},
        matvec {VectorX<Scalar>::Zero(coefficient_vector.rows())} {}
 
 
    /*
     *  PUBLIC METHODS
     */
 
    void addColumnwise(const size_t column, const Scalar value) { sequential_double += value * coefficient_vector(column); }
 
    void addRowwise(const size_t row, const Scalar value) { this->matvec(row) += value * this->nonsequential_double; }
 
    const VectorX<Scalar>& evaluation() const { return this->matvec; }
 
    void increment() {
        this->matvec(this->index) += this->sequential_double;
        this->sequential_double = 0;
        this->index++;
    }
 
    bool isFinished() {
        if (this->index == this->end) {
            this->index = 0;
            return true;
        } else {
            this->nonsequential_double = this->coefficient_vector(this->index);
            return false;
        }
    }
};
 
}  // namespace GQCP