@@ -77,7 +77,7 @@ namespace ReSolve
7777
7878 real_type expected_x[3 ] = {1.0 , -0.5 , 0.25 };
7979
80- real_type tol = 1e-8 ;
80+ real_type tol = 1e-12 ;
8181 for (index_type i = 0 ; i < n; ++i)
8282 {
8383 if (fabs (x->getData (memory::HOST)[i] - expected_x[i]) > tol)
@@ -104,11 +104,17 @@ namespace ReSolve
104104
105105 cholmod_sparse* L = randomSparseLowerTriangular ((size_t ) n);
106106 cholmod_sparse* L_tr = cholmod_transpose (L, 1 , &Common);
107- cholmod_sparse* A_chol = cholmod_ssmult (L, L_tr, 0 , 1 , 0 , &Common);
107+ cholmod_sparse* L_times_L_tr = cholmod_ssmult (L, L_tr, 0 , 1 , 0 , &Common);
108108
109- matrix::Csr* A = new matrix::Csr ((index_type) A_chol->nrow , (index_type) A_chol->ncol , (index_type) A_chol->nzmax );
109+ matrix::Csr* A = new matrix::Csr ((index_type) L_times_L_tr->nrow ,
110+ (index_type) L_times_L_tr->ncol ,
111+ (index_type) L_times_L_tr->nzmax );
110112 A->copyDataFrom (
111- static_cast <int *>(A_chol->p ), static_cast <int *>(A_chol->i ), static_cast <double *>(A_chol->x ), memory::HOST, memspace_);
113+ static_cast <int *>(L_times_L_tr->p ),
114+ static_cast <int *>(L_times_L_tr->i ),
115+ static_cast <double *>(L_times_L_tr->x ),
116+ memory::HOST,
117+ memspace_);
112118
113119 ReSolve::hykkt::CholeskySolver solver (memspace_);
114120
@@ -133,13 +139,11 @@ namespace ReSolve
133139
134140 if (memspace_ == memory::DEVICE)
135141 {
136- // x_expected->syncData(memory::HOST);
137142 x->syncData (memory::HOST);
138143 }
139144
140145 // Verify result
141- // TODO only works with 1e-4 tolerance
142- real_type tol = 1e-4 ;
146+ real_type tol = 1e-12 ;
143147 for (index_type j = 0 ; j < n; ++j)
144148 {
145149 if (fabs (x->getData (memory::HOST)[j] - x_expected->getData (memory::HOST)[j]) > tol)
@@ -155,7 +159,7 @@ namespace ReSolve
155159
156160 cholmod_free_sparse (&L, &Common);
157161 cholmod_free_sparse (&L_tr, &Common);
158- cholmod_free_sparse (&A_chol , &Common);
162+ cholmod_free_sparse (&L_times_L_tr , &Common);
159163
160164 delete A;
161165 delete x_expected;
@@ -212,8 +216,7 @@ namespace ReSolve
212216 }
213217
214218 // Verify result
215- // TODO only works with 1e-4 tolerance
216- real_type tol = 1e-4 ;
219+ real_type tol = 1e-12 ;
217220 for (index_type j = 0 ; j < n; ++j)
218221 {
219222 if (fabs (x->getData (memory::HOST)[j] - x_expected->getData (memory::HOST)[j]) > tol)
@@ -273,7 +276,13 @@ namespace ReSolve
273276 if (i == j || static_cast <double >(rand ()) / RAND_MAX < density)
274277 {
275278 L_i.push_back ((int ) j);
276- L_x.push_back (2.0 * static_cast <double >(rand ()) / RAND_MAX - 1.0 );
279+ // force diagonal entry to be non-zero
280+ double value = 2.0 ;
281+ if (i != j)
282+ {
283+ value = 2.0 * static_cast <double >(rand ()) / RAND_MAX - 1.0 ;
284+ }
285+ L_x.push_back (value);
277286 L_p[i + 1 ]++;
278287 nnz++;
279288 }
0 commit comments