Skip to content
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

CRS matrix Tpetra wrapper type #111

Merged
merged 8 commits into from
Mar 31, 2025
Merged

CRS matrix Tpetra wrapper type #111

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
2fa4e83
ENH: tpetra crs matrix constructor with local dim
Adrian-Diaz Mar 8, 2025
f0bf2c7
ENH: CRS Matrix multiply operator
Adrian-Diaz Mar 10, 2025
47f53f0
WIP: Functor approach to arithmetic operators on MATAR types
Adrian-Diaz Mar 11, 2025
622c07e
ENH: functor style matrix multiply operator for flexibility
Adrian-Diaz Mar 12, 2025
f721968
Merge branch 'Adrian's-Branch' of https://github.com/Adrian-Diaz/MATA…
Adrian-Diaz Mar 12, 2025
b47f87b
ENH: Operators + CRS sorting bug
Adrian-Diaz Mar 25, 2025
0ee3aba
BUG: correct for sorting
Adrian-Diaz Mar 25, 2025
1c5b7f1
ENH: addition and subtract operators for vectors
Adrian-Diaz Mar 31, 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
124 changes: 91 additions & 33 deletions examples/test_tpetra_crs.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,15 +41,18 @@ using namespace mtr; // matar namespace

void TpetraCRSMatrixExample();

void TpetraCRSMatrixMultiplyExample();

int main(int argc, char* argv[])
{
MPI_Init(&argc, &argv);
int process_rank;
MPI_Comm_rank(MPI_COMM_WORLD, &process_rank);
Kokkos::initialize();
{
// Run TpetraFArray 1D example
// Run TpetraCRS example
TpetraCRSMatrixExample();
TpetraCRSMatrixMultiplyExample();
} // end of kokkos scope
Kokkos::finalize();
MPI_Barrier(MPI_COMM_WORLD);
Expand All @@ -64,7 +67,7 @@ void TpetraCRSMatrixExample()
MPI_Comm_rank(MPI_COMM_WORLD, &process_rank);

if(process_rank==0)
printf("\n====================Running TpetraCRSMatrix example====================\n");
printf("\n====================Running TpetraCRSMatrix example with map argument====================\n");

//construct a row map over MPI ranks
long long int n = 100; //global dimension
Expand All @@ -91,41 +94,96 @@ void TpetraCRSMatrixExample()
RaggedRightArrayKokkos<double, Kokkos::LayoutRight> input_values(matrix_strides,"ragged_values");
FOR_ALL(i, 0, nlocal,{
for(int j = 0; j < matrix_strides(i); j++){
input_values(i,j) = 3*(min_global_index+j);
input_values(i,j) = 3*j;
}
});

//NOTE! the constructor resorts the crs graph since Trilinos requires this; the values in each row are likely
//in different places (in the same row) afterwards but the correspondence to global column ID is always preserved
TpetraCRSMatrix<double, Kokkos::LayoutRight> mymatrix(input_pmap, matrix_strides, input_crs, input_values);
//TpetraCRSMatrix<double, Kokkos::LayoutRight> mymatrix(input_pmap, matrix_strides);
mymatrix.print();

// //local size
// int nlocal = myarray.size();

// // set values on host copy of data
// if(process_rank==0)
// printf("Printing host copy of data (should be global ids):\n");
// for (int i = 0; i < nlocal; i++) {
// //set each array element to the corresponding global index
// //we get global indices using a partition map member in the array
// myarray.host(i) = myarray.pmap.getGlobalIndex(i);
// }

// myarray.update_device();

// // Print host copy of data
// myarray.print();
// Kokkos::fence();

// // Manupulate data on device and update host
// FOR_ALL(i, 0, nlocal,{
// myarray(i) = 2*myarray(i);
// });
// myarray.update_host();
// Kokkos::fence();
// if(process_rank==0)
// printf("---Data multiplied by 2 on device---\n");

// // Print host copy of data
// myarray.print();
// Kokkos::fence();
//test case that doesnt pass a map object; reset inputs since the graph and values array were resorted for Trilinos compatibility
if(process_rank==0)
printf("\n====================Running TpetraCRSMatrix example with local dim argument====================\n");

//global indices array
FOR_ALL(i, 0, nlocal,{
for(int j = 0; j < matrix_strides(i); j++){
input_crs(i,j) = j;
}
});

//values array
FOR_ALL(i, 0, nlocal,{
for(int j = 0; j < matrix_strides(i); j++){
input_values(i,j) = 3*j;
}
});
TpetraCRSMatrix<double, Kokkos::LayoutRight> mymatrix2(nlocal, matrix_strides, input_crs, input_values);
//TpetraCRSMatrix<double, Kokkos::LayoutRight> mymatrix(input_pmap, matrix_strides);
mymatrix2.print();
}

