c_matrix_operations.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_matrix_operations.h"
#include "../_definitions/definitions.h"  // USE_CBLAS
 
#if (USE_CBLAS == 1)
    #include "./cblas_interface.h"
#endif
 
 
// ============
// dense matvec
// ============
 
 
template <typename DataType>
void cMatrixOperations<DataType>::dense_matvec(
        const DataType* A,
        const DataType* b,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        const FlagType A_is_row_major,
        DataType* c)
{
    #if (USE_CBLAS == 1)
 
    // Using OpenBlas
    CBLAS_LAYOUT layout;
    if (A_is_row_major)
    {
        layout = CblasRowMajor;
    }
    else
    {
        layout = CblasColMajor;
    }
 
    CBLAS_TRANSPOSE transpose = CblasNoTrans;
    int lda = num_rows;
    int incb = 1;
    int incc = 1;
    DataType alpha = 1.0;
    DataType beta = 0.0;
 
    cblas_interface::xgemv(layout, transpose, num_rows, num_columns, alpha, A,
                           lda, b, incb, beta, c, incc);
 
    #else
 
    // Not using OpenBlas
    LongIndexType j;
    long double sum;
    LongIndexType chunk = 5;
    LongIndexType num_columns_chunked = num_columns - (num_columns % chunk);
 
    // Determine major order of A
    if (A_is_row_major)
    {
        // For row major (C ordering) matrix A
        for (LongIndexType i=0; i < num_rows; ++i)
        {
            sum = 0.0;
            for (j=0; j < num_columns_chunked; j+= chunk)
            {
                sum += A[i*num_columns + j] * b[j] +
                       A[i*num_columns + j+1] * b[j+1] +
                       A[i*num_columns + j+2] * b[j+2] +
                       A[i*num_columns + j+3] * b[j+3] +
                       A[i*num_columns + j+4] * b[j+4];
            }
 
            for (j= num_columns_chunked; j < num_columns; ++j)
            {
                sum += A[i*num_columns + j] * b[j];
            }
 
            c[i] = static_cast<DataType>(sum);
        }
    }
    else
    {
        // For column major (Fortran ordering) matrix A
        for (LongIndexType i=0; i < num_rows; ++i)
        {
            sum = 0.0;
            for (j=0; j < num_columns; ++j)
            {
                sum += A[i + num_rows*j] * b[j];
            }
            c[i] = static_cast<DataType>(sum);
        }
    }
 
    #endif
}
 
 
// =================
// dense matvec plus
// =================
 
 
template <typename DataType>
void cMatrixOperations<DataType>::dense_matvec_plus(
        const DataType* A,
        const DataType* b,
        const DataType alpha,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        const FlagType A_is_row_major,
        DataType* c)
{
    if (alpha == 0.0)
    {
        return;
    }
 
    LongIndexType j;
    long double sum;
    LongIndexType chunk = 5;
    LongIndexType num_columns_chunked = num_columns - (num_columns % chunk);
 
    // Determine major order of A
    if (A_is_row_major)
    {
        // For row major (C ordering) matrix A
        for (LongIndexType i=0; i < num_rows; ++i)
        {
            sum = 0.0;
            for (j=0; j < num_columns_chunked; j+= chunk)
            {
                sum += A[i*num_columns + j] * b[j] +
                       A[i*num_columns + j+1] * b[j+1] +
                       A[i*num_columns + j+2] * b[j+2] +
                       A[i*num_columns + j+3] * b[j+3] +
                       A[i*num_columns + j+4] * b[j+4];
            }
 
            for (j= num_columns_chunked; j < num_columns; ++j)
            {
                sum += A[i*num_columns + j] * b[j];
            }
 
            c[i] += alpha * static_cast<DataType>(sum);
        }
    }
    else
    {
        // For column major (Fortran ordering) matrix A
        for (LongIndexType i=0; i < num_rows; ++i)
        {
            sum = 0.0;
            for (j=0; j < num_columns; ++j)
            {
                sum += A[i + num_rows*j] * b[j];
            }
            c[i] += alpha* static_cast<DataType>(sum);
        }
    }
}
 
 
// =======================
// dense transposed matvec
// =======================
 
 
template <typename DataType>
void cMatrixOperations<DataType>::dense_transposed_matvec(
        const DataType* A,
        const DataType* b,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        const FlagType A_is_row_major,
        DataType* c)
{
    LongIndexType i;
    long double sum;
    LongIndexType chunk = 5;
    LongIndexType num_rows_chunked = num_rows - (num_rows % chunk);
 
    // Determine major order of A
    if (A_is_row_major)
    {
        // For row major (C ordering) matrix A
        for (LongIndexType j=0; j < num_columns; ++j)
        {
            sum = 0.0;
            for (i=0; i < num_rows; ++i)
            {
                sum += A[i*num_columns + j] * b[i];
            }
            c[j] = static_cast<DataType>(sum);
        }
    }
    else
    {
        // For column major (Fortran ordering) matrix A
        for (LongIndexType j=0; j < num_columns; ++j)
        {
            sum = 0.0;
            for (i=0; i < num_rows_chunked; i += chunk)
            {
                sum += A[i + num_rows*j] * b[i] +
                       A[i+1 + num_rows*j] * b[i+1] +
                       A[i+2 + num_rows*j] * b[i+2] +
                       A[i+3 + num_rows*j] * b[i+3] +
                       A[i+4 + num_rows*j] * b[i+4];
            }
 
            for (i=num_rows_chunked; i < num_rows; ++i)
            {
                sum += A[i + num_rows*j] * b[i];
            }
 
            c[j] = static_cast<DataType>(sum);
        }
    }
}
 
 
// ============================
// dense transposed matvec plus
// ============================
 
 
template <typename DataType>
void cMatrixOperations<DataType>::dense_transposed_matvec_plus(
        const DataType* A,
        const DataType* b,
        const DataType alpha,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        const FlagType A_is_row_major,
        DataType* c)
{
    if (alpha == 0.0)
    {
        return;
    }
 
    LongIndexType i;
    long double sum;
    LongIndexType chunk = 5;
    LongIndexType num_rows_chunked = num_rows - (num_rows % chunk);
 
    // Determine major order of A
    if (A_is_row_major)
    {
        // For row major (C ordering) matrix A
        for (LongIndexType j=0; j < num_columns; ++j)
        {
            sum = 0.0;
            for (i=0; i < num_rows; ++i)
            {
                sum += A[i*num_columns + j] * b[i];
            }
            c[j] += alpha * static_cast<DataType>(sum);
        }
    }
    else
    {
        // For column major (Fortran ordering) matrix A
        for (LongIndexType j=0; j < num_columns; ++j)
        {
            sum = 0.0;
            for (i=0; i < num_rows_chunked; i += chunk)
            {
                sum += A[i + num_rows*j] * b[i] +
                       A[i+1 + num_rows*j] * b[i+1] +
                       A[i+2 + num_rows*j] * b[i+2] +
                       A[i+3 + num_rows*j] * b[i+3] +
                       A[i+4 + num_rows*j] * b[i+4];
            }
 
            for (i=num_rows_chunked; i < num_rows; ++i)
            {
                sum += A[i + num_rows*j] * b[i];
            }
 
            c[j] += alpha * static_cast<DataType>(sum);
        }
    }
}
 
 
// ==========
// csr matvec
// ==========
 
 
template <typename DataType>
void cMatrixOperations<DataType>::csr_matvec(
        const DataType* A_data,
        const LongIndexType* A_column_indices,
        const LongIndexType* A_index_pointer,
        const DataType* b,
        const LongIndexType num_rows,
        DataType* c)
{
    LongIndexType index_pointer;
    LongIndexType row;
    LongIndexType column;
    long double sum;
 
    for (row=0; row < num_rows; ++row)
    {
        sum = 0.0;
        for (index_pointer=A_index_pointer[row];
             index_pointer < A_index_pointer[row+1];
             ++index_pointer)
        {
            column = A_column_indices[index_pointer];
            sum += A_data[index_pointer] * b[column];
        }
        c[row] = static_cast<DataType>(sum);
    }
}
 
 
// ===============
// csr matvec plus
// ===============
 
 
template <typename DataType>
void cMatrixOperations<DataType>::csr_matvec_plus(
        const DataType* A_data,
        const LongIndexType* A_column_indices,
        const LongIndexType* A_index_pointer,
        const DataType* b,
        const DataType alpha,
        const LongIndexType num_rows,
        DataType* c)
{
    if (alpha == 0.0)
    {
        return;
    }
 
    LongIndexType index_pointer;
    LongIndexType row;
    LongIndexType column;
    long double sum;
 
    for (row=0; row < num_rows; ++row)
    {
        sum = 0.0;
        for (index_pointer=A_index_pointer[row];
             index_pointer < A_index_pointer[row+1];
             ++index_pointer)
        {
            column = A_column_indices[index_pointer];
            sum += A_data[index_pointer] * b[column];
        }
        c[row] += alpha * static_cast<DataType>(sum);
    }
}
 
 
// =====================
// csr transposed matvec
// =====================
 
 
template <typename DataType>
void cMatrixOperations<DataType>::csr_transposed_matvec(
        const DataType* A_data,
        const LongIndexType* A_column_indices,
        const LongIndexType* A_index_pointer,
        const DataType* b,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        DataType* c)
{
    LongIndexType index_pointer;
    LongIndexType row;
    LongIndexType column;
 
    // Initialize output to zero
    for (column=0; column < num_columns; ++column)
    {
        c[column] = 0.0;
    }
 
    for (row=0; row < num_rows; ++row)
    {
        for (index_pointer=A_index_pointer[row];
             index_pointer < A_index_pointer[row+1];
             ++index_pointer)
        {
            column = A_column_indices[index_pointer];
            c[column] += A_data[index_pointer] * b[row];
        }
    }
}
 
 
// ==========================
// csr transposed matvec plus
// ==========================
 
 
template <typename DataType>
void cMatrixOperations<DataType>::csr_transposed_matvec_plus(
        const DataType* A_data,
        const LongIndexType* A_column_indices,
        const LongIndexType* A_index_pointer,
        const DataType* b,
        const DataType alpha,
        const LongIndexType num_rows,
        DataType* c)
{
    if (alpha == 0.0)
    {
        return;
    }
 
    LongIndexType index_pointer;
    LongIndexType row;
    LongIndexType column;
 
    for (row=0; row < num_rows; ++row)
    {
        for (index_pointer=A_index_pointer[row];
             index_pointer < A_index_pointer[row+1];
             ++index_pointer)
        {
            column = A_column_indices[index_pointer];
            c[column] += alpha * A_data[index_pointer] * b[row];
        }
    }
}
 
 
// ==========
// csc matvec
// ==========
 
 
template <typename DataType>
void cMatrixOperations<DataType>::csc_matvec(
        const DataType* A_data,
        const LongIndexType* A_row_indices,
        const LongIndexType* A_index_pointer,
        const DataType* b,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        DataType* c)
{
    LongIndexType index_pointer;
    LongIndexType row;
    LongIndexType column;
 
    // Initialize output to zero
    for (row=0; row < num_rows; ++row)
    {
        c[row] = 0.0;
    }
 
    for (column=0; column < num_columns; ++column)
    {
        for (index_pointer=A_index_pointer[column];
             index_pointer < A_index_pointer[column+1];
             ++index_pointer)
        {
            row = A_row_indices[index_pointer];
            c[row] += A_data[index_pointer] * b[column];
        }
    }
}
 
 
// ===============
// csc matvec plus
// ===============
 
 
template <typename DataType>
void cMatrixOperations<DataType>::csc_matvec_plus(
        const DataType* A_data,
        const LongIndexType* A_row_indices,
        const LongIndexType* A_index_pointer,
        const DataType* b,
        const DataType alpha,
        const LongIndexType num_columns,
        DataType* c)
{
    if (alpha == 0.0)
    {
        return;
    }
 
    LongIndexType index_pointer;
    LongIndexType row;
    LongIndexType column;
 
    for (column=0; column < num_columns; ++column)
    {
        for (index_pointer=A_index_pointer[column];
             index_pointer < A_index_pointer[column+1];
             ++index_pointer)
        {
            row = A_row_indices[index_pointer];
            c[row] += alpha * A_data[index_pointer] * b[column];
        }
    }
}
 
 
// =====================
// csc transposed matvec
// =====================
 
 
template <typename DataType>
void cMatrixOperations<DataType>::csc_transposed_matvec(
        const DataType* A_data,
        const LongIndexType* A_row_indices,
        const LongIndexType* A_index_pointer,
        const DataType* b,
        const LongIndexType num_columns,
        DataType* c)
{
    LongIndexType index_pointer;
    LongIndexType row;
    LongIndexType column;
    long double sum;
 
    for (column=0; column < num_columns; ++column)
    {
        sum = 0.0;
        for (index_pointer=A_index_pointer[column];
             index_pointer < A_index_pointer[column+1];
             ++index_pointer)
        {
            row = A_row_indices[index_pointer];
            sum += A_data[index_pointer] * b[row];
        }
        c[column] = static_cast<DataType>(sum);
    }
}
 
 
// ==========================
// csc transposed matvec plus
// ==========================
 
 
template <typename DataType>
void cMatrixOperations<DataType>::csc_transposed_matvec_plus(
        const DataType* A_data,
        const LongIndexType* A_row_indices,
        const LongIndexType* A_index_pointer,
        const DataType* b,
        const DataType alpha,
        const LongIndexType num_columns,
        DataType* c)
{
    if (alpha == 0.0)
    {
        return;
    }
 
    LongIndexType index_pointer;
    LongIndexType row;
    LongIndexType column;
    long double sum;
 
    for (column=0; column < num_columns; ++column)
    {
        sum = 0.0;
        for (index_pointer=A_index_pointer[column];
             index_pointer < A_index_pointer[column+1];
             ++index_pointer)
        {
            row = A_row_indices[index_pointer];
            sum += A_data[index_pointer] * b[row];
        }
        c[column] += static_cast<DataType>(alpha * sum);
    }
}
 
 
// ==================
// create band matrix
// ==================
 
 
template <typename DataType>
void cMatrixOperations<DataType>::create_band_matrix(
        const DataType* diagonals,
        const DataType* supdiagonals,
        const IndexType non_zero_size,
        const FlagType tridiagonal,
        DataType** matrix)
{
    for (IndexType j=0; j < non_zero_size; ++j)
    {
        // Diagonals
        matrix[j][j] = diagonals[j];
 
        // Off diagonals
        if (j < non_zero_size-1)
        {
            // Sup-diagonal
            matrix[j][j+1] = supdiagonals[j];
 
            // Sub-diagonal, making symmetric tri-diagonal matrix
            if (tridiagonal)
            {
                matrix[j+1][j] = supdiagonals[j];
            }
        }
    }
}
 
 
// ===============================
// Explicit template instantiation
// ===============================
 
template class cMatrixOperations<float>;
template class cMatrixOperations<double>;
template class cMatrixOperations<long double>;