src/INTERP_KERNEL/InterpKernelMatrix.hxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
   2 //
   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.
   7 //
   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.
  12 //
  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
  16 //
  17 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 #ifndef __INTERPKERNELMATRIX_HXX_
  20 #define __INTERPKERNELMATRIX_HXX__
  21
  22 #include "InterpolationUtils.hxx"
  23
  24 #include <vector>
  25 #include <iostream>
  26 #include <ostream>
  27 #include <istream>
  28 #include <map>
  29
  30 namespace INTERP_KERNEL
  31 {
  32   template<class T, NumberingPolicy type>
  33   class Matrix;
  34
  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);
  39
  40   template<class T, NumberingPolicy type=ALL_C_MODE>
  41   class Matrix
  42   {
  43
  44     class KeyComparator
  45     {
  46     public:
  47       KeyComparator(int val):_val(val) { }
  48       bool operator()(const typename std::pair<int,T>& val) { return val.first==_val; }
  49     protected:
  50       int _val;
  51     };
  52
  53     class Row : public std::vector< typename std::pair<int,T> >
  54     {
  55     public:
  56       Row():std::vector< typename std::pair<int,T> >(){};
  57       Row (const Row& row)
  58       {
  59         this->resize(row.size());
  60         for (int i=0; i<this->size(); i++)
  61           (*this)[i]=row[i];
  62       }
  63       Row& operator= (const Row& row)
  64       {
  65         this->resize(row.size());
  66         for (int i=0; i<this->size(); i++)
  67           (*this)[i]=row[i];
  68         return *this;
  69       }
  70       typename std::vector< std::pair<int,T> >::const_iterator find(int elem) const
  71       {
  72         return std::find_if(std::vector< typename std::pair<int,T> >::begin(),std::vector< typename std::pair<int,T> >::end(),KeyComparator(elem));
  73       }
  74
  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))); }
  76
  77       void insert(const std::pair<int,T>& myPair) { push_back(myPair); }
  78     };
  79
  80   private:
  81     unsigned int _nb_rows;
  82     T* _coeffs;
  83     unsigned int* _cols;
  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);
  88     bool _is_configured;
  89   public:
  90     typedef Row value_type;
  91   public:
  92     Matrix():_coeffs(0), _cols(0), _nb_rows(0), _is_configured(false)
  93     { }
  94     Matrix(int nbrows):_coeffs(0), _cols(0), _is_configured(false)
  95     { _nb_rows=nbrows; }
  96     Matrix(std::vector<std::map<int,T> > & matrix) :
  97       _coeffs(0), _cols(0), _is_configured(false)
  98     {
  99       _nb_rows=matrix.size();
 100       for (int i=0; i<_nb_rows; i++)
 101         {
 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);
 106         }
 107
 108     }
 109     /*!Copy constructor
 110      */
 111     Matrix(const Matrix & m)
 112     {
 113       _is_configured=m._is_configured;
 114       _nb_rows=m._nb_rows;
 115       _auxiliary_matrix=m._auxiliary_matrix;
 116       _ncols_offset=m._ncols_offset;
 117       if (_is_configured)
 118         {
 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));
 124         }
 125     }
 126
 127     ~Matrix()
 128     {
 129       delete[] _coeffs;
 130       delete[] _cols;
 131     }
 132
 133     Matrix& operator=(const Matrix& m)
 134     {
 135       _is_configured=m._is_configured;
 136       _nb_rows=m._nb_rows;
 137       _auxiliary_matrix=m._auxiliary_matrix;
 138       _ncols_offset=m._ncols_offset;
 139       if (_is_configured)
 140         {
 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));
 146         }
 147       return this;
 148     }
 149
 150     /*! declares a method that specifies the number of rows */
 151     void resize(unsigned int nbrows)
 152     {
 153       _nb_rows=nbrows;
 154       _auxiliary_matrix.resize(nbrows);
 155     }
 156
 157     /*! sets (i,j) coefficient to \a value */
 158     void setIJ(int irow, int icol, T value)
 159     {
 160       if (_is_configured)
 161         throw Exception("filling a configured matrix");
 162       if (_auxiliary_matrix.empty())
 163         _auxiliary_matrix.resize(_nb_rows);
 164
 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)
 167           {
 168             _auxiliary_matrix[OTT<int,type>::ind2C(irow)][i].second = value;
 169             return;
 170           }
 171       _auxiliary_matrix[OTT<int,type>::ind2C(irow)].push_back(std::make_pair(icol, value));
 172     }
 173
 174     /*!
 175
 176       Matrix multiplies vector \a input and stores the result in
 177       vector \a output.
 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.
 181     */
 182     void multiply(const T* const input, T* const output)
 183     {
 184       if (!_is_configured)
 185         configure();
 186
 187       for (int i=0; i< _nb_rows; i++)
 188         {
 189           output[i]=0;
 190           for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++) {
 191             int icol = _cols[j];
 192             output[i]+=input[icol]*_coeffs[j];
 193           }
 194         }
 195     }
 196
 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 */
 201
 202     void configure()
 203     {
 204       _ncols_offset.resize(_nb_rows+1);
 205       _ncols_offset[0]=0;
 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++)
 214         {
 215           for (unsigned int j=0; j<_auxiliary_matrix[i].size(); j++)
 216             {
 217               *cols_ptr++ = OTT<int,type>::ind2C(_auxiliary_matrix[i][j].first);
 218               *coeffs_ptr++ = _auxiliary_matrix[i][j].second;
 219             }
 220         }
 221       _auxiliary_matrix.clear();
 222       _is_configured=true;
 223     }
 224
 225     /*!
 226      * 0 <= irow < n
 227      */
 228     Row &operator [] (unsigned int irow)
 229     {
 230       return _auxiliary_matrix[irow];
 231     }
 232
 233   };
 234
 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
 242
 243     for instance, matrix
 244     | 1.0 0.0 0.5 |
 245     | 0.0 1.0 0.0 |
 246     | 0.2 0.0 1.0 |
 247     will be displayed in 0-indexing as
 248     0
 249     3
 250     2 0 1.0 2 0.5
 251     1 1 1.0
 252     2 0 0.2 2 1.0
 253   */
 254
 255   template<class T, NumberingPolicy type>
 256   std::ostream& operator<<(std::ostream& out, const Matrix<T, type>& m)
 257   {
 258     if (m._is_configured)
 259       {
 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++)
 263           {
 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];
 267             out<<std::endl;
 268           }
 269       }
 270     else
 271       {
 272         out << OTT<uint,type>::indFC(0) <<"\n";
 273         out << m._nb_rows <<"\n";
 274         for (uint i=0; i<m._nb_rows; i++)
 275           {
 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;
 280             out <<"\n";
 281           }
 282       }
 283     return out;
 284   }
 285
 286   template<class T, NumberingPolicy type>
 287   std::istream& operator>>(std::istream& in, Matrix<T,type>& m)
 288   {
 289     int index_base_test;
 290     in >> index_base_test;
 291     if (index_base_test!=OTT<uint,type>::indFC(0))
 292       {
 293         std::cerr << "file index is "<<index_base_test<<std::endl;
 294         throw Exception("incompatible indexing reading matrix");
 295       }
 296     in >> m._nb_rows;
 297     m._auxiliary_matrix.resize(m._nb_rows);
 298     for (uint i=0; i<m._nb_rows; i++)
 299       {
 300         uint ncols;
 301         in >> ncols;
 302         m._auxiliary_matrix[i].resize(ncols);
 303         double value;
 304         uint col;
 305         for (uint j=0; j<ncols; j++)
 306           {
 307             in>>col;
 308             in>>value;
 309             m._auxiliary_matrix[i].push_back(std::make_pair(col, value));
 310           }
 311       }
 312     return in;
 313   }
 314 }
 315
 316 #endif