1 // Copyright (C) 2007-2014 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #ifndef __INTERPKERNELMATRIX_HXX_
21 #define __INTERPKERNELMATRIX_HXX__
23 #include "InterpolationUtils.hxx"
31 namespace INTERP_KERNEL
33 template<class T, NumberingPolicy type>
36 template<class U, NumberingPolicy type>
37 std::ostream& operator<<(std::ostream& in, const Matrix<U,type>& m);
38 template<class U, NumberingPolicy type>
39 std::istream& operator>>(std::istream& in, Matrix<U,type>& m);
41 template<class T, NumberingPolicy type=ALL_C_MODE>
48 KeyComparator(int val):_val(val) { }
49 bool operator()(const typename std::pair<int,T>& val) { return val.first==_val; }
54 class Row : public std::vector< std::pair<int,T> >
57 Row():std::vector< std::pair<int,T> >(){};
60 this->resize(row.size());
61 for (int i=0; i<this->size(); i++)
64 Row& operator= (const Row& row)
66 this->resize(row.size());
67 for (int i=0; i<this->size(); i++)
71 typename std::vector< std::pair<int,T> >::const_iterator find(int elem) const
73 return std::find_if(std::vector< std::pair<int,T> >::begin(),std::vector< std::pair<int,T> >::end(),KeyComparator(elem));
76 void erase(int elem) { std::vector< std::pair<int,T> >::erase(std::find_if(std::vector< std::pair<int,T> >::begin(),std::vector< std::pair<int,T> >::end(),KeyComparator(elem))); }
78 void insert(const std::pair<int,T>& myPair) { vector<std::pair<int,T> >::push_back(myPair); }
82 unsigned int _nb_rows;
85 std::vector<unsigned int> _ncols_offset;
86 std::vector< Row > _auxiliary_matrix;
87 friend std::ostream& operator<<<>(std::ostream& in, const Matrix<T,type>& m);
88 friend std::istream& operator>><>(std::istream& in, Matrix<T,type>& m);
91 typedef Row value_type;
93 Matrix():_coeffs(0), _cols(0), _nb_rows(0), _is_configured(false)
95 Matrix(int nbrows):_coeffs(0), _cols(0), _is_configured(false)
97 Matrix(std::vector<std::map<int,T> > & matrix) :
98 _coeffs(0), _cols(0), _is_configured(false)
100 _nb_rows=matrix.size();
101 _auxiliary_matrix.resize(_nb_rows);
102 for (int i=0; i<_nb_rows; i++)
104 typename std::map<int,T>::iterator it;
105 for (it=matrix[i].begin(); it != matrix[i].end(); it++)
106 _auxiliary_matrix[i].push_back(*it);//MN: pq push_back plutot que simple affectation?
111 Matrix(const Matrix & m)
113 _is_configured=m._is_configured;
115 _auxiliary_matrix=m._auxiliary_matrix;
116 _ncols_offset=m._ncols_offset;
119 int size=_ncols_offset[_nb_rows];
120 _coeffs = new double[size];
121 _cols = new unsigned int[size];
122 memcpy(_coeffs, m._coeffs, size*sizeof(double));
123 memcpy(_cols, m._cols, size*sizeof(int));
133 Matrix& operator=(const Matrix& m)
135 _is_configured=m._is_configured;
137 _auxiliary_matrix=m._auxiliary_matrix;
138 _ncols_offset=m._ncols_offset;
141 int size=_ncols_offset[_nb_rows];
142 _coeffs = new double[size];
143 _cols = new unsigned int[size];
144 memcpy(_coeffs, m._coeffs, size*sizeof(double));
145 memcpy(_cols, m._cols, size*sizeof(int));
150 /*! declares a method that specifies the number of rows */
151 void resize(unsigned int nbrows)
154 _auxiliary_matrix.resize(nbrows);
157 /*! sets (i,j) coefficient to \a value */
158 void setIJ(int irow, int icol, T value)
161 throw Exception("filling a configured matrix");
162 if (_auxiliary_matrix.empty())
163 _auxiliary_matrix.resize(_nb_rows);
165 for (unsigned int i=0; i< _auxiliary_matrix[OTT<int,type>::ind2C(irow)].size(); i++)
166 if (_auxiliary_matrix[OTT<int,type>::ind2C(irow)][i].first == icol)
168 _auxiliary_matrix[OTT<int,type>::ind2C(irow)][i].second = value;
171 _auxiliary_matrix[OTT<int,type>::ind2C(irow)].push_back(std::make_pair(icol, value));
175 Matrix multiplies vector \a input and stores the result in
177 The vector pointed by \a input must be dimensioned
178 to the number of columns while the vector pointed by output must be
179 dimensioned to the number of rows.
181 void multiply(const T* const input, T* const output)
186 for (int i=0; i< _nb_rows; i++)
189 for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
192 output[i]+=input[icol]*_coeffs[j];
198 Matrix multiplies vector \a input and stores the result in
200 input and output are supposed to represent the same field
201 discretised on two different on meshes.
202 nb_comp is the number of components of the fields input and output
203 The vector pointed by \a input must be dimensioned
204 to the number of columns times nb_comp while the vector pointed by output must be
205 dimensioned to the number of rows times nb_comp.
207 void multiply(const T* const input, T* const output, int nb_comp)
212 for (int i=0; i< _nb_rows; i++)
214 for(int comp = 0; comp < nb_comp; comp++)
215 output[i*nb_comp+comp]=0.;
216 for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
219 for(int comp = 0; comp < nb_comp; comp++)
220 output[i*nb_comp+comp]+=input[icol*nb_comp+comp]*_coeffs[j];
225 Transpose-multiplies vector \a input and stores the result in
227 nb_cols is the number of columns of the matrix, (it is not an attribute of the class)
228 The vector pointed by \a input must be dimensioned
229 to the number of lines _nb_rows while the vector pointed by output must be
230 dimensioned to the number of columns nb_cols.
232 void transposeMultiply(const T* const input, T* const output, int nb_cols)
237 for (int icol=0; icol< nb_cols; icol++)
239 for (int i=0; i< _nb_rows; i++)
241 for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
244 output[icol]+=input[i]*_coeffs[j];
249 Transpose-multiplies vector \a input and stores the result in
251 input and output are supposed to represent the same field
252 discretised on two different on meshes.
253 nb_comp is the number of components of the fields input and output
254 nb_cols is the number of columns of the matrix, (it is not an attribute of the class)
255 The vector pointed by \a input must be dimensioned
256 to _nb_rows*nb_comp while the vector pointed by output must be
257 dimensioned to nb_cols*nb_comp.
259 void transposeMultiply(const T* const input, T* const output, int nb_cols, int nb_comp)
264 for (int icol=0; icol< nb_cols; icol++)
265 for(int comp = 0; comp < nb_comp; comp++)
266 output[icol*nb_comp+comp]=0.;
268 for (int i=0; i< _nb_rows; i++)
270 for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
273 for(int comp = 0; comp < nb_comp; comp++)
274 output[icol*nb_comp+comp]+=input[i*nb_comp+comp]*_coeffs[j];
279 Sums the coefficients of each column of the matrix
280 nb_cols is the number of columns of the matrix, (it is not an attribute of the class)
281 The vector output must be dimensioned to nb_cols
283 void colSum(std::vector< T >& output, int nb_cols)
287 for (int icol=0; icol< nb_cols; icol++)
289 for (int i=0; i< _nb_rows; i++)
291 for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
294 output[icol]+=_coeffs[j];
300 Sums the coefficients of each row of the matrix
301 The vector output must be dimensioned to _nb_rows
303 void rowSum(std::vector< T >& output)
307 for (int i=0; i< _nb_rows; i++)
310 for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
311 output[i]+=_coeffs[j];
315 /*! This operation freezes the profile of the matrix
316 and puts it under a CSR form so that it becomes
317 efficient both in terms of memory occupation and
318 in terms of multiplication */
322 _ncols_offset.resize(_nb_rows+1);
324 for (unsigned int i=0; i<_nb_rows; i++)
325 _ncols_offset[i+1]=_ncols_offset[i]+_auxiliary_matrix[i].size();
326 int nbcoeffs= _ncols_offset[_nb_rows];
327 _cols=new unsigned int[nbcoeffs];
328 _coeffs=new T[nbcoeffs];
329 unsigned int* cols_ptr=_cols;
330 T* coeffs_ptr=_coeffs;
331 for (unsigned int i=0; i<_nb_rows; i++)
333 for (unsigned int j=0; j<_auxiliary_matrix[i].size(); j++)
335 *cols_ptr++ = OTT<int,type>::ind2C(_auxiliary_matrix[i][j].first);
336 *coeffs_ptr++ = _auxiliary_matrix[i][j].second;
339 _auxiliary_matrix.clear();
346 Row &operator [] (unsigned int irow)
348 return _auxiliary_matrix[irow];
358 /*! output to an ascii file
359 only nonzero elements are written
360 - the first line contains the indexing (0 or 1)
361 - the second line contains the number of rows.
362 - for each row, a line contains:
363 - the number of nonzero coeffs
364 - and for each coeff : icol, value
370 will be displayed in 0-indexing as
378 template<class T, NumberingPolicy type>
379 std::ostream& operator<<(std::ostream& out, const Matrix<T, type>& m)
381 if (m._is_configured)
383 out << OTT<unsigned int,type>::indFC(0) <<std::endl;
384 out << m._nb_rows<<std::endl;
385 for (unsigned int i=0; i<m._nb_rows; i++)
387 out << m._ncols_offset[i+1]-m._ncols_offset[i];
388 for (unsigned int j=m._ncols_offset[i]; j<m._ncols_offset[i+1]; j++)
389 out <<"\t"<< OTT<unsigned int,type>::indFC(m._cols[j]) <<"\t"<<m._coeffs[j];
395 out << OTT<unsigned int,type>::indFC(0) <<"\n";
396 out << m._nb_rows <<"\n";
397 for (unsigned int i=0; i<m._nb_rows; i++)
399 out<< m._auxiliary_matrix[i].size();
400 for (unsigned int j=0; j<m._auxiliary_matrix[i].size(); j++)
401 out << "\t" <<m._auxiliary_matrix[i][j].first <<"\t"
402 <<m._auxiliary_matrix[i][j].second;
409 template<class T, NumberingPolicy type>
410 std::istream& operator>>(std::istream& in, Matrix<T,type>& m)
413 in >> index_base_test;
414 if (index_base_test!=OTT<unsigned int,type>::indFC(0))
416 std::cerr << "file index is "<<index_base_test<<std::endl;
417 throw Exception("incompatible indexing reading matrix");
420 m._auxiliary_matrix.resize(m._nb_rows);
421 for (unsigned int i=0; i<m._nb_rows; i++)
425 m._auxiliary_matrix[i].resize(ncols);
428 for (unsigned int j=0; j<ncols; j++)
432 m._auxiliary_matrix[i].push_back(std::make_pair(col, value));