eigen_sparse_matrix.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 
22 // Local includes
23 #include "libmesh/libmesh_config.h"
24 
25 #ifdef LIBMESH_HAVE_EIGEN
26 
29 #include "libmesh/dense_matrix.h"
30 #include "libmesh/dof_map.h"
32 
33 namespace libMesh
34 {
35 
36 
37 //-----------------------------------------------------------------------
38 // EigenSparseMatrix members
39 template <typename T>
41  const numeric_index_type n_in,
42  const numeric_index_type m_l,
43  const numeric_index_type n_l,
44  const numeric_index_type nnz,
45  const numeric_index_type,
46  const numeric_index_type)
47 {
48  // noz ignored... only used for multiple processors!
49  libmesh_assert_equal_to (m_in, m_l);
50  libmesh_assert_equal_to (n_in, n_l);
51  libmesh_assert_equal_to (m_in, n_in);
52  libmesh_assert_greater (nnz, 0);
53 
54  _mat.resize(m_in, n_in);
55  _mat.reserve(Eigen::Matrix<numeric_index_type, Eigen::Dynamic, 1>::Constant(m_in,nnz));
56 
57  this->_is_initialized = true;
58 }
59 
60 
61 
62 template <typename T>
64 {
65  // Ignore calls on initialized objects
66  if (this->initialized())
67  return;
68 
69  // We need the DofMap for this!
70  libmesh_assert(this->_dof_map);
71 
72  // Clear intialized matrices
73  if (this->initialized())
74  this->clear();
75 
76  const numeric_index_type n_rows = this->_dof_map->n_dofs();
77  const numeric_index_type n_cols = n_rows;
78 
79 #ifndef NDEBUG
80  // The following variables are only used for assertions,
81  // so avoid declaring them when asserts are inactive.
82  const numeric_index_type n_l = this->_dof_map->n_dofs_on_processor(0);
83  const numeric_index_type m_l = n_l;
84 #endif
85 
86  // Laspack Matrices only work for uniprocessor cases
87  libmesh_assert_equal_to (n_rows, n_cols);
88  libmesh_assert_equal_to (m_l, n_rows);
89  libmesh_assert_equal_to (n_l, n_cols);
90 
91  const std::vector<numeric_index_type>& n_nz = this->_dof_map->get_n_nz();
92 
93 #ifndef NDEBUG
94  // The following variables are only used for assertions,
95  // so avoid declaring them when asserts are inactive.
96  const std::vector<numeric_index_type>& n_oz = this->_dof_map->get_n_oz();
97 #endif
98 
99  // Make sure the sparsity pattern isn't empty
100  libmesh_assert_equal_to (n_nz.size(), n_l);
101  libmesh_assert_equal_to (n_oz.size(), n_l);
102 
103  if (n_rows==0)
104  {
105  _mat.resize(0,0);
106  return;
107  }
108 
109  _mat.resize(n_rows,n_cols);
110  _mat.reserve(n_nz);
111 
112  this->_is_initialized = true;
113 
114  libmesh_assert_equal_to (n_rows, this->m());
115  libmesh_assert_equal_to (n_cols, this->n());
116 }
117 
118 
119 
120 template <typename T>
122  const std::vector<numeric_index_type>& rows,
123  const std::vector<numeric_index_type>& cols)
124 
125 {
126  libmesh_assert (this->initialized());
127  unsigned int n_rows = libmesh_cast_int<unsigned int>(rows.size());
128  unsigned int n_cols = libmesh_cast_int<unsigned int>(cols.size());
129  libmesh_assert_equal_to (dm.m(), n_rows);
130  libmesh_assert_equal_to (dm.n(), n_cols);
131 
132 
133  for (unsigned int i=0; i<n_rows; i++)
134  for (unsigned int j=0; j<n_cols; j++)
135  this->add(rows[i],cols[j],dm(i,j));
136 }
137 
138 
139 
140 template <typename T>
142 {
143  EigenSparseVector<T>& dest = libmesh_cast_ref<EigenSparseVector<T>&>(dest_in);
144 
145  dest._vec = _mat.diagonal();
146 }
147 
148 
149 
150 template <typename T>
152 {
153  EigenSparseMatrix<T>& dest = libmesh_cast_ref<EigenSparseMatrix<T>&>(dest_in);
154 
155  dest._mat = _mat.transpose();
156 }
157 
158 
159 
160 template <typename T>
162  SparseMatrix<T>(comm),
163  _closed (false)
164 {
165 }
166 
167 
168 
169 template <typename T>
171 {
172  this->clear ();
173 }
174 
175 
176 
177 template <typename T>
179 {
180  _mat.resize(0,0);
181 
182  _closed = false;
183  this->_is_initialized = false;
184 }
185 
186 
187 
188 template <typename T>
190 {
191  _mat.setZero();
192 }
193 
194 
195 
196 template <typename T>
198 {
199  libmesh_assert (this->initialized());
200 
201  return _mat.rows();
202 }
203 
204 
205 
206 template <typename T>
208 {
209  libmesh_assert (this->initialized());
210 
211  return _mat.cols();
212 }
213 
214 
215 
216 template <typename T>
218 {
219  return 0;
220 }
221 
222 
223 
224 template <typename T>
226 {
227  return this->m();
228 }
229 
230 
231 
232 template <typename T>
234  const numeric_index_type j,
235  const T value)
236 {
237  libmesh_assert (this->initialized());
238  libmesh_assert_less (i, this->m());
239  libmesh_assert_less (j, this->n());
240 
241  _mat.coeffRef(i,j) = value;
242 }
243 
244 
245 
246 template <typename T>
248  const numeric_index_type j,
249  const T value)
250 {
251  libmesh_assert (this->initialized());
252  libmesh_assert_less (i, this->m());
253  libmesh_assert_less (j, this->n());
254 
255  _mat.coeffRef(i,j) += value;
256 }
257 
258 
259 
260 template <typename T>
262  const std::vector<numeric_index_type>& dof_indices)
263 {
264  this->add_matrix (dm, dof_indices, dof_indices);
265 }
266 
267 
268 
269 template <typename T>
270 void EigenSparseMatrix<T>::add (const T a_in, SparseMatrix<T> &X_in)
271 {
272  libmesh_assert (this->initialized());
273  libmesh_assert_equal_to (this->m(), X_in.m());
274  libmesh_assert_equal_to (this->n(), X_in.n());
275 
276  EigenSparseMatrix<T> &X = libmesh_cast_ref<EigenSparseMatrix<T>&> (X_in);
277 
278  _mat += X._mat*a_in;
279 }
280 
281 
282 
283 
284  template <typename T>
286  const numeric_index_type j) const
287 {
288  libmesh_assert (this->initialized());
289  libmesh_assert_less (i, this->m());
290  libmesh_assert_less (j, this->n());
291 
292  return _mat.coeff(i,j);
293 }
294 
295 
296 
297 //------------------------------------------------------------------
298 // Explicit instantiations
299 template class EigenSparseMatrix<Number>;
300 
301 } // namespace libMesh
302 
303 
304 #endif // #ifdef LIBMESH_HAVE_EIGEN

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo