Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
asr.cc
Go to the documentation of this file.
1
21#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
30namespace rheolef {
31
32// ----------------------------------------------------------------------------
33// csr2asr convertion
34// ----------------------------------------------------------------------------
35template<class T, class M, class A>
36void
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// ----------------------------------------------------------------------------
63template <class T, class M, class A>
64void
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// ----------------------------------------------------------------------------
83template <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}
108template <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}
132template <class T, class M, class A>
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}
142template <class T, class M, class A>
145{
146 // send all in a pseudo-sequential matrix on process 0
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}
171template <class T, class M, class A>
174{
176 return put_seq (ods);
177 } else {
178 return put_mpi (ods);
179 }
180}
181template <class T, class M, class A>
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
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
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) \
224template class asr<T,M,A>;
225
227#ifdef _RHEOLEF_HAVE_MPI
229#endif // _RHEOLEF_HAVE_MPI
230
231} // namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
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
pair_set< T, A > row_type
Definition asr.h:53
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
see the distributor page for the full documentation
Definition distributor.h:69
static const size_type decide
Definition distributor.h:83
idiststream: see the diststream page for the full documentation
Definition diststream.h:336
static size_type io_proc()
This routine returns the rank of a process that can perform i/o.
Definition diststream.cc:64
std::istream & is()
Definition diststream.h:400
const communicator & comm() const
Definition diststream.h:356
std::bitset< last > flag_type
Definition iorheo.h:440
static flag_type format_field
Definition iorheo.h:445
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
std::ostream & os()
Definition diststream.h:247
static size_type io_proc()
Definition diststream.cc:79
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
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