Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
sparse_matrix_ez.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: sparse_matrix_ez.templates.h 30040 2013-07-18 17:06:48Z maier @f$
3 //
4 // Copyright (C) 2002 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__sparse_matrix_ez_templates_h
18 #define __deal2__sparse_matrix_ez_templates_h
19 
20 
21 #include <deal.II/lac/sparse_matrix_ez.h>
22 #include <deal.II/lac/vector.h>
23 
24 #include <iostream>
25 #include <iomanip>
26 #include <algorithm>
27 #include <cmath>
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 //---------------------------------------------------------------------------
32 
33 template <typename number>
35 {
36  n_columns = 0;
37 }
38 
39 
40 template <typename number>
42  :
43  Subscriptor (m),
44  n_columns (0)
45 {
46  Assert(m.n() == 0, ExcNotImplemented());
47  Assert(m.m() == 0, ExcNotImplemented());
48 }
49 
50 
51 template <typename number>
53  const size_type n_cols,
54  const size_type default_row_length,
55  const unsigned int default_increment)
56 {
57  reinit(n_rows, n_cols, default_row_length, default_increment);
58 }
59 
60 
61 template <typename number>
63 {}
64 
65 
66 template <typename number>
69 {
71  return *this;
72 }
73 
74 
75 template <typename number>
78 {
80 
81  typename std::vector<Entry>::iterator e = data.begin();
82  const typename std::vector<Entry>::iterator end = data.end();
83 
84  while (e != end)
85  {
86  (e++)->value = 0.;
87  }
88 
89  return *this;
90 }
91 
92 
93 
94 template <typename number>
95 void
97  const size_type n_cols,
98  size_type default_row_length,
99  unsigned int default_increment,
100  size_type reserve)
101 {
102  clear();
103 
104  increment = default_increment;
105 
106  n_columns = n_cols;
107  row_info.resize(n_rows);
108  if (reserve != 0)
109  data.reserve(reserve);
110  data.resize(default_row_length * n_rows);
111 
112  for (size_type i=0; i<n_rows; ++i)
113  row_info[i].start = i * default_row_length;
114 }
115 
116 
117 template <typename number>
118 void
120 {
121  n_columns = 0;
122  row_info.resize(0);
123  data.resize(0);
124 }
125 
126 
127 template <typename number>
128 bool
130 {
131  return ((n_columns == 0) && (row_info.size()==0));
132 }
133 
134 
135 template <typename number>
136 template <typename somenumber>
137 void
139  const Vector<somenumber> &src) const
140 {
141  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
142  Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
143 
144  const size_type end_row = row_info.size();
145  for (size_type row = 0; row < end_row; ++row)
146  {
147  const RowInfo &ri = row_info[row];
148  typename std::vector<Entry>::const_iterator
149  entry = data.begin() + ri.start;
150  double s = 0.;
151  for (unsigned short i=0; i<ri.length; ++i,++entry)
152  {
153  Assert (entry->column != Entry::invalid,
154  ExcInternalError());
155  s += entry->value * src(entry->column);
156  }
157  dst(row) = s;
158  }
159 }
160 
161 
162 template <typename number>
163 number
165 {
166  number sum = 0.;
167  const_iterator start = begin();
168  const_iterator final = end();
169 
170  while (start != final)
171  {
172  const double value = start->value();
173  sum += value*value;
174  ++start;
175  }
176  return std::sqrt(sum);
177 }
178 
179 
180 
181 template <typename number>
182 template <typename somenumber>
183 void
185  const Vector<somenumber> &src) const
186 {
187  dst = 0.;
188  Tvmult_add(dst, src);
189 }
190 
191 
192 template <typename number>
193 template <typename somenumber>
194 void
196  const Vector<somenumber> &src) const
197 {
198  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
199  Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
200 
201  const size_type end_row = row_info.size();
202  for (size_type row = 0; row < end_row; ++row)
203  {
204  const RowInfo &ri = row_info[row];
205  typename std::vector<Entry>::const_iterator
206  entry = data.begin() + ri.start;
207  double s = 0.;
208  for (unsigned short i=0; i<ri.length; ++i,++entry)
209  {
210  Assert (entry->column != Entry::invalid,
211  ExcInternalError());
212  s += entry->value * src(entry->column);
213  }
214  dst(row) += s;
215  }
216 }
217 
218 template <typename number>
219 template <typename somenumber>
220 void
222  const Vector<somenumber> &src) const
223 {
224  Assert(n() == dst.size(), ExcDimensionMismatch(n(),dst.size()));
225  Assert(m() == src.size(), ExcDimensionMismatch(m(),src.size()));
226 
227  const size_type end_row = row_info.size();
228  for (size_type row = 0; row < end_row; ++row)
229  {
230  const RowInfo &ri = row_info[row];
231  typename std::vector<Entry>::const_iterator
232  entry = data.begin() + ri.start;
233  for (unsigned short i=0; i<ri.length; ++i,++entry)
234  {
235  Assert (entry->column != Entry::invalid,
236  ExcInternalError());
237  dst(entry->column) += entry->value * src(row);
238  }
239  }
240 }
241 
242 
243 template <typename number>
244 template <typename somenumber>
245 void
247  const Vector<somenumber> &src,
248  const number om) const
249 {
250  Assert (m() == n(), ExcNotQuadratic());
251  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
252  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
253 
254  somenumber *dst_ptr = dst.begin();
255  const somenumber *src_ptr = src.begin();
256  typename std::vector<RowInfo>::const_iterator ri = row_info.begin();
257  const typename std::vector<RowInfo>::const_iterator end = row_info.end();
258 
259  for (; ri != end; ++dst_ptr, ++src_ptr, ++ri)
260  {
261  Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
262  *dst_ptr = om **src_ptr / data[ri->start + ri->diagonal].value;
263  }
264 }
265 
266 
267 
268 template <typename number>
269 template <typename somenumber>
270 void
272  const Vector<somenumber> &src,
273  const number om) const
274 {
275  Assert (m() == n(), ExcNotQuadratic());
276  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
277  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
278 
279  somenumber *dst_ptr = dst.begin();
280  const somenumber *src_ptr = src.begin();
281  typename std::vector<RowInfo>::const_iterator ri = row_info.begin();
282  const typename std::vector<RowInfo>::const_iterator end = row_info.end();
283 
284  for (; ri != end; ++dst_ptr, ++src_ptr, ++ri)
285  {
286  Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
287  number s = *src_ptr;
288  const size_type end_row = ri->start + ri->diagonal;
289  for (size_type i=ri->start; i<end_row; ++i)
290  s -= data[i].value * dst(data[i].column);
291 
292  *dst_ptr = om * s / data[ri->start + ri->diagonal].value;
293  }
294 }
295 
296 
297 
298 template <typename number>
299 template <typename somenumber>
300 void
302  const Vector<somenumber> &src,
303  const number om) const
304 {
305  Assert (m() == n(), ExcNotQuadratic());
306  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
307  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
308 
309  somenumber *dst_ptr = dst.begin()+dst.size()-1;
310  const somenumber *src_ptr = src.begin()+src.size()-1;
311  typename std::vector<RowInfo>::const_reverse_iterator
312  ri = row_info.rbegin();
313  const typename std::vector<RowInfo>::const_reverse_iterator
314  end = row_info.rend();
315 
316  for (; ri != end; --dst_ptr, --src_ptr, ++ri)
317  {
318  Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
319  number s = *src_ptr;
320  const size_type end_row = ri->start + ri->length;
321  for (size_type i=ri->start+ri->diagonal+1; i<end_row; ++i)
322  s -= data[i].value * dst(data[i].column);
323 
324  *dst_ptr = om * s / data[ri->start + ri->diagonal].value;
325  }
326 }
327 
328 
329 
330 template <typename number>
331 template <typename somenumber>
332 void
334  const Vector<somenumber> &src,
335  const number om,
336  const std::vector<std::size_t> &) const
337 {
338  Assert (m() == n(), ExcNotQuadratic());
339  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
340  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
341 
342  somenumber *dst_ptr = dst.begin();
343  const somenumber *src_ptr = src.begin();
344  typename std::vector<RowInfo>::const_iterator ri;
345  const typename std::vector<RowInfo>::const_iterator end = row_info.end();
346 
347  // Forward
348  for (ri = row_info.begin(); ri != end; ++dst_ptr, ++src_ptr, ++ri)
349  {
350  Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
351  number s = 0;
352  const size_type end_row = ri->start + ri->diagonal;
353  for (size_type i=ri->start; i<end_row; ++i)
354  s += data[i].value * dst(data[i].column);
355 
356  *dst_ptr = *src_ptr - s * om;
357  *dst_ptr /= data[ri->start + ri->diagonal].value;
358  }
359  // Diagonal
360  dst_ptr = dst.begin();
361  for (ri = row_info.begin(); ri != end; ++dst_ptr, ++ri)
362  *dst_ptr *= om*(2.-om) * data[ri->start + ri->diagonal].value;
363 
364  // Backward
365  typename std::vector<RowInfo>::const_reverse_iterator rri;
366  const typename std::vector<RowInfo>::const_reverse_iterator
367  rend = row_info.rend();
368  dst_ptr = dst.begin()+dst.size()-1;
369  for (rri = row_info.rbegin(); rri != rend; --dst_ptr, ++rri)
370  {
371  const size_type end_row = rri->start + rri->length;
372  number s = 0;
373  for (size_type i=rri->start+rri->diagonal+1; i<end_row; ++i)
374  s += data[i].value * dst(data[i].column);
375 
376  *dst_ptr -= s * om;
377  *dst_ptr /= data[rri->start + rri->diagonal].value;
378  }
379 }
380 
381 
382 
383 template <typename number>
384 std::size_t
386 {
387  return
388  sizeof (*this)
389  + sizeof(size_type) * row_info.capacity()
390  + sizeof(typename SparseMatrixEZ<number>::Entry) * data.capacity();
391 }
392 
393 
394 
395 template <typename number>
398 {
399  return row_info[row].length;
400 }
401 
402 
403 
404 template <typename number>
407 {
408  typename std::vector<RowInfo>::const_iterator row = row_info.begin();
409  const typename std::vector<RowInfo>::const_iterator endrow = row_info.end();
410 
411  // Add up entries actually used
412  size_type used = 0;
413  for (; row != endrow ; ++ row)
414  used += row->length;
415  return used;
416 }
417 
418 
419 
420 template <typename number>
421 void
423  size_type &used,
424  size_type &allocated,
425  size_type &reserved,
426  std::vector<size_type> &used_by_line,
427  const bool full) const
428 {
429  typename std::vector<RowInfo>::const_iterator row = row_info.begin();
430  const typename std::vector<RowInfo>::const_iterator endrow = row_info.end();
431 
432  // Add up entries actually used
433  used = 0;
434  size_type max_length = 0;
435  for (; row != endrow ; ++ row)
436  {
437  used += row->length;
438  if (max_length < row->length)
439  max_length = row->length;
440  }
441 
442  // Number of entries allocated is
443  // position of last entry used
444  --row;
445  allocated = row->start + row->length;
446  reserved = data.capacity();
447 
448 
449  if (full)
450  {
451  used_by_line.resize(max_length+1);
452 
453  for (row = row_info.begin() ; row != endrow; ++row)
454  {
455  ++used_by_line[row->length];
456  }
457  }
458 }
459 
460 
461 template <typename number>
462 void
463 SparseMatrixEZ<number>::print (std::ostream &out) const
464 {
465  AssertThrow (out, ExcIO());
466 
467  const_iterator i = begin();
468  const const_iterator e = end();
469  while (i != e)
470  {
471  out << i->row() << '\t'
472  << i->column() << '\t'
473  <<i->value() << std::endl;
474  ++i;
475  }
476 }
477 
478 
479 template <typename number>
480 void
482  const unsigned int precision,
483  const bool scientific,
484  const unsigned int width_,
485  const char *zero_string,
486  const double denominator) const
487 {
488  AssertThrow (out, ExcIO());
489  Assert (m() != 0, ExcNotInitialized());
490  Assert (n() != 0, ExcNotInitialized());
491 
492  unsigned int width = width_;
493 
494  std::ios::fmtflags old_flags = out.flags();
495  unsigned int old_precision = out.precision (precision);
496 
497  if (scientific)
498  {
499  out.setf (std::ios::scientific, std::ios::floatfield);
500  if (!width)
501  width = precision+7;
502  }
503  else
504  {
505  out.setf (std::ios::fixed, std::ios::floatfield);
506  if (!width)
507  width = precision+2;
508  }
509 
510  // TODO: Skip nonexisting entries
511  for (size_type i=0; i<m(); ++i)
512  {
513  for (size_type j=0; j<n(); ++j)
514  {
515  const Entry *entry = locate(i,j);
516  if (entry)
517  out << std::setw(width)
518  << entry->value *denominator << ' ';
519  else
520  out << std::setw(width) << zero_string << ' ';
521  }
522  out << std::endl;
523  };
524 
525  // reset output format
526  out.precision(old_precision);
527  out.flags (old_flags);
528 }
529 
530 
531 template <typename number>
532 void
533 SparseMatrixEZ<number>::block_write (std::ostream &out) const
534 {
535  AssertThrow (out, ExcIO());
536 
537  // first the simple objects,
538  // bracketed in [...]
539  out << '[' << row_info.size() << "]["
540  << n_columns << "]["
541  << data.size() << "]["
542  << increment << "][";
543  // then write out real data
544  typename std::vector<RowInfo>::const_iterator r = row_info.begin();
545  out.write(reinterpret_cast<const char *>(&*r),
546  sizeof(RowInfo) * row_info.size());
547 
548 // Just in case that vector entries are not stored consecutively
549 // const typename std::vector<RowInfo>::const_iterator re = row_info.end();
550 
551 // while (r != re)
552 // {
553 // out.write(reinterpret_cast<const char*>(&*r),
554 // sizeof(RowInfo));
555 // ++r;
556 // }
557 
558  out << "][";
559 
560  typename std::vector<Entry>::const_iterator d = data.begin();
561  out.write(reinterpret_cast<const char *>(&*d),
562  sizeof(Entry) * data.size());
563 
564 // Just in case that vector entries are not stored consecutively
565 // const typename std::vector<Entry>::const_iterator de = data.end();
566 
567 // while (d != de)
568 // {
569 // out.write(reinterpret_cast<const char*>(&*d),
570 // sizeof(Entry));
571 // ++d;
572 // }
573 
574  out << ']';
575 
576  AssertThrow (out, ExcIO());
577 }
578 
579 
580 #define DEAL_II_CHECK_INPUT(in,a,c) \
581  {in >> c; AssertThrow(c == a, \
582  ExcMessage("Unexpected character in input stream"));}
583 
584 template <typename number>
585 void
587 {
588  AssertThrow (in, ExcIO());
589 
590  char c;
591  int n;
592  // first read in simple data
593  DEAL_II_CHECK_INPUT(in,'[',c);
594  in >> n;
595  row_info.resize(n);
596 
597  DEAL_II_CHECK_INPUT(in,']',c);
598  DEAL_II_CHECK_INPUT(in,'[',c);
599  in >> n_columns;
600 
601  DEAL_II_CHECK_INPUT(in,']',c);
602  DEAL_II_CHECK_INPUT(in,'[',c);
603  in >> n;
604  data.resize(n);
605 
606  DEAL_II_CHECK_INPUT(in,']',c);
607  DEAL_II_CHECK_INPUT(in,'[',c);
608  in >> increment;
609 
610  DEAL_II_CHECK_INPUT(in,']',c);
611  DEAL_II_CHECK_INPUT(in,'[',c);
612 
613  // then read data
614  in.read(reinterpret_cast<char *>(&row_info[0]),
615  sizeof(RowInfo) * row_info.size());
616 
617  DEAL_II_CHECK_INPUT(in,']',c);
618  DEAL_II_CHECK_INPUT(in,'[',c);
619 
620  in.read(reinterpret_cast<char *>(&data[0]),
621  sizeof(Entry) * data.size());
622 
623  DEAL_II_CHECK_INPUT(in,']',c);
624 }
625 
626 #undef DEAL_II_CHECK_INPUT
627 
628 
629 
630 DEAL_II_NAMESPACE_CLOSE
631 
632 #endif // __deal2__sparse_matrix_ez_templates_h
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
std::size_t memory_consumption() const
SparseMatrixEZ< number > & operator=(const SparseMatrixEZ< number > &)
void vmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
void vmult_add(Vector< somenumber > &dst, const Vector< somenumber > &src) const
void block_read(std::istream &in)
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
::ExceptionBase & ExcIO()
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
#define Assert(cond, exc)
Definition: exceptions.h:299
size_type n() const
void reinit(const size_type n_rows, const size_type n_columns, size_type default_row_length=0, unsigned int default_increment=1, size_type reserve=0)
void compute_statistics(size_type &used, size_type &allocated, size_type &reserved, std::vector< size_type > &used_by_line, const bool compute_by_line) const
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1., const std::vector< std::size_t > &pos_right_of_diagonal=std::vector< std::size_t >()) const
iterator begin()
::ExceptionBase & ExcInvalidConstructorCall()
void Tvmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
types::global_dof_index size_type
size_type get_row_length(const size_type row) const
size_type n_nonzero_elements() const
std::size_t size() const
void print_formatted(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const unsigned int width=0, const char *zero_string=" ", const double denominator=1.) const
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcNotInitialized()
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void Tvmult_add(Vector< somenumber > &dst, const Vector< somenumber > &src) const
::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
void block_write(std::ostream &out) const
size_type m() const
::ExceptionBase & ExcInternalError()
void print(std::ostream &out) const