c_csr_matrix.cpp

Go to the documentation of this file.

/*
 *  SPDX-FileCopyrightText: Copyright 2021, Siavash Ameli <sameli@berkeley.edu>
 *  SPDX-License-Identifier: BSD-3-Clause
 *  SPDX-FileType: SOURCE
 *
 *  This program is free software: you can redistribute it and/or modify it
 *  under the terms of the license found in the LICENSE.txt file in the root
 *  directory of this source tree.
 */
 
 
// =======
// Headers
// =======
 
#include "./c_csr_matrix.h"
#include <cstddef>  // NULL
#include "../_definitions/definitions.h"  // USE_OPENMP
#if defined(USE_OPENMP) && (USE_OPENMP == 1)
    #include <omp.h>  // omp_in_parallel
#endif
#include "../_c_arithmetics/c_arithmetics.h"  // c_arithmetics
#include "../_c_basic_algebra/c_matrix_operations.h"  // cMatrixOperations
 
 
// =============
// constructor 1
// =============
 
 
template <typename DataType>
cCSRMatrix<DataType>::cCSRMatrix():
    A_data(NULL),
    A_indices(NULL),
    A_index_pointer(NULL)
{
}
 
 
// =============
// constructor 2
// =============
 
 
template <typename DataType>
cCSRMatrix<DataType>::cCSRMatrix(
        const DataType* A_data_,
        const LongIndexType* A_indices_,
        const LongIndexType* A_index_pointer_,
        const LongIndexType num_rows_,
        const LongIndexType num_columns_,
        const FlagType A_is_symmetric_):
 
    // Base class constructors
    cLinearOperatorBase(num_rows_, num_columns_),
    cMatrix<DataType>(A_is_symmetric_),
 
    // Initializer list
    A_data(A_data_),
    A_indices(A_indices_),
    A_index_pointer(A_index_pointer_)
{
}
 
 
// ==========
// destructor
// ==========
 
 
template <typename DataType>
cCSRMatrix<DataType>::~cCSRMatrix()
{
}
 
 
// ==================
// is identity matrix
// ==================
 
 
template <typename DataType>
FlagType cCSRMatrix<DataType>::is_identity_matrix() const
{
    FlagType matrix_is_identity = 1;
    LongIndexType index_pointer;
    LongIndexType column;
    DataType matrix_element;
    const DataType diagonal = 1.0;
    const DataType off_diagonal = 0.0;
 
    // Check matrix element-wise
    #if defined(USE_OPENMP) && (USE_OPENMP == 1)
    #pragma omp parallel for \
        schedule(static) \
        if (!omp_in_parallel()) \
        default(none) \
        shared(matrix_is_identity, diagonal, off_diagonal) \
        private(index_pointer, column, matrix_element)
    #endif
    for (LongIndexType row=0; row < this->num_rows; ++row)
    {
        if (matrix_is_identity)
        {
            for (index_pointer=this->A_index_pointer[row];
                 index_pointer < this->A_index_pointer[row+1];
                 ++index_pointer)
            {
                column = this->A_indices[index_pointer];
 
                if (!((this->A_is_symmetric) && (column >= row)))
                {
                    matrix_element = this->A_data[index_pointer];
 
                    if (((row == column) && \
                         (!c_arithmetics::is_equal(matrix_element,
                                                   diagonal))) || \
                        ((row != column) && \
                         (!c_arithmetics::is_equal(matrix_element,
                                                   off_diagonal))))
                    {
                        #if defined(USE_OPENMP) && (USE_OPENMP == 1)
                        #pragma omp atomic write
                        #endif
                        matrix_is_identity = 0;
 
                        break;
                    }
                }
            }
        }
    }
 
    return matrix_is_identity;
}
 
 
// =======
// get nnz
// =======
 
 
template <typename DataType>
LongIndexType cCSRMatrix<DataType>::get_nnz() const
{
    return this->A_index_pointer[this->num_rows];
}
 
 
// ===
// dot
// ===
 
 
template <typename DataType>
void cCSRMatrix<DataType>::dot(
        const DataType* vector,
        DataType* product)
{
    cMatrixOperations<DataType>::csr_matvec(
            this->A_data,
            this->A_indices,
            this->A_index_pointer,
            vector,
            this->num_rows,
            product);
}
 
 
// ========
// dot plus
// ========
 
 
template <typename DataType>
void cCSRMatrix<DataType>::dot_plus(
        const DataType* vector,
        const DataType alpha,
        DataType* product)
{
    cMatrixOperations<DataType>::csr_matvec_plus(
            this->A_data,
            this->A_indices,
            this->A_index_pointer,
            vector,
            alpha,
            this->num_rows,
            product);
}
 
 
// =============
// transpose dot
// =============
 
 
template <typename DataType>
void cCSRMatrix<DataType>::transpose_dot(
        const DataType* vector,
        DataType* product)
{
    cMatrixOperations<DataType>::csr_transposed_matvec(
            this->A_data,
            this->A_indices,
            this->A_index_pointer,
            this->A_is_symmetric,
            vector,
            this->num_rows,
            this->num_columns,
            product);
}
 
 
// ==================
// transpose dot plus
// ==================
 
 
template <typename DataType>
void cCSRMatrix<DataType>::transpose_dot_plus(
        const DataType* vector,
        const DataType alpha,
        DataType* product)
{
    cMatrixOperations<DataType>::csr_transposed_matvec_plus(
            this->A_data,
            this->A_indices,
            this->A_index_pointer,
            this->A_is_symmetric,
            vector,
            alpha,
            this->num_rows,
            product);
}
 
 
// ===============================
// Explicit template instantiation
// ===============================
 
template class cCSRMatrix<float>;
template class cCSRMatrix<double>;
template class cCSRMatrix<long double>;