Skip to content

Added methods to diagonally scale csr matrices on the left and right. #283

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
shakedregev merged 14 commits into from
May 27, 2025
Merged

Added methods to diagonally scale csr matrices on the left and right. #283

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
810677d
added diagonal scaling, compiles, need to write tests
shakedregev May 21, 2025
a80e3c9
created vector, need to add new verify function
shakedregev May 21, 2025
88854f9
cpu working, cuda not
shakedregev May 21, 2025
7be8ef5
left diagonal scale test passing on cpu and cuda
shakedregev May 22, 2025
60ecd8e
tests pass on hip
shakedregev May 23, 2025
82a0d85
code works, tests pass
shakedregev May 23, 2025
75fedea
corrected name format
shakedregev May 23, 2025
13ae6f5
Revert "corrected name format"
shakedregev May 23, 2025
cfcd70d
fixed naming conflicts
shakedregev May 23, 2025
ce25196
verified on cuda and changed names
shakedregev May 23, 2025
7acca92
made variable names follow conventions
shakedregev May 23, 2025
9682336
Update resolve/matrix/MatrixHandler.cpp
shakedregev May 27, 2025
dd2556a
fixed whitespace inconsistencies
shakedregev May 27, 2025
14af253
fixed whitespace inconsistencies
shakedregev May 27, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
187 changes: 150 additions & 37 deletions resolve/cuda/cudaKernels.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
/**
* @file cudaKernels.cu
* @author Kasia Swirydowicz (kasia.swirydowicz@pnnl.gov)
* @brief
* @brief
* @date 2023-12-08
*
*
*
*
*/

