1 // Copyright (C) 2007-2008 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.
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
19 #ifndef __INTERPKERNELMATRIX_HXX_
20 #define __INTERPKERNELMATRIX_HXX__
22 #include "InterpolationUtils.hxx"
30 namespace INTERP_KERNEL
32 template<class T, NumberingPolicy type>
35 template<class U, NumberingPolicy type>
36 std::ostream& operator<<(std::ostream& in, const Matrix<U,type>& m);
37 template<class U, NumberingPolicy type>
38 std::istream& operator>>(std::istream& in, Matrix<U,type>& m);
40 template<class T, NumberingPolicy type=ALL_C_MODE>
47 KeyComparator(int val):_val(val) { }
48 bool operator()(const typename std::pair<int,T>& val) { return val.first==_val; }
53 class Row : public std::vector< typename std::pair<int,T> >
56 Row():std::vector< typename std::pair<int,T> >(){};
59 this->resize(row.size());
60 for (int i=0; i<this->size(); i++)
63 Row& operator= (const Row& row)
65 this->resize(row.size());
66 for (int i=0; i<this->size(); i++)
70 typename std::vector< std::pair<int,T> >::const_iterator find(int elem) const
72 return std::find_if(std::vector< typename std::pair<int,T> >::begin(),std::vector< typename std::pair<int,T> >::end(),KeyComparator(elem));
75 void erase(int elem) { std::vector< typename std::pair<int,T> >::erase(std::find_if(std::vector< typename std::pair<int,T> >::begin(),std::vector< typename std::pair<int,T> >::end(),KeyComparator(elem))); }
77 void insert(const std::pair<int,T>& myPair) { push_back(myPair); }
81 unsigned int _nb_rows;
84 std::vector<unsigned int> _ncols_offset;
85 std::vector< Row > _auxiliary_matrix;
86 friend std::ostream& operator<<<>(std::ostream& in, const Matrix<T,type>& m);
87 friend std::istream& operator>><>(std::istream& in, Matrix<T,type>& m);
90 typedef Row value_type;
92 Matrix():_coeffs(0), _cols(0), _nb_rows(0), _is_configured(false)
94 Matrix(int nbrows):_coeffs(0), _cols(0), _is_configured(false)
96 Matrix(std::vector<std::map<int,T> > & matrix) :
97 _coeffs(0), _cols(0), _is_configured(false)
99 _nb_rows=matrix.size();
100 for (int i=0; i<_nb_rows; i++)
102 _auxiliary_matrix[i].resize(matrix[i].size());
103 typename std::map<int,T>::iterator it;
104 for (it=matrix[i].begin(); it != matrix[i].end(); it++)
105 _auxiliary_matrix[i].push_back(*it);
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 uint[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 uint[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 (uint 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));
176 Matrix multiplies vector \a input and stores the result in
178 The vector pointed by \a input must be dimensioned
179 to the number of columns while the vector pointed by output must be
180 dimensioned to the number of rows.
182 void multiply(const T* const input, T* const output)
187 for (int i=0; i< _nb_rows; i++)
190 for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++) {
192 output[i]+=input[icol]*_coeffs[j];
197 /*! This operation freezes the profile of the matrix
198 and puts it under a CSR form so that it becomes
199 efficient both in terms of memory occupation and
200 in terms of multiplication */
204 _ncols_offset.resize(_nb_rows+1);
206 for (unsigned int i=0; i<_nb_rows; i++)
207 _ncols_offset[i+1]=_ncols_offset[i]+_auxiliary_matrix[i].size();
208 int nbcoeffs= _ncols_offset[_nb_rows];
209 _cols=new unsigned int[nbcoeffs];
210 _coeffs=new T[nbcoeffs];
211 unsigned int* cols_ptr=_cols;
212 T* coeffs_ptr=_coeffs;
213 for (unsigned int i=0; i<_nb_rows; i++)
215 for (unsigned int j=0; j<_auxiliary_matrix[i].size(); j++)
217 *cols_ptr++ = OTT<int,type>::ind2C(_auxiliary_matrix[i][j].first);
218 *coeffs_ptr++ = _auxiliary_matrix[i][j].second;
221 _auxiliary_matrix.clear();
228 Row &operator [] (unsigned int irow)
230 return _auxiliary_matrix[irow];
235 /*! output to an ascii file
236 only nonzero elements are written
237 - the first line contains the indexing (0 or 1)
238 - the second line contains the number of rows.
239 - for each row, a line contains:
240 - the number of nonzero coeffs
241 - and for each coeff : icol, value
247 will be displayed in 0-indexing as
255 template<class T, NumberingPolicy type>
256 std::ostream& operator<<(std::ostream& out, const Matrix<T, type>& m)
258 if (m._is_configured)
260 out << OTT<uint,type>::indFC(0) <<std::endl;
261 out << m._nb_rows<<std::endl;
262 for (uint i=0; i<m._nb_rows; i++)
264 out << m._ncols_offset[i+1]-m._ncols_offset[i];
265 for (uint j=m._ncols_offset[i]; j<m._ncols_offset[i+1]; j++)
266 out <<"\t"<< OTT<uint,type>::indFC(m._cols[j]) <<"\t"<<m._coeffs[j];
272 out << OTT<uint,type>::indFC(0) <<"\n";
273 out << m._nb_rows <<"\n";
274 for (uint i=0; i<m._nb_rows; i++)
276 out<< m._auxiliary_matrix[i].size();
277 for (uint j=0; j<m._auxiliary_matrix[i].size(); j++)
278 out << "\t" <<m._auxiliary_matrix[i][j].first <<"\t"
279 <<m._auxiliary_matrix[i][j].second;
286 template<class T, NumberingPolicy type>
287 std::istream& operator>>(std::istream& in, Matrix<T,type>& m)
290 in >> index_base_test;
291 if (index_base_test!=OTT<uint,type>::indFC(0))
293 std::cerr << "file index is "<<index_base_test<<std::endl;
294 throw Exception("incompatible indexing reading matrix");
297 m._auxiliary_matrix.resize(m._nb_rows);
298 for (uint i=0; i<m._nb_rows; i++)
302 m._auxiliary_matrix[i].resize(ncols);
305 for (uint j=0; j<ncols; j++)
309 m._auxiliary_matrix[i].push_back(std::make_pair(col, value));