Back to index

salome-med  6.5.0
InterpKernelMatrix.hxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #ifndef __INTERPKERNELMATRIX_HXX_
00021 #define __INTERPKERNELMATRIX_HXX__
00022 
00023 #include "InterpolationUtils.hxx"
00024 
00025 #include <vector>
00026 #include <iostream>
00027 #include <ostream>
00028 #include <istream>
00029 #include <map>
00030 
00031 namespace INTERP_KERNEL
00032 {
00033   template<class T, NumberingPolicy type>
00034   class Matrix;
00035   
00036   template<class U, NumberingPolicy type>
00037   std::ostream& operator<<(std::ostream& in, const Matrix<U,type>& m);
00038   template<class U, NumberingPolicy type>
00039   std::istream& operator>>(std::istream& in, Matrix<U,type>& m);
00040         
00041   template<class T, NumberingPolicy type=ALL_C_MODE>
00042   class Matrix
00043   {
00044 
00045     class KeyComparator
00046     {
00047     public:
00048       KeyComparator(int val):_val(val) { }
00049       bool operator()(const typename std::pair<int,T>& val) { return val.first==_val; }
00050     protected:
00051       int _val;
00052     };
00053 
00054     class Row : public std::vector< std::pair<int,T> >
00055     {
00056     public:
00057       Row():std::vector< std::pair<int,T> >(){};
00058       Row (const Row& row)
00059       {
00060         this->resize(row.size());
00061         for (int i=0; i<this->size(); i++)
00062           (*this)[i]=row[i];
00063       }
00064       Row& operator= (const Row& row)
00065       {
00066         this->resize(row.size());
00067         for (int i=0; i<this->size(); i++)
00068           (*this)[i]=row[i];
00069         return *this;
00070       }
00071       typename std::vector< std::pair<int,T> >::const_iterator find(int elem) const
00072       {
00073         return std::find_if(std::vector< std::pair<int,T> >::begin(),std::vector< std::pair<int,T> >::end(),KeyComparator(elem));
00074       }
00075 
00076       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))); }
00077 
00078       void insert(const std::pair<int,T>& myPair) { vector<std::pair<int,T> >::push_back(myPair); }
00079     };
00080     
00081   private:
00082     unsigned int _nb_rows;
00083     T* _coeffs;
00084     unsigned int* _cols;
00085     std::vector<unsigned int> _ncols_offset;
00086     std::vector< Row > _auxiliary_matrix;
00087     friend std::ostream& operator<<<>(std::ostream& in, const Matrix<T,type>& m);
00088     friend std::istream& operator>><>(std::istream& in, Matrix<T,type>& m);
00089     bool _is_configured;
00090   public:
00091     typedef Row value_type;
00092   public:
00093     Matrix():_coeffs(0), _cols(0), _nb_rows(0), _is_configured(false)
00094     { }
00095     Matrix(int nbrows):_coeffs(0), _cols(0), _is_configured(false)
00096     { _nb_rows=nbrows; }
00097     Matrix(std::vector<std::map<int,T> > & matrix) :
00098       _coeffs(0), _cols(0), _is_configured(false)
00099     {
00100       _nb_rows=matrix.size();
00101       _auxiliary_matrix.resize(_nb_rows);
00102       for (int i=0; i<_nb_rows; i++)
00103         {
00104           typename std::map<int,T>::iterator it;
00105           for (it=matrix[i].begin(); it != matrix[i].end(); it++)
00106             _auxiliary_matrix[i].push_back(*it);//MN: pq push_back plutot que simple affectation?
00107         }      
00108     }
00111     Matrix(const Matrix & m)
00112     {
00113       _is_configured=m._is_configured;
00114       _nb_rows=m._nb_rows;
00115       _auxiliary_matrix=m._auxiliary_matrix;
00116       _ncols_offset=m._ncols_offset;
00117       if (_is_configured)
00118         {
00119           int size=_ncols_offset[_nb_rows];
00120           _coeffs = new double[size];
00121           _cols = new unsigned int[size];
00122           memcpy(_coeffs, m._coeffs, size*sizeof(double));
00123           memcpy(_cols, m._cols, size*sizeof(int));
00124         }
00125     }
00126     
00127     ~Matrix()
00128     {
00129       delete[] _coeffs;
00130       delete[] _cols;
00131     }
00132 
00133     Matrix& operator=(const Matrix& m)
00134     {
00135       _is_configured=m._is_configured;
00136       _nb_rows=m._nb_rows;
00137       _auxiliary_matrix=m._auxiliary_matrix;
00138       _ncols_offset=m._ncols_offset;
00139       if (_is_configured)
00140         {
00141           int size=_ncols_offset[_nb_rows];
00142           _coeffs = new double[size];
00143           _cols = new unsigned int[size];
00144           memcpy(_coeffs, m._coeffs, size*sizeof(double));
00145           memcpy(_cols, m._cols, size*sizeof(int));
00146         }
00147       return this;
00148     }
00149 
00151     void resize(unsigned int nbrows)
00152     {
00153       _nb_rows=nbrows;
00154       _auxiliary_matrix.resize(nbrows);
00155     }
00156     
00158     void setIJ(int irow, int icol, T value)
00159     {
00160       if (_is_configured)
00161         throw Exception("filling a configured matrix");
00162       if (_auxiliary_matrix.empty())
00163         _auxiliary_matrix.resize(_nb_rows);
00164       
00165       for (unsigned int i=0; i< _auxiliary_matrix[OTT<int,type>::ind2C(irow)].size(); i++)
00166         if (_auxiliary_matrix[OTT<int,type>::ind2C(irow)][i].first == icol)
00167           {
00168             _auxiliary_matrix[OTT<int,type>::ind2C(irow)][i].second = value;
00169             return;
00170           }
00171       _auxiliary_matrix[OTT<int,type>::ind2C(irow)].push_back(std::make_pair(icol, value));
00172     }
00173     
00181     void multiply(const T* const input, T* const output)
00182     {
00183       if (!_is_configured)
00184         configure();
00185       
00186       for (int i=0; i< _nb_rows; i++)
00187         {
00188           output[i]=0.;
00189           for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
00190             {
00191               int icol = _cols[j];
00192               output[i]+=input[icol]*_coeffs[j];
00193             }
00194         }
00195     }
00196 
00207     void multiply(const T* const input, T* const output, int nb_comp)
00208     {
00209       if (!_is_configured)
00210         configure();
00211       
00212       for (int i=0; i< _nb_rows; i++)
00213         {
00214           for(int comp = 0; comp < nb_comp; comp++)
00215             output[i*nb_comp+comp]=0.;
00216           for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
00217             {
00218               int icol = _cols[j];
00219               for(int comp = 0; comp < nb_comp; comp++)
00220                 output[i*nb_comp+comp]+=input[icol*nb_comp+comp]*_coeffs[j];
00221             }
00222         }
00223     }   
00232     void transposeMultiply(const T* const input, T* const output, int nb_cols)
00233     {
00234       if (!_is_configured)
00235         configure();
00236       
00237       for (int icol=0; icol< nb_cols; icol++)
00238         output[icol]=0.;
00239       for (int i=0; i< _nb_rows; i++)
00240         {
00241           for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
00242             {
00243               int icol = _cols[j];
00244               output[icol]+=input[i]*_coeffs[j];
00245             }
00246         }
00247     }
00259     void transposeMultiply(const T* const input, T* const output, int nb_cols, int nb_comp)
00260     {
00261       if (!_is_configured)
00262         configure();
00263       
00264       for (int icol=0; icol< nb_cols; icol++)
00265         for(int comp = 0; comp < nb_comp; comp++)
00266           output[icol*nb_comp+comp]=0.;
00267 
00268       for (int i=0; i< _nb_rows; i++)
00269         {
00270           for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
00271             {
00272               int icol = _cols[j];
00273               for(int comp = 0; comp < nb_comp; comp++)
00274                 output[icol*nb_comp+comp]+=input[i*nb_comp+comp]*_coeffs[j];
00275             }
00276         }
00277     }
00278     /*
00279       Sums the coefficients of each column of the matrix
00280       nb_cols is the number of columns of the matrix, (it is not an attribute of the class) 
00281       The vector output must be dimensioned to nb_cols
00282     */
00283     void colSum(std::vector< T >& output, int nb_cols)
00284     {
00285       if (!_is_configured)
00286         configure();
00287       for (int icol=0; icol< nb_cols; icol++)
00288         output[icol]=0.;
00289       for (int i=0; i< _nb_rows; i++)
00290         {
00291           for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
00292             {
00293               int icol = _cols[j];
00294               output[icol]+=_coeffs[j];
00295             }
00296         }
00297     }
00298 
00299     /*
00300       Sums the coefficients of each row of the matrix
00301       The vector output must be dimensioned to _nb_rows
00302     */
00303     void rowSum(std::vector< T >& output)
00304     {
00305       if (!_is_configured)
00306         configure();
00307       for (int i=0; i< _nb_rows; i++)
00308         {
00309           output[i]=0;
00310           for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++) 
00311             output[i]+=_coeffs[j];
00312         }
00313     }
00314 
00320     void configure()
00321     {
00322       _ncols_offset.resize(_nb_rows+1);
00323       _ncols_offset[0]=0;
00324       for (unsigned int i=0; i<_nb_rows; i++)
00325         _ncols_offset[i+1]=_ncols_offset[i]+_auxiliary_matrix[i].size();
00326       int nbcoeffs= _ncols_offset[_nb_rows];
00327       _cols=new unsigned int[nbcoeffs];
00328       _coeffs=new T[nbcoeffs];
00329       unsigned int* cols_ptr=_cols;
00330       T* coeffs_ptr=_coeffs;
00331       for (unsigned int i=0; i<_nb_rows; i++)
00332         {
00333           for (unsigned int j=0; j<_auxiliary_matrix[i].size(); j++)
00334             {
00335               *cols_ptr++ = OTT<int,type>::ind2C(_auxiliary_matrix[i][j].first);
00336               *coeffs_ptr++ = _auxiliary_matrix[i][j].second;
00337             }
00338         }
00339       _auxiliary_matrix.clear();
00340       _is_configured=true;
00341     }
00342 
00346     Row &operator [] (unsigned int irow)
00347     {
00348       return _auxiliary_matrix[irow];
00349     }
00350 
00351     int getNbRows()
00352     {
00353       return _nb_rows;
00354     }
00355     
00356   };
00357   
00378   template<class T, NumberingPolicy type>
00379   std::ostream& operator<<(std::ostream& out, const Matrix<T, type>& m)
00380   {
00381     if (m._is_configured)
00382       {
00383         out << OTT<unsigned int,type>::indFC(0) <<std::endl;
00384         out << m._nb_rows<<std::endl;
00385         for (unsigned int i=0; i<m._nb_rows; i++)
00386           {
00387             out << m._ncols_offset[i+1]-m._ncols_offset[i];
00388             for (unsigned int j=m._ncols_offset[i]; j<m._ncols_offset[i+1]; j++)
00389               out <<"\t"<< OTT<unsigned int,type>::indFC(m._cols[j]) <<"\t"<<m._coeffs[j];
00390             out<<std::endl;
00391           }
00392       }
00393     else
00394       {
00395         out << OTT<unsigned int,type>::indFC(0) <<"\n";
00396         out << m._nb_rows <<"\n";
00397         for (unsigned int i=0; i<m._nb_rows; i++)
00398           {
00399             out<< m._auxiliary_matrix[i].size();
00400             for (unsigned int j=0; j<m._auxiliary_matrix[i].size(); j++)
00401               out << "\t" <<m._auxiliary_matrix[i][j].first <<"\t"
00402                   <<m._auxiliary_matrix[i][j].second;
00403             out <<"\n";
00404           }
00405       }
00406     return out;
00407   }
00408   
00409   template<class T, NumberingPolicy type>
00410   std::istream& operator>>(std::istream& in, Matrix<T,type>& m)
00411   {
00412     int index_base_test;
00413     in >> index_base_test;
00414     if (index_base_test!=OTT<unsigned int,type>::indFC(0))
00415       {
00416         std::cerr << "file index is "<<index_base_test<<std::endl;
00417         throw Exception("incompatible indexing reading matrix");
00418       }
00419     in >> m._nb_rows;
00420     m._auxiliary_matrix.resize(m._nb_rows);
00421     for (unsigned int i=0; i<m._nb_rows; i++)
00422       {
00423         unsigned int ncols;
00424         in >> ncols;
00425         m._auxiliary_matrix[i].resize(ncols);
00426         double value;
00427         unsigned int col;
00428         for (unsigned int j=0; j<ncols; j++)
00429           {
00430             in>>col;
00431             in>>value;
00432             m._auxiliary_matrix[i].push_back(std::make_pair(col, value));
00433           }
00434       }
00435     return in;
00436   }
00437 }
00438 
00439 #endif