Matrix.h

 
// Description: A simple, parallel Matrix class with the utility needed to
//              support the basis generation methods of this library.  A
//              distributed Matrix has its rows distributed across processors.
 
#ifndef included_Matrix_h
#define included_Matrix_h
 
#include "Vector.h"
#include <vector>
#include <complex>
#include <memory>
#include <string>
#include <climits>
 
namespace CAROM {
 
class Matrix
{
public:
    Matrix();
 
    Matrix(
        int num_rows,
        int num_cols,
        bool distributed,
        bool randomized = false);
 
    Matrix(
        double* mat,
        int num_rows,
        int num_cols,
        bool distributed,
        bool copy_data = true);
 
    Matrix(
        const Matrix& other);
 
    ~Matrix();
 
    Matrix&
    operator = (
        const Matrix& rhs);
 
    Matrix&
    operator = (
        const double a);
 
    Matrix&
    operator += (
        const Matrix& rhs);
 
    Matrix&
    operator -= (
        const Matrix& rhs);
 
    void
    setSize(
        int num_rows,
        int num_cols)
    {
        // Check integer multiplication overflow
        if (num_rows > INT_MAX / num_cols)
            CAROM_ERROR("Matrix::setSize- new size exceeds maximum integer value!\n");
 
        int new_size = num_rows*num_cols;
        if (new_size > d_alloc_size) {
            if (!d_owns_data) {
                CAROM_ERROR("Can not reallocate externally owned storage.");
            }
            if (d_mat) {
                delete [] d_mat;
            }
 
            // Allocate new array and initialize all values to zero.
            d_mat = new double [new_size] {0.0};
            d_alloc_size = new_size;
        }
        d_num_rows = num_rows;
        d_num_cols = num_cols;
        if (d_distributed) {
            calculateNumDistributedRows();
        }
    }
 
    bool
    distributed() const
    {
        return d_distributed;
    }
 
    bool balanced() const;
 
    int
    numRows() const
    {
        return d_num_rows;
    }
 
    int
    numDistributedRows() const
    {
        if (!d_distributed) {
            return d_num_rows;
        }
        return d_num_distributed_rows;
    }
 
    int
    numColumns() const
    {
        return d_num_cols;
    }
 
    std::unique_ptr<Matrix>
    getFirstNColumns(int n) const;
 
    void
    getFirstNColumns(
        int n,
        Matrix& result) const;
 
    std::unique_ptr<Matrix>
    mult(
        const Matrix& other) const
    {
        Matrix* result = new Matrix(d_num_rows, other.d_num_cols, d_distributed);
        mult(other, *result);
        return std::unique_ptr<Matrix>(result);
    }
 
    void
    mult(
        const Matrix& other,
        Matrix& result) const;
 
    std::unique_ptr<Vector>
    mult(
        const Vector& other) const
    {
        Vector* result = new Vector(d_num_rows, d_distributed);
        mult(other, *result);
        return std::unique_ptr<Vector>(result);
    }
 
    void
    mult(
        const Vector& other,
        Vector& result) const;
 
    void
    pointwise_mult(
        int this_row,
        const Vector& other,
        Vector& result) const;
 
    void
    pointwise_mult(
        int this_row,
        Vector& other) const;
 
    std::unique_ptr<Matrix>
    elementwise_mult(
        const Matrix& other) const
    {
        Matrix *result = new Matrix(d_num_rows, d_num_cols, d_distributed);
        elementwise_mult(other, *result);
        return std::unique_ptr<Matrix>(result);
    }
 
    void
    elementwise_mult(
        const Matrix& other,
        Matrix& result) const;
 
    std::unique_ptr<Matrix>
    elementwise_square() const
    {
        Matrix *result = new Matrix(d_num_rows, d_num_cols, d_distributed);
        elementwise_square(*result);
        return std::unique_ptr<Matrix>(result);
    }
 
    void
    elementwise_square(
        Matrix& result) const;
 
    void
    multPlus(
        Vector& a,
        const Vector& b,
        double c) const;
 
    std::unique_ptr<Matrix>
    transposeMult(
        const Matrix& other) const
    {
        Matrix* result = new Matrix(d_num_cols, other.d_num_cols, false);
        transposeMult(other, *result);
        return std::unique_ptr<Matrix>(result);
    }
 
    void
    transposeMult(
        const Matrix& other,
        Matrix& result) const;
 
    std::unique_ptr<Vector>
    transposeMult(
        const Vector& other) const
    {
        Vector* result = new Vector(d_num_cols, false);
        transposeMult(other, *result);
        return std::unique_ptr<Vector>(result);
    }
 
    void
    transposeMult(
        const Vector& other,
        Vector& result) const;
 
    std::unique_ptr<Matrix>
    inverse() const
    {
        Matrix* result = new Matrix(d_num_rows, d_num_cols, false);
        inverse(*result);
        return std::unique_ptr<Matrix>(result);
    }
 
