Line data Source code
1 : /*
2 : * Copyright (c) 2014-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 MATRIX_DIAGONAL_H_
34 : #define MATRIX_DIAGONAL_H_
35 :
36 : #include "lib_algebra/operator/interface/linear_iterator.h"
37 :
38 : /**
39 : * D = MatrixDiagonal(mat) creates a LinearOperator which acts like D = diag(mat)
40 : */
41 : namespace ug{
42 : template <typename M, typename X, typename Y = X>
43 : class MatrixDiagonal : public virtual ILinearOperator<X,Y>,
44 : public M
45 : {
46 : public:
47 : // Domain space
48 : typedef X domain_function_type;
49 :
50 : // Range space
51 : typedef Y codomain_function_type;
52 :
53 : // Matrix type
54 : typedef M matrix_type;
55 : typedef MatrixOperator<M, X, Y> mo_type;
56 :
57 : SmartPtr<mo_type> m_mo;
58 :
59 : public:
60 0 : MatrixDiagonal(SmartPtr<mo_type> mo) {
61 0 : m_mo = mo;
62 0 : }
63 :
64 : // Init Operator J(u)
65 0 : virtual void init(const X& u)
66 : {
67 0 : m_mo->init(u);
68 0 : }
69 :
70 : // Init Operator L
71 0 : virtual void init() { m_mo->init(); }
72 :
73 : // Apply Operator f = L*u (e.g. d = J(u)*c in iterative scheme)
74 0 : virtual void apply(Y& f, const X& u)
75 : {
76 0 : matrix_type &A = m_mo->get_matrix();
77 0 : for(size_t i=0; i<A.num_rows(); i++)
78 0 : f[i] = A(i,i)*u[i];
79 0 : }
80 :
81 : // Apply Operator, i.e. f = f - L*u;
82 0 : virtual void apply_sub(Y& f, const X& u)
83 : {
84 0 : matrix_type &A = m_mo->get_matrix();
85 0 : for(size_t i=0; i<A.num_rows(); i++)
86 0 : f[i] -= A(i,i)*u[i];
87 0 : }
88 : };
89 :
90 :
91 : /**
92 : * Dinv = MatrixDiagonalInverse(mat) creates a LinearOperator which acts like Dinv = diag(mat)^{-1}
93 : */
94 : template <typename M, typename X, typename Y = X>
95 : class MatrixDiagonalInverse : public virtual ILinearOperator<X,Y>,
96 : public M
97 : {
98 : public:
99 : // Domain space
100 : typedef X domain_function_type;
101 :
102 : // Range space
103 : typedef Y codomain_function_type;
104 :
105 : // Matrix type
106 : typedef M matrix_type;
107 : typedef MatrixOperator<M, X, Y> mo_type;
108 :
109 : SmartPtr<mo_type> m_mo;
110 :
111 : public:
112 0 : MatrixDiagonalInverse(SmartPtr<mo_type> mo) {
113 0 : m_mo = mo;
114 0 : }
115 :
116 : // Init Operator J(u)
117 0 : virtual void init(const X& u)
118 : {
119 0 : m_mo->init(u);
120 0 : }
121 :
122 : // Init Operator L
123 0 : virtual void init() { m_mo->init(); }
124 :
125 : // Apply Operator f = L*u (e.g. d = J(u)*c in iterative scheme)
126 0 : virtual void apply(Y& f, const X& u)
127 : {
128 0 : matrix_type &A = m_mo->get_matrix();
129 0 : for(size_t i=0; i<A.num_rows(); i++)
130 : InverseMatMult(f[i], 1.0, A(i,i), u[i]);
131 0 : }
132 :
133 : // Apply Operator, i.e. f = f - L*u;
134 0 : virtual void apply_sub(Y& f, const X& u)
135 : {
136 0 : matrix_type &A = m_mo->get_matrix();
137 : typename X::value_type t;
138 0 : for(size_t i=0; i<A.num_rows(); i++)
139 : {
140 : InverseMatMult(t, 1.0, A(i,i), u[i]);
141 0 : f[i] -= t;
142 : }
143 :
144 0 : }
145 : };
146 :
147 : }
148 : #endif /* MATRIX_DIAGONAL_H_ */
|
