ug4
Loading...
Searching...
No Matches
analyzing_solver.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 * Author: Martin Rupp
4 *
5 * This file is part of UG4.
6 *
7 * UG4 is free software: you can redistribute it and/or modify it under the
8 * terms of the GNU Lesser General Public License version 3 (as published by the
9 * Free Software Foundation) with the following additional attribution
10 * requirements (according to LGPL/GPL v3 §7):
11 *
12 * (1) The following notice must be displayed in the Appropriate Legal Notices
13 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 *
15 * (2) The following notice must be displayed at a prominent place in the
16 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 *
18 * (3) The following bibliography is recommended for citation and must be
19 * preserved in all covered files:
20 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 * parallel geometric multigrid solver on hierarchically distributed grids.
22 * Computing and visualization in science 16, 4 (2013), 151-164"
23 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 * flexible software system for simulating pde based models on high performance
25 * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 *
27 * This program is distributed in the hope that it will be useful,
28 * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 * GNU Lesser General Public License for more details.
31 */
32
33#ifndef __H__LIB_ALGEBRA__LAPACK_ANALYZING_SOLVER__
34#define __H__LIB_ALGEBRA__LAPACK_ANALYZING_SOLVER__
35
36#include "../interface/linear_operator_inverse.h"
37#include "../interface/matrix_operator.h"
38#include "common/error.h"
41
42namespace ug{
43
45
46template <typename M, typename X, typename Y = X>
48 : public virtual ILinearOperatorInverse<X,Y>
49{
50 public:
53
56
58 typedef M matrix_type;
59
60 public:
61 virtual bool apply(Y& u, const X& f)
62 {
63 return m_pLinearOperatorInverse->apply(u, f);
64 }
65
66 virtual bool apply_return_defect(Y& u, X& f)
67 {
68 return m_pLinearOperatorInverse->apply_return_defect(u, f);
69 }
70
72 {
73 m_pLinearOperatorInverse = pLinearOperatorInverse;
74 }
75
77 virtual ~AnalyzingSolver() {};
78
80 virtual bool supports_parallel() const
81 {
82 return m_pLinearOperatorInverse->supports_parallel();
83 }
84
85 public:
86
87
88
89
90 void check(const matrix_type &A)
91 {
92 UG_LOG("ANALYZING SOLVER:\n");
93 UG_LOG(" Matrix is of dimension " << A.num_rows() << "\n");
94 if(A.num_rows() != A.num_cols())
95 { UG_LOG(" Matrix is not quadratic???\n"); }
96 if(GetRows(A(0,0)) == 1)
97 { UG_LOG(" Submatrices are DOUBLE.\n")}
98 else
99 {UG_LOG(" Submatrices are Matrices of size " << GetRows(A(0,0)) << " x " << GetCols(A(0,0)) << "\n");}
100
101
103 // check symmetry
104
105
107
108
110 size_t m_size = A.num_rows() * nrOfRows;
111
113 mat.resize_and_clear(m_size, m_size);
114
115 for(size_t r=0; r<A.num_rows(); r++)
116 for(typename matrix_type::const_row_iterator it = A.begin_row(r); it != A.end_row(r); ++it)
117 {
118 size_t rr = r*nrOfRows;
119 size_t cc = it.index()*nrOfRows;
120 for(size_t r2=0; r2<nrOfRows; r2++)
121 for(size_t c2=0; c2<nrOfRows; c2++)
122 {
123 if(BlockRef(it.value(), r2, c2) != 0.0)
124 mat(rr + r2, cc + c2) = BlockRef(it.value(), r2, c2);
125 }
126 }
127 mat.defragment();
128 checksub(mat);
129 }
130 virtual bool init(SmartPtr<ILinearOperator<Y,X> > A, const Y&u)
131 {
132 check(A);
133 return m_pLinearOperatorInverse->init(A, u);
134 }
135
137 {
138 check(A);
139 return m_pLinearOperatorInverse->init(A);
140 }
141
143 {
144 // cast operator
146 A.template cast_dynamic<MatrixOperator<M,Y,X> >();
147
148 // check if correct types are present
149 if(op.invalid())
150 UG_THROW("IMatrixOperatorInverse::init:"
151 " Passed operator is not matrix-based.");
152
153 // forward request
154 check(*op);
155 }
156 virtual const char *name() const { return m_pLinearOperatorInverse->name(); }
157
158 virtual std::string config_string() const
159 {
160 return "AnalyzingSolver " + m_pLinearOperatorInverse->config_string();
161 }
162 private:
164};
165
166
167} // end namespace ug
168
169#endif /* __H__LIB_ALGEBRA__LAPACK_LU_OPERATOR__ */
Definition smart_pointer.h:108
bool invalid() const
returns true if the pointer is invalid, false if not.
Definition smart_pointer.h:212
Definition analyzing_solver.h:49
virtual ~AnalyzingSolver()
virtual destructor
Definition analyzing_solver.h:77
AnalyzingSolver(SmartPtr< ILinearOperatorInverse< X, Y > > pLinearOperatorInverse)
Definition analyzing_solver.h:71
virtual std::string config_string() const
returns information about configuration parameters
Definition analyzing_solver.h:158
virtual bool init(SmartPtr< ILinearOperator< Y, X > > A)
initializes for the inverse for a linear operator
Definition analyzing_solver.h:136
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition analyzing_solver.h:80
virtual bool init(SmartPtr< ILinearOperator< Y, X > > A, const Y &u)
initializes for the inverse for a linearized operator at linearization point u
Definition analyzing_solver.h:130
virtual bool apply_return_defect(Y &u, X &f)
applies inverse operator, i.e. returns u = A^{-1} f and returns defect d := f - A*u
Definition analyzing_solver.h:66
virtual bool apply(Y &u, const X &f)
applies inverse operator, i.e. returns u = A^{-1} f
Definition analyzing_solver.h:61
SmartPtr< ILinearOperatorInverse< X, Y > > m_pLinearOperatorInverse
Definition analyzing_solver.h:163
M matrix_type
Matrix type.
Definition analyzing_solver.h:58
Y codomain_function_type
Range space.
Definition analyzing_solver.h:55
void check(const matrix_type &A)
Definition analyzing_solver.h:90
virtual const char * name() const
returns the name of the operator inverse
Definition analyzing_solver.h:156
void check(SmartPtr< ILinearOperator< Y, X > > A)
Definition analyzing_solver.h:142
X domain_function_type
Domain space.
Definition analyzing_solver.h:52
describes a linear mapping X->Y
Definition linear_operator.h:80
describes an inverse linear mapping X->Y
Definition linear_operator_inverse.h:80
Wrapper for sequential matrices to handle them in parallel.
Definition parallel_matrix.h:65
ParallelMatrix< SparseMatrix< double > > matrix_type
Definition cpu_algebra_types.h:79
#define UG_THROW(msg)
Definition error.h:57
#define UG_LOG(msg)
Definition log.h:367
the ug namespace
size_t GetCols(const T &t)
double & BlockRef(T &vec, size_t i)
Definition blocks.h:66
void checksub(const CPUAlgebra::matrix_type &A)
Definition analyzing_solver.cpp:37
size_t GetRows(const T &t)
Definition communication_policies.h:58