src/INTERP_KERNEL/InterpKernelMatrix.hxx

   1 // Copyright (C) 2007-2020  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, or (at your option) any later version.
   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
  20 #ifndef __INTERPKERNELMATRIX_HXX_
  21 #define __INTERPKERNELMATRIX_HXX__
  22
  23 #include "InterpolationUtils.hxx"
  24
  25 #include <vector>
  26 #include <iostream>
  27 #include <ostream>
  28 #include <istream>
  29 #include <map>
  30
  31 namespace INTERP_KERNEL
  32 {
  33   template<class T, NumberingPolicy type>
  34   class Matrix;
  35
  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);
  40
  41   template<class T, NumberingPolicy type=ALL_C_MODE>
  42   class Matrix
  43   {
  44
  45     class KeyComparator
  46     {
  47     public:
  48       KeyComparator(int val):_val(val) { }
  49       bool operator()(const typename std::pair<int,T>& val) { return val.first==_val; }
  50     protected:
  51       int _val;
  52     };
  53
  54     class Row : public std::vector< std::pair<int,T> >
  55     {
  56     public:
  57       Row():std::vector< std::pair<int,T> >(){};
  58       Row (const Row& row)
  59       {
  60         this->resize(row.size());
  61         for (int i=0; i<this->size(); i++)
  62           (*this)[i]=row[i];
  63       }
  64       Row& operator= (const Row& row)
  65       {
  66         this->resize(row.size());
  67         for (int i=0; i<this->size(); i++)
  68           (*this)[i]=row[i];
  69         return *this;
  70       }
  71       typename std::vector< std::pair<int,T> >::const_iterator find(int elem) const
  72       {
  73         return std::find_if(std::vector< std::pair<int,T> >::begin(),std::vector< std::pair<int,T> >::end(),KeyComparator(elem));
  74       }
  75
  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))); }
  77
  78       void insert(const std::pair<int,T>& myPair) { vector<std::pair<int,T> >::push_back(myPair); }
  79     };
  80
  81   private:
  82     unsigned int _nb_rows;
  83     T* _coeffs;
  84     unsigned int* _cols;
  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);
  89     bool _is_configured;
  90   public:
  91     typedef Row value_type;
  92   public:
  93     Matrix():_coeffs(0), _cols(0), _nb_rows(0), _is_configured(false)
  94     { }
  95     Matrix(int nbrows):_coeffs(0), _cols(0), _is_configured(false)
  96     { _nb_rows=nbrows; }
  97     Matrix(std::vector<std::map<int,T> > & matrix) :
  98       _coeffs(0), _cols(0), _is_configured(false)
  99     {
 100       _nb_rows=matrix.size();
 101       _auxiliary_matrix.resize(_nb_rows);
 102       for (int i=0; i<_nb_rows; i++)
 103         {
 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?
 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 unsigned int[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 unsigned int[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 (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)
 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       Matrix multiplies vector \a input and stores the result in
 176       vector \a output.
 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.
 180     */
 181     void multiply(const T* const input, T* const output)
 182     {
 183       if (!_is_configured)
 184         configure();
 185
 186       for (int i=0; i< _nb_rows; i++)
 187         {
 188           output[i]=0.;
 189           for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
 190             {
 191               int icol = _cols[j];
 192               output[i]+=input[icol]*_coeffs[j];
 193             }
 194         }
 195     }
 196
 197     /*!
 198       Matrix multiplies vector \a input and stores the result in
 199       vector \a output.
 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.
 206     */
 207     void multiply(const T* const input, T* const output, int nb_comp)
 208     {
 209       if (!_is_configured)
 210         configure();
 211
 212       for (int i=0; i< _nb_rows; i++)
 213         {
 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++)
 217             {
 218               int icol = _cols[j];
 219               for(int comp = 0; comp < nb_comp; comp++)
 220                 output[i*nb_comp+comp]+=input[icol*nb_comp+comp]*_coeffs[j];
 221             }
 222         }
 223     }
 224     /*!
 225       Transpose-multiplies vector \a input and stores the result in
 226       vector \a output.
 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.
 231     */
 232     void transposeMultiply(const T* const input, T* const output, int nb_cols)
 233     {
 234       if (!_is_configured)
 235         configure();
 236
 237       for (int icol=0; icol< nb_cols; icol++)
 238         output[icol]=0.;
 239       for (int i=0; i< _nb_rows; i++)
 240         {
 241           for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
 242             {
 243               int icol = _cols[j];
 244               output[icol]+=input[i]*_coeffs[j];
 245             }
 246         }
 247     }
 248     /*!
 249       Transpose-multiplies vector \a input and stores the result in
 250       vector \a output.
 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.
 258     */
 259     void transposeMultiply(const T* const input, T* const output, int nb_cols, int nb_comp)
 260     {
 261       if (!_is_configured)
 262         configure();
 263
 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.;
 267
 268       for (int i=0; i< _nb_rows; i++)
 269         {
 270           for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
 271             {
 272               int icol = _cols[j];
 273               for(int comp = 0; comp < nb_comp; comp++)
 274                 output[icol*nb_comp+comp]+=input[i*nb_comp+comp]*_coeffs[j];
 275             }
 276         }
 277     }
 278     /*
 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
 282     */
 283     void colSum(std::vector< T >& output, int nb_cols)
 284     {
 285       if (!_is_configured)
 286         configure();
 287       for (int icol=0; icol< nb_cols; icol++)
 288         output[icol]=0.;
 289       for (int i=0; i< _nb_rows; i++)
 290         {
 291           for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
 292             {
 293               int icol = _cols[j];
 294               output[icol]+=_coeffs[j];
 295             }
 296         }
 297     }
 298
 299     /*
 300       Sums the coefficients of each row of the matrix
 301       The vector output must be dimensioned to _nb_rows
 302     */
 303     void rowSum(std::vector< T >& output)
 304     {
 305       if (!_is_configured)
 306         configure();
 307       for (int i=0; i< _nb_rows; i++)
 308         {
 309           output[i]=0;
 310           for (unsigned int j=_ncols_offset[i]; j< _ncols_offset[i+1]; j++)
 311             output[i]+=_coeffs[j];
 312         }
 313     }
 314
 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 */
 319
 320     void configure()
 321     {
 322       _ncols_offset.resize(_nb_rows+1);
 323       _ncols_offset[0]=0;
 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++)
 332         {
 333           for (unsigned int j=0; j<_auxiliary_matrix[i].size(); j++)
 334             {
 335               *cols_ptr++ = OTT<int,type>::ind2C(_auxiliary_matrix[i][j].first);
 336               *coeffs_ptr++ = _auxiliary_matrix[i][j].second;
 337             }
 338         }
 339       _auxiliary_matrix.clear();
 340       _is_configured=true;
 341     }
 342
 343     /*!
 344      * 0 <= irow < n
 345      */
 346     Row &operator [] (unsigned int irow)
 347     {
 348       return _auxiliary_matrix[irow];
 349     }
 350
 351     int getNbRows()
 352     {
 353       return _nb_rows;
 354     }
 355
 356   };
 357
 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
 365
 366     for instance, matrix
 367     | 1.0 0.0 0.5 |
 368     | 0.0 1.0 0.0 |
 369     | 0.2 0.0 1.0 |
 370     will be displayed in 0-indexing as
 371     0
 372     3
 373     2 0 1.0 2 0.5
 374     1 1 1.0
 375     2 0 0.2 2 1.0
 376   */
 377
 378   template<class T, NumberingPolicy type>
 379   std::ostream& operator<<(std::ostream& out, const Matrix<T, type>& m)
 380   {
 381     if (m._is_configured)
 382       {
 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++)
 386           {
 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];
 390             out<<std::endl;
 391           }
 392       }
 393     else
 394       {
 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++)
 398           {
 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;
 403             out <<"\n";
 404           }
 405       }
 406     return out;
 407   }
 408
 409   template<class T, NumberingPolicy type>
 410   std::istream& operator>>(std::istream& in, Matrix<T,type>& m)
 411   {
 412     int index_base_test;
 413     in >> index_base_test;
 414     if (index_base_test!=OTT<unsigned int,type>::indFC(0))
 415       {
 416         std::cerr << "file index is "<<index_base_test<<std::endl;
 417         throw Exception("incompatible indexing reading matrix");
 418       }
 419     in >> m._nb_rows;
 420     m._auxiliary_matrix.resize(m._nb_rows);
 421     for (unsigned int i=0; i<m._nb_rows; i++)
 422       {
 423         unsigned int ncols;
 424         in >> ncols;
 425         m._auxiliary_matrix[i].resize(ncols);
 426         double value;
 427         unsigned int col;
 428         for (unsigned int j=0; j<ncols; j++)
 429           {
 430             in>>col;
 431             in>>value;
 432             m._auxiliary_matrix[i].push_back(std::make_pair(col, value));
 433           }
 434       }
 435     return in;
 436   }
 437 }
 438
 439 #endif