17 #ifndef __deal2__sparse_matrix_ez_templates_h
18 #define __deal2__sparse_matrix_ez_templates_h
21 #include <deal.II/lac/sparse_matrix_ez.h>
22 #include <deal.II/lac/vector.h>
33 template <
typename number>
40 template <
typename number>
51 template <
typename number>
55 const unsigned int default_increment)
57 reinit(n_rows, n_cols, default_row_length, default_increment);
61 template <
typename number>
66 template <
typename number>
75 template <
typename number>
81 typename std::vector<Entry>::iterator e = data.begin();
82 const typename std::vector<Entry>::iterator end = data.end();
94 template <
typename number>
99 unsigned int default_increment,
104 increment = default_increment;
107 row_info.resize(n_rows);
109 data.reserve(reserve);
110 data.resize(default_row_length * n_rows);
113 row_info[i].start = i * default_row_length;
117 template <
typename number>
127 template <
typename number>
131 return ((n_columns == 0) && (row_info.size()==0));
135 template <
typename number>
136 template <
typename somenumber>
144 const size_type end_row = row_info.size();
145 for (
size_type row = 0; row < end_row; ++row)
147 const RowInfo &ri = row_info[row];
148 typename std::vector<Entry>::const_iterator
149 entry = data.begin() + ri.
start;
151 for (
unsigned short i=0; i<ri.
length; ++i,++entry)
155 s += entry->value * src(entry->column);
162 template <
typename number>
170 while (start !=
final)
172 const double value = start->
value();
176 return std::sqrt(sum);
181 template <
typename number>
182 template <
typename somenumber>
188 Tvmult_add(dst, src);
192 template <
typename number>
193 template <
typename somenumber>
201 const size_type end_row = row_info.size();
202 for (
size_type row = 0; row < end_row; ++row)
204 const RowInfo &ri = row_info[row];
205 typename std::vector<Entry>::const_iterator
206 entry = data.begin() + ri.
start;
208 for (
unsigned short i=0; i<ri.
length; ++i,++entry)
212 s += entry->value * src(entry->column);
218 template <
typename number>
219 template <
typename somenumber>
227 const size_type end_row = row_info.size();
228 for (
size_type row = 0; row < end_row; ++row)
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)
237 dst(entry->column) += entry->value * src(row);
243 template <
typename number>
244 template <
typename somenumber>
248 const number om)
const
250 Assert (m() == n(), ExcNotQuadratic());
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();
259 for (; ri != end; ++dst_ptr, ++src_ptr, ++ri)
261 Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
262 *dst_ptr = om **src_ptr / data[ri->start + ri->diagonal].value;
268 template <
typename number>
269 template <
typename somenumber>
273 const number om)
const
275 Assert (m() == n(), ExcNotQuadratic());
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();
284 for (; ri != end; ++dst_ptr, ++src_ptr, ++ri)
286 Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
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);
292 *dst_ptr = om * s / data[ri->start + ri->diagonal].value;
298 template <
typename number>
299 template <
typename somenumber>
303 const number om)
const
305 Assert (m() == n(), ExcNotQuadratic());
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();
316 for (; ri != end; --dst_ptr, --src_ptr, ++ri)
318 Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
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);
324 *dst_ptr = om * s / data[ri->start + ri->diagonal].value;
330 template <
typename number>
331 template <
typename somenumber>
336 const std::vector<std::size_t> &)
const
338 Assert (m() == n(), ExcNotQuadratic());
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();
348 for (ri = row_info.begin(); ri != end; ++dst_ptr, ++src_ptr, ++ri)
350 Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
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);
356 *dst_ptr = *src_ptr - s * om;
357 *dst_ptr /= data[ri->start + ri->diagonal].value;
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;
365 typename std::vector<RowInfo>::const_reverse_iterator rri;
366 const typename std::vector<RowInfo>::const_reverse_iterator
367 rend = row_info.rend();
369 for (rri = row_info.rbegin(); rri != rend; --dst_ptr, ++rri)
371 const size_type end_row = rri->start + rri->length;
373 for (
size_type i=rri->start+rri->diagonal+1; i<end_row; ++i)
374 s += data[i].value * dst(data[i].column);
377 *dst_ptr /= data[rri->start + rri->diagonal].value;
383 template <
typename number>
389 +
sizeof(
size_type) * row_info.capacity()
395 template <
typename number>
399 return row_info[row].length;
404 template <
typename number>
408 typename std::vector<RowInfo>::const_iterator row = row_info.begin();
409 const typename std::vector<RowInfo>::const_iterator endrow = row_info.end();
413 for (; row != endrow ; ++ row)
420 template <
typename number>
426 std::vector<size_type> &used_by_line,
427 const bool full)
const
429 typename std::vector<RowInfo>::const_iterator row = row_info.begin();
430 const typename std::vector<RowInfo>::const_iterator endrow = row_info.end();
435 for (; row != endrow ; ++ row)
438 if (max_length < row->length)
439 max_length = row->length;
445 allocated = row->start + row->length;
446 reserved = data.capacity();
451 used_by_line.resize(max_length+1);
453 for (row = row_info.begin() ; row != endrow; ++row)
455 ++used_by_line[row->length];
461 template <
typename number>
471 out << i->
row() <<
'\t'
473 <<i->
value() << std::endl;
479 template <
typename number>
482 const unsigned int precision,
483 const bool scientific,
484 const unsigned int width_,
485 const char *zero_string,
486 const double denominator)
const
492 unsigned int width = width_;
494 std::ios::fmtflags old_flags = out.flags();
495 unsigned int old_precision = out.precision (precision);
499 out.setf (std::ios::scientific, std::ios::floatfield);
505 out.setf (std::ios::fixed, std::ios::floatfield);
515 const Entry *entry = locate(i,j);
517 out << std::setw(width)
518 << entry->
value *denominator <<
' ';
520 out << std::setw(width) << zero_string <<
' ';
526 out.precision(old_precision);
527 out.flags (old_flags);
531 template <
typename number>
539 out <<
'[' << row_info.size() <<
"]["
541 << data.size() <<
"]["
542 << increment <<
"][";
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());
560 typename std::vector<Entry>::const_iterator d = data.begin();
561 out.write(reinterpret_cast<const char *>(&*d),
562 sizeof(
Entry) * data.size());
580 #define DEAL_II_CHECK_INPUT(in,a,c) \
581 {in >> c; AssertThrow(c == a, \
582 ExcMessage("Unexpected character in input stream"));}
584 template <
typename number>
593 DEAL_II_CHECK_INPUT(in,
'[',c);
597 DEAL_II_CHECK_INPUT(in,
']',c);
598 DEAL_II_CHECK_INPUT(in,
'[',c);
601 DEAL_II_CHECK_INPUT(in,
']',c);
602 DEAL_II_CHECK_INPUT(in,
'[',c);
606 DEAL_II_CHECK_INPUT(in,
']',c);
607 DEAL_II_CHECK_INPUT(in,
'[',c);
610 DEAL_II_CHECK_INPUT(in,
']',c);
611 DEAL_II_CHECK_INPUT(in,
'[',c);
614 in.read(reinterpret_cast<char *>(&row_info[0]),
615 sizeof(
RowInfo) * row_info.size());
617 DEAL_II_CHECK_INPUT(in,
']',c);
618 DEAL_II_CHECK_INPUT(in,
'[',c);
620 in.read(reinterpret_cast<char *>(&data[0]),
621 sizeof(
Entry) * data.size());
623 DEAL_II_CHECK_INPUT(in,
']',c);
626 #undef DEAL_II_CHECK_INPUT
630 DEAL_II_NAMESPACE_CLOSE
632 #endif // __deal2__sparse_matrix_ez_templates_h
void block_read(std::istream &in)
types::global_dof_index size_type
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
void vmult_add(Vector< somenumber > &dst, const Vector< somenumber > &src) const
#define AssertThrow(cond, exc)
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
::ExceptionBase & ExcIO()
Iterator is invalid, probably due to an error.
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
#define Assert(cond, exc)
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
void Tvmult_add(Vector< somenumber > &dst, const Vector< somenumber > &src) const
::ExceptionBase & ExcInvalidConstructorCall()
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 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
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
void vmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
size_type get_row_length(const size_type row) const
void Tvmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
void print(std::ostream &out) const
SparseMatrixEZ< number > & operator=(const SparseMatrixEZ< number > &)
size_type n_nonzero_elements() const
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcNotInitialized()
std::size_t memory_consumption() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
::ExceptionBase & ExcInternalError()
void block_write(std::ostream &out) const