src/MEDCoupling/MEDCouplingMatrix.cxx

   1 // Copyright (C) 2007-2016  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 // Author : Anthony Geay
  20
  21 #include "MEDCouplingMatrix.hxx"
  22
  23 #include "InterpKernelMatrixTools.hxx"
  24
  25 #include <sstream>
  26
  27 using namespace MEDCoupling;
  28
  29 DenseMatrix *DenseMatrix::New(int nbRows, int nbCols)
  30 {
  31   return new DenseMatrix(nbRows,nbCols);
  32 }
  33
  34 DenseMatrix *DenseMatrix::New(DataArrayDouble *array, int nbRows, int nbCols)
  35 {
  36   return new DenseMatrix(array,nbRows,nbCols);
  37 }
  38
  39 DenseMatrix *DenseMatrix::deepCopy() const
  40 {
  41   MCAuto<DataArrayDouble> arr(getData()->deepCopy());
  42   MCAuto<DenseMatrix> ret(DenseMatrix::New(arr,getNumberOfRows(),getNumberOfCols()));
  43   return ret.retn();
  44 }
  45
  46 DenseMatrix *DenseMatrix::shallowCpy() const
  47 {
  48   MCAuto<DenseMatrix> ret(DenseMatrix::New(const_cast<DataArrayDouble *>(getData()),getNumberOfRows(),getNumberOfCols()));
  49   return ret.retn();
  50 }
  51
  52 std::size_t DenseMatrix::getHeapMemorySizeWithoutChildren() const
  53 {
  54   return sizeof(DenseMatrix);
  55 }
  56
  57 std::vector<const BigMemoryObject *> DenseMatrix::getDirectChildrenWithNull() const
  58 {
  59   std::vector<const BigMemoryObject *> ret;
  60   ret.push_back((const DataArrayDouble *)_data);
  61   return ret;
  62 }
  63
  64 void DenseMatrix::updateTime() const
  65 {
  66   const DataArrayDouble *pt(_data);
  67   if(pt)
  68     updateTimeWith(*pt);
  69 }
  70
  71 /*!
  72  * This method scratch \a this to use a new input. The shape of \a this can be modified freely without any constraints.
  73  *
  74  * \param [in] array - The array containing data that is expected to be taken as new data.
  75  * \param [in] nbRows - The new number of rows (>0 or -1). If -1, the current number of rows will be taken.
  76  * \param [in] nbCols - The new number of columns (>0 or -1). If -1, the current number of cols will be taken.
  77  *
  78  * \sa reShape
  79  */
  80 void DenseMatrix::reBuild(DataArrayDouble *array, int nbRows, int nbCols)
  81 {
  82   int nbr(getNumberOfRowsExt(nbRows)),nbc(getNumberOfColsExt(nbCols));
  83   CheckArraySizes(array,nbr,nbc);
  84   DataArrayDouble *data(_data);
  85   if(data!=array)
  86     {
  87       _data=array; _data->incrRef();
  88       declareAsNew();
  89     }
  90   if(nbr!=_nb_rows)
  91     {
  92       _nb_rows=nbr;
  93       declareAsNew();
  94     }
  95   if(nbc!=_nb_cols)
  96     {
  97       _nb_cols=nbc;
  98       declareAsNew();
  99     }
 100 }
 101
 102 /*!
 103  * This method does \b not change the content of the data in \a this. It only changes the shape (with a same number of elements in the matrix).
 104  * If the number of elements needs to be changed call reBuild method instead.
 105  *
 106  * \param [in] nbRows - The new number of rows (>0)
 107  * \param [in] nbCols - The new number of columns (>0)
 108  * \throw if the \c nbRows*nbCols is not equal to \c this->getNbOfElems()
 109  * \sa reBuild
 110  */
 111 void DenseMatrix::reShape(int nbRows, int nbCols)
 112 {
 113   if(nbRows<0 || nbCols<0)
 114     throw INTERP_KERNEL::Exception("DenseMatrix::reShape : number of rows and number of cols must be > 0 both !");
 115   if(nbRows*nbCols!=getNbOfElems())
 116     throw INTERP_KERNEL::Exception("DenseMatrix::reShape : This method is designed to change only the shape ! Number of elements must remain the same !");
 117   if(_nb_rows!=nbRows)
 118     {
 119       _nb_rows=nbRows;
 120       declareAsNew();
 121     }
 122   if(_nb_cols!=nbCols)
 123     {
 124       _nb_cols=nbCols;
 125       declareAsNew();
 126     }
 127 }
 128
 129 void DenseMatrix::transpose()
 130 {
 131   const MemArray<double>& mem(getData()->accessToMemArray());
 132   double *pt(mem.toNoInterlace(getNumberOfCols()));
 133   std::copy(pt,pt+getNbOfElems(),getData()->getPointer());//declareAsNew done here automatically by getPointer
 134   free(pt);
 135   std::swap(_nb_rows,_nb_cols);
 136   updateTime();
 137 }
 138
 139 bool DenseMatrix::isEqual(const DenseMatrix& other, double eps) const
 140 {
 141   std::string tmp;
 142   return isEqualIfNotWhy(other,eps,tmp);
 143 }
 144
 145 bool DenseMatrix::isEqualIfNotWhy(const DenseMatrix& other, double eps, std::string& reason) const
 146 {
 147   if(_nb_rows!=other._nb_rows)
 148     {
 149       std::ostringstream oss; oss << "Number of rows differs (" << _nb_rows << "!=" << other._nb_rows << ") !";
 150       reason+=oss.str();
 151       return false;
 152     }
 153   if(_nb_cols!=other._nb_cols)
 154       {
 155         std::ostringstream oss; oss << "Number of cols differs (" << _nb_cols << "!=" << other._nb_cols << ") !";
 156         reason+=oss.str();
 157         return false;
 158       }
 159   std::string tmp1;
 160   if(!_data->isEqualIfNotWhy(*other._data,eps,tmp1))
 161     {
 162       reason+="Data differs : "+tmp1;
 163       return false;
 164     }
 165   return true;
 166 }
 167
 168 DataArrayDouble *DenseMatrix::matVecMult(const DataArrayDouble *vec) const
 169 {
 170   return MatVecMult(this,vec);
 171 }
 172
 173 DataArrayDouble *DenseMatrix::MatVecMult(const DenseMatrix *mat, const DataArrayDouble *vec)
 174 {
 175   if(!mat || !vec)
 176     throw INTERP_KERNEL::Exception("DenseMatrix::MatVecMult : input matrix or vec is NULL !");
 177   vec->checkAllocated();
 178   if(vec->getNumberOfComponents()!=1)
 179     throw INTERP_KERNEL::Exception("DenseMatrix::MatVecMult : input vector must have only one component !");
 180   if(vec->getNumberOfTuples()!=mat->getNumberOfCols())
 181     throw INTERP_KERNEL::Exception("DenseMatrix::MatVecMult : Number of columns of this must be equal to number of tuples of vec !");
 182   MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(mat->getNumberOfRows(),1);
 183   INTERP_KERNEL::matrixProduct(mat->getData()->begin(),mat->getNumberOfRows(),mat->getNumberOfCols(),vec->begin(),vec->getNumberOfTuples(),1,ret->getPointer());
 184   return ret.retn();
 185 }
 186
 187 DenseMatrix *DenseMatrix::Add(const DenseMatrix *a1, const DenseMatrix *a2)
 188 {
 189   if(!a1 || !a2)
 190     throw INTERP_KERNEL::Exception("DenseMatrix::Add : input matrices must be not NULL !");
 191   CheckSameSize(a1,a2);
 192   MCAuto<DataArrayDouble> data(DataArrayDouble::Add(a1->getData(),a2->getData()));
 193   MCAuto<DenseMatrix> ret(DenseMatrix::New(data,a1->getNumberOfRows(),a1->getNumberOfCols()));
 194   return ret.retn();
 195 }
 196
 197 void DenseMatrix::addEqual(const DenseMatrix *other)
 198 {
 199   if(!other)
 200     throw INTERP_KERNEL::Exception("DenseMatrix::addEqual : other must be not NULL !");
 201   CheckSameSize(this,other);
 202   getData()->addEqual(other->getData());
 203 }
 204
 205 DenseMatrix *DenseMatrix::Substract(const DenseMatrix *a1, const DenseMatrix *a2)
 206 {
 207   if(!a1 || !a2)
 208     throw INTERP_KERNEL::Exception("DenseMatrix::Substract : input matrices must be not NULL !");
 209   CheckSameSize(a1,a2);
 210   MCAuto<DataArrayDouble> data(DataArrayDouble::Substract(a1->getData(),a2->getData()));
 211   MCAuto<DenseMatrix> ret(DenseMatrix::New(data,a1->getNumberOfRows(),a1->getNumberOfCols()));
 212   return ret.retn();
 213 }
 214
 215 void DenseMatrix::substractEqual(const DenseMatrix *other)
 216 {
 217   if(!other)
 218     throw INTERP_KERNEL::Exception("DenseMatrix::substractEqual : other must be not NULL !");
 219   CheckSameSize(this,other);
 220   getData()->substractEqual(other->getData());
 221 }
 222
 223 DenseMatrix *DenseMatrix::Multiply(const DenseMatrix *a1, const DenseMatrix *a2)
 224 {
 225   if(!a1 || !a2)
 226     throw INTERP_KERNEL::Exception("DenseMatrix::Multiply : input matrices must be not NULL !");
 227   CheckCompatibleSizeForMul(a1,a2);
 228   int nbr(a1->getNumberOfRows()),nbc(a2->getNumberOfCols());
 229   MCAuto<DataArrayDouble> data(DataArrayDouble::New()); data->alloc(nbr*nbc,1);
 230   MCAuto<DenseMatrix> ret(DenseMatrix::New(data,a1->getNumberOfRows(),a2->getNumberOfCols()));
 231   INTERP_KERNEL::matrixProduct(a1->getData()->begin(),a1->getNumberOfRows(),a1->getNumberOfCols(),a2->getData()->begin(),a2->getNumberOfRows(),a2->getNumberOfCols(),data->getPointer());
 232   return ret.retn();
 233 }
 234
 235 DenseMatrix *DenseMatrix::Multiply(const DenseMatrix *a1, const DataArrayDouble *a2)
 236 {
 237   if(!a1 || !a2 || !a2->isAllocated())
 238     throw INTERP_KERNEL::Exception("DenseMatrix::Multiply #2 : input matrices must be not NULL and a2 allocated !");
 239   if(a2->getNumberOfComponents()!=1)
 240     throw INTERP_KERNEL::Exception("DenseMatrix::Multiply #2 : The 2nd member must have exactly one component !");
 241   MCAuto<DenseMatrix> a2Bis(DenseMatrix::New(const_cast<DataArrayDouble *>(a2),a2->getNumberOfTuples(),1));
 242   return DenseMatrix::Multiply(a1,a2Bis);
 243 }
 244
 245 DenseMatrix::~DenseMatrix()
 246 {
 247 }
 248
 249 DenseMatrix::DenseMatrix(int nbRows, int nbCols):_nb_rows(nbRows),_nb_cols(nbCols),_data(DataArrayDouble::New())
 250 {
 251   if(_nb_rows<0 || _nb_cols<0)
 252     throw INTERP_KERNEL::Exception("constructor of DenseMatrix : number of rows and number of cols must be > 0 both !");
 253   int nbOfTuples(_nb_rows*_nb_cols);
 254   _data->alloc(nbOfTuples,1);
 255 }
 256
 257 DenseMatrix::DenseMatrix(DataArrayDouble *array, int nbRows, int nbCols):_nb_rows(nbRows),_nb_cols(nbCols)
 258 {
 259   CheckArraySizes(array,_nb_rows,_nb_cols);
 260   _data=array; _data->incrRef();
 261 }
 262
 263 int DenseMatrix::getNumberOfRowsExt(int nbRows) const
 264 {
 265   if(nbRows<-1)
 266     throw INTERP_KERNEL::Exception("DenseMatrix::getNumberOfRowsExt : invalid input must be >= -1 !");
 267   if(nbRows==-1)
 268     return _nb_rows;
 269   else
 270     return nbRows;
 271 }
 272
 273 int DenseMatrix::getNumberOfColsExt(int nbCols) const
 274 {
 275   if(nbCols<-1)
 276     throw INTERP_KERNEL::Exception("DenseMatrix::getNumberOfColsExt : invalid input must be >= -1 !");
 277   if(nbCols==-1)
 278     return _nb_cols;
 279   else
 280     return nbCols;
 281 }
 282
 283 void DenseMatrix::checkValidData() const
 284 {
 285   if(!getData())
 286     throw INTERP_KERNEL::Exception("DenseMatrix::checkValidData : data is NULL !");
 287   if(!getData()->isAllocated())
 288     throw INTERP_KERNEL::Exception("DenseMatrix::checkValidData : data is not allocated !");
 289   if(getData()->getNumberOfComponents()!=1)
 290     throw INTERP_KERNEL::Exception("DenseMatrix::checkValidData : data has not 1 component !");
 291 }
 292
 293 void DenseMatrix::CheckArraySizes(DataArrayDouble *array, int nbRows, int nbCols)
 294 {
 295   if(nbRows<0 || nbCols<0)
 296     throw INTERP_KERNEL::Exception("constructor #2 of DenseMatrix : number of rows and number of cols must be > 0 both !");
 297   if(!array || !array->isAllocated())
 298     throw INTERP_KERNEL::Exception("constructor #2 of DenseMatrix : input array is empty or not allocated !");
 299   if(array->getNumberOfComponents()!=1)
 300     throw INTERP_KERNEL::Exception("constructor #2 of DenseMatrix : input array must have exactly one component !");
 301   std::size_t nbr((std::size_t)nbRows),nbc((std::size_t)nbCols);
 302   if(nbr*nbc!=array->getNbOfElems())
 303     throw INTERP_KERNEL::Exception("constructor #2 of DenseMatrix : the number of elems in input array is not equal to the product of nbRows and nbCols !");
 304 }
 305
 306 void DenseMatrix::CheckSameSize(const DenseMatrix *a1, const DenseMatrix *a2)
 307 {
 308   if(!a1 || !a2)
 309     throw INTERP_KERNEL::Exception("DenseMatrix::CheckSameSize : a1 or a2 is NULL !");
 310   a1->checkValidData(); a2->checkValidData();
 311   if(a1->getNumberOfRows()!=a2->getNumberOfRows())
 312     throw INTERP_KERNEL::Exception("DenseMatrix::CheckSameSize : number of rows mismatches !");
 313   if(a1->getNumberOfCols()!=a2->getNumberOfCols())
 314     throw INTERP_KERNEL::Exception("DenseMatrix::CheckSameSize : number of columns mismatches !");
 315
 316 }
 317
 318 void DenseMatrix::CheckCompatibleSizeForMul(const DenseMatrix *a1, const DenseMatrix *a2)
 319 {
 320   if(!a1 || !a2)
 321       throw INTERP_KERNEL::Exception("DenseMatrix::CheckCompatibleSizeForMul : a1 or a2 is NULL !");
 322   a1->checkValidData(); a2->checkValidData();
 323   if(a1->getNumberOfCols()!=a2->getNumberOfRows())
 324     throw INTERP_KERNEL::Exception("DenseMatrix::CheckCompatibleSizeForMul : number of cols of a1 must be equal to number of rows of a2 !");
 325 }
 326