Skip to content

Commit 0ac0955

Browse files
peleshpelesh
committed
Fix bugs, add functionality tests.
1 parent 03d3a66 commit 0ac0955

File tree

4 files changed

+378
-5
lines changed

4 files changed

+378
-5
lines changed

‎examples/r_randGMRES.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -89,8 +89,8 @@ int main(int argc, char *argv[])
8989
matrix_handler->setValuesChanged(true, "hip");
9090

9191
Rf->setup(A);
92-
FGMRES->setRestart(800);
93-
FGMRES->setMaxit(800);
92+
FGMRES->setRestart(150);
93+
FGMRES->setMaxit(2500);
9494
FGMRES->setTol(1e-12);
9595
FGMRES->setup(A);
9696
GS->setup(FGMRES->getKrand(), FGMRES->getRestart());

‎resolve/LinSolverIterativeRandFGMRES.cpp

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -313,7 +313,7 @@ namespace ReSolve
313313
// residual norm estimate
314314
rnorm = fabs(h_rs_[i + 1]);
315315

316-
io::Logger::misc() << "it: "<<it<< "norm of residual "
316+
io::Logger::misc() << "it: "<<it<< " --> norm of the residual "
317317
<< std::scientific << std::setprecision(16)
318318
<< rnorm << "\n";
319319
// check convergence
@@ -371,6 +371,9 @@ namespace ReSolve
371371

372372
rand_manager_->reset();
373373

374+
if (rand_method_ == cs) {
375+
mem_.setZeroArrayOnDevice(d_S_->getData(memory::DEVICE), d_S_->getSize() * d_S_->getNumVectors());
376+
}
374377
vec_v->setData( d_V_->getVectorData(0, memory::DEVICE), memory::DEVICE);
375378
vec_s->setData( d_S_->getVectorData(0, memory::DEVICE), memory::DEVICE);
376379
rand_manager_->Theta(vec_v, vec_s);
@@ -389,7 +392,7 @@ namespace ReSolve
389392
// rnorm = ||V_0||
390393
rnorm = sqrt(rnorm);
391394

392-
io::Logger::misc() << "End of cycle, COMPUTER norm of residual "
395+
io::Logger::misc() << "End of cycle, COMPUTED norm of residual "
393396
<< std::scientific << std::setprecision(16)
394397
<< rnorm << "\n";
395398

@@ -496,4 +499,4 @@ namespace ReSolve
496499
{
497500
return fgmres_iters_;
498501
}
499-
}//namespace
502+
} // namespace ReSolve

‎tests/functionality/testRandGMRES_Cuda.cpp

