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 <omp.h>  // omp_in_parallel
#include "../_cu_definitions/cu_types.h" // __nv_fp8_e5m2, __nv_fp8_e4m3,
                                         // __half, __nv_bfloat16
#include "../_cu_arithmetics/cu_arithmetics.h"  // cu_arithmetics
#include "./cublas_api.h"  // cublas_api
#include "./cusparse_api.h"  // cusparse_api
#include "../_definitions/definitions.h"  // LARGE_ARRAY_SIZE
#include <stdexcept>  // std::invalid_argument
 
 
// ============
// dense matvec
// ============
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::dense_matvec(
        cublasHandle_t cublas_handle,
        const DataType* RESTRICT A,
        const DataType* RESTRICT b,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        const FlagType A_is_row_major,
        DataType* RESTRICT c)
{
    cublasOperation_t trans;
    int m;
    int n;
    int lda;
    DataType alpha = cu_arithmetics::cast<float, DataType>(1.0f);
    DataType beta = cu_arithmetics::cast<float, DataType>(0.0f);
    int incb = 1;
    int incc = 1;
 
    // Since cublas accepts column major (Fortran) ordering, use transpose for
    // row_major matrix.
    if (A_is_row_major)
    {
        // A is row-major, not compatible with cublas. Use transpose instead.
        trans = CUBLAS_OP_T;
        m = num_columns;
        n = num_rows;
    }
    else
    {
        // A is column-major, compatible with cublas.
        trans = CUBLAS_OP_N;
        m = num_rows;
        n = num_columns;
    }
 
    lda = m;
 
    // Calling cublas
    cublasStatus_t status = cublas_api::cublasXgemv<DataType>(
            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* RESTRICT A,
        const DataType* RESTRICT b,
        const DataType alpha,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        const FlagType A_is_row_major,
        DataType* RESTRICT c)
{
    DataType zero = cu_arithmetics::cast<float, DataType>(0.0f);
    if (cu_arithmetics::is_equal(alpha, zero))
    {
        return;
    }
    
    cublasOperation_t trans;
    int m;
    int n;
    int lda;
    DataType beta = cu_arithmetics::cast<float, DataType>(1.0f);
    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_api::cublasXgemv<DataType>(
            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* RESTRICT A,
        const DataType* RESTRICT b,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        const FlagType A_is_row_major,
        DataType* RESTRICT c)
{
    cublasOperation_t trans;
    int m;
    int n;
    int lda;
    DataType alpha = cu_arithmetics::cast<float, DataType>(1.0f);
    DataType beta = cu_arithmetics::cast<float, DataType>(0.0f);
    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_api::cublasXgemv<DataType>(
            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* RESTRICT A,
        const DataType* RESTRICT b,
        const DataType alpha,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        const FlagType A_is_row_major,
        DataType* RESTRICT c)
{
    DataType zero = cu_arithmetics::cast<float, DataType>(0.0f);
    if (cu_arithmetics::is_equal(alpha, zero))
    {
        return;
    }
 
    cublasOperation_t trans;
    int m;
    int n;
    int lda;
    DataType beta = cu_arithmetics::cast<float, DataType>(0.0f);
    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_api::cublasXgemv<DataType>(
            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* RESTRICT A_data,
        const LongIndexType* RESTRICT A_column_indices,
        const LongIndexType* RESTRICT A_index_pointer,
        const DataType* RESTRICT b,
        const LongIndexType num_rows,
        DataType* RESTRICT c)
{
    throw std::runtime_error("Function not implemented.");
}
 
 
// ===============
// csr matvec plus
// ===============
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::csr_matvec_plus(
        cusparseHandle_t cusparse_handle,
        const DataType* RESTRICT A_data,
        const LongIndexType* RESTRICT A_column_indices,
        const LongIndexType* RESTRICT A_index_pointer,
        const DataType* RESTRICT b,
        const DataType alpha,
        const LongIndexType num_rows,
        DataType* RESTRICT c)
{
    throw std::runtime_error("Function not implemented.");
}
 
 
// =====================
// csr transposed matvec
// =====================
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::csr_transposed_matvec(
        cusparseHandle_t cusparse_handle,
        const DataType* RESTRICT A_data,
        const LongIndexType* RESTRICT A_column_indices,
        const LongIndexType* RESTRICT A_index_pointer,
        const DataType* RESTRICT b,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        DataType* RESTRICT c)
{
    throw std::runtime_error("Function not implemented.");
}
 
 
// ==========================
// csr transposed matvec plus
// ==========================
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::csr_transposed_matvec_plus(
        cusparseHandle_t cusparse_handle,
        const DataType* RESTRICT A_data,
        const LongIndexType* RESTRICT A_column_indices,
        const LongIndexType* RESTRICT A_index_pointer,
        const DataType* RESTRICT b,
        const DataType alpha,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        DataType* RESTRICT c)
{
    throw std::runtime_error("Function not implemented.");
}
 
 
// ==========
// csc matvec
// ==========
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::csc_matvec(
        cusparseHandle_t cusparse_handle,
        const DataType* RESTRICT A_data,
        const LongIndexType* RESTRICT A_row_indices,
        const LongIndexType* RESTRICT A_index_pointer,
        const DataType* RESTRICT b,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        DataType* RESTRICT c)
{
    throw std::runtime_error("Function not implemented.");
}
 
 
// ===============
// csc matvec plus
// ===============
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::csc_matvec_plus(
        cusparseHandle_t cusparse_handle,
        const DataType* RESTRICT A_data,
        const LongIndexType* RESTRICT A_row_indices,
        const LongIndexType* RESTRICT A_index_pointer,
        const DataType* RESTRICT b,
        const DataType alpha,
        const LongIndexType num_rows,
        const LongIndexType num_columns,
        DataType* RESTRICT c)
{
    throw std::runtime_error("Function not implemented.");
}
 
 
// =====================
// csc transposed matvec
// =====================
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::csc_transposed_matvec(
        cusparseHandle_t cusparse_handle,
        const DataType* RESTRICT A_data,
        const LongIndexType* RESTRICT A_row_indices,
        const LongIndexType* RESTRICT A_index_pointer,
        const DataType* RESTRICT b,
        const LongIndexType num_columns,
        DataType* RESTRICT c)
{
    throw std::runtime_error("Function not implemented.");
}
 
 
// ==========================
// csc transposed matvec plus
// ==========================
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::csc_transposed_matvec_plus(
        cusparseHandle_t cusparse_handle,
        const DataType* RESTRICT A_data,
        const LongIndexType* RESTRICT A_row_indices,
        const LongIndexType* RESTRICT A_index_pointer,
        const DataType* RESTRICT b,
        const DataType alpha,
        const LongIndexType num_columns,
        DataType* RESTRICT c)
{
    throw std::runtime_error("Function not implemented.");
}
 
 
// ==================
// create band matrix
// ==================
 
 
template <typename DataType>
void cuMatrixOperations<DataType>::create_band_matrix(
        cusparseHandle_t cublas_handle,
        const DataType* RESTRICT diagonals,
        const DataType* RESTRICT supdiagonals,
        const IndexType non_zero_size,
        const FlagType tridiagonal,
        DataType** RESTRICT matrix)
{
    throw std::runtime_error("Function not implemented.");
}
 
 
// ===============================
// Explicit template instantiation
// ===============================
 
#if defined(USE_CUDA_FP8_E5M2) && (USE_CUDA_FP8_E5M2 == 1)
    template class cuMatrixOperations<__nv_fp8_e5m2>;
#endif
 
#if defined(USE_CUDA_FP8_E4M3) && (USE_CUDA_FP8_E4M3 == 1)
    template class cuMatrixOperations<__nv_fp8_e4m3>;
#endif
 
#if defined(USE_CUDA_FP16) && (USE_CUDA_FP16 == 1)
    template class cuMatrixOperations<__half>;
#endif
 
#if defined(USE_CUDA_BF16) && (USE_CUDA_BF16 == 1)
    template class cuMatrixOperations<__nv_bfloat16>;
#endif
 
#if defined(USE_CUDA_FP32) && (USE_CUDA_FP32 == 1)
    template class cuMatrixOperations<float>;
#endif
 
#if defined(USE_CUDA_FP64) && (USE_CUDA_FP64 == 1)
    template class cuMatrixOperations<double>;
#endif