1 // Copyright (C) 2007-2021 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay
21 #include "MEDCouplingMatrix.hxx"
23 #include "InterpKernelMatrixTools.hxx"
27 using namespace MEDCoupling;
29 DenseMatrix *DenseMatrix::New(mcIdType nbRows, mcIdType nbCols)
31 return new DenseMatrix(nbRows,nbCols);
34 DenseMatrix *DenseMatrix::New(DataArrayDouble *array, mcIdType nbRows, mcIdType nbCols)
36 return new DenseMatrix(array,nbRows,nbCols);
39 DenseMatrix *DenseMatrix::deepCopy() const
41 MCAuto<DataArrayDouble> arr(getData()->deepCopy());
42 MCAuto<DenseMatrix> ret(DenseMatrix::New(arr,getNumberOfRows(),getNumberOfCols()));
46 DenseMatrix *DenseMatrix::shallowCpy() const
48 MCAuto<DenseMatrix> ret(DenseMatrix::New(const_cast<DataArrayDouble *>(getData()),getNumberOfRows(),getNumberOfCols()));
52 std::size_t DenseMatrix::getHeapMemorySizeWithoutChildren() const
54 return sizeof(DenseMatrix);
57 std::vector<const BigMemoryObject *> DenseMatrix::getDirectChildrenWithNull() const
59 std::vector<const BigMemoryObject *> ret;
60 ret.push_back((const DataArrayDouble *)_data);
64 void DenseMatrix::updateTime() const
66 const DataArrayDouble *pt(_data);
72 * This method scratch \a this to use a new input. The shape of \a this can be modified freely without any constraints.
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.
80 void DenseMatrix::reBuild(DataArrayDouble *array, mcIdType nbRows, mcIdType nbCols)
82 mcIdType nbr(getNumberOfRowsExt(nbRows)),nbc(getNumberOfColsExt(nbCols));
83 CheckArraySizes(array,nbr,nbc);
84 DataArrayDouble *data(_data);
87 _data=array; _data->incrRef();
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.
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()
111 void DenseMatrix::reShape(mcIdType nbRows, mcIdType nbCols)
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 !");
129 void DenseMatrix::transpose()
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
135 std::swap(_nb_rows,_nb_cols);
139 bool DenseMatrix::isEqual(const DenseMatrix& other, double eps) const
142 return isEqualIfNotWhy(other,eps,tmp);
145 bool DenseMatrix::isEqualIfNotWhy(const DenseMatrix& other, double eps, std::string& reason) const
147 if(_nb_rows!=other._nb_rows)
149 std::ostringstream oss; oss << "Number of rows differs (" << _nb_rows << "!=" << other._nb_rows << ") !";
153 if(_nb_cols!=other._nb_cols)
155 std::ostringstream oss; oss << "Number of cols differs (" << _nb_cols << "!=" << other._nb_cols << ") !";
160 if(!_data->isEqualIfNotWhy(*other._data,eps,tmp1))
162 reason+="Data differs : "+tmp1;
168 DataArrayDouble *DenseMatrix::matVecMult(const DataArrayDouble *vec) const
170 return MatVecMult(this,vec);
173 DataArrayDouble *DenseMatrix::MatVecMult(const DenseMatrix *mat, const DataArrayDouble *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());
187 DenseMatrix *DenseMatrix::Add(const DenseMatrix *a1, const DenseMatrix *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()));
197 void DenseMatrix::addEqual(const DenseMatrix *other)
200 throw INTERP_KERNEL::Exception("DenseMatrix::addEqual : other must be not NULL !");
201 CheckSameSize(this,other);
202 getData()->addEqual(other->getData());
205 DenseMatrix *DenseMatrix::Substract(const DenseMatrix *a1, const DenseMatrix *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()));
215 void DenseMatrix::substractEqual(const DenseMatrix *other)
218 throw INTERP_KERNEL::Exception("DenseMatrix::substractEqual : other must be not NULL !");
219 CheckSameSize(this,other);
220 getData()->substractEqual(other->getData());
223 DenseMatrix *DenseMatrix::Multiply(const DenseMatrix *a1, const DenseMatrix *a2)
226 throw INTERP_KERNEL::Exception("DenseMatrix::Multiply : input matrices must be not NULL !");
227 CheckCompatibleSizeForMul(a1,a2);
228 mcIdType 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());
235 DenseMatrix *DenseMatrix::Multiply(const DenseMatrix *a1, const DataArrayDouble *a2)
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);
245 DenseMatrix::~DenseMatrix()
249 DenseMatrix::DenseMatrix(mcIdType nbRows, mcIdType nbCols):_nb_rows(nbRows),_nb_cols(nbCols),_data(DataArrayDouble::New())
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 mcIdType nbOfTuples(_nb_rows*_nb_cols);
254 _data->alloc(nbOfTuples,1);
257 DenseMatrix::DenseMatrix(DataArrayDouble *array, mcIdType nbRows, mcIdType nbCols):_nb_rows(nbRows),_nb_cols(nbCols)
259 CheckArraySizes(array,_nb_rows,_nb_cols);
260 _data=array; _data->incrRef();
263 mcIdType DenseMatrix::getNumberOfRowsExt(mcIdType nbRows) const
266 throw INTERP_KERNEL::Exception("DenseMatrix::getNumberOfRowsExt : invalid input must be >= -1 !");
273 mcIdType DenseMatrix::getNumberOfColsExt(mcIdType nbCols) const
276 throw INTERP_KERNEL::Exception("DenseMatrix::getNumberOfColsExt : invalid input must be >= -1 !");
283 void DenseMatrix::checkValidData() const
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 !");
293 void DenseMatrix::CheckArraySizes(DataArrayDouble *array, mcIdType nbRows, mcIdType nbCols)
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 if(nbRows*nbCols!=array->getNbOfElems())
302 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 !");
305 void DenseMatrix::CheckSameSize(const DenseMatrix *a1, const DenseMatrix *a2)
308 throw INTERP_KERNEL::Exception("DenseMatrix::CheckSameSize : a1 or a2 is NULL !");
309 a1->checkValidData(); a2->checkValidData();
310 if(a1->getNumberOfRows()!=a2->getNumberOfRows())
311 throw INTERP_KERNEL::Exception("DenseMatrix::CheckSameSize : number of rows mismatches !");
312 if(a1->getNumberOfCols()!=a2->getNumberOfCols())
313 throw INTERP_KERNEL::Exception("DenseMatrix::CheckSameSize : number of columns mismatches !");
317 void DenseMatrix::CheckCompatibleSizeForMul(const DenseMatrix *a1, const DenseMatrix *a2)
320 throw INTERP_KERNEL::Exception("DenseMatrix::CheckCompatibleSizeForMul : a1 or a2 is NULL !");
321 a1->checkValidData(); a2->checkValidData();
322 if(a1->getNumberOfCols()!=a2->getNumberOfRows())
323 throw INTERP_KERNEL::Exception("DenseMatrix::CheckCompatibleSizeForMul : number of cols of a1 must be equal to number of rows of a2 !");