Line data Source code
1 : /*
2 : * Copyright (c) 2010-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 : /*
34 : * blockMatrix.h
35 : *
36 : * Header File for general block matrix / double accessing
37 : * i.e. setAt(mat, i, j, d) -> mat(i,j) = d
38 : * and setAt(f, 0, 0, d) -> f = d.
39 : * This means doubles and block matrices can be accessed
40 : * by the same methods.
41 : */
42 :
43 : // todo: also with float / complex<float> / complex<double>
44 :
45 : #ifndef __H__UG__SMALL_ALGEBRA__DOUBLE__
46 : #define __H__UG__SMALL_ALGEBRA__DOUBLE__
47 :
48 : #include "blocks.h"
49 : #include "common/common.h"
50 :
51 : namespace ug{
52 :
53 :
54 : //////////////////////////////////////////////////////
55 : template<typename T> number BlockNorm(const T &t);
56 : template <>
57 : inline number BlockNorm(const number &a)
58 : {
59 0 : return a>0 ? a : -a;
60 : }
61 :
62 : template<typename T> number BlockNorm2(const T &t);
63 : template <>
64 : inline number BlockNorm2(const number &a)
65 : {
66 0 : return a*a;
67 : }
68 :
69 : template<typename T> number BlockMaxNorm(const T &t);
70 : template <>
71 : inline number BlockMaxNorm(const number &a)
72 : {
73 0 : return a>0 ? a : -a;
74 : }
75 :
76 : //////////////////////////////////////////////////////
77 : // get/set specialization for numbers
78 :
79 : template<> inline number &BlockRef(number &m, size_t i)
80 : {
81 : UG_ASSERT(i == 0, "block is number, doesnt have component (" << i << ").");
82 : return m;
83 : }
84 : template<> inline const number &BlockRef(const number &m, size_t i)
85 : {
86 : UG_ASSERT(i == 0, "block is number, doesnt have component (" << i << ").");
87 : return m;
88 : }
89 :
90 : template<> inline number &BlockRef(number &m, size_t i, size_t j)
91 : {
92 : UG_ASSERT(i == 0 && j == 0, "block is number, doesnt have component (" << i << ", " << j << ").");
93 : return m;
94 : }
95 : template<> inline const number &BlockRef(const number &m, size_t i, size_t j)
96 : {
97 : UG_ASSERT(i == 0 && j == 0, "block is number, doesnt have component (" << i << ", " << j << ").");
98 : return m;
99 : }
100 :
101 : //////////////////////////////////////////////////////
102 : // algebra stuff to avoid temporary variables
103 :
104 : inline void AssignMult(number &dest, const number &b, const number &vec)
105 : {
106 0 : dest = b*vec;
107 : }
108 : // dest += vec*b
109 : inline void AddMult(number &dest, const number &b, const number &vec)
110 : {
111 0 : dest += b*vec;
112 0 : }
113 :
114 :
115 : // dest -= vec*b
116 : inline void SubMult(number &dest, const number &b, const number &vec)
117 : {
118 : dest -= b*vec;
119 : }
120 :
121 :
122 : //////////////////////////////////////////////////////
123 : //setSize(t, a, b) for numbers
124 : template<>
125 : inline void SetSize(number &d, size_t a)
126 : {
127 : UG_ASSERT(a == 1, "block is number, cannot change size to " << a << ".");
128 : return;
129 : }
130 :
131 : template<>
132 : inline void SetSize(number &d, size_t a, size_t b)
133 : {
134 : UG_ASSERT(a == 1 && b == 1, "block is number, cannot change size to (" << a << ", " << b << ").");
135 : return;
136 : }
137 :
138 : template<>
139 : inline size_t GetSize(const number &t)
140 : {
141 : return 1;
142 : }
143 :
144 : template<>
145 : inline size_t GetRows(const number &t)
146 : {
147 : return 1;
148 : }
149 :
150 : template<>
151 : inline size_t GetCols(const number &t)
152 : {
153 : return 1;
154 : }
155 : ///////////////////////////////////////////////////////////////////
156 :
157 : inline bool InverseMatMult(number &dest, const double &beta, const number &mat, const number &vec)
158 : {
159 0 : dest = beta*vec/mat;
160 : return true;
161 : }
162 :
163 : ///////////////////////////////////////////////////////////////////
164 : // traits: information for numbers
165 :
166 :
167 : template<>
168 : struct block_traits<number>
169 : {
170 : typedef number vec_type;
171 : typedef number inverse_type;
172 :
173 : enum { is_static = true};
174 : enum { static_num_rows = 1};
175 : enum { static_num_cols = 1};
176 : enum { static_size = 1 };
177 : enum { depth = 0 };
178 : };
179 :
180 : template<> struct block_multiply_traits<number, number>
181 : {
182 : typedef number ReturnType;
183 : };
184 :
185 : template<typename T>
186 : struct block_multiply_traits<number, T>
187 : {
188 : typedef T ReturnType;
189 : };
190 :
191 : template<typename T>
192 : struct block_multiply_traits<T, number>
193 : {
194 : typedef T ReturnType;
195 : };
196 :
197 : inline bool GetInverse(number &inv, const number &m)
198 : {
199 0 : inv = 1.0/m;
200 : return (m != 0.0);
201 : }
202 :
203 : inline bool Invert(number &m)
204 : {
205 : bool b = (m != 0.0);
206 : m = 1/m;
207 : return b;
208 : }
209 :
210 : } // namespace ug
211 :
212 : #endif
213 :
214 :
215 :
|
