]> SALOME platform Git repositories - modules/med.git/blob - medtool/src/MEDCoupling/MEDCouplingMatrix.cxx
Salome HOME
Move medtool folder to MED base repository
[modules/med.git] / 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