Go to the documentation of this file.00001
00002 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00003 #pragma implementation
00004 #endif
00005
00006 #include "vnl_sparse_matrix_linear_system.h"
00007 #include <vcl_cassert.h>
00008 #include <vnl/vnl_copy.h>
00009
00010 VCL_DEFINE_SPECIALIZATION
00011 void vnl_sparse_matrix_linear_system<double>::get_rhs(vnl_vector<double>& b) const
00012 {
00013 b = b_;
00014 }
00015
00016 VCL_DEFINE_SPECIALIZATION
00017 void vnl_sparse_matrix_linear_system<double>::transpose_multiply(vnl_vector<double> const& b, vnl_vector<double> & x) const
00018 {
00019 A_.pre_mult(b,x);
00020 }
00021
00022 VCL_DEFINE_SPECIALIZATION
00023 void vnl_sparse_matrix_linear_system<float>::get_rhs(vnl_vector<double>& b) const
00024 {
00025 vnl_copy(b_, b);
00026 }
00027
00028 VCL_DEFINE_SPECIALIZATION
00029 void vnl_sparse_matrix_linear_system<float>::transpose_multiply(vnl_vector<double> const& b, vnl_vector<double> & x) const
00030 {
00031 static vnl_vector<float> x_float;
00032 static vnl_vector<float> b_float;
00033
00034 if (x_float.size() != x.size()) x_float = vnl_vector<float> (x.size());
00035 if (b_float.size() != b.size()) b_float = vnl_vector<float> (b.size());
00036
00037 vnl_copy(b, b_float);
00038 A_.pre_mult(b_float,x_float);
00039 vnl_copy(x_float, x);
00040 }
00041
00042 VCL_DEFINE_SPECIALIZATION
00043 void vnl_sparse_matrix_linear_system<double>::multiply(vnl_vector<double> const& x, vnl_vector<double> & b) const
00044 {
00045 A_.mult(x,b);
00046 }
00047
00048
00049 VCL_DEFINE_SPECIALIZATION
00050 void vnl_sparse_matrix_linear_system<float>::multiply(vnl_vector<double> const& x, vnl_vector<double> & b) const
00051 {
00052 static vnl_vector<float> x_float;
00053 static vnl_vector<float> b_float;
00054
00055 if (x_float.size() != x.size()) x_float = vnl_vector<float> (x.size());
00056 if (b_float.size() != b.size()) b_float = vnl_vector<float> (b.size());
00057
00058 vnl_copy(x, x_float);
00059 A_.mult(x_float,b_float);
00060 vnl_copy(b_float, b);
00061 }
00062
00063
00064 template<class T>
00065 void vnl_sparse_matrix_linear_system<T>::apply_preconditioner(vnl_vector<double> const& x, vnl_vector<double> & px) const
00066 {
00067 assert(x.size() == px.size());
00068
00069 if (jacobi_precond_.size() == 0) {
00070 vnl_vector<T> tmp(get_number_of_unknowns());
00071 A_.diag_AtA(tmp);
00072 const_cast<vnl_vector<double> &>(jacobi_precond_) = vnl_vector<double> (tmp.size());
00073 for (unsigned int i=0; i < tmp.size(); ++i)
00074 const_cast<vnl_vector<double> &>(jacobi_precond_)[i] = 1.0 / double(tmp[i]);
00075 }
00076
00077 px = dot_product(x,jacobi_precond_);
00078 }
00079
00080 template class vnl_sparse_matrix_linear_system<double>;
00081 template class vnl_sparse_matrix_linear_system<float>;
00082
