00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__sparse_matrix_ez_h
00013 #define __deal2__sparse_matrix_ez_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/subscriptor.h>
00018 #include <deal.II/base/smartpointer.h>
00019 #include <deal.II/lac/exceptions.h>
00020
00021 #include <vector>
00022
00023 DEAL_II_NAMESPACE_OPEN
00024
00025 template<typename number> class Vector;
00026 template<typename number> class FullMatrix;
00027
00095 template <typename number>
00096 class SparseMatrixEZ : public Subscriptor
00097 {
00098 public:
00104 struct Entry
00105 {
00111 Entry();
00112
00117 Entry (const unsigned int column,
00118 const number& value);
00119
00123 unsigned int column;
00124
00128 number value;
00129
00133 static const unsigned int invalid = numbers::invalid_unsigned_int;
00134 };
00135
00142 struct RowInfo
00143 {
00147 RowInfo (unsigned int start = Entry::invalid);
00148
00153 unsigned int start;
00158 unsigned short length;
00164 unsigned short diagonal;
00168 static const unsigned short
00169 invalid_diagonal = static_cast<unsigned short>(-1);
00170 };
00171
00172 public:
00173
00177 class const_iterator
00178 {
00179 private:
00183 class Accessor
00184 {
00185 public:
00192 Accessor (const SparseMatrixEZ<number> *matrix,
00193 const unsigned int row,
00194 const unsigned short index);
00195
00201 unsigned int row() const;
00202
00208 unsigned short index() const;
00209
00215 unsigned int column() const;
00216
00220 number value() const;
00221
00222 protected:
00226 const SparseMatrixEZ<number>* matrix;
00227
00231 unsigned int a_row;
00232
00236 unsigned short a_index;
00237
00242 friend class const_iterator;
00243 };
00244
00245 public:
00249 const_iterator(const SparseMatrixEZ<number> *matrix,
00250 const unsigned int row,
00251 const unsigned short index);
00252
00258 const_iterator& operator++ ();
00259
00265 const_iterator& operator++ (int);
00266
00270 const Accessor& operator* () const;
00271
00275 const Accessor* operator-> () const;
00276
00283 bool operator == (const const_iterator&) const;
00287 bool operator != (const const_iterator&) const;
00288
00298 bool operator < (const const_iterator&) const;
00299
00300 private:
00305 Accessor accessor;
00306
00322 #ifdef DEAL_II_NESTED_CLASS_FRIEND_BUG
00323 template <typename> friend class SparseMatrixEZ;
00324 #endif
00325 };
00326
00331 typedef number value_type;
00332
00342 SparseMatrixEZ ();
00343
00354 SparseMatrixEZ (const SparseMatrixEZ &);
00355
00369 explicit SparseMatrixEZ (const unsigned int n_rows,
00370 const unsigned int n_columns,
00371 const unsigned int default_row_length = 0,
00372 const unsigned int default_increment = 1);
00373
00379 ~SparseMatrixEZ ();
00380
00385 SparseMatrixEZ<number>& operator = (const SparseMatrixEZ<number> &);
00386
00403 SparseMatrixEZ<number>& operator = (const double d);
00404
00420 void reinit (const unsigned int n_rows,
00421 const unsigned int n_columns,
00422 unsigned int default_row_length = 0,
00423 unsigned int default_increment = 1,
00424 unsigned int reserve = 0);
00425
00433 void clear ();
00435
00444 bool empty () const;
00445
00452 unsigned int m () const;
00453
00460 unsigned int n () const;
00461
00466 unsigned int get_row_length (const unsigned int row) const;
00467
00472 unsigned int n_nonzero_elements () const;
00473
00479 std::size_t memory_consumption () const;
00480
00490 template <class STREAM>
00491 void print_statistics (STREAM& s, bool full = false);
00492
00506 void compute_statistics (unsigned int& used,
00507 unsigned int& allocated,
00508 unsigned int& reserved,
00509 std::vector<unsigned int>& used_by_line,
00510 const bool compute_by_line) const;
00512
00525 void set (const unsigned int i, const unsigned int j,
00526 const number value);
00527
00537 void add (const unsigned int i,
00538 const unsigned int j,
00539 const number value);
00540
00568 template <typename number2>
00569 void add (const std::vector<unsigned int> &indices,
00570 const FullMatrix<number2> &full_matrix,
00571 const bool elide_zero_values = true);
00572
00580 template <typename number2>
00581 void add (const std::vector<unsigned int> &row_indices,
00582 const std::vector<unsigned int> &col_indices,
00583 const FullMatrix<number2> &full_matrix,
00584 const bool elide_zero_values = true);
00585
00603 template <typename number2>
00604 void add (const unsigned int row,
00605 const std::vector<unsigned int> &col_indices,
00606 const std::vector<number2> &values,
00607 const bool elide_zero_values = true);
00608
00626 template <typename number2>
00627 void add (const unsigned int row,
00628 const unsigned int n_cols,
00629 const unsigned int *col_indices,
00630 const number2 *values,
00631 const bool elide_zero_values = true,
00632 const bool col_indices_are_sorted = false);
00633
00660 template <class MATRIX>
00661 SparseMatrixEZ<number> &
00662 copy_from (const MATRIX &source);
00663
00675 template <class MATRIX>
00676 void add (const number factor,
00677 const MATRIX &matrix);
00679
00701 number operator () (const unsigned int i,
00702 const unsigned int j) const;
00703
00709 number el (const unsigned int i,
00710 const unsigned int j) const;
00712
00721 template <typename somenumber>
00722 void vmult (Vector<somenumber> &dst,
00723 const Vector<somenumber> &src) const;
00724
00733 template <typename somenumber>
00734 void Tvmult (Vector<somenumber> &dst,
00735 const Vector<somenumber> &src) const;
00736
00743 template <typename somenumber>
00744 void vmult_add (Vector<somenumber> &dst,
00745 const Vector<somenumber> &src) const;
00746
00755 template <typename somenumber>
00756 void Tvmult_add (Vector<somenumber> &dst,
00757 const Vector<somenumber> &src) const;
00759
00766 number l2_norm () const;
00768
00782 template <typename somenumber>
00783 void precondition_Jacobi (Vector<somenumber> &dst,
00784 const Vector<somenumber> &src,
00785 const number omega = 1.) const;
00786
00791 template <typename somenumber>
00792 void precondition_SSOR (Vector<somenumber> &dst,
00793 const Vector<somenumber> &src,
00794 const number om = 1.,
00795 const std::vector<unsigned int> &pos_right_of_diagonal = std::vector<unsigned int>()) const;
00796
00802 template <typename somenumber>
00803 void precondition_SOR (Vector<somenumber> &dst,
00804 const Vector<somenumber> &src,
00805 const number om = 1.) const;
00806
00812 template <typename somenumber>
00813 void precondition_TSOR (Vector<somenumber> &dst,
00814 const Vector<somenumber> &src,
00815 const number om = 1.) const;
00816
00831 template <class MATRIXA, class MATRIXB>
00832 void conjugate_add (const MATRIXA& A,
00833 const MATRIXB& B,
00834 const bool transpose = false);
00836
00844 const_iterator begin () const;
00845
00849 const_iterator end () const;
00850
00858 const_iterator begin (const unsigned int r) const;
00859
00865 const_iterator end (const unsigned int r) const;
00867
00878 void print (std::ostream &out) const;
00879
00919 void print_formatted (std::ostream &out,
00920 const unsigned int precision = 3,
00921 const bool scientific = true,
00922 const unsigned int width = 0,
00923 const char *zero_string = " ",
00924 const double denominator = 1.) const;
00925
00933 void block_write (std::ostream &out) const;
00934
00952 void block_read (std::istream &in);
00954
00961 DeclException0(ExcNoDiagonal);
00962
00966 DeclException2 (ExcInvalidEntry,
00967 int, int,
00968 << "The entry with index (" << arg1 << ',' << arg2
00969 << ") does not exist.");
00970
00971 DeclException2(ExcEntryAllocationFailure,
00972 int, int,
00973 << "An entry with index (" << arg1 << ',' << arg2
00974 << ") cannot be allocated.");
00976 private:
00983 const Entry* locate (const unsigned int row,
00984 const unsigned int col) const;
00985
00992 Entry* locate (const unsigned int row,
00993 const unsigned int col);
00994
00998 Entry* allocate (const unsigned int row,
00999 const unsigned int col);
01000
01010 template <typename somenumber>
01011 void threaded_vmult (Vector<somenumber> &dst,
01012 const Vector<somenumber> &src,
01013 const unsigned int begin_row,
01014 const unsigned int end_row) const;
01015
01027 template <typename somenumber>
01028 void threaded_matrix_norm_square (const Vector<somenumber> &v,
01029 const unsigned int begin_row,
01030 const unsigned int end_row,
01031 somenumber *partial_sum) const;
01032
01044 template <typename somenumber>
01045 void threaded_matrix_scalar_product (const Vector<somenumber> &u,
01046 const Vector<somenumber> &v,
01047 const unsigned int begin_row,
01048 const unsigned int end_row,
01049 somenumber *partial_sum) const;
01050
01056 unsigned int n_columns;
01057
01061 std::vector<RowInfo> row_info;
01062
01066 std::vector<Entry> data;
01067
01071 unsigned int increment;
01072
01080 #ifdef DEAL_II_NESTED_CLASS_FRIEND_BUG
01081 friend class const_iterator;
01082 friend class const_iterator::Accessor;
01083 #endif
01084 };
01085
01087
01088
01089 template <typename number>
01090 inline
01091 SparseMatrixEZ<number>::Entry::Entry(unsigned int column,
01092 const number& value)
01093 :
01094 column(column),
01095 value(value)
01096 {}
01097
01098
01099
01100 template <typename number>
01101 inline
01102 SparseMatrixEZ<number>::Entry::Entry()
01103 :
01104 column(invalid),
01105 value(0)
01106 {}
01107
01108
01109 template <typename number>
01110 inline
01111 SparseMatrixEZ<number>::RowInfo::RowInfo(unsigned int start)
01112 :
01113 start(start),
01114 length(0),
01115 diagonal(invalid_diagonal)
01116 {}
01117
01118
01119
01120 template <typename number>
01121 inline
01122 SparseMatrixEZ<number>::const_iterator::Accessor::
01123 Accessor (const SparseMatrixEZ<number> *matrix,
01124 const unsigned int r,
01125 const unsigned short i)
01126 :
01127 matrix(matrix),
01128 a_row(r),
01129 a_index(i)
01130 {}
01131
01132
01133 template <typename number>
01134 inline
01135 unsigned int
01136 SparseMatrixEZ<number>::const_iterator::Accessor::row() const
01137 {
01138 return a_row;
01139 }
01140
01141
01142 template <typename number>
01143 inline
01144 unsigned int
01145 SparseMatrixEZ<number>::const_iterator::Accessor::column() const
01146 {
01147 return matrix->data[matrix->row_info[a_row].start+a_index].column;
01148 }
01149
01150
01151 template <typename number>
01152 inline
01153 unsigned short
01154 SparseMatrixEZ<number>::const_iterator::Accessor::index() const
01155 {
01156 return a_index;
01157 }
01158
01159
01160
01161 template <typename number>
01162 inline
01163 number
01164 SparseMatrixEZ<number>::const_iterator::Accessor::value() const
01165 {
01166 return matrix->data[matrix->row_info[a_row].start+a_index].value;
01167 }
01168
01169
01170 template <typename number>
01171 inline
01172 SparseMatrixEZ<number>::const_iterator::
01173 const_iterator(const SparseMatrixEZ<number> *matrix,
01174 const unsigned int r,
01175 const unsigned short i)
01176 :
01177 accessor(matrix, r, i)
01178 {
01179
01180 if (r==accessor.matrix->m() && i==0) return;
01181
01182
01183
01184
01185
01186
01187
01188
01189 if (accessor.a_index >= accessor.matrix->row_info[accessor.a_row].length)
01190 {
01191 do
01192 {
01193 ++accessor.a_row;
01194 }
01195
01196
01197
01198
01199 while (accessor.a_row < accessor.matrix->m()
01200 && accessor.matrix->row_info[accessor.a_row].length == 0);
01201 }
01202 }
01203
01204
01205 template <typename number>
01206 inline
01207 typename SparseMatrixEZ<number>::const_iterator&
01208 SparseMatrixEZ<number>::const_iterator::operator++ ()
01209 {
01210 Assert (accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
01211
01212
01213 ++(accessor.a_index);
01214
01215
01216
01217 if (accessor.a_index >= accessor.matrix->row_info[accessor.a_row].length)
01218 {
01219 accessor.a_index = 0;
01220
01221
01222 do
01223 {
01224 ++accessor.a_row;
01225 }
01226 while (accessor.a_row < accessor.matrix->m()
01227 && accessor.matrix->row_info[accessor.a_row].length == 0);
01228 }
01229 return *this;
01230 }
01231
01232
01233 template <typename number>
01234 inline
01235 const typename SparseMatrixEZ<number>::const_iterator::Accessor&
01236 SparseMatrixEZ<number>::const_iterator::operator* () const
01237 {
01238 return accessor;
01239 }
01240
01241
01242 template <typename number>
01243 inline
01244 const typename SparseMatrixEZ<number>::const_iterator::Accessor*
01245 SparseMatrixEZ<number>::const_iterator::operator-> () const
01246 {
01247 return &accessor;
01248 }
01249
01250
01251 template <typename number>
01252 inline
01253 bool
01254 SparseMatrixEZ<number>::const_iterator::operator == (
01255 const const_iterator& other) const
01256 {
01257 return (accessor.row() == other.accessor.row() &&
01258 accessor.index() == other.accessor.index());
01259 }
01260
01261
01262 template <typename number>
01263 inline
01264 bool
01265 SparseMatrixEZ<number>::const_iterator::
01266 operator != (const const_iterator& other) const
01267 {
01268 return ! (*this == other);
01269 }
01270
01271
01272 template <typename number>
01273 inline
01274 bool
01275 SparseMatrixEZ<number>::const_iterator::
01276 operator < (const const_iterator& other) const
01277 {
01278 return (accessor.row() < other.accessor.row() ||
01279 (accessor.row() == other.accessor.row() &&
01280 accessor.index() < other.accessor.index()));
01281 }
01282
01283
01284
01285 template <typename number>
01286 inline
01287 unsigned int SparseMatrixEZ<number>::m () const
01288 {
01289 return row_info.size();
01290 }
01291
01292
01293 template <typename number>
01294 inline
01295 unsigned int SparseMatrixEZ<number>::n () const
01296 {
01297 return n_columns;
01298 }
01299
01300
01301 template <typename number>
01302 inline
01303 typename SparseMatrixEZ<number>::Entry*
01304 SparseMatrixEZ<number>::locate (const unsigned int row,
01305 const unsigned int col)
01306 {
01307 Assert (row<m(), ExcIndexRange(row,0,m()));
01308 Assert (col<n(), ExcIndexRange(col,0,n()));
01309
01310 const RowInfo& r = row_info[row];
01311 const unsigned int end = r.start + r.length;
01312 for (unsigned int i=r.start;i<end;++i)
01313 {
01314 Entry * const entry = &data[i];
01315 if (entry->column == col)
01316 return entry;
01317 if (entry->column == Entry::invalid)
01318 return 0;
01319 }
01320 return 0;
01321 }
01322
01323
01324
01325 template <typename number>
01326 inline
01327 const typename SparseMatrixEZ<number>::Entry*
01328 SparseMatrixEZ<number>::locate (const unsigned int row,
01329 const unsigned int col) const
01330 {
01331 SparseMatrixEZ<number>* t = const_cast<SparseMatrixEZ<number>*> (this);
01332 return t->locate(row,col);
01333 }
01334
01335
01336 template <typename number>
01337 inline
01338 typename SparseMatrixEZ<number>::Entry*
01339 SparseMatrixEZ<number>::allocate (const unsigned int row,
01340 const unsigned int col)
01341 {
01342 Assert (row<m(), ExcIndexRange(row,0,m()));
01343 Assert (col<n(), ExcIndexRange(col,0,n()));
01344
01345 RowInfo& r = row_info[row];
01346 const unsigned int end = r.start + r.length;
01347
01348 unsigned int i = r.start;
01349
01350
01351
01352 if (r.diagonal != RowInfo::invalid_diagonal && col >= row)
01353 i += r.diagonal;
01354
01355 while (i<end && data[i].column < col) ++i;
01356
01357
01358 if (i != end && data[i].column == col)
01359 return &data[i];
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369 if (row != row_info.size()-1)
01370 {
01371 if (end >= row_info[row+1].start)
01372 {
01373
01374 Assert(increment!=0,ExcEntryAllocationFailure(row,col));
01375
01376
01377 data.insert(data.begin()+end, increment, Entry());
01378
01379
01380 for (unsigned int rn=row+1;rn<row_info.size();++rn)
01381 row_info[rn].start += increment;
01382 }
01383 }
01384 else
01385 {
01386 if (end >= data.size())
01387 {
01388
01389
01390
01391 data.push_back(Entry());
01392 }
01393 }
01394
01395 Entry* entry = &data[i];
01396
01397 Entry temp = *entry;
01398
01399
01400
01401
01402 entry->column = col;
01403 entry->value = 0;
01404
01405 ++r.length;
01406 if (col == row)
01407 r.diagonal = i - r.start;
01408 else
01409 if (col<row && r.diagonal!= RowInfo::invalid_diagonal)
01410 ++r.diagonal;
01411
01412 if (i == end)
01413 return entry;
01414
01415
01416
01417 for (unsigned int j = i+1; j < end; ++j)
01418 {
01419
01420
01421 Assert (data[j].column != Entry::invalid, ExcInternalError());
01422
01423
01424 std::swap (data[j], temp);
01425 }
01426 Assert (data[end].column == Entry::invalid, ExcInternalError());
01427
01428 data[end] = temp;
01429
01430 return entry;
01431 }
01432
01433
01434
01435 template <typename number>
01436 inline
01437 void SparseMatrixEZ<number>::set (const unsigned int i,
01438 const unsigned int j,
01439 const number value)
01440 {
01441
01442 Assert (numbers::is_finite(value), ExcNumberNotFinite());
01443
01444 Assert (i<m(), ExcIndexRange(i,0,m()));
01445 Assert (j<n(), ExcIndexRange(j,0,n()));
01446
01447 if (value == 0.)
01448 {
01449 Entry* entry = locate(i,j);
01450 if (entry != 0)
01451 entry->value = 0.;
01452 }
01453 else
01454 {
01455 Entry* entry = allocate(i,j);
01456 entry->value = value;
01457 }
01458 }
01459
01460
01461
01462 template <typename number>
01463 inline
01464 void SparseMatrixEZ<number>::add (const unsigned int i,
01465 const unsigned int j,
01466 const number value)
01467 {
01468
01469 Assert (numbers::is_finite(value), ExcNumberNotFinite());
01470
01471 Assert (i<m(), ExcIndexRange(i,0,m()));
01472 Assert (j<n(), ExcIndexRange(j,0,n()));
01473
01474
01475 if (value == 0)
01476 return;
01477
01478 Entry* entry = allocate(i,j);
01479 entry->value += value;
01480 }
01481
01482
01483 template <typename number>
01484 template <typename number2>
01485 void SparseMatrixEZ<number>::add (const std::vector<unsigned int> &indices,
01486 const FullMatrix<number2> &full_matrix,
01487 const bool elide_zero_values)
01488 {
01489
01490 for (unsigned int i=0; i<indices.size(); ++i)
01491 for (unsigned int j=0; j<indices.size(); ++j)
01492 if ((full_matrix(i,j) != 0) || (elide_zero_values == false))
01493 add (indices[i], indices[j], full_matrix(i,j));
01494 }
01495
01496
01497
01498 template <typename number>
01499 template <typename number2>
01500 void SparseMatrixEZ<number>::add (const std::vector<unsigned int> &row_indices,
01501 const std::vector<unsigned int> &col_indices,
01502 const FullMatrix<number2> &full_matrix,
01503 const bool elide_zero_values)
01504 {
01505
01506 for (unsigned int i=0; i<row_indices.size(); ++i)
01507 for (unsigned int j=0; j<col_indices.size(); ++j)
01508 if ((full_matrix(i,j) != 0) || (elide_zero_values == false))
01509 add (row_indices[i], col_indices[j], full_matrix(i,j));
01510 }
01511
01512
01513
01514
01515 template <typename number>
01516 template <typename number2>
01517 void SparseMatrixEZ<number>::add (const unsigned int row,
01518 const std::vector<unsigned int> &col_indices,
01519 const std::vector<number2> &values,
01520 const bool elide_zero_values)
01521 {
01522
01523 for (unsigned int j=0; j<col_indices.size(); ++j)
01524 if ((values[j] != 0) || (elide_zero_values == false))
01525 add (row, col_indices[j], values[j]);
01526 }
01527
01528
01529
01530 template <typename number>
01531 template <typename number2>
01532 void SparseMatrixEZ<number>::add (const unsigned int row,
01533 const unsigned int n_cols,
01534 const unsigned int *col_indices,
01535 const number2 *values,
01536 const bool elide_zero_values,
01537 const bool )
01538 {
01539
01540 for (unsigned int j=0; j<n_cols; ++j)
01541 if ((values[j] != 0) || (elide_zero_values == false))
01542 add (row, col_indices[j], values[j]);
01543 }
01544
01545
01546
01547
01548 template <typename number>
01549 inline
01550 number SparseMatrixEZ<number>::el (const unsigned int i,
01551 const unsigned int j) const
01552 {
01553 const Entry* entry = locate(i,j);
01554 if (entry)
01555 return entry->value;
01556 return 0.;
01557 }
01558
01559
01560
01561 template <typename number>
01562 inline
01563 number SparseMatrixEZ<number>::operator() (const unsigned int i,
01564 const unsigned int j) const
01565 {
01566 const Entry* entry = locate(i,j);
01567 if (entry)
01568 return entry->value;
01569 Assert(false, ExcInvalidEntry(i,j));
01570 return 0.;
01571 }
01572
01573
01574 template <typename number>
01575 inline
01576 typename SparseMatrixEZ<number>::const_iterator
01577 SparseMatrixEZ<number>::begin () const
01578 {
01579 const_iterator result(this, 0, 0);
01580 return result;
01581 }
01582
01583 template <typename number>
01584 inline
01585 typename SparseMatrixEZ<number>::const_iterator
01586 SparseMatrixEZ<number>::end () const
01587 {
01588 return const_iterator(this, m(), 0);
01589 }
01590
01591 template <typename number>
01592 inline
01593 typename SparseMatrixEZ<number>::const_iterator
01594 SparseMatrixEZ<number>::begin (const unsigned int r) const
01595 {
01596 Assert (r<m(), ExcIndexRange(r,0,m()));
01597 const_iterator result (this, r, 0);
01598 return result;
01599 }
01600
01601 template <typename number>
01602 inline
01603 typename SparseMatrixEZ<number>::const_iterator
01604 SparseMatrixEZ<number>::end (const unsigned int r) const
01605 {
01606 Assert (r<m(), ExcIndexRange(r,0,m()));
01607 const_iterator result(this, r+1, 0);
01608 return result;
01609 }
01610
01611 template<typename number>
01612 template <class MATRIX>
01613 inline
01614 SparseMatrixEZ<number>&
01615 SparseMatrixEZ<number>::copy_from (const MATRIX& M)
01616 {
01617 reinit(M.m(), M.n());
01618
01619 typename MATRIX::const_iterator start = M.begin();
01620 const typename MATRIX::const_iterator final = M.end();
01621
01622 while (start != final)
01623 {
01624 if (start->value() != 0.)
01625 set(start->row(), start->column(), start->value());
01626 ++start;
01627 }
01628 return *this;
01629 }
01630
01631 template<typename number>
01632 template <class MATRIX>
01633 inline
01634 void
01635 SparseMatrixEZ<number>::add (const number factor,
01636 const MATRIX& M)
01637 {
01638 Assert (M.m() == m(), ExcDimensionMismatch(M.m(), m()));
01639 Assert (M.n() == n(), ExcDimensionMismatch(M.n(), n()));
01640
01641 if (factor == 0.)
01642 return;
01643
01644 typename MATRIX::const_iterator start = M.begin();
01645 const typename MATRIX::const_iterator final = M.end();
01646
01647 while (start != final)
01648 {
01649 if (start->value() != 0.)
01650 add(start->row(), start->column(), factor * start->value());
01651 ++start;
01652 }
01653 }
01654
01655
01656
01657 template<typename number>
01658 template <class MATRIXA, class MATRIXB>
01659 inline void
01660 SparseMatrixEZ<number>::conjugate_add (const MATRIXA& A,
01661 const MATRIXB& B,
01662 const bool transpose)
01663 {
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680 typename MATRIXB::const_iterator b1 = B.begin();
01681 const typename MATRIXB::const_iterator b_final = B.end();
01682 if (transpose)
01683 while (b1 != b_final)
01684 {
01685 const unsigned int i = b1->column();
01686 const unsigned int k = b1->row();
01687 typename MATRIXB::const_iterator b2 = B.begin();
01688 while (b2 != b_final)
01689 {
01690 const unsigned int j = b2->column();
01691 const unsigned int l = b2->row();
01692
01693 const typename MATRIXA::value_type a = A.el(k,l);
01694
01695 if (a != 0.)
01696 add (i, j, a * b1->value() * b2->value());
01697 ++b2;
01698 }
01699 ++b1;
01700 }
01701 else
01702 {
01703
01704
01705
01706
01707 std::vector<unsigned int> minrow(B.n(), B.m());
01708 std::vector<unsigned int> maxrow(B.n(), 0);
01709 while (b1 != b_final)
01710 {
01711 const unsigned int r = b1->row();
01712 if (r < minrow[b1->column()])
01713 minrow[b1->column()] = r;
01714 if (r > maxrow[b1->column()])
01715 maxrow[b1->column()] = r;
01716 ++b1;
01717 }
01718
01719 typename MATRIXA::const_iterator ai = A.begin();
01720 const typename MATRIXA::const_iterator ae = A.end();
01721
01722 while (ai != ae)
01723 {
01724 const typename MATRIXA::value_type a = ai->value();
01725
01726
01727 if (a == 0.) continue;
01728
01729
01730
01731
01732
01733 b1 = B.begin(minrow[ai->row()]);
01734 const typename MATRIXB::const_iterator
01735 be1 = B.end(maxrow[ai->row()]);
01736 const typename MATRIXB::const_iterator
01737 be2 = B.end(maxrow[ai->column()]);
01738
01739 while (b1 != be1)
01740 {
01741 const double b1v = b1->value();
01742
01743
01744
01745
01746 if (b1->column() == ai->row() && (b1v != 0.))
01747 {
01748 const unsigned int i = b1->row();
01749
01750 typename MATRIXB::const_iterator
01751 b2 = B.begin(minrow[ai->column()]);
01752 while (b2 != be2)
01753 {
01754 if (b2->column() == ai->column())
01755 {
01756 const unsigned int j = b2->row();
01757 add (i, j, a * b1v * b2->value());
01758 }
01759 ++b2;
01760 }
01761 }
01762 ++b1;
01763 }
01764 ++ai;
01765 }
01766 }
01767 }
01768
01769
01770 template <typename number>
01771 template <class STREAM>
01772 inline
01773 void
01774 SparseMatrixEZ<number>::print_statistics(STREAM& out, bool full)
01775 {
01776 unsigned int used;
01777 unsigned int allocated;
01778 unsigned int reserved;
01779 std::vector<unsigned int> used_by_line;
01780
01781 compute_statistics (used, allocated, reserved, used_by_line, full);
01782
01783 out << "SparseMatrixEZ:used entries:" << used << std::endl
01784 << "SparseMatrixEZ:allocated entries:" << allocated << std::endl
01785 << "SparseMatrixEZ:reserved entries:" << reserved << std::endl;
01786
01787 if (full)
01788 {
01789 for (unsigned int i=0; i< used_by_line.size();++i)
01790 if (used_by_line[i] != 0)
01791 out << "SparseMatrixEZ:entries\t" << i
01792 << "\trows\t" << used_by_line[i]
01793 << std::endl;
01794
01795 }
01796 }
01797
01798
01799 DEAL_II_NAMESPACE_CLOSE
01800
01801 #endif
01802
01803