src/medtool/src/MEDCoupling/MEDCouplingMatrix.cxx

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