ug4
Loading...
Searching...
No Matches
mapsparsematrix_impl.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2012-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__UG__MAP_ALGEBRA__SPARSEMATRIX_IMPL__
34#define __H__UG__MAP_ALGEBRA__SPARSEMATRIX_IMPL__
35
36// creation etc
37
38#include <fstream>
39#include <cstring>
40
43#include "mapsparsematrix.h"
44#include <vector>
45#include <algorithm>
46
47
48
49namespace ug{
50
51template<typename T>
55
56template<typename T>
57bool MapSparseMatrix<T>::resize(size_t newRows, size_t newCols)
58{
59 PROFILE_BEGIN_GROUP(MapSparseMatrix_resize, "algebra MapSparseMatrix");
60 data.resize(newRows);
61 cols = newCols;
62 return true;
63}
64
65
66template<typename T>
68{
69 PROFILE_BEGIN_GROUP(MapSparseMatrix_set_as_transpose_of, "algebra MapSparseMatrix");
70 resize(B.num_cols(), B.num_rows());
71
72 for(size_t r=0; r<B.num_rows(); r++)
73 for(const_row_iterator it = B.begin_row(r); it != B.end_row(r); ++it)
74 operator()(it.index(), r) = scale*it.value();
75 return true;
76}
77
78template<typename T>
79template<typename vector_t>
80inline void MapSparseMatrix<T>::mat_mult_add_row(size_t row, typename vector_t::value_type &dest, double alpha, const vector_t &v) const
81{
82 for(const_row_iterator conn = begin_row(row); conn != end_row(row); ++conn)
83 MatMultAdd(dest, 1.0, dest, alpha, conn.value(), v[conn.index()]);
84}
85
86
87// calculate dest = alpha1*v1 + beta1*A*w1 (A = this matrix)
88template<typename T>
89template<typename vector_t>
90bool MapSparseMatrix<T>::axpy(vector_t &dest,
91 const number &alpha1, const vector_t &v1,
92 const number &beta1, const vector_t &w1) const
93{
94 PROFILE_BEGIN_GROUP(SparseMatrix_axpy, "algebra MapSparseMatrix");
95 if(alpha1 == 0.0)
96 {
97 for(size_t i=0; i < num_rows(); i++)
98 {
100 if(conn == end_row(i)) continue;
101 MatMult(dest[i], beta1, conn.value(), w1[conn.index()]);
102 for(++conn; conn != end_row(i); ++conn)
103 // res[i] += conn.value() * x[conn.index()];
104 MatMultAdd(dest[i], 1.0, dest[i], beta1, conn.value(), w1[conn.index()]);
105 }
106 }
107 else if(&dest == &v1)
108 {
109 if(alpha1 != 1.0) {
110 for(size_t i=0; i < num_rows(); i++)
111 {
112 dest[i] *= alpha1;
113 mat_mult_add_row(i, dest[i], beta1, w1);
114 }
115 }
116 else
117 for(size_t i=0; i < num_rows(); i++)
118 mat_mult_add_row(i, dest[i], beta1, w1);
119
120 }
121 else
122 {
123 for(size_t i=0; i < num_rows(); i++)
124 {
125 VecScaleAssign(dest[i], alpha1, v1[i]);
126 mat_mult_add_row(i, dest[i], beta1, w1);
127 }
128 }
129 return true;
130}
131
132// calculate dest = alpha1*v1 + beta1*A^T*w1 (A = this matrix)
133template<typename T>
134template<typename vector_t>
136 const number &alpha1, const vector_t &v1,
137 const number &beta1, const vector_t &w1) const
138{
139 PROFILE_BEGIN_GROUP(MapSparseMatrix_axpy_transposed, "algebra MapSparseMatrix");
140
141 if(&dest == &v1) {
142 if(alpha1 == 0.0)
143 dest.set(0.0);
144 else if(alpha1 != 1.0)
145 dest *= alpha1;
146 }
147 else if(alpha1 == 0.0)
148 dest.set(0.0);
149 else
150 VecScaleAssign(dest, alpha1, v1);
151
152 for(size_t i=0; i<num_rows(); i++)
153 {
154 for(const_row_iterator conn = begin_row(i); conn != end_row(i); ++conn)
155 {
156 if(conn.value() != 0.0)
157 // dest[conn.index()] += beta1 * conn.value() * w1[i];
158 MatMultTransposedAdd(dest[conn.index()], 1.0, dest[conn.index()], beta1, conn.value(), w1[i]);
159 }
160 }
161 return true;
162}
163
164template<typename T>
166{
167 PROFILE_BEGIN_GROUP(MapSparseMatrix_set, "algebra MapSparseMatrix");
168
169 for(size_t row=0; row<num_rows(); row++)
170 for(row_iterator it = begin_row(row); it != end_row(row); ++it)
171 {
172 if(it.index() == row)
173 it.value() = a;
174 else
175 it.value() = 0.0;
176 }
177 return true;
178}
179
180
181template<typename T>
182inline bool MapSparseMatrix<T>::is_isolated(size_t i) const
183{
184 UG_ASSERT(i < num_rows() && i >= 0, *this << ": " << i << " out of bounds.");
185
186 for(const_row_iterator it = begin_row(i); it != end_row(i); ++it)
187 if(it.index() != i && it.value() != 0.0)
188 return false;
189 return true;
190}
191
192
193template<typename T>
194void MapSparseMatrix<T>::set_matrix_row(size_t row, connection *c, size_t nr)
195{
196 PROFILE_BEGIN_GROUP(MapSparseMatrix_set_matrix_row, "algebra MapSparseMatrix");
197 for(size_t i=0; i<nr; i++)
198 operator()(row, c[i].iIndex) = c[i].dValue;
199}
200
201template<typename T>
202void MapSparseMatrix<T>::add_matrix_row(size_t row, connection *c, size_t nr)
203{
204 PROFILE_BEGIN_GROUP(MapSparseMatrix_add_matrix_row, "algebra MapSparseMatrix");
205 for(size_t i=0; i<nr; i++)
206 operator()(row, c[i].iIndex) += c[i].dValue;
207}
208
209
210template<typename T>
212{
213 PROFILE_BEGIN_GROUP(MapSparseMatrix_set_as_copy_of, "algebra MapSparseMatrix");
214 resize(B.num_rows(), B.num_cols());
215 for(size_t i=0; i < B.num_rows(); i++)
216 for(const_row_iterator it = B.begin_row(i); it != B.end_row(i); ++it)
217 operator()(i, it.index()) = scale*it.value();
218 return true;
219}
220
221
222
223template<typename T>
225{
226 PROFILE_BEGIN_GROUP(MapSparseMatrix_scale, "algebra MapSparseMatrix");
227 for(size_t i=0; i < num_rows(); i++)
228 for(row_iterator it = begin_row(i); it != end_row(i); ++it)
229 it.value() *= d;
230 return true;
231}
232
233
234
235
236
237
238
239} // namespace ug
240
241#endif
242
Definition mapsparsematrix.h:334
size_t index() const
Definition mapsparsematrix.h:339
const value_type & value()
Definition mapsparsematrix.h:338
Definition mapsparsematrix.h:326
sparse matrix for big, variable sparse matrices.
Definition mapsparsematrix.h:54
bool axpy_transposed(vector_t &dest, const number &alpha1, const vector_t &v1, const number &beta1, const vector_t &w1) const
calculate dest = alpha1*v1 + beta1*A^T*w1 (A = this matrix)
Definition mapsparsematrix_impl.h:135
bool set_as_copy_of(const MapSparseMatrix< value_type > &B, double scale=1.0)
create/recreate this as a copy of SparseMatrix B
Definition mapsparsematrix_impl.h:211
void mat_mult_add_row(size_t row, typename vector_t::value_type &dest, double alpha, const vector_t &v) const
calculates dest += alpha * A[row, .] v;
Definition mapsparsematrix_impl.h:80
void set(const M &mat)
set local matrix
bool resize(size_t newRows, size_t newCols)
resizes the SparseMatrix
Definition mapsparsematrix_impl.h:57
size_t num_cols() const
returns the number of cols
Definition mapsparsematrix.h:290
size_t num_rows() const
returns number of rows
Definition mapsparsematrix.h:287
bool is_isolated(size_t i) const
check for isolated condition of an index
Definition mapsparsematrix_impl.h:182
row_iterator begin_row(size_t r)
Definition mapsparsematrix.h:343
MapSparseMatrix()
constructor for empty SparseMatrix
Definition mapsparsematrix_impl.h:52
bool axpy(vector_t &dest, const number &alpha1, const vector_t &v1, const number &beta1, const vector_t &w1) const
calculate dest = alpha1*v1 + beta1*A*w1 (A = this matrix)
Definition mapsparsematrix_impl.h:90
row_iterator end_row(size_t r)
Definition mapsparsematrix.h:344
void set_matrix_row(size_t row, connection *c, size_t nr)
Definition mapsparsematrix_impl.h:194
void add_matrix_row(size_t row, connection *c, size_t nr)
Definition mapsparsematrix_impl.h:202
bool set_as_transpose_of(const SparseMatrix< value_type > &B, double scale=1.0)
write in a empty SparseMatrix (this) the transpose SparseMatrix of B.
Definition mapsparsematrix_impl.h:67
bool scale(double d)
Definition mapsparsematrix_impl.h:224
sparse matrix for big, variable sparse matrices.
Definition sparsematrix.h:99
size_t num_rows() const
returns number of rows
Definition sparsematrix.h:342
row_iterator end_row(size_t r)
Definition sparsematrix.h:485
size_t num_cols() const
returns the number of cols
Definition sparsematrix.h:345
row_iterator begin_row(size_t r)
Definition sparsematrix.h:484
#define UG_ASSERT(expr, msg)
Definition assert.h:70
double number
Definition types.h:124
the ug namespace
bool MatMultTransposedAdd(vector_t &dest, const number &alpha1, const vector_t &v1, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = alpha1*v1 + beta1 * A1 *w1;
Definition operations_mat.h:121
bool MatMultAdd(vector_t &dest, const number &alpha1, const vector_t &v1, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = alpha1*v1 + beta1 * A1 *w1;
Definition operations_mat.h:68
bool MatMult(vector_t &dest, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = beta1 * A1;
Definition operations_mat.h:59
void VecScaleAssign(double &dest, double alpha1, const double &v1)
calculates dest = alpha1*v1. for doubles
Definition operations_vec.h:49
#define PROFILE_BEGIN_GROUP(name, groups)
Definition profiler.h:255
bool resize(size_t newRows, size_t newCols)
size_t num_rows() const
Definition sparsematrix_interface.h:38
const_row_iterator end_row(size_t row) const
const_row_iterator begin_row(size_t row) const
Definition mapsparsematrix.h:64
size_t iIndex
Definition mapsparsematrix.h:65