    void
    inverse(
        Matrix& result) const;
 
    void
    inverse();
 
    std::unique_ptr<Vector>
    getColumn(int column) const
    {
        Vector* result = new Vector(d_num_rows, true);
        getColumn(column, *result);
        return std::unique_ptr<Vector>(result);
    }
 
    void
    getColumn(int column,
              Vector& result) const;
 
    void transpose();
 
    void transposePseudoinverse();
 
    std::unique_ptr<Matrix>
    qr_factorize() const;
 
    void qr_factorize(std::vector<std::unique_ptr<Matrix>> & QR) const;
 
    void
    qrcp_pivots_transpose(int* row_pivot,
                          int* row_pivot_owner,
                          int  pivots_requested) const;
 
    void
    orthogonalize(bool double_pass = false, double zero_tol = 1.0e-15);
 
    void
    orthogonalize_last(int ncols = -1, bool double_pass = false,
                       double zero_tol = 1.0e-15);
 
    void
    rescale_rows_max();
 
    void
    rescale_cols_max();
 
    const double&
    item(int row,
         int col) const
    {
        CAROM_ASSERT((0 <= row) && (row < numRows()));
        CAROM_ASSERT((0 <= col) && (col < numColumns()));
        return d_mat[row*d_num_cols+col];
    }
 
    double&
    item(int row,
         int col)
    {
        CAROM_ASSERT((0 <= row) && (row < numRows()));
        CAROM_ASSERT((0 <= col) && (col < numColumns()));
        return d_mat[row*d_num_cols+col];
    }
 
    const double& operator() (int row, int col) const
    {
        return item(row, col);
    }
 
    double& operator() (int row, int col)
    {
        return item(row, col);
    }
 
    void print(const char * prefix) const;
 
    void write(const std::string& base_file_name) const;
 
    void read(const std::string& base_file_name);
 
    void local_read(const std::string& base_file_name, int rank);
 
    double *getData() const
    {
        return d_mat;
    }
 
    void distribute(const int &local_num_rows);
 
    void gather();
 
private:
    void
    calculateNumDistributedRows();
 
    void
    qrcp_pivots_transpose_serial(int* row_pivot,
                                 int* row_pivot_owner,
                                 int  pivots_requested) const;
 
    void
    qrcp_pivots_transpose_distributed(int* row_pivot,
                                      int* row_pivot_owner,
                                      int  pivots_requested) const;
 
    void
    qrcp_pivots_transpose_distributed_scalapack(int* row_pivot,
            int* row_pivot_owner,
            int  pivots_requested) const;
 
    void
    qrcp_pivots_transpose_distributed_elemental(int* row_pivot,
            int* row_pivot_owner,
            int  pivots_requested) const;
 
    void
    qrcp_pivots_transpose_distributed_elemental_balanced
    (int* row_pivot, int* row_pivot_owner, int pivots_requested) const;
 
    void
    qrcp_pivots_transpose_distributed_elemental_unbalanced
    (int* row_pivot, int* row_pivot_owner, int pivots_requested) const;
 
    void qr_factorize(bool computeR, std::vector<Matrix*> & QR) const;
 
    double* d_mat;
 
    int d_num_rows;
 
    int d_num_distributed_rows;
 
    int d_num_cols;
 
    int d_alloc_size;
 
    bool d_distributed;
 
    int d_num_procs;
 
    int d_rank;
 
    bool d_owns_data;
};
 
// NOTE(goxberry@gmail.com; oxberry1@llnl.gov): Passing by value is
// supposed to be more idiomatic C++11; consider passing by value
// instead of passing by reference.
Matrix outerProduct(const Vector &v, const Vector &w);
 
Matrix DiagonalMatrixFactory(const Vector &v);
 
Matrix IdentityMatrixFactory(const Vector &v);
 
struct EigenPair {
 
    Matrix* ev;
 
    std::vector<double> eigs;
};
 
struct ComplexEigenPair {
 
    Matrix* ev_real;
 
    Matrix* ev_imaginary;
 
    std::vector<std::complex<double>> eigs;
};
 
struct SerialSVDDecomposition
{
    Matrix* U;
 
    Vector* S;
 
    Matrix* V;
};
 
void SerialSVD(Matrix* A,
               Matrix* U,
               Vector* S,
               Matrix* V);
 
SerialSVDDecomposition SerialSVD(const Matrix & A);
 
EigenPair SymmetricRightEigenSolve(const Matrix & A);
 
ComplexEigenPair NonSymmetricRightEigenSolve(const Matrix & A);
 
std::unique_ptr<Matrix> SpaceTimeProduct(const CAROM::Matrix & As,
        const CAROM::Matrix & At, const CAROM::Matrix & Bs,
        const CAROM::Matrix & Bt,
        const std::vector<double> *tscale=NULL,
        const bool At0at0=false, const bool Bt0at0=false,
        const bool lagB=false, const bool skip0=false);
}
 
#endif