1 // This is core/vnl/vnl_sparse_matrix_linear_system.cxx
2
3 #include <cassert>
4 #include "vnl_sparse_matrix_linear_system.h"
5 #include <vnl/vnl_copy.h>
6
7 template <>
get_rhs(vnl_vector<double> & b) const8 void vnl_sparse_matrix_linear_system<double>::get_rhs(vnl_vector<double>& b) const
9 {
10 b = b_;
11 }
12
13 template <>
transpose_multiply(vnl_vector<double> const & b,vnl_vector<double> & x) const14 void vnl_sparse_matrix_linear_system<double>::transpose_multiply(vnl_vector<double> const& b, vnl_vector<double> & x) const
15 {
16 A_.pre_mult(b,x);
17 }
18
19 template <>
get_rhs(vnl_vector<double> & b) const20 void vnl_sparse_matrix_linear_system<float>::get_rhs(vnl_vector<double>& b) const
21 {
22 vnl_copy(b_, b);
23 }
24
25 template <>
transpose_multiply(vnl_vector<double> const & b,vnl_vector<double> & x) const26 void vnl_sparse_matrix_linear_system<float>::transpose_multiply(vnl_vector<double> const& b, vnl_vector<double> & x) const
27 {
28 static vnl_vector<float> x_float;
29 static vnl_vector<float> b_float;
30
31 if (x_float.size() != x.size()) x_float = vnl_vector<float> (x.size());
32 if (b_float.size() != b.size()) b_float = vnl_vector<float> (b.size());
33
34 vnl_copy(b, b_float);
35 A_.pre_mult(b_float,x_float);
36 vnl_copy(x_float, x);
37 }
38
39 template <>
multiply(vnl_vector<double> const & x,vnl_vector<double> & b) const40 void vnl_sparse_matrix_linear_system<double>::multiply(vnl_vector<double> const& x, vnl_vector<double> & b) const
41 {
42 A_.mult(x,b);
43 }
44
45
46 template <>
multiply(vnl_vector<double> const & x,vnl_vector<double> & b) const47 void vnl_sparse_matrix_linear_system<float>::multiply(vnl_vector<double> const& x, vnl_vector<double> & b) const
48 {
49 static vnl_vector<float> x_float;
50 static vnl_vector<float> b_float;
51
52 if (x_float.size() != x.size()) x_float = vnl_vector<float> (x.size());
53 if (b_float.size() != b.size()) b_float = vnl_vector<float> (b.size());
54
55 vnl_copy(x, x_float);
56 A_.mult(x_float,b_float);
57 vnl_copy(b_float, b);
58 }
59
60
61 template<class T>
apply_preconditioner(vnl_vector<double> const & x,vnl_vector<double> & px) const62 void vnl_sparse_matrix_linear_system<T>::apply_preconditioner(vnl_vector<double> const& x, vnl_vector<double> & px) const
63 {
64 assert(x.size() == px.size());
65
66 if (jacobi_precond_.size() == 0) {
67 vnl_vector<T> tmp(get_number_of_unknowns());
68 A_.diag_AtA(tmp);
69 const_cast<vnl_vector<double> &>(jacobi_precond_) = vnl_vector<double> (tmp.size());
70 for (unsigned int i=0; i < tmp.size(); ++i)
71 const_cast<vnl_vector<double> &>(jacobi_precond_)[i] = 1.0 / double(tmp[i]);
72 }
73
74 px = dot_product(x,jacobi_precond_);
75 }
76
77 template class vnl_sparse_matrix_linear_system<double>;
78 template class vnl_sparse_matrix_linear_system<float>;
79