Rheolef  7.1
an efficient C++ finite element environment
asr.cc
Go to the documentation of this file.
1 #ifdef USE_NEW_ASR
22 #define asr asr
23 #endif // USE_NEW_ASR
24 
25 #include "rheolef/asr.h"
26 #include "rheolef/csr.h"
27 #include "rheolef/iorheo.h"
28 #include "rheolef/mm_io.h"
29 
30 namespace rheolef {
31 
32 // ----------------------------------------------------------------------------
33 // csr2asr convertion
34 // ----------------------------------------------------------------------------
35 template<class T, class M, class A>
36 void
38 {
40  dia_ia = a.begin(),
41  ext_ia = a.ext_begin();
42  size_type first_dis_i = a.row_ownership().first_index();
43  size_type first_dis_j = a.col_ownership().first_index();
44  for (size_type i = 0, n = a.nrow(); i < n; i++) {
45  for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; ++p) {
46  size_type dis_j = first_dis_j + (*p).first;
47  const T& value = (*p).second;
48  semi_dis_entry (i, dis_j) += value;
49  }
50  if (is_sequential<M>::value || a.ext_nnz() == 0) continue;
51  size_type dis_i = first_dis_i + i;
52  for (typename csr<T,M>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++) {
53  size_type dis_j = a.jext2dis_j ((*p).first);
54  const T& value = (*p).second;
55  semi_dis_entry (i, dis_j) += value;
56  }
57  }
58  dis_entry_assembly();
59 }
60 // ----------------------------------------------------------------------------
61 // after a asr::dis_entry_assembly, may recompute nnz and dis_nnz
62 // ----------------------------------------------------------------------------
63 template <class T, class M, class A>
64 void
66 {
67  _nnz = 0;
68  for (typename base::const_iterator
69  iter = base::begin(),
70  last = base::end(); iter != last; ++iter) {
71  _nnz += (*iter).size();
72  }
73  _dis_nnz = _nnz;
74 #ifdef _RHEOLEF_HAVE_MPI
76  _dis_nnz = mpi::all_reduce (comm(), _nnz, std::plus<size_type>());
77  }
78 #endif // _RHEOLEF_HAVE_MPI
79 }
80 // ----------------------------------------------------------------------------
81 // io
82 // ----------------------------------------------------------------------------
83 template <class T, class M, class A>
86 {
87  std::ostream& os = ods.os();
88  std::string name = iorheo::getbasename(ods.os());
89  if (name == "") name = "a";
90  os << std::setprecision (std::numeric_limits<T>::digits10)
91  << name << " = spalloc(" << nrow() << "," << ncol() << "," << nnz() << ");" << std::endl;
92  if (_nnz == 0) return ods;
93  size_type i = 0;
94  for (typename base::const_iterator
95  iter = base::begin(),
96  last = base::end(); iter != last; ++iter, ++i) {
97  const row_type& row = *iter;
98  for (typename row_type::const_iterator
99  row_iter = row.begin(),
100  row_last = row.end(); row_iter != row_last; ++row_iter) {
101  os << name << "(" << first_dis_i+i+1 << ","
102  << (*row_iter).first+1 << ") = "
103  << (*row_iter).second << ";" << std::endl;
104  }
105  }
106  return ods;
107 }
108 template <class T, class M, class A>
111 {
112  std::ostream& os = ods.os();
113  os << std::setprecision(std::numeric_limits<T>::digits10)
114  << "%%MatrixMarket matrix coordinate real general" << std::endl
115  << nrow() << " " << ncol() << " " << nnz() << std::endl;
116  if (_nnz == 0) return ods;
117  size_type i = 0;
118  for (typename base::const_iterator
119  iter = base::begin(),
120  last = base::end(); iter != last; ++iter, ++i) {
121  const row_type& row = *iter;
122  for (typename row_type::const_iterator
123  row_iter = row.begin(),
124  row_last = row.end(); row_iter != row_last; ++row_iter) {
125  os << first_dis_i+i+1 << " "
126  << (*row_iter).first+1 << " "
127  << (*row_iter).second << std::endl;
128  }
129  }
130  return ods;
131 }
132 template <class T, class M, class A>
134 asr<T,M,A>::put_seq (odiststream& ods, size_type first_dis_i) const
135 {
136  iorheo::flag_type format = iorheo::flags(ods.os()) & iorheo::format_field;
137  if (format [iorheo::matlab] || format [iorheo::sparse_matlab])
138  { return put_seq_sparse_matlab (ods,first_dis_i); }
139  // default is matrix market format
140  return put_seq_matrix_market (ods,first_dis_i);
141 }
142 template <class T, class M, class A>
145 {
146  // send all in a pseudo-sequential matrix on process 0
147  size_type io_proc = odiststream::io_proc();
148  size_type my_proc = comm().rank();
149  distributor io_row_ownership (dis_nrow(), comm(), (my_proc == io_proc ? dis_nrow() : 0));
150  distributor io_col_ownership (dis_ncol(), comm(), (my_proc == io_proc ? dis_ncol() : 0));
151  asr<T> a (io_row_ownership, io_col_ownership);
152  size_type first_dis_i = row_ownership().first_index();
153  size_type i = 0;
154  for (typename base::const_iterator
155  iter = base::begin(),
156  last = base::end(); iter != last; ++iter, ++i) {
157  const row_type& row = *iter;
158  for (typename row_type::const_iterator
159  row_iter = row.begin(),
160  row_last = row.end(); row_iter != row_last; ++row_iter) {
161  a.dis_entry (first_dis_i + i, (*row_iter).first) += (*row_iter).second;
162  }
163  }
164  a.dis_entry_assembly();
165  if (my_proc == io_proc) {
166  // print sequential matrix on io_proc
167  a.put_seq (ods);
168  }
169  return ods;
170 }
171 template <class T, class M, class A>
174 {
176  return put_seq (ods);
177  } else {
178  return put_mpi (ods);
179  }
180 }
181 template <class T, class M, class A>
182 idiststream&
183 asr<T,M,A>::get (idiststream& ips)
184 {
185  // matrix market format:
186  // %%mm_header
187  // nrow ncol nnz
188  // {i, j, aij}*
189  // we suppose nrow ncol already readed
190  // and the matrix has good dimensions
191  struct matrix_market mm = read_matrix_market_header (ips);
192  check_macro (mm.format != matrix_market::hermitian, "Hermitian matrix not yet supported");
193  size_type dis_nrow, dis_ncol, dis_nnz;
194  ips >> dis_nrow >> dis_ncol >> dis_nnz;
195  communicator comm = ips.comm();
196  distributor row_ownership (dis_nrow, comm, distributor::decide);
197  distributor col_ownership (dis_ncol, comm, distributor::decide);
198  resize (row_ownership, col_ownership);
199  size_type my_proc = comm.rank();
200  if (my_proc == idiststream::io_proc()) {
201  std::istream& is = ips.is();
202  for (size_type p = 0; p < dis_nnz; p++) {
203  size_type dis_i, dis_j;
204  T value;
205  is >> dis_i >> dis_j >> value;
206  dis_i--; dis_j--; // 1..N convention -> 0..N-1
207  dis_entry(dis_i,dis_j) += value;
208  if (mm.format == matrix_market::general || dis_i == dis_j) continue;
209  // when symmetries, add also the corespondings coefs, for simplicity
210  if (mm.format == matrix_market::symmetric) {
211  dis_entry(dis_j,dis_i) += value;
212  } else if (mm.format == matrix_market::skew_symmetric) {
213  dis_entry(dis_j,dis_i) += -value;
214  }
215  }
216  }
217  dis_entry_assembly();
218  return ips;
219 }
220 // ----------------------------------------------------------------------------
221 // instanciation in library
222 // ----------------------------------------------------------------------------
223 #define _RHEOLEF_instanciation(T,M,A) \
224 template class asr<T,M,A>;
225 
227 #ifdef _RHEOLEF_HAVE_MPI
229 #endif // _RHEOLEF_HAVE_MPI
230 
231 } // namespace rheolef
see the Float page for the full documentation
odiststream & put_seq_sparse_matlab(odiststream &ops, size_type first_dis_i=0) const
Definition: asr.cc:85
idiststream & get(idiststream &ips)
Definition: asr.cc:183
odiststream & put_mpi(odiststream &ops) const
Definition: asr.cc:144
odiststream & put(odiststream &ops) const
Definition: asr.cc:173
base::size_type size_type
Definition: asr.h:55
void build_from_csr(const csr_rep< T, M > &)
Definition: asr.cc:37
void _recompute_nnz()
Definition: asr.cc:65
odiststream & put_seq(odiststream &ops, size_type first_dis_i=0) const
Definition: asr.cc:134
odiststream & put_seq_matrix_market(odiststream &ops, size_type first_dis_i=0) const
Definition: asr.cc:110
see the csr page for the full documentation
Definition: csr.h:317
rep::base::const_iterator const_iterator
Definition: disarray.h:465
see the distributor page for the full documentation
Definition: distributor.h:62
static const size_type decide
Definition: distributor.h:76
odiststream: see the diststream page for the full documentation
Definition: diststream.h:126
std::ostream & os()
Definition: diststream.h:236
static size_type io_proc()
Definition: diststream.cc:78
base::const_iterator const_iterator
Definition: pair_set.h:84
rheolef::std value
Expr1::float_type T
Definition: field_expr.h:261
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format matlab
This file is part of Rheolef.
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
struct matrix_market read_matrix_market_header(idiststream &ips)
Definition: mm_io.cc:30
Definition: sphere.icc:25
static const format_type general
Definition: mm_io.h:32
static const format_type symmetric
Definition: mm_io.h:33
static const format_type skew_symmetric
Definition: mm_io.h:34
static const format_type hermitian
Definition: mm_io.h:35
format_type format
Definition: mm_io.h:37