Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
matrix_block.h
1 // ---------------------------------------------------------------------
2 // @f$Id: matrix_block.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 2007 - 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__matrix_block_h
18 #define __deal2__matrix_block_h
19 
20 #include <deal.II/base/config.h>
21 #include <deal.II/base/named_data.h>
22 #include <deal.II/base/smartpointer.h>
23 #include <deal.II/base/std_cxx1x/shared_ptr.h>
24 #include <deal.II/base/memory_consumption.h>
25 #include <deal.II/base/mg_level_object.h>
26 #include <deal.II/lac/block_indices.h>
27 #include <deal.II/lac/block_sparsity_pattern.h>
28 #include <deal.II/lac/sparse_matrix.h>
29 #include <deal.II/lac/full_matrix.h>
30 
31 
32 DEAL_II_NAMESPACE_OPEN
33 
34 template <class MATRIX> class MatrixBlock;
35 
36 namespace internal
37 {
38  template <class MATRIX>
39  void
40  reinit(MatrixBlock<MATRIX> &v, const BlockSparsityPattern &p);
41 
42  template <typename number>
43  void
45 }
46 
104 template <class MATRIX>
105 class MatrixBlock
106  : public Subscriptor
107 {
108 public:
113 
118  MatrixBlock();
119 
124 
132 
144  void reinit(const BlockSparsityPattern &sparsity);
145 
146  operator MATRIX &();
147  operator const MATRIX &() const;
148 
156  void add (const size_type i,
157  const size_type j,
158  const typename MATRIX::value_type value);
159 
190  template <typename number>
191  void add (const std::vector<size_type> &indices,
192  const FullMatrix<number> &full_matrix,
193  const bool elide_zero_values = true);
194 
224  template <typename number>
225  void add (const std::vector<size_type> &row_indices,
226  const std::vector<size_type> &col_indices,
227  const FullMatrix<number> &full_matrix,
228  const bool elide_zero_values = true);
229 
264  template <typename number>
265  void add (const size_type row_index,
266  const std::vector<size_type> &col_indices,
267  const std::vector<number> &values,
268  const bool elide_zero_values = true);
269 
287  template <typename number>
288  void add (const size_type row,
289  const size_type n_cols,
290  const size_type *col_indices,
291  const number *values,
292  const bool elide_zero_values = true,
293  const bool col_indices_are_sorted = false);
294 
303  template<class VECTOR>
304  void vmult (VECTOR &w, const VECTOR &v) const;
305 
314  template<class VECTOR>
315  void vmult_add (VECTOR &w, const VECTOR &v) const;
316 
325  template<class VECTOR>
326  void Tvmult (VECTOR &w, const VECTOR &v) const;
327 
336  template<class VECTOR>
337  void Tvmult_add (VECTOR &w, const VECTOR &v) const;
338 
342  std::size_t memory_consumption () const;
343 
351  DeclException2(ExcBlockIndexMismatch, size_type, size_type,
352  << "Block index " << arg1 << " does not match " << arg2);
353 
367 
372 private:
390 
391  friend void internal::reinit<>(MatrixBlock<MATRIX> &, const BlockSparsityPattern &);
392 };
393 
394 
404 template <class MATRIX>
406  :
407  private NamedData<std_cxx1x::shared_ptr<MatrixBlock<MATRIX> > >
408 {
409 public:
414 
419 
425  void add(size_type row, size_type column, const std::string &name);
426 
434  void reinit(const BlockSparsityPattern &sparsity);
435 
453  void clear (bool really_clean = false);
454 
458  std::size_t memory_consumption () const;
459 
464  const value_type &block(size_type i) const;
465 
471 
477  MATRIX &matrix(size_type i);
478 
486 };
487 
488 
498 template <class MATRIX>
500  : public Subscriptor
501 {
502 public:
507 
527  MGMatrixBlockVector(const bool edge_matrices = false,
528  const bool edge_flux_matrices = false);
529 
533  unsigned int size () const;
534 
543  void add(size_type row, size_type column, const std::string &name);
544 
567  void reinit_edge(const MGLevelObject<BlockSparsityPattern> &sparsity);
580 
598  void clear (bool really_clean = false);
599 
604  const value_type &block(size_type i) const;
605 
611 
616  const value_type &block_in(size_type i) const;
617 
623 
628  const value_type &block_out(size_type i) const;
629 
635 
640  const value_type &block_up(size_type i) const;
641 
647 
652  const value_type &block_down(size_type i) const;
653 
659 
663  std::size_t memory_consumption () const;
664 private:
667 
669  const bool edge_matrices;
670 
672  const bool edge_flux_matrices;
673 
684 };
685 
686 
687 //----------------------------------------------------------------------//
688 
689 namespace internal
690 {
691  template <class MATRIX>
692  void
693  reinit(MatrixBlock<MATRIX> &v, const BlockSparsityPattern &p)
694  {
697  }
698 
699 
700  template <typename number>
701  void
702  reinit(MatrixBlock<::SparseMatrix<number> > &v, const BlockSparsityPattern &p)
703  {
706  v.matrix.reinit(p.block(v.row, v.column));
707  }
708 }
709 
710 
711 template <class MATRIX>
712 inline
714  :
715  row(deal_II_numbers::invalid_size_type),
716  column(deal_II_numbers::invalid_size_type)
717 {}
718 
719 
720 template <class MATRIX>
721 inline
723  :
724  Subscriptor(),
725  row(M.row),
726  column(M.column),
727  matrix(M.matrix),
728  row_indices(M.row_indices),
729  column_indices(M.column_indices)
730 {}
731 
732 
733 template <class MATRIX>
734 inline
736  :
737  row(i), column(j)
738 {}
739 
740 
741 template <class MATRIX>
742 inline
743 void
745 {
746  internal::reinit(*this, sparsity);
747 }
748 
749 
750 template <class MATRIX>
751 inline
753 {
754  return matrix;
755 }
756 
757 
758 template <class MATRIX>
759 inline
761 {
762  return matrix;
763 }
764 
765 
766 template <class MATRIX>
767 inline void
769  const size_type gi,
770  const size_type gj,
771  const typename MATRIX::value_type value)
772 {
773  Assert(row_indices.size() != 0, ExcNotInitialized());
774  Assert(column_indices.size() != 0, ExcNotInitialized());
775 
776  const std::pair<unsigned int, size_type> bi
777  = row_indices.global_to_local(gi);
778  const std::pair<unsigned int, size_type> bj
779  = column_indices.global_to_local(gj);
780 
781  Assert (bi.first == row, ExcBlockIndexMismatch(bi.first, row));
782  Assert (bj.first == column, ExcBlockIndexMismatch(bj.first, column));
783 
784  matrix.add(bi.second, bj.second, value);
785 }
786 
787 
788 template <class MATRIX>
789 template <typename number>
790 inline
791 void
792 MatrixBlock<MATRIX>::add (const std::vector<size_type> &r_indices,
793  const std::vector<size_type> &c_indices,
794  const FullMatrix<number> &values,
795  const bool elide_zero_values)
796 {
797  Assert(row_indices.size() != 0, ExcNotInitialized());
798  Assert(column_indices.size() != 0, ExcNotInitialized());
799 
800  AssertDimension (r_indices.size(), values.m());
801  AssertDimension (c_indices.size(), values.n());
802 
803  for (size_type i=0; i<row_indices.size(); ++i)
804  add (r_indices[i], c_indices.size(), &c_indices[0], &values(i,0),
805  elide_zero_values);
806 }
807 
808 
809 template <class MATRIX>
810 template <typename number>
811 inline
812 void
814  const size_type n_cols,
815  const size_type *col_indices,
816  const number *values,
817  const bool,
818  const bool)
819 {
820  Assert(row_indices.size() != 0, ExcNotInitialized());
821  Assert(column_indices.size() != 0, ExcNotInitialized());
822 
823  const std::pair<unsigned int, size_type> bi
824  = row_indices.global_to_local(b_row);
825 
826  // In debug mode, we check whether
827  // all indices are in the correct
828  // block.
829 
830  // Actually, for the time being, we
831  // leave it at this. While it may
832  // not be the most efficient way,
833  // it is at least thread safe.
834 //#ifdef DEBUG
835  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
836 
837  for (size_type j=0; j<n_cols; ++j)
838  {
839  const std::pair<unsigned int, size_type> bj
840  = column_indices.global_to_local(col_indices[j]);
841  Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
842 
843  matrix.add(bi.second, bj.second, values[j]);
844  }
845 //#endif
846 }
847 
848 
849 template <class MATRIX>
850 template <typename number>
851 inline
852 void
853 MatrixBlock<MATRIX>::add (const std::vector<size_type> &indices,
854  const FullMatrix<number> &values,
855  const bool elide_zero_values)
856 {
857  Assert(row_indices.size() != 0, ExcNotInitialized());
858  Assert(column_indices.size() != 0, ExcNotInitialized());
859 
860  AssertDimension (indices.size(), values.m());
861  Assert (values.n() == values.m(), ExcNotQuadratic());
862 
863  for (size_type i=0; i<indices.size(); ++i)
864  add (indices[i], indices.size(), &indices[0], &values(i,0),
865  elide_zero_values);
866 }
867 
868 
869 
870 template <class MATRIX>
871 template <typename number>
872 inline
873 void
875  const std::vector<size_type> &col_indices,
876  const std::vector<number> &values,
877  const bool elide_zero_values)
878 {
879  Assert(row_indices.size() != 0, ExcNotInitialized());
880  Assert(column_indices.size() != 0, ExcNotInitialized());
881 
882  AssertDimension (col_indices.size(), values.size());
883  add (row, col_indices.size(), &col_indices[0], &values[0],
884  elide_zero_values);
885 }
886 
887 
888 template <class MATRIX>
889 template <class VECTOR>
890 inline
891 void
892 MatrixBlock<MATRIX>::vmult (VECTOR &w, const VECTOR &v) const
893 {
894  matrix.vmult(w,v);
895 }
896 
897 
898 template <class MATRIX>
899 template <class VECTOR>
900 inline
901 void
902 MatrixBlock<MATRIX>::vmult_add (VECTOR &w, const VECTOR &v) const
903 {
904  matrix.vmult_add(w,v);
905 }
906 
907 
908 template <class MATRIX>
909 template <class VECTOR>
910 inline
911 void
912 MatrixBlock<MATRIX>::Tvmult (VECTOR &w, const VECTOR &v) const
913 {
914  matrix.Tvmult(w,v);
915 }
916 
917 
918 template <class MATRIX>
919 template <class VECTOR>
920 inline
921 void
922 MatrixBlock<MATRIX>::Tvmult_add (VECTOR &w, const VECTOR &v) const
923 {
924  matrix.Tvmult_add(w,v);
925 }
926 
927 
928 template <class MATRIX>
929 inline
930 std::size_t
932 {
933  return (sizeof(*this)
935  - sizeof(matrix));
936 }
937 
938 //----------------------------------------------------------------------//
939 
940 template <class MATRIX>
941 inline void
943  size_type row, size_type column,
944  const std::string &name)
945 {
946  std_cxx1x::shared_ptr<value_type> p(new value_type(row, column));
948 }
949 
950 
951 template <class MATRIX>
952 inline void
954 {
955  for (size_type i=0; i<this->size(); ++i)
956  {
957  block(i).reinit(sparsity);
958  }
959 }
960 
961 
962 template <class MATRIX>
963 inline void
965 {
966  if (really_clean)
967  {
968  Assert(false, ExcNotImplemented());
969  }
970  else
971  {
972  for (size_type i=0; i<this->size(); ++i)
973  matrix(i).clear();
974  }
975 }
976 
977 
978 
979 template <class MATRIX>
980 inline const MatrixBlock<MATRIX> &
982 {
983  return *this->read(i);
984 }
985 
986 
987 template <class MATRIX>
988 inline MatrixBlock<MATRIX> &
990 {
991  return *(*this)(i);
992 }
993 
994 
995 template <class MATRIX>
996 inline MATRIX &
998 {
999  return (*this)(i)->matrix;
1000 }
1001 
1002 
1003 
1004 //----------------------------------------------------------------------//
1005 
1006 template <class MATRIX>
1007 inline
1009  const bool e, const bool f)
1010  :
1011  edge_matrices(e),
1012  edge_flux_matrices(f)
1013 {}
1014 
1015 
1016 template <class MATRIX>
1017 inline
1018 unsigned int
1020 {
1021  return matrices.size();
1022 }
1023 
1024 
1025 template <class MATRIX>
1026 inline void
1028  size_type row, size_type column,
1029  const std::string &name)
1030 {
1032  p[0].row = row;
1033  p[0].column = column;
1034 
1035  matrices.add(p, name);
1036  if (edge_matrices)
1037  {
1038  matrices_in.add(p, name);
1039  matrices_out.add(p, name);
1040  }
1041  if (edge_flux_matrices)
1042  {
1043  flux_matrices_up.add(p, name);
1044  flux_matrices_down.add(p, name);
1045  }
1046 }
1047 
1048 
1049 template <class MATRIX>
1050 inline const MGLevelObject<MatrixBlock<MATRIX> > &
1052 {
1053  return matrices.read(i);
1054 }
1055 
1056 
1057 template <class MATRIX>
1060 {
1061  return matrices(i);
1062 }
1063 
1064 
1065 template <class MATRIX>
1066 inline const MGLevelObject<MatrixBlock<MATRIX> > &
1068 {
1069  return matrices_in.read(i);
1070 }
1071 
1072 
1073 template <class MATRIX>
1076 {
1077  return matrices_in(i);
1078 }
1079 
1080 
1081 template <class MATRIX>
1082 inline const MGLevelObject<MatrixBlock<MATRIX> > &
1084 {
1085  return matrices_out.read(i);
1086 }
1087 
1088 
1089 template <class MATRIX>
1092 {
1093  return matrices_out(i);
1094 }
1095 
1096 
1097 template <class MATRIX>
1098 inline const MGLevelObject<MatrixBlock<MATRIX> > &
1100 {
1101  return flux_matrices_up.read(i);
1102 }
1103 
1104 
1105 template <class MATRIX>
1108 {
1109  return flux_matrices_up(i);
1110 }
1111 
1112 
1113 template <class MATRIX>
1114 inline const MGLevelObject<MatrixBlock<MATRIX> > &
1116 {
1117  return flux_matrices_down.read(i);
1118 }
1119 
1120 
1121 template <class MATRIX>
1124 {
1125  return flux_matrices_down(i);
1126 }
1127 
1128 
1129 template <class MATRIX>
1130 inline void
1132 {
1133  for (size_type i=0; i<this->size(); ++i)
1134  {
1135  MGLevelObject<MatrixBlock<MATRIX> > &o = block(i);
1136  const size_type row = o[o.min_level()].row;
1137  const size_type col = o[o.min_level()].column;
1138 
1139  o.resize(sparsity.min_level(), sparsity.max_level());
1140  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1141  {
1142  o[level].row = row;
1143  o[level].column = col;
1144  internal::reinit(o[level], sparsity[level]);
1145  }
1146  }
1147 }
1148 
1149 
1150 template <class MATRIX>
1151 inline void
1153 {
1154  for (size_type i=0; i<this->size(); ++i)
1155  {
1156  MGLevelObject<MatrixBlock<MATRIX> > &o = block(i);
1157  const size_type row = o[o.min_level()].row;
1158  const size_type col = o[o.min_level()].column;
1159 
1160  block_in(i).resize(sparsity.min_level(), sparsity.max_level());
1161  block_out(i).resize(sparsity.min_level(), sparsity.max_level());
1162  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1163  {
1164  block_in(i)[level].row = row;
1165  block_in(i)[level].column = col;
1166  internal::reinit(block_in(i)[level], sparsity[level]);
1167  block_out(i)[level].row = row;
1168  block_out(i)[level].column = col;
1169  internal::reinit(block_out(i)[level], sparsity[level]);
1170  }
1171  }
1172 }
1173 
1174 
1175 template <class MATRIX>
1176 inline void
1178 {
1179  for (size_type i=0; i<this->size(); ++i)
1180  {
1181  MGLevelObject<MatrixBlock<MATRIX> > &o = block(i);
1182  const size_type row = o[o.min_level()].row;
1183  const size_type col = o[o.min_level()].column;
1184 
1185  block_up(i).resize(sparsity.min_level(), sparsity.max_level());
1186  block_down(i).resize(sparsity.min_level(), sparsity.max_level());
1187  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1188  {
1189  block_up(i)[level].row = row;
1190  block_up(i)[level].column = col;
1191  internal::reinit(block_up(i)[level], sparsity[level]);
1192  block_down(i)[level].row = row;
1193  block_down(i)[level].column = col;
1194  internal::reinit(block_down(i)[level], sparsity[level]);
1195  }
1196 
1197  }
1198 }
1199 
1200 
1201 template <class MATRIX>
1202 inline void
1204 {
1205  for (size_type i=0; i<mo.size(); ++i)
1206  {
1207  MGLevelObject<MatrixBlock<MATRIX> > &o = mo(i);
1208  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1209  o[level].matrix.clear();
1210  }
1211 }
1212 
1213 
1214 template <class MATRIX>
1215 inline void
1217 {
1218  if (really_clean)
1219  {
1220  Assert(false, ExcNotImplemented());
1221  }
1222  else
1223  {
1224  clear_object(matrices);
1225  clear_object(matrices_in);
1226  clear_object(matrices_out);
1227  clear_object(flux_matrices_up);
1228  clear_object(flux_matrices_down);
1229  }
1230 }
1231 
1232 
1233 
1234 
1235 
1236 DEAL_II_NAMESPACE_CLOSE
1237 
1238 #endif
MatrixBlock< MATRIX > value_type
Definition: matrix_block.h:418
void unsubscribe(const char *identifier=0) const
BlockIndices column_indices
Definition: matrix_block.h:389
const value_type & block_down(size_type i) const
const value_type & block_up(size_type i) const
std::size_t memory_consumption() const
void clear_object(NamedData< MGLevelObject< MatrixBlock< MATRIX > > > &)
Clear one of the matrix objects.
NamedData< MGLevelObject< MatrixBlock< MATRIX > > > matrices
The level matrices.
Definition: matrix_block.h:675
void subscribe(const char *identifier=0) const
const value_type & block_in(size_type i) const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
const std::string & name(unsigned int i) const
Name of object at index.
void Tvmult(VECTOR &w, const VECTOR &v) const
Definition: matrix_block.h:912
void Tvmult_add(VECTOR &w, const VECTOR &v) const
Definition: matrix_block.h:922
unsigned int size() const
NamedData< MGLevelObject< MatrixBlock< MATRIX > > > flux_matrices_down
The DG flux from a level to the lower level.
Definition: matrix_block.h:681
void vmult_add(VECTOR &w, const VECTOR &v) const
Definition: matrix_block.h:902
std::size_t memory_consumption() const
Definition: matrix_block.h:931
Auxiliary class aiding in the handling of block structures like in BlockVector or FESystem...
Definition: block_indices.h:55
const value_type & block(size_type i) const
Definition: matrix_block.h:981
void add(const size_type i, const size_type j, const typename MATRIX::value_type value)
Definition: matrix_block.h:768
void add(size_type row, size_type column, const std::string &name)
const bool edge_matrices
Flag for storing matrices_in and matrices_out.
Definition: matrix_block.h:669
MGMatrixBlockVector(const bool edge_matrices=false, const bool edge_flux_matrices=false)
types::global_dof_index size_type
Definition: matrix_block.h:112
SparsityPatternBase & block(const size_type row, const size_type column)
void reinit(const BlockSparsityPattern &sparsity)
Definition: matrix_block.h:953
types::global_dof_index size_type
Definition: matrix_block.h:413
types::global_dof_index size_type
Definition: matrix_block.h:506
void clear(bool really_clean=false)
Definition: matrix_block.h:964
size_type n() const
void reinit_matrix(const MGLevelObject< BlockSparsityPattern > &sparsity)
unsigned int global_dof_index
Definition: types.h:100
#define Assert(cond, exc)
Definition: exceptions.h:299
std::size_t memory_consumption(const T &t)
unsigned int max_level() const
void vmult(VECTOR &w, const VECTOR &v) const
Definition: matrix_block.h:892
const BlockIndices & get_column_indices() const
void reinit(const BlockSparsityPattern &sparsity)
Definition: matrix_block.h:744
unsigned int size() const
Number of stored data objects.
NamedData< MGLevelObject< MatrixBlock< MATRIX > > > matrices_out
The matrix from the refinement edge to the interior of a level.
Definition: matrix_block.h:679
void reinit_edge(const MGLevelObject< BlockSparsityPattern > &sparsity)
size_type m() const
DeclException2(ExcBlockIndexMismatch, size_type, size_type,<< "Block index "<< arg1<< " does not match "<< arg2)
void clear(bool really_clean=false)
size_type row
Definition: matrix_block.h:359
BlockIndices row_indices
Definition: matrix_block.h:380
const BlockIndices & get_row_indices() const
NamedData< MGLevelObject< MatrixBlock< MATRIX > > > flux_matrices_up
The DG flux from the lower level to a level.
Definition: matrix_block.h:683
void reinit_edge_flux(const MGLevelObject< BlockSparsityPattern > &sparsity)
const value_type & block(size_type i) const
unsigned int min_level() const
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcNotInitialized()
const value_type & block_out(size_type i) const
std::size_t memory_consumption() const
NamedData< MGLevelObject< MatrixBlock< MATRIX > > > matrices_in
The matrix from the interior of a level to the refinement edge.
Definition: matrix_block.h:677
const bool edge_flux_matrices
Flag for storing flux_matrices_up and flux_matrices_down.
Definition: matrix_block.h:672
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel)
void add(size_type row, size_type column, const std::string &name)
Definition: matrix_block.h:942
MATRIX matrix
Definition: matrix_block.h:371
MATRIX & matrix(size_type i)
Definition: matrix_block.h:997
MGLevelObject< MatrixBlock< MATRIX > > value_type
Definition: matrix_block.h:511
size_type column
Definition: matrix_block.h:366