cu_matrix_operations.cu

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 "./cu_matrix_operations.h"
#include <cassert>  // assert
#include "./cublas_interface.h"  // cublas_interface
#include "./cusparse_interface.h"  // cusparse_interface
 
 
// ============
// dense matvec
// ============
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::dense_matvec(
        cublasHandle_t cublas_handle,
        const DataType* A,
        const DataType* b,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        const FlagType A_is_row_major,
        DataType* c)
{
    cublasOperation_t trans;
    int m;
    int n;
    int lda;
    DataType alpha = 1.0;
    DataType beta = 0.0;
    int incb = 1;
    int incc = 1;
 
    // Since cublas accepts column major (Fortran) ordering, use transpose for
    // row_major matrix.
    if (A_is_row_major)
    {
        trans = CUBLAS_OP_T;
        m = num_columns;
        n = num_rows;
    }
    else
    {
        trans = CUBLAS_OP_N;
        m = num_rows;
        n = num_columns;
    }
 
    lda = m;
 
    // Calling cublas
    cublasStatus_t status = cublas_interface::cublasXgemv(cublas_handle, trans,
                                                          m, n, &alpha, A, lda,
                                                          b, incb, &beta, c,
                                                          incc);
    assert(status == CUBLAS_STATUS_SUCCESS);
}
 
 
// =================
// dense matvec plus
// =================
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::dense_matvec_plus(
        cublasHandle_t cublas_handle,
        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)
{
    cublasOperation_t trans;
    int m;
    int n;
    int lda;
    DataType beta = 1.0;
    int incb = 1;
    int incc = 1;
 
    // Since cublas accepts column major (Fortran) ordering, use transpose for
    // row_major matrix.
    if (A_is_row_major)
    {
        trans = CUBLAS_OP_T;
        m = num_columns;
        n = num_rows;
    }
    else
    {
        trans = CUBLAS_OP_N;
        m = num_rows;
        n = num_columns;
    }
 
    lda = m;
 
    // Calling cublas
    cublasStatus_t status = cublas_interface::cublasXgemv(cublas_handle, trans,
                                                          m, n, &alpha, A, lda,
                                                          b, incb, &beta, c,
                                                          incc);
    assert(status == CUBLAS_STATUS_SUCCESS);
}
 
 
// =======================
// dense transposed matvec
// =======================
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::dense_transposed_matvec(
        cublasHandle_t cublas_handle,
        const DataType* A,
        const DataType* b,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        const FlagType A_is_row_major,
        DataType* c)
{
    cublasOperation_t trans;
    int m;
    int n;
    int lda;
    DataType alpha = 1.0;
    DataType beta = 0.0;
    int incb = 1;
    int incc = 1;
 
    // Since cublas accepts column major (Fortran) ordering, use non-transpose
    // for row_major matrix.
    if (A_is_row_major)
    {
        trans = CUBLAS_OP_N;
        m = num_columns;
        n = num_rows;
    }
    else
    {
        trans = CUBLAS_OP_T;
        m = num_rows;
        n = num_columns;
    }
 
    lda = m;
 
    // Calling cublas
    cublasStatus_t status = cublas_interface::cublasXgemv(cublas_handle, trans,
                                                          m, n, &alpha, A, lda,
                                                          b, incb, &beta, c,
                                                          incc);
    assert(status == CUBLAS_STATUS_SUCCESS);
}
 
 
// ============================
// dense transposed matvec plus
// ============================
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::dense_transposed_matvec_plus(
        cublasHandle_t cublas_handle,
        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;
    }
 
    cublasOperation_t trans;
    int m;
    int n;
    int lda;
    DataType beta = 0.0;
    int incb = 1;
    int incc = 1;
 
    // Since cublas accepts column major (Fortran) ordering, use non-transpose
    // for row_major matrix.
    if (A_is_row_major)
    {
        trans = CUBLAS_OP_N;
        m = num_columns;
        n = num_rows;
    }
    else
    {
        trans = CUBLAS_OP_T;
        m = num_rows;
        n = num_columns;
    }
 
    lda = m;
 
    // Calling cublas
    cublasStatus_t status = cublas_interface::cublasXgemv(cublas_handle, trans,
                                                          m, n, &alpha, A, lda,
                                                          b, incb, &beta, c,
                                                          incc);
    assert(status == CUBLAS_STATUS_SUCCESS);
}
 
 
// ==========
// csr matvec
// ==========
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::csr_matvec(
        cusparseHandle_t cusparse_handle,
        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;
    DataType 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] = sum;
    }
}
 
 
// ===============
// csr matvec plus
// ===============
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::csr_matvec_plus(
        cusparseHandle_t cusparse_handle,
        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;
    DataType 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 * sum;
    }
}
 
 
// =====================
// csr transposed matvec
// =====================
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::csr_transposed_matvec(
        cusparseHandle_t cusparse_handle,
        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 cuMatrixOperations<DataType>::csr_transposed_matvec_plus(
        cusparseHandle_t cusparse_handle,
        const DataType* A_data,
        const LongIndexType* A_column_indices,
        const LongIndexType* A_index_pointer,
        const DataType* b,
        const DataType alpha,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        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 cuMatrixOperations<DataType>::csc_matvec(
        cusparseHandle_t cusparse_handle,
        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 cuMatrixOperations<DataType>::csc_matvec_plus(
        cusparseHandle_t cusparse_handle,
        const DataType* A_data,
        const LongIndexType* A_row_indices,
        const LongIndexType* A_index_pointer,
        const DataType* b,
        const DataType alpha,
        const LongIndexType num_rows,
        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 cuMatrixOperations<DataType>::csc_transposed_matvec(
        cusparseHandle_t cusparse_handle,
        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;
    DataType 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] = sum;
    }
}
 
 
// ==========================
// csc transposed matvec plus
// ==========================
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::csc_transposed_matvec_plus(
        cusparseHandle_t cusparse_handle,
        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;
    DataType 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] += alpha * sum;
    }
}
 
 
// ==================
// create band matrix
// ==================
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::create_band_matrix(
        cusparseHandle_t cusparse_handle,
        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 cuMatrixOperations<float>;
template class cuMatrixOperations<double>;