Lines changed: 185 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,185 @@
1+
#include <string>
2+
#include <iostream>
3+
#include <iomanip>
4+
5+
#include <resolve/matrix/Coo.hpp>
6+
#include <resolve/matrix/Csr.hpp>
7+
#include <resolve/matrix/Csc.hpp>
8+
#include <resolve/vector/Vector.hpp>
9+
#include <resolve/matrix/io.hpp>
10+
#include <resolve/matrix/MatrixHandler.hpp>
11+
#include <resolve/vector/VectorHandler.hpp>
12+
#include <resolve/LinSolverDirectKLU.hpp>
13+
#include <resolve/LinSolverDirectCuSparseILU0.hpp>
14+
#include <resolve/LinSolverIterativeRandFGMRES.hpp>
15+
#include <resolve/workspace/LinAlgWorkspace.hpp>
16+
17+
using namespace ReSolve::constants;
18+
19+
int main(int argc, char *argv[])
20+
{
21+
// Use the same data types as those you specified in ReSolve build.
22+
using real_type = ReSolve::real_type;
23+
using vector_type = ReSolve::vector::Vector;
24+
25+
26+
//we want error sum to be 0 at the end
27+
//that means PASS.
28+
//otheriwse it is a FAIL.
29+
int error_sum = 0;
30+
int status;
31+
const std::string data_path = (argc == 2) ? argv[1] : "./";
32+
33+
34+
std::string matrixFileName = data_path + "data/SiO2.mtx";
35+
std::string rhsFileName = data_path + "data/SiO2_rhs.mtx";
36+
37+
38+
ReSolve::matrix::Coo* A_coo;
39+
ReSolve::matrix::Csr* A;
40+
41+
ReSolve::LinAlgWorkspaceCUDA* workspace_CUDA = new ReSolve::LinAlgWorkspaceCUDA();
42+
workspace_CUDA->initializeHandles();
43+
ReSolve::MatrixHandler* matrix_handler = new ReSolve::MatrixHandler(workspace_CUDA);
44+
ReSolve::VectorHandler* vector_handler = new ReSolve::VectorHandler(workspace_CUDA);
45+
real_type* rhs = nullptr;
46+
real_type* x = nullptr;
47+
48+
vector_type* vec_rhs;
49+
vector_type* vec_x;
50+
51+
ReSolve::GramSchmidt* GS = new ReSolve::GramSchmidt(vector_handler, ReSolve::GramSchmidt::cgs2);
52+
53+
ReSolve::LinSolverDirectCuSparseILU0* Rf = new ReSolve::LinSolverDirectCuSparseILU0(workspace_CUDA);
54+
ReSolve::LinSolverIterativeRandFGMRES* FGMRES = new ReSolve::LinSolverIterativeRandFGMRES(matrix_handler, vector_handler,ReSolve::LinSolverIterativeRandFGMRES::cs, GS, "cuda");
55+
56+
std::ifstream mat_file(matrixFileName);
57+
if(!mat_file.is_open())
58+
{
59+
std::cout << "Failed to open file " << matrixFileName << "\n";
60+
return -1;
61+
}
62+
63+
std::ifstream rhs_file(rhsFileName);
64+
if(!rhs_file.is_open())
65+
{
66+
std::cout << "Failed to open file " << rhsFileName << "\n";
67+
return -1;
68+
}
69+
A_coo = ReSolve::io::readMatrixFromFile(mat_file);
70+
A = new ReSolve::matrix::Csr(A_coo->getNumRows(),
71+
A_coo->getNumColumns(),
72+
A_coo->getNnz(),
73+
A_coo->symmetric(),
74+
A_coo->expanded());
75+
76+
rhs = ReSolve::io::readRhsFromFile(rhs_file);
77+
78+
x = new real_type[A->getNumRows()];
79+
vec_rhs = new vector_type(A->getNumRows());
80+
vec_x = new vector_type(A->getNumRows());
81+
vec_x->allocate(ReSolve::memory::HOST);
82+
83+
//iinit guess is 0
84+
vec_x->allocate(ReSolve::memory::DEVICE);
85+
vec_x->setToZero(ReSolve::memory::DEVICE);
86+
87+
mat_file.close();
88+
rhs_file.close();
89+
90+
matrix_handler->coo2csr(A_coo,A, "cuda");
91+
vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
92+
93+
real_type norm_b;
94+
matrix_handler->setValuesChanged(true, "cuda");
95+
96+
status = Rf->setup(A);
97+
error_sum += status;
98+
99+
FGMRES->setRestart(150);
100+
FGMRES->setMaxit(2500);
101+
FGMRES->setTol(1e-12);
102+
FGMRES->setup(A);
103+
104+
status = GS->setup(FGMRES->getKrand(), FGMRES->getRestart());
105+
error_sum += status;
106+
107+
//matrix_handler->setValuesChanged(true, "cuda");
108+
status = FGMRES->resetMatrix(A);
109+
error_sum += status;
110+
111+
status = FGMRES->setupPreconditioner("LU", Rf);
112+
error_sum += status;
113+
114+
FGMRES->setFlexible(1);
115+
116+
vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
117+
FGMRES->solve(vec_rhs, vec_x);
118+
119+
norm_b = vector_handler->dot(vec_rhs, vec_rhs, "cuda");
120+
norm_b = std::sqrt(norm_b);
121+
real_type final_norm_first = FGMRES->getFinalResidualNorm();
122+
std::cout << "Randomized FGMRES results (first run): \n"
123+
<< "\t Sketching method: : CountSketch\n"
124+
<< "\t Initial residual norm: ||b-Ax_0||_2 : "
125+
<< std::scientific << std::setprecision(16)
126+
<< FGMRES->getInitResidualNorm()<<" \n"
127+
<< "\t Initial relative residual norm: ||b-Ax_0||_2/||b||_2 : "
128+
<< FGMRES->getInitResidualNorm()/norm_b<<" \n"
129+
<< "\t Final residual norm: ||b-Ax||_2 : "
130+
<< FGMRES->getFinalResidualNorm() <<" \n"
131+
<< "\t Final relative residual norm: ||b-Ax||_2/||b||_2 : "
132+
<< FGMRES->getFinalResidualNorm()/norm_b <<" \n"
133+
<< "\t Number of iterations : " << FGMRES->getNumIter() << "\n";
134+
135+
delete FGMRES;
136+
delete GS;
137+
GS = new ReSolve::GramSchmidt(vector_handler, ReSolve::GramSchmidt::cgs2);
138+
FGMRES = new ReSolve::LinSolverIterativeRandFGMRES(matrix_handler, vector_handler,ReSolve::LinSolverIterativeRandFGMRES::fwht, GS, "cuda");
139+
140+
141+
FGMRES->setRestart(150);
142+
FGMRES->setMaxit(2500);
143+
FGMRES->setTol(1e-12);
144+
FGMRES->setup(A);
145+
146+
status = GS->setup(FGMRES->getKrand(), FGMRES->getRestart());
147+
error_sum += status;
148+
149+
status = FGMRES->setupPreconditioner("LU", Rf);
150+
error_sum += status;
151+
152+
vec_x->setToZero(ReSolve::memory::DEVICE);
153+
FGMRES->solve(vec_rhs, vec_x);
154+
155+
156+
std::cout << "Randomized FGMRES results (second run): \n"
157+
<< "\t Sketching method: : FWHT\n"
158+
<< "\t Initial residual norm: ||b-Ax_0||_2 : "
159+
<< std::scientific << std::setprecision(16)
160+
<< FGMRES->getInitResidualNorm()<<" \n"
161+
<< "\t Initial relative residual norm: ||b-Ax_0||_2/||b||_2 : "
162+
<< FGMRES->getInitResidualNorm()/norm_b<<" \n"
163+
<< "\t Final residual norm: ||b-Ax||_2 : "
164+
<< FGMRES->getFinalResidualNorm() <<" \n"
165+
<< "\t Final relative residual norm: ||b-Ax||_2/||b||_2 : "
166+
<< FGMRES->getFinalResidualNorm()/norm_b <<" \n"
167+
<< "\t Number of iterations : " << FGMRES->getNumIter() << "\n";
168+
169+
if ((error_sum == 0) && (final_norm_first/norm_b < 1e-11) && (FGMRES->getFinalResidualNorm()/norm_b < 1e-11 )) {
170+
std::cout<<"Test 5 (randomized GMRES) PASSED"<<std::endl<<std::endl;;
171+
} else {
172+
std::cout<<"Test 5 (randomized GMRES) FAILED, error sum: "<<error_sum<<std::endl<<std::endl;;
173+
}
174+
delete A;
175+
delete A_coo;
176+
delete Rf;
177+
delete [] x;
178+
delete [] rhs;
179+
delete vec_x;
180+
delete workspace_CUDA;
181+
delete matrix_handler;
182+
delete vector_handler;
183+
184+
return error_sum;
185+
}

0 commit comments

Comments
 (0)