solver(2rheolef) [debian man page]
solver(2rheolef) rheolef-6.1 solver(2rheolef) NAME
solver - direct or interative solver interface DESCRIPTION
The class implements a matrix factorization: LU factorization for an unsymmetric matrix and Choleski fatorisation for a symmetric one. Let a be a square invertible matrix in csr format (see csr(2)). csr<Float> a; We get the factorization by: solver<Float> sa (a); Each call to the direct solver for a*x = b writes either: vec<Float> x = sa.solve(b); When the matrix is modified in a computation loop but conserves its sparsity pattern, an efficient re-factorization writes: sa.update_values (new_a); x = sa.solve(b); This approach skip the long step of the symbolic factization step. ITERATIVE SOLVER
The factorization can also be incomplete, i.e. a pseudo-inverse, suitable for preconditionning iterative methods. In that case, the sa.solve(b) call runs a conjugate gradient when the matrix is symmetric, or a generalized minimum residual algorithm when the matrix is unsymmetric. AUTOMATIC CHOICE AND CUSTOMIZATION
The symmetry of the matrix is tested via the a.is_symmetric() property (see csr(2)) while the choice between direct or iterative solver is switched from the a.pattern_dimension() value. When the pattern is 3D, an iterative method is faster and less memory consuming. Otherwhise, for 1D or 2D problems, the direct method is prefered. These default choices can be supersetted by using explicit options: solver_option_type opt; opt.iterative = true; solver<Float> sa (a, opt); See the solver.h header for the complete list of available options. IMPLEMENTATION NOTE
The implementation bases on the pastix library. IMPLEMENTATION
template <class T, class M = rheo_default_memory_model> class solver_basic : public smart_pointer<solver_rep<T,M> > { public: // typedefs: typedef solver_rep<T,M> rep; typedef smart_pointer<rep> base; // allocator: solver_basic (); explicit solver_basic (const csr<T,M>& a, const solver_option_type& opt = solver_option_type()); void update_values (const csr<T,M>& a); // accessors: vec<T,M> trans_solve (const vec<T,M>& b) const; vec<T,M> solve (const vec<T,M>& b) const; }; // factorizations: template <class T, class M> solver_basic<T,M> ldlt(const csr<T,M>& a, const solver_option_type& opt = solver_option_type()); template <class T, class M> solver_basic<T,M> lu (const csr<T,M>& a, const solver_option_type& opt = solver_option_type()); template <class T, class M> solver_basic<T,M> ic0 (const csr<T,M>& a, const solver_option_type& opt = solver_option_type()); template <class T, class M> solver_basic<T,M> ilu0(const csr<T,M>& a, const solver_option_type& opt = solver_option_type()); typedef solver_basic<Float> solver; SEE ALSO
csr(2), csr(2) rheolef-6.1 rheolef-6.1 solver(2rheolef)
Check Out this Related Man Page
tensor(2rheolef) rheolef-6.1 tensor(2rheolef) NAME
tensor - a N*N tensor, N=1,2,3 SYNOPSYS
The tensor class defines a 3*3 tensor, as the value of a tensorial valued field. Basic algebra with scalars, vectors of R^3 (i.e. the point class) and tensor objects are supported. IMPLEMENTATION
template<class T> class tensor_basic { public: typedef size_t size_type; typedef T element_type; // allocators: tensor_basic (const T& init_val = 0); tensor_basic (T x[3][3]); tensor_basic (const tensor_basic<T>& a); // affectation: tensor_basic<T>& operator = (const tensor_basic<T>& a); tensor_basic<T>& operator = (const T& val); // modifiers: void fill (const T& init_val); void reset (); void set_row (const point_basic<T>& r, size_t i, size_t d = 3); void set_column (const point_basic<T>& c, size_t j, size_t d = 3); // accessors: T& operator()(size_type i, size_type j); T operator()(size_type i, size_type j) const; point_basic<T> row(size_type i) const; point_basic<T> col(size_type i) const; size_t nrow() const; // = 3, for template matrix compatibility size_t ncol() const; // inputs/outputs: std::ostream& put (std::ostream& s, size_type d = 3) const; std::istream& get (std::istream&); // algebra: bool operator== (const tensor_basic<T>&) const; bool operator!= (const tensor_basic<T>& b) const { return ! operator== (b); } template <class U> friend tensor_basic<U> operator- (const tensor_basic<U>&); template <class U> friend tensor_basic<U> operator+ (const tensor_basic<U>&, const tensor_basic<U>&); template <class U> friend tensor_basic<U> operator- (const tensor_basic<U>&, const tensor_basic<U>&); template <class U> friend tensor_basic<U> operator* (int k, const tensor_basic<U>& a); template <class U> friend tensor_basic<U> operator* (const U& k, const tensor_basic<U>& a); template <class U> friend tensor_basic<U> operator* (const tensor_basic<U>& a, int k); template <class U> friend tensor_basic<U> operator* (const tensor_basic<U>& a, const U& k); template <class U> friend tensor_basic<U> operator/ (const tensor_basic<U>& a, int k); template <class U> friend tensor_basic<U> operator/ (const tensor_basic<U>& a, const U& k); template <class U> friend point_basic<U> operator* (const tensor_basic<U>&, const point_basic<U>&); template <class U> friend point_basic<U> operator* (const point_basic<U>& yt, const tensor_basic<U>& a); point_basic<T> trans_mult (const point_basic<T>& x) const; template <class U> friend tensor_basic<U> trans (const tensor_basic<U>& a, size_t d = 3); template <class U> friend tensor_basic<U> operator* (const tensor_basic<U>& a, const tensor_basic<U>& b); template <class U> friend void prod (const tensor_basic<U>& a, const tensor_basic<U>& b, tensor_basic<U>& result, size_t di=3, size_t dj=3, size_t dk=3); template <class U> friend tensor_basic<U> inv (const tensor_basic<U>& a, size_t d = 3); template <class U> friend tensor_basic<U> diag (const point_basic<U>& d); template <class U> friend tensor_basic<U> identity (size_t d=3); template <class U> friend tensor_basic<U> dyadic (const point_basic<U>& u, const point_basic<U>& v, size_t d=3); // metric and geometric transformations: template <class U> friend U dotdot (const tensor_basic<U>&, const tensor_basic<U>&); template <class U> friend U norm2 (const tensor_basic<U>& a) { return dotdot(a,a); } template <class U> friend U dist2 (const tensor_basic<U>& a, const tensor_basic<U>& b) { return norm2(a-b); } template <class U> friend U norm (const tensor_basic<U>& a) { return ::sqrt(norm2(a)); } template <class U> friend U dist (const tensor_basic<U>& a, const tensor_basic<U>& b) { return norm(a-b); } T determinant (size_type d = 3) const; template <class U> friend U determinant (const tensor_basic<U>& A, size_t d = 3); template <class U> friend bool invert_3x3 (const tensor_basic<U>& A, tensor_basic<U>& result); // spectral: // eigenvalues & eigenvectors: // a = q*d*q^T // a may be symmetric // where q=(q1,q2,q3) are eigenvectors in rows (othonormal matrix) // and d=(d1,d2,d3) are eigenvalues, sorted in decreasing order d1 >= d2 >= d3 // return d point_basic<T> eig (tensor_basic<T>& q, size_t dim = 3) const; point_basic<T> eig (size_t dim = 3) const; // singular value decomposition: // a = u*s*v^T // a can be unsymmetric // where u=(u1,u2,u3) are left pseudo-eigenvectors in rows (othonormal matrix) // v=(v1,v2,v3) are right pseudo-eigenvectors in rows (othonormal matrix) // and s=(s1,s2,s3) are eigenvalues, sorted in decreasing order s1 >= s2 >= s3 // return s point_basic<T> svd (tensor_basic<T>& u, tensor_basic<T>& v, size_t dim = 3) const; // data: T _x[3][3]; }; typedef tensor_basic<Float> tensor; // inputs/outputs: template<class T> inline std::istream& operator>> (std::istream& in, tensor_basic<T>& a) { return a.get (in); } template<class T> inline std::ostream& operator<< (std::ostream& out, const tensor_basic<T>& a) { return a.put (out); } // t = a otimes b template<class T> tensor_basic<T> otimes (const point_basic<T>& a, const point_basic<T>& b, size_t na = 3); // t += a otimes b template<class T> void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na = 3); template<class T> void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na, size_t nb); rheolef-6.1 rheolef-6.1 tensor(2rheolef)