void TpetraCRSMatrixMultiplyExample()
{
int process_rank;
MPI_Comm_rank(MPI_COMM_WORLD, &process_rank);
if(process_rank==0)
printf("\n====================Running Matrix Multiply example====================\n");

int nx = 100; //global dimension
int ny = 5;

//distributed (first dimension gets distributed) dual array with layout left
TpetraDFArray<double> myarray(nx, ny);

//local size
int nxlocal = myarray.dims(0);

// set values on host copy of data
for (int i = 0; i < nxlocal; i++) {
for (int j = 0; j < ny; j++){
//set each array element to a computed global degree of freedom index
//we get global indices for dim0 using a partition map member in the array
myarray.host(i,j) = ny*myarray.pmap.getGlobalIndex(i) + j;
}
}
myarray.update_device();

//construct matrix
//construct strides, index graph, and values arrays
DCArrayKokkos<size_t, Kokkos::LayoutRight> matrix_strides(nxlocal, "matrix_strides");
//set strides; map is contiguous so Trilinos leaves device view of map empty (BE WARNED)
const long long int min_global_index = myarray.pmap.getMinGlobalIndex();
FOR_ALL(i, 0, nxlocal,{
matrix_strides(i) = (min_global_index+i) + 1;
});

//global indices array
RaggedRightArrayKokkos<long long int, Kokkos::LayoutRight> input_crs(matrix_strides,"graph_indices");
FOR_ALL(i, 0, nxlocal,{
for(int j = 0; j < matrix_strides(i); j++){
input_crs(i,j) = j;
}
});

//values array
RaggedRightArrayKokkos<double, Kokkos::LayoutRight> input_values(matrix_strides,"ragged_values");
FOR_ALL(i, 0, nxlocal,{
for(int j = 0; j < matrix_strides(i); j++){
input_values(i,j) = 3*j;
}
});

TpetraCRSMatrix<double, Kokkos::LayoutRight> mymatrix(nxlocal, matrix_strides, input_crs, input_values);

if(process_rank==0)
printf("multiplication result:\n");

//perform multiplication; currently the operator can allocate the result vector for you if it wasnt already
TpetraDFArray<double> result = mymatrix*myarray;
result.update_host();
result.print();
}
54 changes: 54 additions & 0 deletions examples/test_tpetra_farray.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ using namespace mtr; // matar namespace
void TpetraFArrayOneDimensionExample();
void TpetraFArrayTwoDimensionExample();
void TpetraFArraySevenDimensionExample();
void TpetraFArrayAddExample();

int main(int argc, char* argv[])
{
Expand All @@ -58,6 +59,9 @@ int main(int argc, char* argv[])

// Run TpetraFArray 7D example
TpetraFArraySevenDimensionExample();

// Run TpetraFArray Addition example
TpetraFArrayAddExample();
} // end of kokkos scope
Kokkos::finalize();
MPI_Barrier(MPI_COMM_WORLD);
Expand Down Expand Up @@ -232,3 +236,53 @@ void TpetraFArraySevenDimensionExample()
myarray.print();
Kokkos::fence();
}

void TpetraFArrayAddExample()
{
int process_rank;
MPI_Comm_rank(MPI_COMM_WORLD, &process_rank);
if(process_rank==0)
printf("\n====================Running 2D TpetraFarray example====================\n");

int nx = 20; //global dimension
int ny = 5;

//distributed (first dimension gets distributed) dual array with layout left
TpetraDFArray<double> myarray(nx, ny);

//local size
int nxlocal = myarray.dims(0);

// set values on host copy of data
if(process_rank==0)
printf("Printing host copy of data (should be global ids):\n");
for (int i = 0; i < nxlocal; i++) {
for (int j = 0; j < ny; j++){
//set each array element to a computed global degree of freedom index
//we get global indices for dim0 using a partition map member in the array
myarray.host(i,j) = ny*myarray.pmap.getGlobalIndex(i) + j;
}
}

myarray.update_device();

// Print host copy of data
//myarray.print();
Kokkos::fence();

TpetraDFArray<double> myarray2(nx, ny);
// Manupulate data on device and update host
FOR_ALL(i, 0, nxlocal,
j, 0, ny,{
myarray2(i,j) = 2*myarray(i,j);
});
myarray2.update_host();
Kokkos::fence();
if(process_rank==0)
printf("---Data after additition on device---\n");
TpetraDFArray<double> myresult = myarray + myarray2;
// Print host copy of data
myresult.update_host();
myresult.print();
Kokkos::fence();
}
17 changes: 17 additions & 0 deletions src/include/host_types.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5388,6 +5388,23 @@ CSCArray<T>::~CSCArray() {}

// End of CSCArray

////////////////////////////////////////////////
// OperatorFunctor: Base class for math functors.
////////////////////////////////////////////////

class OperatorFunctor {

public:

OperatorFunctor(){}

// Method that update device view
virtual void apply_function(void* Y) const {}

// Deconstructor
~OperatorFunctor (){}
}; // End of TpetraCRSMatrix


//=======================================================================
// end of standard MATAR data-types
Expand Down
Loading