PetscMatTools.hpp

/*

Copyright (C) University of Oxford, 2005-2011

University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.

This file is part of Chaste.

Chaste 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 2.1 of the License, or
(at your option) any later version.

Chaste 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. The offer of Chaste under the terms of the
License is subject to the License being interpreted in accordance with
English Law and subject to any action against the University of Oxford
being under the jurisdiction of the English Courts.

You should have received a copy of the GNU Lesser General Public License
along with Chaste. If not, see <http://www.gnu.org/licenses/>.

*/

#ifndef _PETSCMATTOOLS_HPP_
#define _PETSCMATTOOLS_HPP_

#include "UblasMatrixInclude.hpp" // needs to be 'first'
#include "PetscTools.hpp"
#include <vector>
#include <petscvec.h>
#include <petscmat.h>

class PetscMatTools
{
public:

    static void SetElement(Mat matrix, PetscInt row, PetscInt col, double value);

    static void AddToElement(Mat matrix, PetscInt row, PetscInt col, double value);

    static void Finalise(Mat matrix);

    static void SwitchWriteMode(Mat matrix);

    static void Display(Mat matrix);

    static void SetRow(Mat matrix, PetscInt row, double value);

    static void ZeroRowsWithValueOnDiagonal(Mat matrix, std::vector<unsigned>& rRows, double diagonalValue);

    static void ZeroRowsAndColumnsWithValueOnDiagonal(Mat matrix, std::vector<unsigned>& rRowColIndices, double diagonalValue);

    static void ZeroColumn(Mat matrix, PetscInt col);

    static void Zero(Mat matrix);

    static unsigned GetSize(Mat matrix);

    static void GetOwnershipRange(Mat matrix, PetscInt& lo, PetscInt& hi);

    static double GetElement(Mat matrix, PetscInt row, PetscInt col);

    static void SetOption(Mat matrix, MatOption option);

    static Vec GetMatrixRowDistributed(Mat matrix, unsigned rowIndex);

    static bool CheckEquality(const Mat mat1, const Mat mat2, double tol=1e-10);

    static bool CheckSymmetry(const Mat matrix, double tol=1e-10);

    template<size_t MATRIX_SIZE>
    static void AddMultipleValues(Mat matrix, unsigned* matrixRowAndColIndices, c_matrix<double, MATRIX_SIZE, MATRIX_SIZE>& rSmallMatrix)
    {
        PetscInt matrix_row_indices[MATRIX_SIZE];
        PetscInt num_rows_owned = 0;
        PetscInt global_row;
        PetscInt lo, hi;
        GetOwnershipRange(matrix, lo, hi);

        for (unsigned row = 0; row<MATRIX_SIZE; row++)
        {
            global_row = matrixRowAndColIndices[row];
            if (global_row >= lo && global_row < hi)
            {
                matrix_row_indices[num_rows_owned++] = global_row;
            }
        }

        if ( num_rows_owned == MATRIX_SIZE)
        {
            MatSetValues(matrix,
                         num_rows_owned,
                         matrix_row_indices,
                         MATRIX_SIZE,
                         (PetscInt*) matrixRowAndColIndices,
                         rSmallMatrix.data(),
                         ADD_VALUES);
        }
        else
        {
            // We need continuous data, if some of the rows do not belong to the processor their values
            // are not passed to MatSetValues
            double values[MATRIX_SIZE*MATRIX_SIZE];
            unsigned num_values_owned = 0;
            for (unsigned row = 0; row < MATRIX_SIZE; row++)
            {
                global_row = matrixRowAndColIndices[row];
                if (global_row >= lo && global_row < hi)
                {
                    for (unsigned col = 0; col < MATRIX_SIZE; col++)
                    {
                        values[num_values_owned++] = rSmallMatrix(row, col);
                    }
                }
            }

            MatSetValues(matrix,
                         num_rows_owned,
                         matrix_row_indices,
                         MATRIX_SIZE,
                         (PetscInt*) matrixRowAndColIndices,
                         values,
                         ADD_VALUES);
        }
    }
};

#endif //_PETSCMATTOOLS_HPP_