diff --git a/README.md b/README.md index a837deea..1293fcaa 100644 --- a/README.md +++ b/README.md @@ -91,6 +91,7 @@ Primary authors of this project are Kasia Świrydowicz and Slaven Peles. ReSolve project would not be possible without significant contributions from (in alphabetic order): - Maksudul Alam +- Kaleb Brunhoeber - Ryan Danehy - Nicholson Koukpaizan - Jaelyn Litzinger diff --git a/examples/r_KLU_GLU.cpp b/examples/r_KLU_GLU.cpp index a050e324..3e1ffffe 100644 --- a/examples/r_KLU_GLU.cpp +++ b/examples/r_KLU_GLU.cpp @@ -35,7 +35,6 @@ int main(int argc, char *argv[]) std::string matrixFileNameFull; std::string rhsFileNameFull; - ReSolve::matrix::Coo* A_coo; ReSolve::matrix::Csr* A; ReSolve::LinAlgWorkspaceCUDA* workspace_CUDA = new ReSolve::LinAlgWorkspaceCUDA; workspace_CUDA->initializeHandles(); @@ -82,15 +81,11 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileNameFull << "\n"; return -1; } + bool is_expand_symmetric = true; if (i == 0) { - A_coo = ReSolve::io::readMatrixFromFile(mat_file); - A = new ReSolve::matrix::Csr(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); @@ -98,23 +93,22 @@ int main(int argc, char *argv[]) vec_x->allocate(ReSolve::memory::DEVICE); vec_r = new vector_type(A->getNumRows()); } else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout<<"Finished reading the matrix and rhs, size: "<getNumRows()<<" x "<getNumColumns()<< ", nnz: "<< A->getNnz()<< ", symmetric? "<symmetric()<< ", Expanded? "<expanded()<updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); } - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<initializeHandles(); ReSolve::MatrixHandler* matrix_handler = new ReSolve::MatrixHandler(workspace_CUDA); @@ -84,15 +83,11 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileNameFull << "\n"; return -1; } + bool is_expand_symmetric = true; if (i == 0) { - A_coo = ReSolve::io::readMatrixFromFile(mat_file); - A = new ReSolve::matrix::Csr(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); @@ -101,29 +96,28 @@ int main(int argc, char *argv[]) vec_r = new vector_type(A->getNumRows()); } else { if (i==1) { - A_exp_coo = ReSolve::io::readMatrixFromFile(mat_file); + A_exp = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); } else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_exp_coo); + ReSolve::io::updateMatrixFromFile(mat_file, A_exp); } std::cout<<"Updating values of A_coo!"<updateValues(A_exp_coo->getValues(ReSolve::memory::HOST), ReSolve::memory::HOST, ReSolve::memory::HOST); - //ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + A->updateValues(A_exp->getValues(ReSolve::memory::HOST), ReSolve::memory::HOST, ReSolve::memory::HOST); + //ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout<<"Finished reading the matrix and rhs, size: "<getNumRows()<<" x "<getNumColumns()<< ", nnz: "<< A->getNnz()<< ", symmetric? "<symmetric()<< ", Expanded? "<expanded()<updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); } - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); vec_r = new vector_type(A->getNumRows()); } else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout << "Finished reading the matrix and rhs, size: " << A->getNumRows() << " x " << A->getNumColumns() @@ -105,16 +100,14 @@ int main(int argc, char *argv[]) mat_file.close(); rhs_file.close(); - //Now convert to CSR. + // Update data. if (i < 2) { - A->updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - A->updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); } - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + bool is_expand_symmetric = true; + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); @@ -81,11 +76,9 @@ int main(int argc, char *argv[]) mat_file.close(); rhs_file.close(); - //Now convert to CSR. - A->updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); - std::cout << "COO to CSR completed. Expanded NNZ: " << A->getNnzExpanded() << std::endl; + std::cout << "COO to CSR completed. Expanded NNZ: " << A->getNnz() << std::endl; //Now call direct solver int status; KLU->setup(A); diff --git a/examples/r_KLU_cusolverrf_redo_factorization.cpp b/examples/r_KLU_cusolverrf_redo_factorization.cpp index 3f53a020..78ef3e99 100644 --- a/examples/r_KLU_cusolverrf_redo_factorization.cpp +++ b/examples/r_KLU_cusolverrf_redo_factorization.cpp @@ -35,7 +35,6 @@ int main(int argc, char *argv[] ) std::string matrixFileNameFull; std::string rhsFileNameFull; - ReSolve::matrix::Coo* A_coo; ReSolve::matrix::Csr* A; ReSolve::LinAlgWorkspaceCUDA* workspace_CUDA = new ReSolve::LinAlgWorkspaceCUDA; @@ -93,37 +92,32 @@ int main(int argc, char *argv[] ) std::cout << "Failed to open file " << rhsFileNameFull << "\n"; return -1; } + bool is_expand_symmetric = true; if (i == 0) { - A_coo = ReSolve::io::readMatrixFromFile(mat_file); - A = new ReSolve::matrix::Csr(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); vec_r = new vector_type(A->getNumRows()); } else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout<<"Finished reading the matrix and rhs, size: "<getNumRows()<<" x "<getNumColumns()<< ", nnz: "<< A->getNnz()<< ", symmetric? "<symmetric()<< ", Expanded? "<expanded()<coo2csr(A_coo, A, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - matrix_handler->coo2csr(A_coo, A, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); } - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<setup(A); @@ -215,7 +209,6 @@ int main(int argc, char *argv[] ) //now DELETE delete A; - delete A_coo; delete KLU; delete Rf; delete [] x; diff --git a/examples/r_KLU_rf.cpp b/examples/r_KLU_rf.cpp index 376e4be2..71c4b752 100644 --- a/examples/r_KLU_rf.cpp +++ b/examples/r_KLU_rf.cpp @@ -35,7 +35,6 @@ int main(int argc, char *argv[] ) std::string matrixFileNameFull; std::string rhsFileNameFull; - ReSolve::matrix::Coo* A_coo; ReSolve::matrix::Csr* A; ReSolve::LinAlgWorkspaceCUDA* workspace_CUDA = new ReSolve::LinAlgWorkspaceCUDA; @@ -82,38 +81,33 @@ int main(int argc, char *argv[] ) std::cout << "Failed to open file " << rhsFileNameFull << "\n"; return -1; } + bool is_expand_symmetric = true; if (i == 0) { - A_coo = ReSolve::io::readMatrixFromFile(mat_file); - A = new ReSolve::matrix::Csr(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); vec_r = new vector_type(A->getNumRows()); } else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout<<"Finished reading the matrix and rhs, size: "<getNumRows()<<" x "<getNumColumns()<< ", nnz: "<< A->getNnz()<< ", symmetric? "<symmetric()<< ", Expanded? "<expanded()<updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); } - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<initializeHandles(); @@ -85,44 +84,38 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileNameFull << "\n"; return -1; } + bool is_expand_symmetric = true; if (i == 0) { - A_coo = ReSolve::io::readMatrixFromFile(mat_file); - A = new ReSolve::matrix::Csr(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); vec_x->allocate(ReSolve::memory::HOST);//for KLU vec_x->allocate(ReSolve::memory::DEVICE); vec_r = new vector_type(A->getNumRows()); - } - else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + } else { + ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout<<"Finished reading the matrix and rhs, size: "<getNumRows()<<" x "<getNumColumns()<< ", nnz: "<< A->getNnz()<< ", symmetric? "<symmetric()<< ", Expanded? "<expanded()<updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); } - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<setup(A); matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE); status = KLU->analyze(); diff --git a/examples/r_KLU_rf_FGMRES_reuse_factorization.cpp b/examples/r_KLU_rf_FGMRES_reuse_factorization.cpp index f3715121..93f811e6 100644 --- a/examples/r_KLU_rf_FGMRES_reuse_factorization.cpp +++ b/examples/r_KLU_rf_FGMRES_reuse_factorization.cpp @@ -37,7 +37,6 @@ int main(int argc, char *argv[]) std::string matrixFileNameFull; std::string rhsFileNameFull; - ReSolve::matrix::Coo* A_coo; ReSolve::matrix::Csr* A; ReSolve::LinAlgWorkspaceCUDA* workspace_CUDA = new ReSolve::LinAlgWorkspaceCUDA; @@ -87,14 +86,10 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileNameFull << "\n"; return -1; } + bool is_expand_symmetric = true; if (i == 0) { - A_coo = ReSolve::io::readMatrixFromFile(mat_file); - A = new ReSolve::matrix::Csr(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); @@ -103,23 +98,22 @@ int main(int argc, char *argv[]) vec_r = new vector_type(A->getNumRows()); } else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout<<"Finished reading the matrix and rhs, size: "<getNumRows()<<" x "<getNumColumns()<< ", nnz: "<< A->getNnz()<< ", symmetric? "<symmetric()<< ", Expanded? "<expanded()<updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); } - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<initializeHandles(); @@ -90,15 +89,11 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileNameFull << "\n"; return -1; } + bool is_expand_symmetric = true; if (i == 0) { - A_coo = ReSolve::io::readMatrixFromFile(mat_file); - A = new ReSolve::matrix::Csr(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); @@ -106,8 +101,8 @@ int main(int argc, char *argv[]) vec_x->allocate(ReSolve::memory::DEVICE); vec_r = new vector_type(A->getNumRows()); } else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout << "Finished reading the matrix and rhs, size: " << A->getNumRows() << " x "<< A->getNumColumns() << ", nnz: " << A->getNnz() @@ -117,18 +112,17 @@ int main(int argc, char *argv[]) rhs_file.close(); RESOLVE_RANGE_POP("Matrix Read"); - //Now convert to CSR. + // Update host and device data. RESOLVE_RANGE_PUSH("Convert to CSR"); if (i < 2) { - A->updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); } RESOLVE_RANGE_POP("Convert to CSR"); - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); vec_r = new vector_type(A->getNumRows()); } else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout<<"Finished reading the matrix and rhs, size: "<getNumRows()<<" x "<getNumColumns()<< ", nnz: "<< A->getNnz()<< ", symmetric? "<symmetric()<< ", Expanded? "<expanded()<updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); } - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); vec_r = new vector_type(A->getNumRows()); } else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout<<"Finished reading the matrix and rhs, size: "<getNumRows()<<" x "<getNumColumns()<< ", nnz: "<< A->getNnz()<< ", symmetric? "<symmetric()<< ", Expanded? "<expanded()<updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); } - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); vec_r = new vector_type(A->getNumRows()); } else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout << "Finished reading the matrix and rhs, size: " << A->getNumRows() << " x " << A->getNumColumns() @@ -105,16 +100,14 @@ int main(int argc, char *argv[]) mat_file.close(); rhs_file.close(); - //Now convert to CSR. + // Update data. if (i < 2) { - A->updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - A->updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); } - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<setMatrix(A); int status; diff --git a/examples/r_SysSolverCuda.cpp b/examples/r_SysSolverCuda.cpp index abffc586..ab7ea303 100644 --- a/examples/r_SysSolverCuda.cpp +++ b/examples/r_SysSolverCuda.cpp @@ -35,7 +35,6 @@ int main(int argc, char *argv[]) std::string matrixFileNameFull; std::string rhsFileNameFull; - ReSolve::matrix::Coo* A_coo; ReSolve::matrix::Csr* A; ReSolve::LinAlgWorkspaceCUDA* workspace_CUDA = new ReSolve::LinAlgWorkspaceCUDA; workspace_CUDA->initializeHandles(); @@ -78,23 +77,19 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileNameFull << "\n"; return -1; } + bool is_expand_symmetric = true; if (i == 0) { - A_coo = ReSolve::io::readMatrixFromFile(mat_file); - A = new ReSolve::matrix::Csr(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); vec_x->allocate(ReSolve::memory::HOST);//for KLU vec_x->allocate(ReSolve::memory::DEVICE); } else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout << "Finished reading the matrix and rhs, size: " << A->getNumRows() << " x " << A->getNumColumns() @@ -104,16 +99,15 @@ int main(int argc, char *argv[]) mat_file.close(); rhs_file.close(); - //Now convert to CSR. + // Update host and device data. if (i < 1) { - A->updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); } - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<setMatrix(A); int status = -1; diff --git a/examples/r_SysSolverHip.cpp b/examples/r_SysSolverHip.cpp index 9a1be1fc..fa720a68 100644 --- a/examples/r_SysSolverHip.cpp +++ b/examples/r_SysSolverHip.cpp @@ -36,7 +36,6 @@ int main(int argc, char *argv[] ) std::string matrixFileNameFull; std::string rhsFileNameFull; - ReSolve::matrix::Coo* A_coo; ReSolve::matrix::Csr* A; ReSolve::LinAlgWorkspaceHIP* workspace_HIP = new ReSolve::LinAlgWorkspaceHIP(); @@ -80,21 +79,17 @@ int main(int argc, char *argv[] ) std::cout << "Failed to open file " << rhsFileNameFull << "\n"; return -1; } + bool is_expand_symmetric = true; if (i == 0) { - A_coo = ReSolve::io::readMatrixFromFile(mat_file); - A = new ReSolve::matrix::Csr(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); } else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout << "Finished reading the matrix and rhs, size: " << A->getNumRows() << " x " << A->getNumColumns() @@ -104,16 +99,15 @@ int main(int argc, char *argv[] ) mat_file.close(); rhs_file.close(); - //Now convert to CSR. + // Update host and device data. if (i < 2) { - A->updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); } - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<setMatrix(A); int status = -1; @@ -144,7 +138,6 @@ int main(int argc, char *argv[] ) //now DELETE delete A; - delete A_coo; delete [] x; delete [] rhs; delete vec_rhs; diff --git a/examples/r_SysSolverHipRefine.cpp b/examples/r_SysSolverHipRefine.cpp index 288abc96..d340b51c 100644 --- a/examples/r_SysSolverHipRefine.cpp +++ b/examples/r_SysSolverHipRefine.cpp @@ -37,7 +37,6 @@ int main(int argc, char *argv[]) std::string matrixFileNameFull; std::string rhsFileNameFull; - ReSolve::matrix::Coo* A_coo; ReSolve::matrix::Csr* A; ReSolve::LinAlgWorkspaceHIP* workspace_HIP = new ReSolve::LinAlgWorkspaceHIP(); workspace_HIP->initializeHandles(); @@ -82,23 +81,19 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileNameFull << "\n"; return -1; } + bool is_expand_symmetric = true; if (i == 0) { - A_coo = ReSolve::io::readMatrixFromFile(mat_file); - A = new ReSolve::matrix::Csr(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); vec_x->allocate(ReSolve::memory::HOST);//for KLU vec_x->allocate(ReSolve::memory::DEVICE); } else { - ReSolve::io::readAndUpdateMatrix(mat_file, A_coo); - ReSolve::io::readAndUpdateRhs(rhs_file, &rhs); + ReSolve::io::updateMatrixFromFile(mat_file, A); + ReSolve::io::updateArrayFromFile(rhs_file, &rhs); } std::cout << "Finished reading the matrix and rhs, size: " << A->getNumRows() << " x " << A->getNumColumns() @@ -108,16 +103,15 @@ int main(int argc, char *argv[]) mat_file.close(); rhs_file.close(); - //Now convert to CSR. + // Update host and device data. if (i < 2) { - A->updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); } else { - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); } - std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<getNnz()<setMatrix(A); int status; if (i < 2) { @@ -169,7 +163,6 @@ int main(int argc, char *argv[]) } // for (int i = 0; i < numSystems; ++i) delete A; - delete A_coo; delete solver; delete [] x; delete [] rhs; diff --git a/examples/r_randGMRES.cpp b/examples/r_randGMRES.cpp index 0bfec15d..2a7fb3f8 100644 --- a/examples/r_randGMRES.cpp +++ b/examples/r_randGMRES.cpp @@ -26,7 +26,6 @@ int main(int argc, char *argv[]) std::string rhsFileName = argv[2]; - ReSolve::matrix::Coo* A_coo; ReSolve::matrix::Csr* A; ReSolve::LinAlgWorkspaceHIP* workspace_HIP = new ReSolve::LinAlgWorkspaceHIP(); workspace_HIP->initializeHandles(); @@ -61,14 +60,10 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileName << "\n"; return -1; } - A_coo = ReSolve::io::readMatrixFromFile(mat_file); - A = new ReSolve::matrix::Csr(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + bool is_expand_symmetric = true; + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); @@ -81,7 +76,7 @@ int main(int argc, char *argv[]) mat_file.close(); rhs_file.close(); - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); //Now call direct solver real_type norm_b; @@ -113,7 +108,6 @@ int main(int argc, char *argv[]) delete A; - delete A_coo; delete Rf; delete [] x; delete [] rhs; diff --git a/examples/r_randGMRES_CUDA.cpp b/examples/r_randGMRES_CUDA.cpp index dda08d3c..9c60915a 100644 --- a/examples/r_randGMRES_CUDA.cpp +++ b/examples/r_randGMRES_CUDA.cpp @@ -25,8 +25,6 @@ int main(int argc, char *argv[]) std::string matrixFileName = argv[1]; std::string rhsFileName = argv[2]; - - ReSolve::matrix::Coo* A_coo; ReSolve::matrix::Csr* A; ReSolve::LinAlgWorkspaceCUDA* workspace_CUDA = new ReSolve::LinAlgWorkspaceCUDA(); workspace_CUDA->initializeHandles(); @@ -61,14 +59,10 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileName << "\n"; return -1; } - A_coo = ReSolve::io::readMatrixFromFile(mat_file); - A = new ReSolve::matrix::Csr(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); - - rhs = ReSolve::io::readRhsFromFile(rhs_file); + bool is_expand_symmetric = true; + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); + + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); @@ -81,7 +75,7 @@ int main(int argc, char *argv[]) mat_file.close(); rhs_file.close(); - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); + A->copyData(ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); //Now call direct solver real_type norm_b; @@ -113,7 +107,6 @@ int main(int argc, char *argv[]) delete A; - delete A_coo; delete Rf; delete [] x; delete [] rhs; diff --git a/examples/r_randGMRES_cpu.cpp b/examples/r_randGMRES_cpu.cpp index 332f563c..73edc20b 100644 --- a/examples/r_randGMRES_cpu.cpp +++ b/examples/r_randGMRES_cpu.cpp @@ -25,8 +25,6 @@ int main(int argc, char *argv[]) std::string matrixFileName = argv[1]; std::string rhsFileName = argv[2]; - - ReSolve::matrix::Coo* A_coo; ReSolve::matrix::Csr* A; ReSolve::LinAlgWorkspaceCpu* workspace_Cpu = new ReSolve::LinAlgWorkspaceCpu(); workspace_Cpu->initializeHandles(); @@ -61,14 +59,10 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileName << "\n"; return -1; } - A_coo = ReSolve::io::readMatrixFromFile(mat_file); - A = new ReSolve::matrix::Csr(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()); + bool is_expand_symmetric = true; + A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric); - rhs = ReSolve::io::readRhsFromFile(rhs_file); + rhs = ReSolve::io::createArrayFromFile(rhs_file); x = new real_type[A->getNumRows()]; vec_rhs = new vector_type(A->getNumRows()); vec_x = new vector_type(A->getNumRows()); @@ -81,7 +75,6 @@ int main(int argc, char *argv[]) mat_file.close(); rhs_file.close(); - A->updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); //Now call direct solver real_type norm_b; @@ -112,7 +105,6 @@ int main(int argc, char *argv[]) delete A; - delete A_coo; delete Rf; delete [] x; delete [] rhs; diff --git a/resolve/LinSolverDirectCuSolverGLU.cpp b/resolve/LinSolverDirectCuSolverGLU.cpp index 6de6c61c..eb91fc76 100644 --- a/resolve/LinSolverDirectCuSolverGLU.cpp +++ b/resolve/LinSolverDirectCuSolverGLU.cpp @@ -40,7 +40,7 @@ namespace ReSolve handle_cusolversp_ = workspaceCUDA->getCusolverSpHandle(); A_ = (matrix::Csr*) A; index_type n = A_->getNumRows(); - index_type nnz = A_->getNnzExpanded(); + index_type nnz = A_->getNnz(); //create combined factor addFactors(L,U); @@ -160,7 +160,7 @@ namespace ReSolve status_cusolver_ = cusolverSpDgluReset(handle_cusolversp_, A_->getNumRows(), /* A is original matrix */ - A_->getNnzExpanded(), + A_->getNnz(), descr_A_, A_->getValues( memory::DEVICE), A_->getRowData(memory::DEVICE), diff --git a/resolve/LinSolverDirectCuSolverRf.cpp b/resolve/LinSolverDirectCuSolverRf.cpp index e08c5844..3b176d31 100644 --- a/resolve/LinSolverDirectCuSolverRf.cpp +++ b/resolve/LinSolverDirectCuSolverRf.cpp @@ -57,7 +57,7 @@ namespace ReSolve status_cusolverrf_ = cusolverRfSetResetValuesFastMode(handle_cusolverrf_, CUSOLVERRF_RESET_VALUES_FAST_MODE_ON); error_sum += status_cusolverrf_; status_cusolverrf_ = cusolverRfSetupDevice(n, - A_->getNnzExpanded(), + A_->getNnz(), A_->getRowData(memory::DEVICE), A_->getColData(memory::DEVICE), A_->getValues( memory::DEVICE), @@ -101,7 +101,7 @@ namespace ReSolve { int error_sum = 0; status_cusolverrf_ = cusolverRfResetValues(A_->getNumRows(), - A_->getNnzExpanded(), + A_->getNnz(), A_->getRowData(memory::DEVICE), A_->getColData(memory::DEVICE), A_->getValues( memory::DEVICE), diff --git a/resolve/LinSolverDirectCuSparseILU0.cpp b/resolve/LinSolverDirectCuSparseILU0.cpp index 23fdd2d5..0b3b9666 100644 --- a/resolve/LinSolverDirectCuSparseILU0.cpp +++ b/resolve/LinSolverDirectCuSparseILU0.cpp @@ -29,7 +29,7 @@ namespace ReSolve this->A_ = (matrix::Csr*) A; index_type n = A_->getNumRows(); - index_type nnz = A_->getNnzExpanded(); + index_type nnz = A_->getNnz(); mem_.allocateArrayOnDevice(&d_ILU_vals_,nnz); //copy A values to a buffer first mem_.copyArrayDeviceToDevice(d_ILU_vals_, A_->getValues(ReSolve::memory::DEVICE), nnz); @@ -216,7 +216,7 @@ namespace ReSolve int error_sum = 0; this->A_ = A; index_type n = A_->getNumRows(); - index_type nnz = A_->getNnzExpanded(); + index_type nnz = A_->getNnz(); mem_.copyArrayDeviceToDevice(d_ILU_vals_, A_->getValues(ReSolve::memory::DEVICE), nnz); status_cusparse_ = cusparseDcsrilu02(workspace_->getCusparseHandle(), diff --git a/resolve/LinSolverDirectRocSolverRf.cpp b/resolve/LinSolverDirectRocSolverRf.cpp index 1fbc8775..b7af7428 100644 --- a/resolve/LinSolverDirectRocSolverRf.cpp +++ b/resolve/LinSolverDirectRocSolverRf.cpp @@ -60,11 +60,11 @@ namespace ReSolve status_rocblas_ = rocsolver_dcsrrf_analysis(workspace_->getRocblasHandle(), n, 1, - A_->getNnzExpanded(), + A_->getNnz(), A_->getRowData(ReSolve::memory::DEVICE), //kRowPtr_, A_->getColData(ReSolve::memory::DEVICE), //jCol_, A_->getValues(ReSolve::memory::DEVICE), //vals_, - M_->getNnzExpanded(), + M_->getNnz(), M_->getRowData(ReSolve::memory::DEVICE), M_->getColData(ReSolve::memory::DEVICE), M_->getValues(ReSolve::memory::DEVICE), //vals_, @@ -111,7 +111,7 @@ namespace ReSolve status_rocblas_ = rocsolver_dcsrrf_splitlu(workspace_->getRocblasHandle(), n, - M_->getNnzExpanded(), + M_->getNnz(), M_->getRowData(ReSolve::memory::DEVICE), M_->getColData(ReSolve::memory::DEVICE), M_->getValues(ReSolve::memory::DEVICE), //vals_, @@ -198,11 +198,11 @@ namespace ReSolve mem_.deviceSynchronize(); status_rocblas_ = rocsolver_dcsrrf_refactlu(workspace_->getRocblasHandle(), A_->getNumRows(), - A_->getNnzExpanded(), + A_->getNnz(), A_->getRowData(ReSolve::memory::DEVICE), //kRowPtr_, A_->getColData(ReSolve::memory::DEVICE), //jCol_, A_->getValues(ReSolve::memory::DEVICE), //vals_, - M_->getNnzExpanded(), + M_->getNnz(), M_->getRowData(ReSolve::memory::DEVICE), M_->getColData(ReSolve::memory::DEVICE), M_->getValues(ReSolve::memory::DEVICE), //OUTPUT, @@ -217,7 +217,7 @@ namespace ReSolve //split M, fill L and U with correct values status_rocblas_ = rocsolver_dcsrrf_splitlu(workspace_->getRocblasHandle(), A_->getNumRows(), - M_->getNnzExpanded(), + M_->getNnz(), M_->getRowData(ReSolve::memory::DEVICE), M_->getColData(ReSolve::memory::DEVICE), M_->getValues(ReSolve::memory::DEVICE), //vals_, diff --git a/resolve/LinSolverDirectRocSparseILU0.cpp b/resolve/LinSolverDirectRocSparseILU0.cpp index 876c7625..b996f47d 100644 --- a/resolve/LinSolverDirectRocSparseILU0.cpp +++ b/resolve/LinSolverDirectRocSparseILU0.cpp @@ -29,7 +29,7 @@ namespace ReSolve this->A_ = (matrix::Csr*) A; index_type n = A_->getNumRows(); - index_type nnz = A_->getNnzExpanded(); + index_type nnz = A_->getNnz(); mem_.allocateArrayOnDevice(&d_ILU_vals_,nnz); //copy A values to a buffer first mem_.copyArrayDeviceToDevice(d_ILU_vals_, A_->getValues(ReSolve::memory::DEVICE), nnz); @@ -190,7 +190,7 @@ namespace ReSolve int error_sum = 0; this->A_ = A; index_type n = A_->getNumRows(); - index_type nnz = A_->getNnzExpanded(); + index_type nnz = A_->getNnz(); mem_.copyArrayDeviceToDevice(d_ILU_vals_, A_->getValues(ReSolve::memory::DEVICE), nnz); status_rocsparse_ = rocsparse_dcsrilu0(workspace_->getRocsparseHandle(), @@ -215,7 +215,7 @@ namespace ReSolve status_rocsparse_ = rocsparse_dcsrsv_solve(workspace_->getRocsparseHandle(), rocsparse_operation_none, A_->getNumRows(), - A_->getNnzExpanded(), + A_->getNnz(), &(constants::ONE), descr_L_, d_ILU_vals_, //vals_, @@ -231,7 +231,7 @@ namespace ReSolve status_rocsparse_ = rocsparse_dcsrsv_solve(workspace_->getRocsparseHandle(), rocsparse_operation_none, A_->getNumRows(), - A_->getNnzExpanded(), + A_->getNnz(), &(constants::ONE), descr_U_, d_ILU_vals_, //vals_, @@ -257,7 +257,7 @@ namespace ReSolve status_rocsparse_ = rocsparse_dcsrsv_solve(workspace_->getRocsparseHandle(), rocsparse_operation_none, A_->getNumRows(), - A_->getNnzExpanded(), + A_->getNnz(), &(constants::ONE), descr_L_, d_ILU_vals_, //vals_, @@ -273,7 +273,7 @@ namespace ReSolve status_rocsparse_ = rocsparse_dcsrsv_solve(workspace_->getRocsparseHandle(), rocsparse_operation_none, A_->getNumRows(), - A_->getNnzExpanded(), + A_->getNnz(), &(constants::ONE), descr_U_, d_ILU_vals_, //vals_, diff --git a/resolve/LinSolverDirectSerialILU0.cpp b/resolve/LinSolverDirectSerialILU0.cpp index 33cec473..57d58f82 100644 --- a/resolve/LinSolverDirectSerialILU0.cpp +++ b/resolve/LinSolverDirectSerialILU0.cpp @@ -33,7 +33,7 @@ namespace ReSolve { this->A_ = (matrix::Csr*) A; index_type n = A_->getNumRows(); - index_type nnz = A_->getNnzExpanded(); + index_type nnz = A_->getNnz(); h_ILU_vals_ = new real_type[nnz]; h_aux1_ = new real_type[n]; diff --git a/resolve/matrix/CMakeLists.txt b/resolve/matrix/CMakeLists.txt index 6a92975e..ad38dd5d 100644 --- a/resolve/matrix/CMakeLists.txt +++ b/resolve/matrix/CMakeLists.txt @@ -15,7 +15,6 @@ set(Matrix_SRC Coo.cpp MatrixHandler.cpp MatrixHandlerCpu.cpp - Utilities.cpp ) # C++ code that depends on CUDA SDK libraries @@ -36,7 +35,6 @@ set(Matrix_HEADER_INSTALL Csr.hpp Csc.hpp MatrixHandler.hpp - Utilities.hpp ) # Build shared library ReSolve::matrix diff --git a/resolve/matrix/Coo.cpp b/resolve/matrix/Coo.cpp index a0ea657c..c4622116 100644 --- a/resolve/matrix/Coo.cpp +++ b/resolve/matrix/Coo.cpp @@ -26,6 +26,101 @@ namespace ReSolve { } + /** + * @brief Hijacking constructor + */ + matrix::Coo::Coo(index_type n, + index_type m, + index_type nnz, + bool symmetric, + bool expanded, + index_type** rows, + index_type** cols, + real_type** vals, + memory::MemorySpace memspaceSrc, + memory::MemorySpace memspaceDst) + : Sparse(n, m, nnz, symmetric, expanded) + { + using namespace memory; + int control = -1; + if ((memspaceSrc == memory::HOST) && (memspaceDst == memory::HOST)) { control = 0;} + if ((memspaceSrc == memory::HOST) && (memspaceDst == memory::DEVICE)){ control = 1;} + if ((memspaceSrc == memory::DEVICE) && (memspaceDst == memory::HOST)) { control = 2;} + if ((memspaceSrc == memory::DEVICE) && (memspaceDst == memory::DEVICE)){ control = 3;} + + switch (control) + { + case 0: // cpu->cpu + // Set host data + h_row_data_ = *rows; + h_col_data_ = *cols; + h_val_data_ = *vals; + h_data_updated_ = true; + owns_cpu_vals_ = true; + owns_cpu_data_ = true; + // Make sure there is no device data. + if (d_row_data_ || d_col_data_ || d_val_data_) { + out::error() << "Device data unexpectedly allocated. " + << "Possible bug in matrix::Sparse class.\n"; + } + // Hijack data from the source + *rows = nullptr; + *cols = nullptr; + *vals = nullptr; + break; + case 2: // gpu->cpu + // Set device data and copy it to host + d_row_data_ = *rows; + d_col_data_ = *cols; + d_val_data_ = *vals; + d_data_updated_ = true; + owns_gpu_vals_ = true; + owns_gpu_data_ = true; + copyData(memspaceDst); + // Hijack data from the source + *rows = nullptr; + *cols = nullptr; + *vals = nullptr; + break; + case 1: // cpu->gpu + // Set host data and copy it to device + h_row_data_ = *rows; + h_col_data_ = *cols; + h_val_data_ = *vals; + h_data_updated_ = true; + owns_cpu_vals_ = true; + owns_cpu_data_ = true; + copyData(memspaceDst); + // Hijack data from the source + *rows = nullptr; + *cols = nullptr; + *vals = nullptr; + break; + case 3: // gpu->gpu + // Set device data + d_row_data_ = *rows; + d_col_data_ = *cols; + d_val_data_ = *vals; + d_data_updated_ = true; + owns_gpu_vals_ = true; + owns_gpu_data_ = true; + // Make sure there is no device data. + if (h_row_data_ || h_col_data_ || h_val_data_) { + out::error() << "Host data unexpectedly allocated. " + << "Possible bug in matrix::Sparse class.\n"; + } + // Hijack data from the source + *rows = nullptr; + *cols = nullptr; + *vals = nullptr; + break; + default: + out::error() << "Coo constructor failed! " + << "Possible bug in memory spaces setting.\n"; + break; + } + } + matrix::Coo::~Coo() { } @@ -81,7 +176,6 @@ namespace ReSolve //four cases (for now) index_type nnz_current = nnz_; - if (is_expanded_) {nnz_current = nnz_expanded_;} setNotUpdated(); int control=-1; if ((memspaceIn == memory::HOST) && (memspaceOut == memory::HOST)){ control = 0;} @@ -163,7 +257,6 @@ namespace ReSolve int matrix::Coo::allocateMatrixData(memory::MemorySpace memspace) { index_type nnz_current = nnz_; - if (is_expanded_) {nnz_current = nnz_expanded_;} destroyMatrixData(memspace);//just in case if (memspace == memory::HOST) { @@ -194,9 +287,6 @@ namespace ReSolve using namespace ReSolve::memory; index_type nnz_current = nnz_; - if (is_expanded_) { - nnz_current = nnz_expanded_; - } switch (memspaceOut) { case HOST: @@ -249,12 +339,12 @@ namespace ReSolve * * @param out - Output stream where the matrix data is printed */ - void matrix::Coo::print(std::ostream& out) + void matrix::Coo::print(std::ostream& out, index_type indexing_base) { out << std::scientific << std::setprecision(std::numeric_limits::digits10); for(int i = 0; i < nnz_; ++i) { - out << h_row_data_[i] << " " - << h_col_data_[i] << " " + out << h_row_data_[i] + indexing_base << " " + << h_col_data_[i] + indexing_base << " " << h_val_data_[i] << "\n"; } } diff --git a/resolve/matrix/Coo.hpp b/resolve/matrix/Coo.hpp index 987368f8..b12b4227 100644 --- a/resolve/matrix/Coo.hpp +++ b/resolve/matrix/Coo.hpp @@ -13,6 +13,16 @@ namespace ReSolve { namespace matrix { index_type nnz, bool symmetric, bool expanded); + Coo(index_type n, + index_type m, + index_type nnz, + bool symmetric, + bool expanded, + index_type** rows, + index_type** cols, + real_type** vals, + memory::MemorySpace memspaceSrc, + memory::MemorySpace memspaceDst); ~Coo(); virtual index_type* getRowData(memory::MemorySpace memspace); @@ -24,7 +34,7 @@ namespace ReSolve { namespace matrix { virtual int allocateMatrixData(memory::MemorySpace memspace); - virtual void print(std::ostream& file_out = std::cout); + virtual void print(std::ostream& file_out = std::cout, index_type indexing_base = 0); virtual int copyData(memory::MemorySpace memspaceOut); }; diff --git a/resolve/matrix/Csc.cpp b/resolve/matrix/Csc.cpp index c25bdc09..9287c847 100644 --- a/resolve/matrix/Csc.cpp +++ b/resolve/matrix/Csc.cpp @@ -77,9 +77,7 @@ namespace ReSolve memory::MemorySpace memspaceOut) { index_type nnz_current = nnz_; - if (is_expanded_) { - nnz_current = nnz_expanded_; - } + //four cases (for now) int control = -1; setNotUpdated(); @@ -168,7 +166,6 @@ namespace ReSolve int matrix::Csc::allocateMatrixData(memory::MemorySpace memspace) { index_type nnz_current = nnz_; - if (is_expanded_) {nnz_current = nnz_expanded_;} destroyMatrixData(memspace);//just in case if (memspace == memory::HOST) { @@ -199,9 +196,6 @@ namespace ReSolve using namespace ReSolve::memory; index_type nnz_current = nnz_; - if (is_expanded_) { - nnz_current = nnz_expanded_; - } switch(memspaceOut) { case HOST: @@ -254,13 +248,13 @@ namespace ReSolve * * @param out - Output stream where the matrix data is printed */ - void matrix::Csc::print(std::ostream& out) + void matrix::Csc::print(std::ostream& out, index_type indexing_base) { out << std::scientific << std::setprecision(std::numeric_limits::digits10); for(index_type i = 0; i < m_; ++i) { for (index_type j = h_col_data_[i]; j < h_col_data_[i+1]; ++j) { - out << h_row_data_[j] << " " - << i << " " + out << h_row_data_[j] + indexing_base << " " + << i + indexing_base << " " << h_val_data_[j] << "\n"; } } diff --git a/resolve/matrix/Csc.hpp b/resolve/matrix/Csc.hpp index 824072ae..3395ad6c 100644 --- a/resolve/matrix/Csc.hpp +++ b/resolve/matrix/Csc.hpp @@ -33,7 +33,7 @@ namespace ReSolve { namespace matrix { virtual int allocateMatrixData(memory::MemorySpace memspace); - virtual void print(std::ostream& file_out = std::cout); + virtual void print(std::ostream& file_out = std::cout, index_type indexing_base = 0); virtual int copyData(memory::MemorySpace memspaceOut); }; diff --git a/resolve/matrix/Csr.cpp b/resolve/matrix/Csr.cpp index 268898bc..15c786c2 100644 --- a/resolve/matrix/Csr.cpp +++ b/resolve/matrix/Csr.cpp @@ -6,7 +6,6 @@ #include "Csr.hpp" #include "Coo.hpp" #include -#include namespace ReSolve { @@ -28,16 +27,6 @@ namespace ReSolve { } - matrix::Csr::Csr(matrix::Coo* A_coo, memory::MemorySpace memspace) - : Sparse(A_coo->getNumRows(), - A_coo->getNumColumns(), - A_coo->getNnz(), - A_coo->symmetric(), - A_coo->expanded()) - { - matrix::coo2csr(A_coo, this, memspace); - } - /** * @brief Hijacking constructor * @@ -201,7 +190,6 @@ namespace ReSolve { //four cases (for now) index_type nnz_current = nnz_; - if (is_expanded_) {nnz_current = nnz_expanded_;} setNotUpdated(); int control = -1; if ((memspaceIn == memory::HOST) && (memspaceOut == memory::HOST)) { control = 0;} @@ -285,7 +273,6 @@ namespace ReSolve int matrix::Csr::allocateMatrixData(memory::MemorySpace memspace) { index_type nnz_current = nnz_; - if (is_expanded_) {nnz_current = nnz_expanded_;} destroyMatrixData(memspace);//just in case if (memspace == memory::HOST) { @@ -316,9 +303,6 @@ namespace ReSolve using namespace ReSolve::memory; index_type nnz_current = nnz_; - if (is_expanded_) { - nnz_current = nnz_expanded_; - } switch (memspaceOut) { case HOST: @@ -367,29 +351,19 @@ namespace ReSolve } // switch } - int matrix::Csr::updateFromCoo(matrix::Coo* A_coo, memory::MemorySpace memspaceOut) - { - assert(n_ == A_coo->getNumRows()); - assert(m_ == A_coo->getNumColumns()); - assert(nnz_ == A_coo->getNnz()); - assert(is_symmetric_ == A_coo->symmetric()); // <- Do we need to check for this? - - return matrix::coo2csr(A_coo, this, memspaceOut); - } - /** * @brief Prints matrix data. * * @param out - Output stream where the matrix data is printed */ - void matrix::Csr::print(std::ostream& out) + void matrix::Csr::print(std::ostream& out, index_type indexing_base) { out << std::scientific << std::setprecision(std::numeric_limits::digits10); for(index_type i = 0; i < n_; ++i) { for (index_type j = h_row_data_[i]; j < h_row_data_[i+1]; ++j) { - out << i << " " - << h_col_data_[j] << " " + out << i + indexing_base << " " + << h_col_data_[j] + indexing_base << " " << h_val_data_[j] << "\n"; } } diff --git a/resolve/matrix/Csr.hpp b/resolve/matrix/Csr.hpp index a3896a1d..2f597cd5 100644 --- a/resolve/matrix/Csr.hpp +++ b/resolve/matrix/Csr.hpp @@ -30,8 +30,6 @@ namespace ReSolve { namespace matrix { memory::MemorySpace memspaceSrc, memory::MemorySpace memspaceDst); - Csr(matrix::Coo* mat, memory::MemorySpace memspace); - ~Csr(); virtual index_type* getRowData(memory::MemorySpace memspace); @@ -43,12 +41,10 @@ namespace ReSolve { namespace matrix { virtual int allocateMatrixData(memory::MemorySpace memspace); - virtual void print(std::ostream& file_out = std::cout); + virtual void print(std::ostream& file_out = std::cout, index_type indexing_base = 0); virtual int copyData(memory::MemorySpace memspaceOut); - int updateFromCoo(matrix::Coo* mat, memory::MemorySpace memspaceOut); - }; }} // namespace ReSolve::matrix diff --git a/resolve/matrix/MatrixHandler.cpp b/resolve/matrix/MatrixHandler.cpp index b61a93b5..6682a315 100644 --- a/resolve/matrix/MatrixHandler.cpp +++ b/resolve/matrix/MatrixHandler.cpp @@ -6,7 +6,6 @@ #include #include #include -#include #include "MatrixHandler.hpp" #include "MatrixHandlerCpu.hpp" @@ -109,16 +108,6 @@ namespace ReSolve { } } - /** - * @brief Converts COO to CSR matrix format. - * - * Conversion takes place on CPU, and then CSR matrix is copied to `memspace`. - */ - int MatrixHandler::coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr, memory::MemorySpace memspace) - { - return matrix::coo2csr(A_coo, A_csr, memspace); - } - /** * @brief Matrix vector product: result = alpha * A * x + beta * result * diff --git a/resolve/matrix/MatrixHandler.hpp b/resolve/matrix/MatrixHandler.hpp index b1de27e7..f34ec32c 100644 --- a/resolve/matrix/MatrixHandler.hpp +++ b/resolve/matrix/MatrixHandler.hpp @@ -53,7 +53,6 @@ namespace ReSolve { ~MatrixHandler(); int csc2csr(matrix::Csc* A_csc, matrix::Csr* A_csr, memory::MemorySpace memspace); - int coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr, memory::MemorySpace memspace); /// Should compute vec_result := alpha*A*vec_x + beta*vec_result, but at least on cpu alpha and beta are flipped int matvec(matrix::Sparse* A, diff --git a/resolve/matrix/MatrixHandlerCuda.cpp b/resolve/matrix/MatrixHandlerCuda.cpp index 7f6746a7..915fdde4 100644 --- a/resolve/matrix/MatrixHandlerCuda.cpp +++ b/resolve/matrix/MatrixHandlerCuda.cpp @@ -57,7 +57,7 @@ namespace ReSolve { status = cusparseCreateCsr(&matA, A->getNumRows(), A->getNumColumns(), - A->getNnzExpanded(), + A->getNnz(), A->getRowData(memory::DEVICE), A->getColData(memory::DEVICE), A->getValues( memory::DEVICE), @@ -138,7 +138,7 @@ namespace ReSolve { } matrix_row_sums(A->getNumRows(), - A->getNnzExpanded(), + A->getNnz(), A->getRowData(memory::DEVICE), A->getValues(memory::DEVICE), d_r); diff --git a/resolve/matrix/MatrixHandlerHip.cpp b/resolve/matrix/MatrixHandlerHip.cpp index 44d9d8be..6cb56a41 100644 --- a/resolve/matrix/MatrixHandlerHip.cpp +++ b/resolve/matrix/MatrixHandlerHip.cpp @@ -59,7 +59,7 @@ namespace ReSolve { rocsparse_operation_none, A->getNumRows(), A->getNumColumns(), - A->getNnzExpanded(), + A->getNnz(), descrA, A->getValues( memory::DEVICE), A->getRowData(memory::DEVICE), @@ -77,7 +77,7 @@ namespace ReSolve { rocsparse_operation_none, A->getNumRows(), A->getNumColumns(), - A->getNnzExpanded(), + A->getNnz(), alpha, descrA, A->getValues( memory::DEVICE), @@ -127,7 +127,7 @@ namespace ReSolve { mem_.deviceSynchronize(); matrix_row_sums(A->getNumRows(), - A->getNnzExpanded(), + A->getNnz(), A->getRowData(memory::DEVICE), A->getValues(memory::DEVICE), d_r); diff --git a/resolve/matrix/Sparse.cpp b/resolve/matrix/Sparse.cpp index c795f0c9..7a59c553 100644 --- a/resolve/matrix/Sparse.cpp +++ b/resolve/matrix/Sparse.cpp @@ -30,7 +30,6 @@ namespace ReSolve { { this->is_symmetric_ = false; this->is_expanded_ = true; //default is a normal non-symmetric fully expanded matrix - this->nnz_expanded_ = nnz; setNotUpdated(); @@ -70,11 +69,6 @@ namespace ReSolve { is_symmetric_{symmetric}, is_expanded_{expanded} { - if (is_expanded_) { - this->nnz_expanded_ = nnz_; - } else { - this->nnz_expanded_ = 0; - } setNotUpdated(); //set everything to nullptr @@ -141,16 +135,6 @@ namespace ReSolve { return this->nnz_; } - /** - * @brief get number of non-zeros in expanded matrix. - * - * @return number of non-zeros in expanded matrix. - */ - index_type matrix::Sparse::getNnzExpanded() - { - return this->nnz_expanded_; - } - /** * @brief check if matrix is symmetric. * @@ -191,16 +175,6 @@ namespace ReSolve { this->is_expanded_ = expanded; } - /** - * @brief Set number of non-zeros in expanded matrix. - * - * @param[in] nnz_expanded_new - new number of non-zeros in expanded matrix - */ - void matrix::Sparse::setNnzExpanded(index_type nnz_expanded_new) - { - this->nnz_expanded_ = nnz_expanded_new; - } - /** * @brief Set number of non-zeros. * @@ -355,7 +329,6 @@ namespace ReSolve { { index_type nnz_current = nnz_; - if (is_expanded_) {nnz_current = nnz_expanded_;} //four cases (for now) setNotUpdated(); int control=-1; diff --git a/resolve/matrix/Sparse.hpp b/resolve/matrix/Sparse.hpp index ed5a66a5..d971f24f 100644 --- a/resolve/matrix/Sparse.hpp +++ b/resolve/matrix/Sparse.hpp @@ -32,13 +32,11 @@ namespace ReSolve { namespace matrix { index_type getNumRows(); index_type getNumColumns(); index_type getNnz(); - index_type getNnzExpanded(); bool symmetric(); bool expanded(); void setSymmetric(bool symmetric); void setExpanded(bool expanded); - void setNnzExpanded(index_type nnz_expanded_new); void setNnz(index_type nnz_new); // for resetting when removing duplicates int setUpdated(memory::MemorySpace what); @@ -54,7 +52,7 @@ namespace ReSolve { namespace matrix { int destroyMatrixData(memory::MemorySpace memspace); - virtual void print(std::ostream& file_out) = 0; + virtual void print(std::ostream& file_out, index_type indexing_base) = 0; virtual int copyData(memory::MemorySpace memspaceOut) = 0; @@ -71,7 +69,6 @@ namespace ReSolve { namespace matrix { index_type n_{0}; ///< number of rows index_type m_{0}; ///< number of columns index_type nnz_{0}; ///< number of non-zeros - index_type nnz_expanded_{0}; ///< number of non-zeros in an expanded matrix (for symmetric matrices) bool is_symmetric_{false}; ///< symmetry flag bool is_expanded_{false}; ///< "expanded" flag diff --git a/resolve/matrix/Utilities.cpp b/resolve/matrix/Utilities.cpp deleted file mode 100644 index bd5a6f9b..00000000 --- a/resolve/matrix/Utilities.cpp +++ /dev/null @@ -1,206 +0,0 @@ -#include -#include -#include -#include -#include -#include "Utilities.hpp" - -namespace ReSolve -{ - /// @brief Helper class for COO matrix sorting - class IndexValuePair - { - public: - IndexValuePair() : idx_(0), value_(0.0) - {} - ~IndexValuePair() - {} - void setIdx (index_type new_idx) - { - idx_ = new_idx; - } - void setValue (real_type new_value) - { - value_ = new_value; - } - - index_type getIdx() - { - return idx_; - } - real_type getValue() - { - return value_; - } - - bool operator < (const IndexValuePair& str) const - { - return (idx_ < str.idx_); - } - - private: - index_type idx_; - real_type value_; - }; - - using out = io::Logger; - namespace matrix - { - /** - * @brief Creates a CSR from a COO matrix. - * - * @param[in] A_coo - * @param[out] A_csr - * @return int - Error code, 0 if successful. - * - * @pre `A_coo` is a valid sparse matrix in unordered COO format. - * Duplicates are allowed. Up-to-date values and indices must be - * on host. - * - * @post `A_csr` is representing the same matrix as `A_coo` but in - * _general_ CSR format. `A_csr` is allocated and stored on host. - * - * @invariant `A_coo` is not changed. - */ - int coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr, memory::MemorySpace memspace) - { - //count nnzs first - index_type nnz_unpacked = 0; - const index_type nnz = A_coo->getNnz(); - const index_type n = A_coo->getNumRows(); - const bool symmetric = A_coo->symmetric(); - const bool expanded = A_coo->expanded(); - - index_type* nnz_counts = new index_type[n]; ///< Number of elements per row - std::fill_n(nnz_counts, n, 0); - const index_type* coo_rows = A_coo->getRowData(memory::HOST); - const index_type* coo_cols = A_coo->getColData(memory::HOST); - const real_type* coo_vals = A_coo->getValues( memory::HOST); - - index_type* diag_control = new index_type[n]; // for DEDUPLICATION of the diagonal - std::fill_n(diag_control, n, 0); - index_type nnz_unpacked_no_duplicates = 0; - index_type nnz_no_duplicates = nnz; - - // Count matrix elements - for (index_type i = 0; i < nnz; ++i) { - nnz_counts[coo_rows[i]]++; - nnz_unpacked++; - nnz_unpacked_no_duplicates++; - // Count elements added after expansion - if ((coo_rows[i] != coo_cols[i]) && (symmetric) && (!expanded)) { - nnz_counts[coo_cols[i]]++; - nnz_unpacked++; - nnz_unpacked_no_duplicates++; - } - // Bookkeeping of diagonal elements that were counted - if (coo_rows[i] == coo_cols[i]) { - if (diag_control[coo_rows[i]] > 0) { - // duplicate - nnz_unpacked_no_duplicates--; - nnz_no_duplicates--; - } - diag_control[coo_rows[i]]++; - } - } - A_csr->setExpanded(true); - A_csr->setNnzExpanded(nnz_unpacked_no_duplicates); - - // Allocate matrix format conversion workspace - index_type* csr_ia = new index_type[n+1]; - std::fill_n(csr_ia, n + 1, 0); - index_type* csr_ja = new index_type[nnz_unpacked]; - real_type* csr_a = new real_type[nnz_unpacked]; - index_type* nnz_shifts = new index_type[n]; - std::fill_n(nnz_shifts, n , 0); - - IndexValuePair* tmp = new IndexValuePair[nnz_unpacked]; - - // Set CSR row pointers - csr_ia[0] = 0; - for (index_type i = 1; i < n + 1; ++i) { - csr_ia[i] = csr_ia[i - 1] + nnz_counts[i - 1] - (diag_control[i-1] - 1); - } - - int r, start; - - for (index_type i = 0; i < nnz; ++i) { - //which row - r = coo_rows[i]; - start = csr_ia[r]; - - if ((start + nnz_shifts[r]) > nnz_unpacked) { - out::warning() << "index out of bounds (case 1) start: " << start - << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; - } - if ((r == coo_cols[i]) && (diag_control[r] > 1)) {//diagonal, and there are duplicates - bool already_there = false; - for (index_type j = start; j < start + nnz_shifts[r]; ++j) - { - index_type c = tmp[j].getIdx(); - if (c == r) { - real_type val = tmp[j].getValue(); - val += coo_vals[i]; - tmp[j].setValue(val); - already_there = true; - out::warning() << " duplicate found, row " << c - << " adding in place " << j << " current value: " << val << std::endl; - } - } - if (!already_there) { // first time this duplicates appears - - tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]); - tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); - - nnz_shifts[r]++; - } - } else { // non-diagonal - tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]); - tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); - nnz_shifts[r]++; - - if ((coo_rows[i] != coo_cols[i]) && (symmetric == 1)) { - r = coo_cols[i]; - start = csr_ia[r]; - - if ((start + nnz_shifts[r]) > nnz_unpacked) { - out::warning() << "index out of bounds (case 2) start: " << start - << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; - } - tmp[start + nnz_shifts[r]].setIdx(coo_rows[i]); - tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); - nnz_shifts[r]++; - } - } - } - - //now sort whatever is inside rows - for (int i = 0; i < n; ++i) { - //now sorting (and adding 1) - int colStart = csr_ia[i]; - int colEnd = csr_ia[i + 1]; - int length = colEnd - colStart; - std::sort(&tmp[colStart], &tmp[colStart] + length); - } - - for (index_type i = 0; i < nnz_unpacked; ++i) - { - csr_ja[i] = tmp[i].getIdx(); - csr_a[i] = tmp[i].getValue(); - } - A_csr->setNnz(nnz_no_duplicates); - A_csr->updateData(csr_ia, csr_ja, csr_a, memory::HOST, memspace); - - delete [] nnz_counts; - delete [] tmp; - delete [] nnz_shifts; - delete [] csr_ia; - delete [] csr_ja; - delete [] csr_a; - delete [] diag_control; - - return 0; - - } - } -} \ No newline at end of file diff --git a/resolve/matrix/Utilities.hpp b/resolve/matrix/Utilities.hpp deleted file mode 100644 index ecaee79e..00000000 --- a/resolve/matrix/Utilities.hpp +++ /dev/null @@ -1,16 +0,0 @@ -#pragma once - -#include - -namespace ReSolve -{ - namespace matrix - { - // Forward declarations - class Coo; - class Csr; - - /// @brief Converts symmetric or general COO to general CSR matrix - int coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr, memory::MemorySpace memspace); - } -} diff --git a/resolve/matrix/io.cpp b/resolve/matrix/io.cpp index d4cdd90e..17eaa014 100644 --- a/resolve/matrix/io.cpp +++ b/resolve/matrix/io.cpp @@ -3,198 +3,752 @@ #include #include #include +#include #include #include #include +#include #include "io.hpp" -namespace ReSolve { namespace io { +namespace ReSolve +{ + + /** + * @class MatrixElementTriplet + * + * @brief Helper class for COO matrix sorting. + * + * Contains triplet of row index, column index and the value of a sparse + * matrix element, as well as methods and operator overloads for its + * management. + * + * The entire code is in this file. Its scope is to support matrix file I/O + * only. + * + */ + class MatrixElementTriplet + { + public: + /// Default constructor initializes all to zero. + MatrixElementTriplet() : rowidx_(0), colidx_(0), value_(0.0) + {} + + /// Constructor that initializes row and column indices and the element value. + MatrixElementTriplet(index_type i, index_type j, real_type v) : rowidx_(i), colidx_(j), value_(v) + {} + + ~MatrixElementTriplet() = default; + + /// Set the row and column indices and the element value. + void set(index_type rowidx, index_type colidx, real_type value) + { + rowidx_ = rowidx; + colidx_ = colidx; + value_ = value; + } + + index_type getRowIdx() const + { + return rowidx_; + } + index_type getColIdx() const + { + return colidx_; + } + real_type getValue() const + { + return value_; + } + + /** + * @brief Overload of `<` operator + * + * Ensures that matrix elements stored in MatrixElementTriplet will be + * sorted by their indices in a row-major order. + * + */ + bool operator < (const MatrixElementTriplet& t) const + { + if (rowidx_ < t.rowidx_) + return true; + + if ((rowidx_ == t.rowidx_) && (colidx_ < t.colidx_)) + return true; + + return false; + } + + /** + * @brief Overload of `==` operator. + * + * This overload is used to indicate when two different instances of + * MatrixElementTriplet correspond to the same matrix element. + */ + bool operator == (const MatrixElementTriplet& str) const + { + return (rowidx_ == str.rowidx_) && (colidx_ == str.colidx_); + } + + /** + * @brief Overload of `+=` operator. + * + * @param t - Triplet to be added in place. + * @return MatrixElementTriplet& - reference to `*this`. + * + * This overload is used to merge duplicates in sparse matrix in COO + * format. It will return error and leave `*this` unchanged if the + * argument corresponds to an element with different row or column + * indices. + */ + MatrixElementTriplet& operator += (const MatrixElementTriplet t) + { + if ((rowidx_ != t.rowidx_) || (colidx_ != t.colidx_)) { + io::Logger::error() << "Adding values into non-matching triplet.\n"; + return *this; + } + value_ += t.value_; + return *this; + } + + /// Utility to print indices (0 index base). + void print() const + { + // Add 1 to indices to restore indexing from MM format + std::cout << rowidx_ << " " << colidx_ << " " << value_ << "\n"; + } - matrix::Coo* readMatrixFromFile(std::istream& file) + private: + index_type rowidx_{0}; + index_type colidx_{0}; + real_type value_{0.0}; + }; + + + namespace io { - if(!file) { - Logger::error() << "Empty input to readMatrixFromFile function ... \n" << std::endl; - return nullptr; - } - std::stringstream ss; - std::string line; - index_type i = 0; - index_type m, n, nnz; - bool symmetric = false; - bool expanded = true; - std::getline(file, line); - //symmetric? - size_t found = line.find("symmetric"); - if (found != std::string::npos) { - symmetric = true; - expanded = false; - } - while (line.at(0) == '%') { - std::getline(file, line); - // std::cout<> n >> m >> nnz; - //create matrix object - matrix::Coo* A = new matrix::Coo(n, m, nnz,symmetric, expanded ); - //create coo arrays - index_type* coo_rows = new index_type[nnz]; - index_type* coo_cols = new index_type[nnz]; - real_type* coo_vals = new real_type[nnz]; - i = 0; - index_type a, b; - real_type c; - while (file >> a >> b >> c) { - coo_rows[i] = a - 1; - coo_cols[i] = b - 1; - coo_vals[i] = c; - i++; + // Static helper functionsdeclarations + static void createMatrixFromFileAsList(std::istream& file, + bool is_expand_symmetric, + std::list& tmp, + index_type& n, + index_type& m, + index_type& nnz, + bool& symmetric, + bool& expanded); + static void createMatrixFromFileAsList(std::istream& file, + matrix::Sparse* A, + std::list& tmp); + // static void print_list(std::list& l); + static int loadToList(std::istream& file, bool is_expand_symmetric, std::list& tmp); + static int removeDuplicates(std::list& tmp); + static int copyListToCoo(const std::list& tmp, matrix::Coo* A); + static int copyListToCsr(const std::list& tmp, matrix::Csr* A); + + + /** + * @brief Create a COO matrix and populate it with data from Matrix Market + * file. + * + * @param file - input Matrix Market file + * @param is_expand_symmetric - whether to expand symmetric matrix to general format + * @return matrix::Coo* - pointer to COO matrix + * + * @pre file is a valid std::istream with Matrix Market data. + * @pre input data is in valid Matrix Market format. + * @post Valid COO matrix sorted in row major order and without duplicates + * is created. + */ + matrix::Coo* createCooFromFile(std::istream& file, bool is_expand_symmetric) + { + if(!file) { + Logger::error() << "Empty input to createCooFromFile function ... \n" << std::endl; + return nullptr; + } + + index_type m = 0, n = 0, nnz = 0; + bool symmetric = false; + bool expanded = true; + + std::list tmp; + + createMatrixFromFileAsList(file, is_expand_symmetric, tmp, n, m, nnz, symmetric, expanded); + + // Create matrix + matrix::Coo* B = new matrix::Coo(n, m, nnz, symmetric, expanded); + B->allocateMatrixData(memory::HOST); + + copyListToCoo(tmp, B); + + return B; } - A->setMatrixData(coo_rows, coo_cols, coo_vals, memory::HOST); - return A; - } + /** + * @brief + * + * @param file - input Matrix Market file + * @param is_expand_symmetric - whether to expand symmetric matrix to general format + * @return matrix::Csr* - pointer to COO matrix + * + * @pre file is a valid std::istream with Matrix Market data. + * @pre input data is in valid Matrix Market format. + * @post Valid CSR matrix sorted in row major order and without duplicates + * is created. + */ + matrix::Csr* createCsrFromFile(std::istream& file, bool is_expand_symmetric) + { + if(!file) { + Logger::error() << "Empty input to createCooFromFile function ... \n" << std::endl; + return nullptr; + } - real_type* readRhsFromFile(std::istream& file) - { - if(!file) { - Logger::error() << "Empty input to " << __func__ << " function ... \n" << std::endl; - return nullptr; + index_type m = 0, n = 0, nnz = 0; + bool symmetric = false; + bool expanded = true; + + std::list tmp; + + createMatrixFromFileAsList(file, is_expand_symmetric, tmp, n, m, nnz, symmetric, expanded); + + // Create matrix + matrix::Csr* A = new matrix::Csr(n, m, nnz, symmetric, expanded); + A->allocateMatrixData(memory::HOST); + + copyListToCsr(tmp, A); + + return A; } - std::stringstream ss; - std::string line; - index_type i = 0; - index_type n, m; + /** + * @brief Imports vector data from a Matrix Market file. + * + * @param file - std::istream to Matrix Market file (data). + * @return real_type* - pointer to array with (dense) vector entries. + * + * @pre `file` is a valid std::istream with Matrix Market data. + * @pre Input data is in valid Matrix Market format. + * @post A raw array with vector data is created. + * + */ + real_type* createArrayFromFile(std::istream& file) + { + if(!file) { + Logger::error() << "Empty input to " << __func__ << " function ... \n" << std::endl; + return nullptr; + } + + std::stringstream ss; + std::string line; + index_type i = 0; + index_type n, m; - std::getline(file, line); - while (line.at(0) == '%') { - std::getline(file, line); - // std::cout << line << std::endl; + std::getline(file, line); + while (line.at(0) == '%') { + std::getline(file, line); + // std::cout << line << std::endl; + } + ss << line; + ss >> n >> m ; + + real_type* vec = new real_type[n]; + real_type a; + while (file >> a) { + vec[i] = a; + i++; + } + return vec; } - ss << line; - ss >> n >> m ; - - real_type* vec = new real_type[n]; - real_type a; - while (file >> a){ - vec[i] = a; - i++; + + /** + * @brief Reads data from a Matrix Market file and updates COO matrix A. + * + * Compute complexity of this function is O(NNZ*log(NNZ)). There is an + * overload of this function that generates a CSR matrix. + * + * @param file - std::istream to Matrix Market file (data). + * @param A - output COO matrix. + * + * @pre `file` is a valid std::istream with Matrix Market data. + * @pre Input data is in valid Matrix Market format. + * @pre Size of matrix stored in the Matrix Market file matches the size of A. + * @post Valid COO matrix sorted in row major order and without duplicates + * is created. + */ + void updateMatrixFromFile(std::istream& file, matrix::Coo* A) + { + if(!file) { + Logger::error() << "Empty input to createCooFromFile function ..." << std::endl; + return; + } + + std::list tmp; + + createMatrixFromFileAsList(file, A, tmp); + + // Populate COO matrix. Complexity O(NNZ) + copyListToCoo(tmp, A); } - return vec; - } - void readAndUpdateMatrix(std::istream& file, matrix::Coo* A) - { - if(!file) { - Logger::error() << "Empty input to readMatrixFromFile function ..." << std::endl; - return; + /** + * @brief Reads data from a Matrix Market file and updates CSR matrix A. + * + * Compute complexity of this function is O(NNZ*log(NNZ)). There is an + * overload of this function that generates a COO matrix. + * + * @param file - std::istream to Matrix Market file (data). + * @param A - output CSR matrix. + * + * @pre `file` is a valid std::istream with Matrix Market data. + * @pre Input data is in valid Matrix Market format. + * @pre Size of matrix stored in the Matrix Market file matches the size of A. + * @post Valid CSR matrix sorted in row major order and without duplicates + * is created. + */ + void updateMatrixFromFile(std::istream& file, matrix::Csr* A) + { + if(!file) { + Logger::error() << "Empty input to updateMatrixFromFile function ..." << std::endl; + return; + } + + std::list tmp; + + createMatrixFromFileAsList(file, A, tmp); + + // Populate COO matrix. Complexity O(NNZ) + copyListToCsr(tmp, A); } - std::stringstream ss; - A->setExpanded(false); - std::string line; - index_type i = 0; - index_type m, n, nnz; - std::getline(file, line); - while (line.at(0) == '%') { - std::getline(file, line); - // std::cout << line << std::endl; + /** + * @brief Reads data from a Matrix Market file and updates array p_rhs. + * + * @param file - std::istream to Matrix Market file (data). + * @param p_rhs - pointer to a pointer to a raw array with vector data. + * + * @todo The righ-hand-side should be of vector type, not a raw array. With + * current implementation it is impossible to verify if the sufficient + * space is allocated to store all the data from the input file. Risk of + * writing past the end of the array is high. + */ + void updateArrayFromFile(std::istream& file, real_type** p_rhs) + { + if (!file) { + Logger::error() << "Empty input to updateArrayFromFile function ..." << std::endl; + return; + } + + real_type* rhs = *p_rhs; + std::stringstream ss; + std::string line; + index_type n, m; + + std::getline(file, line); + while (line.at(0) == '%') { + std::getline(file, line); + // std::cout<> n >> m ; + + if (rhs == nullptr) { + // std::cout << "Allocating array of size " << n << "\n"; + rhs = new real_type[n]; + } + real_type a; + index_type i = 0; + while (file >> a) { + rhs[i] = a; + // std::cout << i << ": " << a << "\n"; + i++; + } } - ss << line; - ss >> n >> m >> nnz; - if ((A->getNumRows() != n) || (A->getNumColumns() != m) || (A->getNnz() < nnz)) { - Logger::error() << "Wrong matrix size: " << A->getNumRows() - << "x" << A->getNumColumns() - << ", NNZ: " << A->getNnz() - << " Cannot update! \n "; - return; + /** + * @brief Writes matrix A to a file in Matrix Market format. + * + * @param A - input matrix. + * @param file_out - std::ostream to output file. + * @return int - 0 if successful, error code otherwise. + * + * @pre `A` is a valid sparse matrix. + * @post Valid Matrix Marked data is written to std::ostream. + * @invariant Matrix `A` elements are unchanged. + */ + int writeMatrixToFile(matrix::Sparse* A, std::ostream& file_out) + { + if (A == nullptr) { + Logger::error() << "Matrix pointer is NULL!\n"; + return -1; + } + + if (A->symmetric() && !A->expanded()) { + file_out << "%%MatrixMarket matrix coordinate real symmetric\n"; + } else { + file_out << "%%MatrixMarket matrix coordinate real general\n"; + } + file_out << "% Generated by Re::Solve \n"; + file_out << A->getNumRows() << " " + << A->getNumColumns() << " " + << A->getNnz() << "\n"; + + index_type indexing_base = 1; + A->print(file_out, indexing_base); + return 0; } - A->setNnz(nnz); - //create coo arrays - index_type* coo_rows = A->getRowData(memory::HOST); - index_type* coo_cols = A->getColData(memory::HOST); - real_type* coo_vals = A->getValues( memory::HOST); - i = 0; - index_type a, b; - real_type c; - while (file >> a >> b >> c) { - coo_rows[i] = a - 1; - coo_cols[i] = b - 1; - coo_vals[i] = c; - i++; + + /** + * @brief Writes vector data to a file in Matrix Market format. + * + * @param vec_x - Input vector. + * @param file_out - std::ostream to output file. + * @return int - 0 if successful, error code otherwise. + * + * @pre `vec_x` is a valid vector. + * @post Valid Matrix Market data is written to std::ostream. + * @invariant Elements of `vec_x` are unchanged. + */ + int writeVectorToFile(vector_type* vec_x, std::ostream& file_out) + { + real_type* x_data = vec_x->getData(memory::HOST); + + file_out << "%%MatrixMarket matrix array real general \n"; + file_out << "% Generated by Re::Solve \n"; + file_out << vec_x->getSize() << " " << 1 << "\n"; + for (int i = 0; i < vec_x->getSize(); ++i) { + file_out << std::setprecision(32) << std::scientific << x_data[i] << "\n"; + } + + return 0; } - } - void readAndUpdateRhs(std::istream& file, real_type** p_rhs) - { - if (!file) { - Logger::error() << "Empty input to readAndUpdateRhs function ..." << std::endl; - return; + + // + // Static helper functions + // + + /** + * @brief Reads Matrix Market data from std::istream and stores it into + * std::list. + * + * @param[in] file - std::istream to Matrix Market file (data). + * @param[in] is_expand_symmetric - whether to expand symmetric matrix to general format + * @param[out] tmp - std::list where to store matrix data + * @param[out] n - number of rows as read from Matrix Market file + * @param[out] m - number of columns as read from Matrix Market file + * @param[out] nnz - calculated number of matrix nonzeros + * @param[out] symmetric - if matrix is symmetric + * @param[out] expanded - if symmetric matrix is expanded to general format + * + * @pre `file` is a valid std::istream with Matrix Market data. + * @pre Input data is in valid Matrix Market format. + * @pre `tmp` is an empty list! + * @post Matrix size is set to values read from the Matrix Market file. + * @post `nnz` is number of nonzeros in the matrix after duplicates are + * removed and the matrix is (optionally) expanded. + * @post `symmetric` and `expanded` flags are set based on data read into + * `tmp` list. + * @post `tmp` list is overwritten with matrix elements read from the input + * stream. + */ + static void createMatrixFromFileAsList(std::istream& file, + bool is_expand_symmetric, + std::list& tmp, + index_type& n, + index_type& m, + index_type& nnz, + bool& symmetric, + bool& expanded) + { + std::stringstream ss; + std::string line; + m = 0; + n = 0; + nnz = 0; + symmetric = false; + expanded = true; + + // Parse header and check if matrix is symmetric + while (std::getline(file, line)) + { + if (line.at(0) != '%') + break; + if (line.find("symmetric") != std::string::npos) { + symmetric = true; + expanded = is_expand_symmetric; + } + } + + // Read the first line with matrix sizes + ss << line; + ss >> n >> m >> nnz; + + // Store COO data in the temporary workspace. Complexity O(NNZ) + loadToList(file, symmetric && is_expand_symmetric, tmp); + + // Sort tmp. Complexity O(NNZ*log(NNZ)). + tmp.sort(); + + // Deduplicate tmp. Complexity O(NNZ). + removeDuplicates(tmp); + + // Set nnz without duplicates and for possibly expanded matrix. + nnz = static_cast(tmp.size()); } - real_type* rhs = *p_rhs; - std::stringstream ss; - std::string line; - index_type n, m; + /** + * @brief + * + * @param[in] file - std::istream to Matrix Market file (data). + * @param[in] A - sparse matrix to be updated + * @param[out] tmp - temporary list with matrix entries + * + * @pre `file` is a valid std::istream with Matrix Market data. + * @pre Input data is in valid Matrix Market format. + * @pre `tmp` is an empty list! + * @pre Matrix size of `A` must match matrix size read from the Matrix Market + * file. + * @pre Number of nonzeros in `A` must not be smaller than the number of + * zeros read from the Matrix Market file. + * @post `tmp` list is overwritten with matrix elements read from the input + * stream. + * @invariant Elements of `A` are unchanged in this function but they are + * expected to be overwritten with values in `tmp` later in the code. + */ + static void createMatrixFromFileAsList(std::istream& file, + matrix::Sparse* A, + std::list& tmp) + { + std::stringstream ss; + std::string line; + // Default is a general matrix + bool symmetric = false; + + // Default is not to expand symmetric matrix + bool is_expand_symmetric = false; + + // Parse header and check if matrix is symmetric + std::getline(file, line); + if (line.find("symmetric") != std::string::npos) { + symmetric = true; + } + if (symmetric != A->symmetric()) { + Logger::error() << "In function updateMatrixFromFile:" + << "Source data does not match the symmetry of destination matrix.\n"; + } + // If the destination matrix is symmetric and expanded, then expand data. + if (A->symmetric()) { + is_expand_symmetric = A->expanded(); + } + + // Skip the header comments + while (line.at(0) == '%') { + std::getline(file, line); + // std::cout << line << std::endl; + } + + // Read the first line with matrix sizes + index_type m, n, nnz; + ss << line; + ss >> n >> m >> nnz; + + // Make sure input data matches matrix A size + if ((A->getNumRows() != n) || (A->getNumColumns() != m)) { + Logger::error() << "Wrong matrix size: " << A->getNumRows() + << "x" << A->getNumColumns() + << ". Cannot update! \n "; + return; + } + + // Store COO data in the temporary workspace. Complexity O(NNZ) + loadToList(file, is_expand_symmetric, tmp); + + // Sort tmp. Complexity O(NNZ*log(NNZ)) + tmp.sort(); + + // Remove duplicates in `tmp` list. Complexity O(NNZ) + removeDuplicates(tmp); - std::getline(file, line); - while (line.at(0) == '%') { - std::getline(file, line); - // std::cout<> n >> m ; - - if (rhs == nullptr) { - // std::cout << "Allocating array of size " << n << "\n"; - rhs = new real_type[n]; - } - real_type a; - index_type i = 0; - while (file >> a) { - rhs[i] = a; - // std::cout << i << ": " << a << "\n"; - i++; + + + // Commented out; needed for debugging only. + // void print_list(std::list& l) + // { + // // Print out the list + // std::cout << "tmp list:\n"; + // for (MatrixElementTriplet& n : l) + // n.print(); + // std::cout << "\n"; + // } + + /** + * @brief Loads data from Matrix Market file to a std::list. + * + * @param[in] file - std::istream to Matrix Market file (data). + * @param[in] is_expand_symmetric - whether to expand symmetric matrix. + * @param[out] tmp - temporary list with matrix entries + * @return int - 0 if successful, error code otherwise. + * + * @pre `file` is a valid std::istream with Matrix Market data. + * @pre Input data is in valid Matrix Market format. + * @pre `tmp` is an empty list! + * @post `tmp` list is overwritten with matrix elements read from the input + * stream. + */ + int loadToList(std::istream& file, bool is_expand_symmetric, std::list& tmp) + { + index_type i, j; + real_type v; + + // If the `tmp` list is not empty, clear it. + if (tmp.size() != 0) { + tmp.clear(); + } + + while (file >> i >> j >> v) { + MatrixElementTriplet triplet(i - 1, j - 1, v); + tmp.push_back(std::move(triplet)); + if (is_expand_symmetric) { + if (i != j) { + MatrixElementTriplet triplet(j - 1, i - 1, v); + tmp.push_back(std::move(triplet)); + } + } + } + + return 0; } - } - int writeMatrixToFile(matrix::Sparse* A, std::ostream& file_out) - { - if (A == nullptr) { - Logger::error() << "Matrix pointer is NULL!\n"; - return -1; + /** + * @brief Removes duplicates from `tmp` list. + * + * @param[in,out] tmp - List with matrix entries. + * @return int - 0 if successful, error code otherwise. + * + * @pre `tmp` contains matrix elements. + * @post Duplicates in `tmp` are added in place to the first instance + * of that matrix element. + */ + int removeDuplicates(std::list& tmp) + { + std::list::iterator it = tmp.begin(); + while (it != tmp.end()) + { + std::list::iterator it_tmp = it; + it++; + if (*it == *it_tmp) { + *it += *it_tmp; + tmp.erase(it_tmp); + } + } + + return 0; } - if (A->symmetric() && !A->expanded()) { - file_out << "%%MatrixMarket matrix coordinate real symmetric\n"; - } else { - file_out << "%%MatrixMarket matrix coordinate real general\n"; + /** + * @brief Writes data from the std::list to COO matrix. + * + * @param[in] tmp - List with matrix entries + * @param[out] A - Output COO matrix + * @return int - 0 if successful, error code otherwise. + * + * @pre `tmp` contains matrix elements sorted in row-major order and + * without duplicates. + * @pre Number of `tmp` elements is not larger than number of nonzeros + * in `A`. + * @post Matrix data in `A` is updated with data from `tmp`. + */ + int copyListToCoo(const std::list& tmp, matrix::Coo* A) + { + index_type* coo_rows = A->getRowData(memory::HOST); + index_type* coo_cols = A->getColData(memory::HOST); + real_type* coo_vals = A->getValues( memory::HOST); + + index_type nnz = static_cast(tmp.size()); + if (A->getNnz() < nnz) { + Logger::error() << "Too many NNZs: " << nnz + << ". Cannot update! \n "; + return 1; + } + A->setNnz(nnz); + + index_type element_counter = 0; + std::list::const_iterator it = tmp.begin(); + while (it != tmp.end()) + { + coo_rows[element_counter] = it->getRowIdx(); + coo_cols[element_counter] = it->getColIdx(); + coo_vals[element_counter] = it->getValue(); + it++; + element_counter++; + } + + // We updated matrix values outside Matrix API. We need to note that. + A->setUpdated(memory::HOST); + + return 0; } - file_out << "% Generated by Re::Solve \n"; - file_out << A->getNumRows() << " " - << A->getNumColumns() << " " - << A->getNnz() << "\n"; - A->print(file_out); - return 0; - } - - int writeVectorToFile(vector_type* vec_x, std::ostream& file_out) - { - real_type* x_data = vec_x->getData(memory::HOST); - // std::ofstream file_out (filename, std::ofstream::out); - file_out << "%%MatrixMarket matrix array real general \n"; - file_out << "% ID: XXX \n"; - file_out << vec_x->getSize() << " " << 1 << "\n"; - for (int i = 0; i < vec_x->getSize(); ++i) { - file_out << std::setprecision(32) << std::scientific << x_data[i] << "\n"; + + + /** + * @brief Writes data from the std::list to CSR matrix. + * + * @param[in] tmp - List with matrix entries + * @param[out] A - Output CSR matrix + * @return int - 0 if successful, error code otherwise. + * + * @pre `tmp` contains matrix elements sorted in row-major order and + * without duplicates. + * @pre Number of `tmp` elements is not larger than number of nonzeros + * in `A`. + * @post Matrix data in `A` is updated with data from `tmp`. + */ + int copyListToCsr(const std::list& tmp, matrix::Csr* A) + { + index_type* csr_rows = A->getRowData(memory::HOST); + index_type* csr_cols = A->getColData(memory::HOST); + real_type* csr_vals = A->getValues( memory::HOST); + + // Set number of nonzeros + index_type nnz = static_cast(tmp.size()); + if (A->getNnz() < nnz) { + Logger::error() << "Too many NNZs: " << nnz + << ". Cannot update! \n "; + return 1; + } + A->setNnz(nnz); + + // Set all iterators + index_type column_index_counter = 0; + index_type row_pointer_counter = 0; + std::list::const_iterator it = tmp.begin(); + + // Set first row pointer to zero + csr_rows[0] = 0; + csr_cols[0] = it->getColIdx(); + csr_vals[0] = it->getValue(); + + for (index_type i = 1; i < nnz; ++i) { + std::list::const_iterator it_tmp = it; + it++; + if (it->getRowIdx() != it_tmp->getRowIdx()) { + row_pointer_counter++; + csr_rows[row_pointer_counter] = i; + } + column_index_counter++; + csr_cols[column_index_counter] = it->getColIdx(); + csr_vals[column_index_counter] = it->getValue(); + } + row_pointer_counter++; + csr_rows[row_pointer_counter] = nnz; + + // We updated matrix values outside Matrix API. We need to note that. + A->setUpdated(memory::HOST); + + return 0; } - // file_out.close(); - return 0; - } -}} // ReSolve::io + } // namespace io +} // namespace ReSolve diff --git a/resolve/matrix/io.hpp b/resolve/matrix/io.hpp index 41486ece..aac6ad3d 100644 --- a/resolve/matrix/io.hpp +++ b/resolve/matrix/io.hpp @@ -7,15 +7,18 @@ namespace ReSolve { namespace vector { namespace ReSolve { namespace matrix { class Sparse; class Coo; + class Csr; }} namespace ReSolve { namespace io { using vector_type = vector::Vector; - matrix::Coo* readMatrixFromFile(std::istream& file); - void readAndUpdateMatrix(std::istream& file, matrix::Coo* A); - real_type* readRhsFromFile(std::istream& file); - void readAndUpdateRhs(std::istream& file, real_type** rhs); + matrix::Coo* createCooFromFile(std::istream& file, bool is_expand_symmetric = true); + matrix::Csr* createCsrFromFile(std::istream& file, bool is_expand_symmetric = true); + void updateMatrixFromFile(std::istream& file, matrix::Coo* A); + void updateMatrixFromFile(std::istream& file, matrix::Csr* A); + real_type* createArrayFromFile(std::istream& file); + void updateArrayFromFile(std::istream& file, real_type** rhs); int writeMatrixToFile(matrix::Sparse* A, std::ostream& file_out); int writeVectorToFile(vector_type* vec_x, std::ostream &file_out); diff --git a/tests/functionality/testKLU.cpp b/tests/functionality/testKLU.cpp index c84f442f..67d6fe83 100644 --- a/tests/functionality/testKLU.cpp +++ b/tests/functionality/testKLU.cpp @@ -14,6 +14,7 @@ //author: KS //functionality test to check whether KLU works correctly. using namespace ReSolve::constants; +using namespace ReSolve::colors; int main(int argc, char *argv[]) { @@ -49,8 +50,7 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << matrixFileName1 << "\n"; return -1; } - ReSolve::matrix::Coo* A_coo = ReSolve::io::readMatrixFromFile(mat1); - ReSolve::matrix::Csr* A = new ReSolve::matrix::Csr(A_coo, ReSolve::memory::HOST); + ReSolve::matrix::Csr* A = ReSolve::io::createCsrFromFile(mat1); mat1.close(); // Read first rhs vector @@ -60,7 +60,7 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileName1 << "\n"; return -1; } - real_type* rhs = ReSolve::io::readRhsFromFile(rhs1_file); + real_type* rhs = ReSolve::io::createArrayFromFile(rhs1_file); real_type* x = new real_type[A->getNumRows()]; vector_type* vec_rhs = new vector_type(A->getNumRows()); vector_type* vec_x = new vector_type(A->getNumRows()); @@ -145,7 +145,7 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << matrixFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateMatrix(mat2, A_coo); + ReSolve::io::updateMatrixFromFile(mat2, A); mat2.close(); // Load the second rhs vector @@ -155,10 +155,9 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateRhs(rhs2_file, &rhs); + ReSolve::io::updateArrayFromFile(rhs2_file, &rhs); rhs2_file.close(); - A->updateFromCoo(A_coo, ReSolve::memory::HOST); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); // and solve it too @@ -197,14 +196,19 @@ int main(int argc, char *argv[]) std::cout<<"\t ||x-x_true||_2/||x_true||_2 : "< 1e-16 ) || (normRmatrix2/normB2 > 1e-16)) { std::cout << "Result inaccurate!\n"; error_sum++; } if (error_sum == 0) { - std::cout<<"Test KLU with KLU refactorization PASSED"<getNumRows()]; vector_type* vec_rhs = new vector_type(A->getNumRows()); vector_type* vec_x = new vector_type(A->getNumRows()); @@ -164,7 +164,7 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << matrixFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateMatrix(mat2, A_coo); + ReSolve::io::updateMatrixFromFile(mat2, A); mat2.close(); // Load the second rhs vector @@ -174,10 +174,9 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateRhs(rhs2_file, &rhs); + ReSolve::io::updateArrayFromFile(rhs2_file, &rhs); rhs2_file.close(); - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); status = GLU->refactorize(); @@ -216,14 +215,19 @@ int main(int argc, char *argv[]) std::cout<<"\t ||x-x_true||_2/||x_true||_2 : "< 1e-16 ) || (normRmatrix2/normB2 > 1e-16)) { std::cout << "Result inaccurate!\n"; error_sum++; } if (error_sum == 0) { - std::cout<<"Test KLU with cuSolverGLU refactorization PASSED"<getNumRows()]; vector_type* vec_rhs = new vector_type(A->getNumRows()); vector_type* vec_x = new vector_type(A->getNumRows()); @@ -163,7 +163,7 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << matrixFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateMatrix(mat2, A_coo); + ReSolve::io::updateMatrixFromFile(mat2, A); mat2.close(); // Load the second rhs vector @@ -173,10 +173,9 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateRhs(rhs2_file, &rhs); + ReSolve::io::updateArrayFromFile(rhs2_file, &rhs); rhs2_file.close(); - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); status = Rf->refactorize(); @@ -214,14 +213,19 @@ int main(int argc, char *argv[]) std::cout<<"\t ||x-x_true||_2/||x_true||_2 : "< 1e-16 ) || (normRmatrix2/normB2 > 1e-16)) { std::cout << "Result inaccurate!\n"; error_sum++; } if (error_sum == 0) { - std::cout<<"Test KLU with cuSolverRf refactorization PASSED"<getNumRows()]; vector_type* vec_rhs = new vector_type(A->getNumRows()); vector_type* vec_x = new vector_type(A->getNumRows()); @@ -197,7 +197,7 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << matrixFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateMatrix(mat2, A_coo); + ReSolve::io::updateMatrixFromFile(mat2, A); mat2.close(); // Load the second rhs vector @@ -207,10 +207,9 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateRhs(rhs2_file, &rhs); + ReSolve::io::updateArrayFromFile(rhs2_file, &rhs); rhs2_file.close(); - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); Rf->setNumericalProperties(1e-14, 1e-1); @@ -262,14 +261,19 @@ int main(int argc, char *argv[]) std::cout<<"\t IR starting res. norm : "<getInitResidualNorm()<<" "<getFinalResidualNorm()<<" "< 1e-12 ) || (normRmatrix2/normB2 > 1e-15)) { std::cout << "Result inaccurate!\n"; error_sum++; } if (error_sum == 0) { - std::cout<<"Test KLU with cuSolverRf refactorization + IR PASSED"< #include #include @@ -11,10 +17,9 @@ #include #include #include -//author: KS -//functionality test to check whether rocsolver_rf works correctly. using namespace ReSolve::constants; +using namespace ReSolve::colors; int main(int argc, char *argv[]) { @@ -55,8 +60,7 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << matrixFileName1 << "\n"; return -1; } - ReSolve::matrix::Coo* A_coo = ReSolve::io::readMatrixFromFile(mat1); - ReSolve::matrix::Csr* A = new ReSolve::matrix::Csr(A_coo, ReSolve::memory::HOST); + ReSolve::matrix::Csr* A = ReSolve::io::createCsrFromFile(mat1); mat1.close(); // Read first rhs vector @@ -66,7 +70,7 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileName1 << "\n"; return -1; } - real_type* rhs = ReSolve::io::readRhsFromFile(rhs1_file); + real_type* rhs = ReSolve::io::createArrayFromFile(rhs1_file); real_type* x = new real_type[A->getNumRows()]; vector_type* vec_rhs = new vector_type(A->getNumRows()); vector_type* vec_x = new vector_type(A->getNumRows()); @@ -168,7 +172,7 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << matrixFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateMatrix(mat2, A_coo); + ReSolve::io::updateMatrixFromFile(mat2, A); mat2.close(); // Load the second rhs vector @@ -178,10 +182,9 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateRhs(rhs2_file, &rhs); + ReSolve::io::updateArrayFromFile(rhs2_file, &rhs); rhs2_file.close(); - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); // this hangs up @@ -221,17 +224,19 @@ int main(int argc, char *argv[]) std::cout<<"\t ||x-x_true||_2/||x_true||_2 : "< 1e-16 ) || (normRmatrix2/normB2 > 1e-16)) { std::cout << "Result inaccurate!\n"; error_sum++; } if (error_sum == 0) { - std::cout<<"Test KLU with rocSolverRf refactorization PASSED"<getNumRows()]; vector_type* vec_rhs = new vector_type(A->getNumRows()); vector_type* vec_x = new vector_type(A->getNumRows()); @@ -117,7 +117,7 @@ int main(int argc, char *argv[]) vec_diff = new vector_type(A->getNumRows()); real_type* x_data = new real_type[A->getNumRows()]; - for (int i=0; igetNumRows(); ++i){ + for (int i = 0; i < A->getNumRows(); ++i) { x_data[i] = 1.0; } @@ -192,7 +192,7 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << matrixFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateMatrix(mat2, A_coo); + ReSolve::io::updateMatrixFromFile(mat2, A); mat2.close(); // Load the second rhs vector @@ -202,10 +202,9 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateRhs(rhs2_file, &rhs); + ReSolve::io::updateArrayFromFile(rhs2_file, &rhs); rhs2_file.close(); - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); status = Rf->refactorize(); @@ -256,15 +255,20 @@ int main(int argc, char *argv[]) std::cout<<"\t IR starting res. norm : "<getInitResidualNorm()<<" "<getFinalResidualNorm()<<" (tol 1e-14)"< 1e-12 ) || (normRmatrix2/normB2 > 1e-9)) { std::cout << "Result inaccurate!\n"; error_sum++; } if (error_sum == 0) { - std::cout<<"Test KLU with rocsolverrf refactorization + IR PASSED"<getNumRows()]; vector_type* vec_rhs = new vector_type(A->getNumRows()); vector_type* vec_x = new vector_type(A->getNumRows()); @@ -151,7 +151,7 @@ int main(int argc, char *argv[]) // Set the reference solution vector (all ones) on both, CPU and GPU real_type* x_data_ref = new real_type[A->getNumRows()]; - for (int i=0; igetNumRows(); ++i) { + for (int i = 0; i < A->getNumRows(); ++i) { x_data_ref[i] = 1.0; } vec_test->update(x_data_ref, ReSolve::memory::HOST, ReSolve::memory::HOST); @@ -201,7 +201,7 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << matrixFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateMatrix(mat2, A_coo); + ReSolve::io::updateMatrixFromFile(mat2, A); mat2.close(); // Load the second rhs vector @@ -210,11 +210,10 @@ int main(int argc, char *argv[]) std::cout << "Failed to open file " << rhsFileName2 << "\n"; return -1; } - ReSolve::io::readAndUpdateRhs(rhs2_file, &rhs); + ReSolve::io::updateArrayFromFile(rhs2_file, &rhs); rhs2_file.close(); // Update system values - A->updateFromCoo(A_coo, ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); // Refactorize and solve @@ -284,19 +283,23 @@ int main(int argc, char *argv[]) std::cout << "\t ||x-x_true||_2/||x_true||_2 : " << normDiffMatrix2/normXtrue << " (relative solution error)\n"; std::cout << "\t ||b-A*x_true||_2 (control) : " << exactSol_normRmatrix2 << " (residual norm with exact solution)\n\n"; + if (!std::isfinite(normRmatrix1/normB1) || !std::isfinite(normRmatrix2/normB2)) { + std::cout << "Result is not a finite number!\n"; + error_sum++; + } if ((normRmatrix1/normB1 > 1e-16 ) || (normRmatrix2/normB2 > 1e-16)) { std::cout << "Result inaccurate!\n"; error_sum++; } if (error_sum == 0) { - std::cout<<"Test KLU with cuSolverGLU refactorization PASSED"<updateFromCoo(A_coo, ReSolve::memory::DEVICE); vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); status = solver.refactorize(); @@ -317,18 +316,21 @@ int main(int argc, char *argv[]) std::cout << "\t IR starting res. norm : " << init_rnorm << "\n"; std::cout << "\t IR final res. norm : " << final_rnorm << " (tol " << std::setprecision(2) << tol << ")\n\n"; + if (!std::isfinite(normRmatrix1/normB1) || !std::isfinite(normRmatrix2/normB2)) { + std::cout << "Result is not a finite number!\n"; + error_sum++; + } if ((normRmatrix1/normB1 > 1e-12 ) || (normRmatrix2/normB2 > 1e-15)) { std::cout << "Result inaccurate!\n"; error_sum++; } if (error_sum == 0) { - std::cout<<"Test KLU with rocsolverrf refactorization + IR PASSED"< #include - +#include #include //author: RD @@ -18,16 +18,18 @@ */ int main() { + using namespace ReSolve::colors; std::string answer("0.99.1"); std::string versionstr; ReSolve::VersionGetVersionStr(versionstr); std::cout << "ReSolveVersionGetVersionStr Test: " << versionstr << std::endl << std::endl; if (versionstr != answer) { - std::cout << "ReSolve version set incorrectly. Test FAILED\n"; + std::cout << "ReSolve version set incorrectly. " + << "Test " << RED << "FAILED" << CLEAR << "\n"; return 1; } - std::cout << "ReSolve version test PASSED\n"; + std::cout << "ReSolve version test " << GREEN << "PASSED" << CLEAR << "\n"; return 0; } \ No newline at end of file diff --git a/tests/unit/matrix/MatrixIoTests.hpp b/tests/unit/matrix/MatrixIoTests.hpp index 34b82f6e..4be9c3b8 100644 --- a/tests/unit/matrix/MatrixIoTests.hpp +++ b/tests/unit/matrix/MatrixIoTests.hpp @@ -5,6 +5,7 @@ #include #include #include +#include namespace ReSolve { namespace tests { @@ -19,9 +20,9 @@ class MatrixIoTests : TestBase { TestStatus status; - // Read string into istream and status it to `readMatrixFromFile` function. + // Read string into istream and status it to `createCooFromFile` function. std::istringstream file(general_coo_matrix_file_); - ReSolve::matrix::Coo* A = ReSolve::io::readMatrixFromFile(file); + ReSolve::matrix::Coo* A = ReSolve::io::createCooFromFile(file); // Check if the matrix data was correctly loaded status = true; @@ -48,8 +49,9 @@ class MatrixIoTests : TestBase delete A; A = nullptr; - std::istringstream file2(symmetric_coo_matrix_file_); - A = ReSolve::io::readMatrixFromFile(file2); + bool is_expand_symmetric = false; + std::istringstream file2(symmetric_duplicates_coo_matrix_file_); + A = ReSolve::io::createCooFromFile(file2, is_expand_symmetric); nnz_answer = static_cast(symmetric_coo_matrix_vals_.size()); if (A->getNnz() != nnz_answer) { @@ -69,6 +71,68 @@ class MatrixIoTests : TestBase status *= verifyAnswer(*A, symmetric_coo_matrix_rows_, symmetric_coo_matrix_cols_, symmetric_coo_matrix_vals_); + delete A; + A = nullptr; + + is_expand_symmetric = true; + std::istringstream file3(symmetric_duplicates_coo_matrix_file_); + A = ReSolve::io::createCooFromFile(file3, is_expand_symmetric); + + nnz_answer = static_cast(symmetric_expanded_coo_matrix_vals_.size()); + if (A->getNnz() != nnz_answer) { + std::cout << "Incorrect NNZ read from the file ...\n"; + status = false; + } + + if (!A->symmetric()) { + std::cout << "Incorrect matrix type, matrix is symmetric ...\n"; + status = false; + } + + if (!A->expanded()) { + std::cout << "Incorrect matrix type, matrix is expanded ...\n"; + status = false; + } + + status *= verifyAnswer(*A, symmetric_expanded_coo_matrix_rows_, symmetric_expanded_coo_matrix_cols_, symmetric_expanded_coo_matrix_vals_); + + delete A; + A = nullptr; + + return status.report(__func__); + } + + + TestOutcome csrMatrixImport() + { + TestStatus status; + status = true; + + bool is_expand_symmetric = true; + std::istringstream file(symmetric_duplicates_coo_matrix_file_); + ReSolve::matrix::Csr* B = ReSolve::io::createCsrFromFile(file, is_expand_symmetric); + + index_type nnz_answer = static_cast(symmetric_expanded_csr_matrix_vals_.size()); + if (B->getNnz() != nnz_answer) { + std::cout << "Incorrect NNZ read from the file ...\n"; + status = false; + } + + if (!B->symmetric()) { + std::cout << "Incorrect matrix type, matrix is symmetric ...\n"; + status = false; + } + + if (!B->expanded()) { + std::cout << "Incorrect matrix type, matrix is expanded ...\n"; + status = false; + } + + status *= verifyAnswer(*B, symmetric_expanded_csr_matrix_rows_, symmetric_expanded_csr_matrix_cols_, symmetric_expanded_csr_matrix_vals_); + + delete B; + B = nullptr; + return status.report(__func__); } @@ -78,7 +142,7 @@ class MatrixIoTests : TestBase TestStatus status; status = true; - // Read string into istream and status it to `readMatrixFromFile` function. + // Read string into istream and status it to `createCooFromFile` function. std::ostringstream buffer; // Deep copy constant test vectors with matrix data to nonconstant ones @@ -105,32 +169,62 @@ class MatrixIoTests : TestBase // Write the matrix to an ostream ReSolve::io::writeMatrixToFile(&A, buffer); status *= (buffer.str() == resolve_general_coo_matrix_file_); - buffer.str(""); + + return status.report(__func__); + } + + TestOutcome csrMatrixExport() + { + TestStatus status; + status = true; + + // Read string into istream and status it to `createCooFromFile` function. + std::ostringstream buffer; + + // Deep copy constant test vectors with matrix data to nonconstant ones + std::vector rows = general_csr_matrix_rows_; + std::vector cols = general_csr_matrix_cols_; + std::vector vals = general_csr_matrix_vals_; + + // Get number of matrix rows + const index_type N = static_cast(rows.size()) - 1; + + // Get number of matrix columns + const index_type M = 1 + *(std::max_element(cols.begin(), cols.end())); + + // Get number of nonzeros + const index_type NNZ = static_cast(vals.size()); // Create the test CSR matrix - ReSolve::matrix::Csr B(&A, ReSolve::memory::HOST); + ReSolve::matrix::Csr A(N, M, NNZ, false, false); + A.setMatrixData(&rows[0], + &cols[0], + &vals[0], + memory::HOST); // Write the matrix to an ostream - ReSolve::io::writeMatrixToFile(&B, buffer); + ReSolve::io::writeMatrixToFile(&A, buffer); status *= (buffer.str() == resolve_row_sorted_general_coo_matrix_file_); return status.report(__func__); } - TestOutcome cooMatrixReadAndUpdate() { TestStatus status; + bool is_symmetric = true; + bool is_expanded = false; + // Create a 5x5 COO matrix with 10 nonzeros - ReSolve::matrix::Coo A(5, 5, 10); + ReSolve::matrix::Coo A(5, 5, 10, is_symmetric, is_expanded); A.allocateMatrixData(memory::HOST); - // Read string into istream and status it to `readMatrixFromFile` function. + // Read string into istream and status it to `createCooFromFile` function. std::istringstream file2(symmetric_coo_matrix_file_); // Update matrix A with data from the matrix market file - ReSolve::io::readAndUpdateMatrix(file2, &A); + ReSolve::io::updateMatrixFromFile(file2, &A); // Check if the matrix data was correctly loaded status = true; @@ -143,19 +237,69 @@ class MatrixIoTests : TestBase status *= verifyAnswer(A, symmetric_coo_matrix_rows_, symmetric_coo_matrix_cols_, symmetric_coo_matrix_vals_); - // Read string into istream and status it to `readMatrixFromFile` function. - std::istringstream file(general_coo_matrix_file_); + // Read next matrix market file into istream. This matrix has duplicates + // that need to be merged and number of nonzeros needs to be recalculated + // accordingly. + std::istringstream file(symmetric_duplicates_coo_matrix_file_); // Update matrix A with data from the matrix market file - ReSolve::io::readAndUpdateMatrix(file, &A); + ReSolve::io::updateMatrixFromFile(file, &A); - nnz_answer = static_cast(general_coo_matrix_vals_.size()); if (A.getNnz() != nnz_answer) { std::cout << "Incorrect NNZ read from the file ...\n"; status = false; } - status *= verifyAnswer(A, general_coo_matrix_rows_, general_coo_matrix_cols_, general_coo_matrix_vals_); + status *= verifyAnswer(A, symmetric_coo_matrix_rows_, symmetric_coo_matrix_cols_, symmetric_coo_matrix_vals_); + + // Create a 5x5 COO matrix with 13 nonzeros + is_expanded = true; + ReSolve::matrix::Coo B(5, 5, 13, is_symmetric, is_expanded); + B.allocateMatrixData(memory::HOST); + + // Read in symmetric matrix data + std::istringstream file3(symmetric_duplicates_coo_matrix_file_); + + // Update matrix B with data from the matrix market file + ReSolve::io::updateMatrixFromFile(file3, &B); + + nnz_answer = static_cast(symmetric_expanded_coo_matrix_vals_.size()); + if (B.getNnz() != nnz_answer) { + std::cout << "Incorrect NNZ read from the file ...\n"; + status = false; + } + + status *= verifyAnswer(B, symmetric_expanded_coo_matrix_rows_, symmetric_expanded_coo_matrix_cols_, symmetric_expanded_coo_matrix_vals_); + + return status.report(__func__); + } + + + TestOutcome csrMatrixReadAndUpdate() + { + TestStatus status; + status = true; + + bool is_symmetric = true; + bool is_expanded = true; + + ReSolve::matrix::Csr A(5, 5, 13, is_symmetric, is_expanded); + A.allocateMatrixData(memory::HOST); + + // Read in symmetric matrix data + std::istringstream file(symmetric_duplicates_coo_matrix_file_); + + // Update matrix B with data from the matrix market file + ReSolve::io::updateMatrixFromFile(file, &A); + + index_type nnz_answer = static_cast(symmetric_expanded_csr_matrix_vals_.size()); + if (A.getNnz() != nnz_answer) { + std::cout << "Incorrect NNZ read from the file ...\n"; + std::cout << A.getNnz() << " ?= " << nnz_answer << "\n"; + status = false; + } + + status *= verifyAnswer(A, symmetric_expanded_csr_matrix_rows_, symmetric_expanded_csr_matrix_cols_, symmetric_expanded_csr_matrix_vals_); return status.report(__func__); } @@ -164,11 +308,11 @@ class MatrixIoTests : TestBase { TestStatus status; - // Read string into istream and status it to `readMatrixFromFile` function. + // Read string into istream and status it to `createCooFromFile` function. std::istringstream file(general_vector_file_); // Create rhs vector and load its data from the input file - real_type* rhs = ReSolve::io::readRhsFromFile(file); + real_type* rhs = ReSolve::io::createArrayFromFile(file); // Check if the matrix data was correctly loaded status = true; @@ -190,14 +334,14 @@ class MatrixIoTests : TestBase { TestStatus status; - // Read string into istream and status it to `readMatrixFromFile` function. + // Read string into istream and status it to `createCooFromFile` function. std::istringstream file(general_vector_file_); - // For now let's test only the case when `readAndUpdateRhs` does not allocate rhs + // For now let's test only the case when `updateArrayFromFile` does not allocate rhs real_type* rhs = new real_type[5]; //nullptr; // Update matrix A with data from the matrix market file - ReSolve::io::readAndUpdateRhs(file, &rhs); + ReSolve::io::updateArrayFromFile(file, &rhs); // Check if the matrix data was correctly loaded status = true; @@ -233,6 +377,29 @@ class MatrixIoTests : TestBase return true; } + bool verifyAnswer(/* const */ ReSolve::matrix::Csr& answer, + const std::vector& row_data, + const std::vector& col_data, + const std::vector& val_data) + { + for (size_t i = 0; i < val_data.size(); ++i) { + if ((answer.getColData(memory::HOST)[i] != col_data[i]) || + (!isEqual(answer.getValues(memory::HOST)[i], val_data[i]))) { + std::cout << "Incorrect matrix value at storage element " << i << ".\n"; + return false; + } + } + + for (size_t i = 0; i < row_data.size(); ++i) { + if(answer.getRowData(memory::HOST)[i] != row_data[i]) { + std::cout << "Incorrect row pointer value at storage element " << i << ".\n"; + return false; + } + } + + return true; + } + private: // // Test examples @@ -241,7 +408,8 @@ class MatrixIoTests : TestBase /// String pretending to be matrix market file /// Same stored in file `matrix_general_coo_ordered.mtx` const std::string general_coo_matrix_file_ = -R"(% This ASCII file represents a sparse MxN matrix with L +R"(%%MatrixMarket matrix coordinate real general +% This ASCII file represents a sparse MxN matrix with L % nonzeros in the following Matrix Market format: % % +----------------------------------------------+ @@ -277,14 +445,14 @@ R"(% This ASCII file represents a sparse MxN matrix with L R"(%%MatrixMarket matrix coordinate real general % Generated by Re::Solve 5 5 8 -0 0 1.000000000000000e+00 -1 1 1.050000000000000e+01 -2 2 1.500000000000000e-02 -0 3 6.000000000000000e+00 -3 1 2.505000000000000e+02 -3 3 -2.800000000000000e+02 -3 4 3.332000000000000e+01 -4 4 1.200000000000000e+01 +1 1 1.000000000000000e+00 +1 4 6.000000000000000e+00 +2 2 1.050000000000000e+01 +3 3 1.500000000000000e-02 +4 2 2.505000000000000e+02 +4 4 -2.800000000000000e+02 +4 5 3.332000000000000e+01 +5 5 1.200000000000000e+01 )"; /// String pretending to be matrix market file. @@ -293,30 +461,68 @@ R"(%%MatrixMarket matrix coordinate real general R"(%%MatrixMarket matrix coordinate real general % Generated by Re::Solve 5 5 8 -0 0 1.000000000000000e+00 -0 3 6.000000000000000e+00 -1 1 1.050000000000000e+01 -2 2 1.500000000000000e-02 -3 1 2.505000000000000e+02 -3 3 -2.800000000000000e+02 -3 4 3.332000000000000e+01 -4 4 1.200000000000000e+01 +1 1 1.000000000000000e+00 +1 4 6.000000000000000e+00 +2 2 1.050000000000000e+01 +3 3 1.500000000000000e-02 +4 2 2.505000000000000e+02 +4 4 -2.800000000000000e+02 +4 5 3.332000000000000e+01 +5 5 1.200000000000000e+01 )"; /// Matching COO matrix data as it is supposed to be read from the file - const std::vector general_coo_matrix_rows_ = {0,1,2,0,3,3,3,4}; - const std::vector general_coo_matrix_cols_ = {0,1,2,3,1,3,4,4}; - const std::vector general_coo_matrix_vals_ = { 1.000e+00, - 1.050e+01, - 1.500e-02, - 6.000e+00, - 2.505e+02, - -2.800e+02, - 3.332e+01, - 1.200e+01 }; + const std::vector general_coo_matrix_rows_ = {0,0,1,2,3,3,3,4}; + const std::vector general_coo_matrix_cols_ = {0,3,1,2,1,3,4,4}; + const std::vector general_coo_matrix_vals_ = { 1.000e+00, + 6.000e+00, + 1.050e+01, + 1.500e-02, + 2.505e+02, + -2.800e+02, + 3.332e+01, + 1.200e+01 }; + + /// Matching CSR matrix data as it is supposed to be read from the file + const std::vector general_csr_matrix_rows_ = {0,2,3,4,7,8}; + const std::vector general_csr_matrix_cols_ = {0,3,1,2,1,3,4,4}; + const std::vector general_csr_matrix_vals_ = { 1.000e+00, + 6.000e+00, + 1.050e+01, + 1.500e-02, + 2.505e+02, + -2.800e+02, + 3.332e+01, + 1.200e+01 }; + // + // [11 15] + // [ 22 23 24 ] + // A = [ 33 35] + // [ 44 ] + // [ 55] + // const std::string symmetric_coo_matrix_file_ = R"(%%MatrixMarket matrix coordinate real symmetric +% This ASCII file represents a sparse MxN matrix with L +% nonzeros in the following Matrix Market format: +% +% +------------------------------------------------+ +% |%%MatrixMarket matrix coordinate real symmetric | <--- header line +% |% | <--+ +% |% comments | |-- 0 or more comment lines +% |% | <--+ +% | M N L | <--- rows, columns, entries +% | I1 J1 A(I1, J1) | <--+ +% | I2 J2 A(I2, J2) | | +% | I3 J3 A(I3, J3) | |-- L lines +% | . . . | | +% | IL JL A(IL, JL) | <--+ +% +------------------------------------------------+ +% +% Indices are 1-based, i.e. A(1,1) is the first element. +% +%================================================================================= % 5 5 9 1 1 11.0 @@ -331,10 +537,38 @@ R"(%%MatrixMarket matrix coordinate real symmetric )"; + // + // [11 15] + // [ 22 23 24 ] + // A = [ 33 35] + // [ 44 ] + // [ 55] + // + // A(1,1), A(5,5) and A(2,4) are stored in duplicate entries. + const std::string symmetric_duplicates_coo_matrix_file_ = +R"(%%MatrixMarket matrix coordinate real symmetric +% + 5 5 13 + 5 5 50.0 + 1 1 10.0 + 1 5 15.0 + 2 2 22.0 + 2 3 23.0 + 2 4 20.0 + 3 3 33.0 + 3 5 35.0 + 4 4 44.0 + 5 5 5.0 + 2 4 2.0 + 2 4 2.0 + 1 1 1.0 + )"; + + /// Matching COO matrix data as it is supposed to be read from the file const std::vector symmetric_coo_matrix_rows_ = {0,0,1,1,1,2,2,3,4}; const std::vector symmetric_coo_matrix_cols_ = {0,4,1,2,3,2,4,3,4}; - const std::vector symmetric_coo_matrix_vals_ = { 11.0, + const std::vector symmetric_coo_matrix_vals_ = {11.0, 15.0, 22.0, 23.0, @@ -342,8 +576,59 @@ R"(%%MatrixMarket matrix coordinate real symmetric 33.0, 35.0, 44.0, - 55.0 }; - + 55.0}; + + // Matching expanded COO matrix data as it is supposed to be created + // + // [11 15] + // [ 22 23 24 ] + // A = [ 23 33 35] + // [ 24 44 ] + // [15 35 55] + // + // Symmetric matrix stored in COO general format + // + const std::vector symmetric_expanded_coo_matrix_rows_ = {0,0,1,1,1,2,2,2,3,3,4,4,4}; + const std::vector symmetric_expanded_coo_matrix_cols_ = {0,4,1,2,3,1,2,4,1,3,0,2,4}; + const std::vector symmetric_expanded_coo_matrix_vals_ = { 11.0, + 15.0, + 22.0, + 23.0, + 24.0, + 23.0, + 33.0, + 35.0, + 24.0, + 44.0, + 15.0, + 35.0, + 55.0 }; + + // Matching expanded CSR matrix data as it is supposed to be converted + // + // [11 15] + // [ 22 23 24 ] + // A = [ 23 33 35] + // [ 24 44 ] + // [15 35 55] + // + // Symmetric matrix in CSR general format + // + const std::vector symmetric_expanded_csr_matrix_rows_ = {0,2,5,8,10,13}; + const std::vector symmetric_expanded_csr_matrix_cols_ = {0,4,1,2,3,1,2,4,1,3,0,2,4}; + const std::vector symmetric_expanded_csr_matrix_vals_ = { 11.0, + 15.0, + 22.0, + 23.0, + 24.0, + 23.0, + 33.0, + 35.0, + 24.0, + 44.0, + 15.0, + 35.0, + 55.0 }; const std::string general_vector_file_ = R"(% This ASCII file represents a sparse MxN matrix with L diff --git a/tests/unit/matrix/runMatrixIoTests.cpp b/tests/unit/matrix/runMatrixIoTests.cpp index 2c14b4f8..485deff0 100644 --- a/tests/unit/matrix/runMatrixIoTests.cpp +++ b/tests/unit/matrix/runMatrixIoTests.cpp @@ -12,8 +12,11 @@ int main(int, char**) ReSolve::tests::TestingResults result; result += test.cooMatrixImport(); + result += test.csrMatrixImport(); result += test.cooMatrixExport(); + result += test.csrMatrixExport(); result += test.cooMatrixReadAndUpdate(); + result += test.csrMatrixReadAndUpdate(); result += test.rhsVectorReadFromFile(); result += test.rhsVectorReadAndUpdate();