#include "cudaKernels.h"
Expand All @@ -17,22 +17,22 @@ namespace ReSolve {
/**
* @brief Computes v^T * [u1 u2] where v is n x k multivector
* and u1 and u2 are n x 1 vectors.
*
*
* @tparam Tv5 - Size of shared memory
*
*
* @param[in] u1 - (n x 1) vector
* @param[in] u2 - (n x 1) vector
* @param[in] v - (n x k) multivector
* @param[out] result - (k x 2) multivector
* @param[in] k - dimension of the subspace
* @param[in] N - size of vectors u1, u2
*/
template <size_t Tv5 = 1024>
__global__ void MassIPTwoVec(const real_type* __restrict__ u1,
const real_type* __restrict__ u2,
const real_type* __restrict__ v,
template <size_t Tv5 = 1024>
__global__ void MassIPTwoVec(const real_type* __restrict__ u1,
const real_type* __restrict__ u2,
const real_type* __restrict__ v,
real_type* result,
const index_type k,
const index_type k,
const index_type N)
{
index_type t = threadIdx.x;
Expand Down Expand Up @@ -119,16 +119,16 @@ namespace ReSolve {

/**
* @brief AXPY y = y - x*alpha where alpha is [k x 1], and x is [N x k] needed in 1 and 2 synch GMRES
*
* @tparam Tmaxk
*
*
* @tparam Tmaxk
*
* @param[in] N - number of rows in x
* @param[in] k - number of columns in x
* @param[in] x_data - double array, size [N x k]
* @param[out] y_data - double array, size [N x 1]
* @param[in] alpha - doble array, size [k x 1]
*/
template <size_t Tmaxk = 1024>
template <size_t Tmaxk = 1024>
__global__ void massAxpy3(index_type N,
index_type k,
const real_type* x_data,
Expand All @@ -155,17 +155,17 @@ namespace ReSolve {

/**
* @brief Pass through matrix rows and sum each as \sum_{j=1}^m abs(a_{ij})
*
*
* @param[in] n - number of rows in the matrix.
* @param[in] nnz - number of non-zero values in the matrix
* @param[in] a_ia - row pointers (CSR storage)
* @param[in] a_val - values (CSR storage)
* @param[out] result - array size [n x 1] containing sums of values in each row.
*/
__global__ void matrixInfNormPart1(const index_type n,
const index_type nnz,
__global__ void matrixInfNormPart1(const index_type n,
const index_type nnz,
const index_type* a_ia,
const real_type* a_val,
const real_type* a_val,
real_type* result)
{
index_type idx = blockIdx.x * blockDim.x + threadIdx.x;
Expand All @@ -178,45 +178,109 @@ namespace ReSolve {
idx += (blockDim.x * gridDim.x);
}
}
/**
* @brief Scales a csr matrix on the left by a diagonal matrix
*
* @param[in] n - number of rows in the matrix
* @param[in] a_row_ptr - row pointers (CSR storage)
* @param[in, out] a_val - values (CSR storage). Changes in place.
* @param[in] d_val - diagonal values
*
* @todo Decide how to allow user to configure grid and block sizes.
*/
__global__ void leftDiagScale(index_type n,
const index_type* a_row_ptr,
real_type* a_val,
const real_type* d_val)
{
// Get row index from thread and block indices
Comment on lines +190 to +196
Copy link
Copy Markdown
Collaborator

@pelesh pelesh May 25, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Either fix these or wait for #261 to be merged.

index_type row = blockIdx.x * blockDim.x + threadIdx.x;

} // namespace kernels
// Check if the thread's row is within matrix bounds
if (row < n) {
// Get the start and end positions for this row in the CSR format
index_type row_start = a_row_ptr[row];
index_type row_end = a_row_ptr[row + 1];

// Get the scaling factor for this row from the diagonal matrix
real_type scale = d_val[row];

// Scale all non-zero elements in this row
for (index_type i = row_start; i < row_end; i++) {
a_val[i] *= scale;
}
}
}

/**
* @brief Scales a csr matrix on the right by a diagonal matrix
*
* @param[in] n - number of rows in the matrix
* @param[in] a_row_ptr - row pointers (CSR storage)
* @param[in] a_col_ind - column indices (CSR storage)
* @param[in, out] a_val - values (CSR storage). Changes in place.
* @param[in] d_val - diagonal values
*
* @todo Decide how to allow user to configure grid and block sizes.
*/
__global__ void rightDiagScale(index_type n,
const index_type* a_row_ptr,
const index_type* a_col_ind,
real_type* a_val,
const real_type* d_val)
{
// Get row index from thread and block indices
index_type row = blockIdx.x * blockDim.x + threadIdx.x;

// Check if the thread's row is within matrix bounds
if (row < n) {
// Get the start and end positions for this row in the CSR format
index_type row_start = a_row_ptr[row];
index_type row_end = a_row_ptr[row + 1];

// Scale all non-zero elements in this row
for (index_type i = row_start; i < row_end; i++) {
a_val[i] *= d_val[a_col_ind[i]];
}
}
}
} // namespace kernels

//
// Kernel wrappers
//

/**
* @brief Computes result = mvec^T * [vec1 vec2]
*
*
* @param n - size of vectors vec1, vec2
* @param i -
* @param vec1 - (n x 1) vector
* @param i -
* @param vec1 - (n x 1) vector
* @param vec2 - (n x 1) vector
* @param mvec - (n x (i+1)) multivector
* @param result - ((i+1) x 2) multivector
*
*
* @todo Input data should be `const`.
* @todo Is it coincidence that the block size is equal to the default
* value of Tv5?
* @todo Should we use dynamic shared memory here instead?
*/
void mass_inner_product_two_vectors(index_type n,
index_type i,
const real_type* vec1,
const real_type* vec2,
const real_type* mvec,
void mass_inner_product_two_vectors(index_type n,
index_type i,
const real_type* vec1,
const real_type* vec2,
const real_type* mvec,
real_type* result)
{
kernels::MassIPTwoVec<<<i, 1024>>>(vec1, vec2, mvec, result, i, n);
}

/**
* @brief Computes y := y - x*alpha
*
*
* @param[in] n - vector size
* @param[in] i - number of vectors in the multivector
* @param[in] x - (n x (i+1)) multivector
* @param[in] x - (n x (i+1)) multivector
* @param[out] y - (n x (i+1)) multivector
* @param[in] alpha - ((i+1) x 1) vector
*/
Expand All @@ -226,20 +290,69 @@ namespace ReSolve {
}

/**
* @brief
*
* @brief Wrapper that scales a csr matrix on the left by a diagonal matrix
*
* @param[in] n - number of rows in the matrix
* @param[in] a_row_ptr - row pointers (CSR storage)
* @param[in, out] a_val - values (CSR storage). Changes in place.
* @param[in] d_val - diagonal values
*
* @todo Decide how to allow user to configure grid and block sizes.
*/
void leftDiagScale(index_type n,
const index_type* a_row_ptr,
real_type* a_val,
const real_type* d_val)
{

// Define block size and number of blocks
const int blockSize = 1;
int numBlocks = (n + blockSize - 1) / blockSize;
// Launch the kernel
kernels::leftDiagScale<<<numBlocks, blockSize>>>(n, a_row_ptr, a_val, d_val);
}

/**
* @brief Wrapper that scales a csr matrix on the right by a diagonal matrix
*
* @param[in] n - number of rows in the matrix
* @param[in] a_row_ptr - row pointers (CSR storage)
* @param[in] a_col_ind - column indices (CSR storage)
* @param[in, out] a_val - values (CSR storage). Changes in place.
* @param[in] d_val - diagonal values
*
* @todo Decide how to allow user to configure grid and block sizes.
*/
void rightDiagScale(index_type n,
const index_type* a_row_ptr,
const index_type* a_col_ind,
real_type* a_val,
const real_type* d_val)
{
// Define block size and number of blocks
const int blockSize = 256;
int numBlocks = (n + blockSize - 1) / blockSize;
// Launch the kernel
kernels::rightDiagScale<<<numBlocks, blockSize>>>(n, a_row_ptr, a_col_ind, a_val, d_val);
}



/**
* @brief
*
* @param[in] n -
* @param[in] nnz -
* @param[in] a_ia -
* @param[in] a_ia -
* @param[in] a_val -
* @param[out] result -
*
*
* @todo Decide how to allow user to configure grid and block sizes.
*/
void matrix_row_sums(index_type n,
index_type nnz,
void matrix_row_sums(index_type n,
index_type nnz,
const index_type* a_ia,
const real_type* a_val,
const real_type* a_val,
real_type* result)
{
kernels::matrixInfNormPart1<<<1000, 1024>>>(n, nnz, a_ia, a_val, result);
Expand Down
33 changes: 22 additions & 11 deletions resolve/cuda/cudaKernels.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,32 +1,43 @@
/**
* @file cudaKernels.h
* @author Kasia Swirydowicz (kasia.swirydowicz@pnnl.gov)
*
*
* @brief Contains prototypes of CUDA kernels.
*
*
* @note These kernels will be used in CUDA specific code, only.
*
*
*/
#pragma once

#include <resolve/Common.hpp>

namespace ReSolve {

void mass_inner_product_two_vectors(index_type n,
index_type i,
const real_type* vec1,
const real_type* vec2,
const real_type* mvec,
void mass_inner_product_two_vectors(index_type n,
index_type i,
const real_type* vec1,
const real_type* vec2,
const real_type* mvec,
real_type* result);

void mass_axpy(index_type n, index_type i, const real_type* x, real_type* y, const real_type* alpha);

void leftDiagScale(index_type n,
const index_type* a_row_ptr,
real_type* a_val,
const real_type* diag);

void rightDiagScale(index_type n,
const index_type* a_row_ptr,
const index_type* a_col_idx,
real_type* a_val,
const real_type* diag);

//needed for matrix inf nrm
void matrix_row_sums(index_type n,
index_type nnz,
void matrix_row_sums(index_type n,
index_type nnz,
const index_type* a_ia,
const real_type* a_val,
const real_type* a_val,
real_type* result);

} // namespace ReSolve
Loading