From 6f4704cbb6e4b264c28bd3c3eb693fbfe77c83dd Mon Sep 17 00:00:00 2001 From: geay Date: Fri, 16 May 2014 13:51:51 +0200 Subject: [PATCH] First draft for dense matrix --- src/MEDCoupling/CMakeLists.txt | 1 + src/MEDCoupling/MEDCouplingFieldDouble.cxx | 14 +- src/MEDCoupling/MEDCouplingFieldDouble.hxx | 1 + src/MEDCoupling/MEDCouplingMatrix.cxx | 326 ++++++++++++++++++ src/MEDCoupling/MEDCouplingMatrix.hxx | 86 +++++ src/MEDCoupling/MEDCouplingMemArray.cxx | 4 +- .../MEDCouplingTimeDiscretization.cxx | 11 + .../MEDCouplingTimeDiscretization.hxx | 1 + src/MEDCoupling_Swig/MEDCoupling.i | 6 + src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 65 +++- src/MEDCoupling_Swig/MEDCouplingCommon.i | 98 ++++++ src/MEDCoupling_Swig/MEDCouplingFinalize.i | 5 + src/MEDCoupling_Swig/MEDCouplingRemapper.i | 6 + src/MEDLoader/Swig/MEDLoader.i | 6 + 14 files changed, 625 insertions(+), 5 deletions(-) create mode 100644 src/MEDCoupling/MEDCouplingMatrix.cxx create mode 100644 src/MEDCoupling/MEDCouplingMatrix.hxx diff --git a/src/MEDCoupling/CMakeLists.txt b/src/MEDCoupling/CMakeLists.txt index 6ac3fbd86..5a3aef6b7 100644 --- a/src/MEDCoupling/CMakeLists.txt +++ b/src/MEDCoupling/CMakeLists.txt @@ -56,6 +56,7 @@ SET(medcoupling_SOURCES MEDCouplingDefinitionTime.cxx MEDCouplingFieldOverTime.cxx MEDCouplingCartesianAMRMesh.cxx + MEDCouplingMatrix.cxx ) SET(medcouplingremapper_SOURCES diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx index 3c2e2eb8c..cfa138f3a 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx @@ -1445,7 +1445,7 @@ void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double time, dou } /*! - * Apply a liner function to a given component of \a this field, so that + * Apply a linear function to a given component of \a this field, so that * a component value (x) becomes \f$ a * x + b \f$. * \param [in] a - the first coefficient of the function. * \param [in] b - the second coefficient of the function. @@ -1457,6 +1457,18 @@ void MEDCouplingFieldDouble::applyLin(double a, double b, int compoId) _time_discr->applyLin(a,b,compoId); } +/*! + * Apply a linear function to all components of \a this field, so that + * values (x) becomes \f$ a * x + b \f$. + * \param [in] a - the first coefficient of the function. + * \param [in] b - the second coefficient of the function. + * \throw If the data array(s) is(are) not set. + */ +void MEDCouplingFieldDouble::applyLin(double a, double b) +{ + _time_discr->applyLin(a,b); +} + /*! * This method sets \a this to a uniform scalar field with one component. * All tuples will have the same value 'value'. diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.hxx b/src/MEDCoupling/MEDCouplingFieldDouble.hxx index 9d0394a20..dec30c8bf 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.hxx @@ -115,6 +115,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void getValueOn(const double *spaceLoc, double time, double *res) const; MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const double *spaceLoc, int nbOfPoints) const; MEDCOUPLING_EXPORT void applyLin(double a, double b, int compoId); + MEDCOUPLING_EXPORT void applyLin(double a, double b); MEDCOUPLING_EXPORT MEDCouplingFieldDouble &operator=(double value); MEDCOUPLING_EXPORT void fillFromAnalytic(int nbOfComp, FunctionToEvaluate func); MEDCOUPLING_EXPORT void fillFromAnalytic(int nbOfComp, const std::string& func); diff --git a/src/MEDCoupling/MEDCouplingMatrix.cxx b/src/MEDCoupling/MEDCouplingMatrix.cxx new file mode 100644 index 000000000..7500d226e --- /dev/null +++ b/src/MEDCoupling/MEDCouplingMatrix.cxx @@ -0,0 +1,326 @@ +// Copyright (C) 2007-2014 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay + +#include "MEDCouplingMatrix.hxx" + +#include "InterpKernelMatrixTools.hxx" + +using namespace ParaMEDMEM; + +DenseMatrix *DenseMatrix::New(int nbRows, int nbCols) +{ + return new DenseMatrix(nbRows,nbCols); +} + +DenseMatrix *DenseMatrix::New(DataArrayDouble *array, int nbRows, int nbCols) +{ + return new DenseMatrix(array,nbRows,nbCols); +} + +DenseMatrix *DenseMatrix::deepCpy() const +{ + MEDCouplingAutoRefCountObjectPtr arr(getData()->deepCpy()); + MEDCouplingAutoRefCountObjectPtr ret(DenseMatrix::New(arr,getNumberOfRows(),getNumberOfCols())); + return ret.retn(); +} + +DenseMatrix *DenseMatrix::shallowCpy() const +{ + MEDCouplingAutoRefCountObjectPtr ret(DenseMatrix::New(const_cast(getData()),getNumberOfRows(),getNumberOfCols())); + return ret.retn(); +} + +std::size_t DenseMatrix::getHeapMemorySizeWithoutChildren() const +{ + return sizeof(DenseMatrix); +} + +std::vector DenseMatrix::getDirectChildren() const +{ + std::vector ret; + const DataArrayDouble *pt(_data); + if(pt) + ret.push_back(pt); + return ret; +} + +void DenseMatrix::updateTime() const +{ + const DataArrayDouble *pt(_data); + if(pt) + updateTimeWith(*pt); +} + +/*! + * This method scratch \a this to use a new input. The shape of \a this can be modified freely without any constraints. + * + * \param [in] array - The array containing data that is expected to be taken as new data. + * \param [in] nbRows - The new number of rows (>0 or -1). If -1, the current number of rows will be taken. + * \param [in] nbCols - The new number of columns (>0 or -1). If -1, the current number of cols will be taken. + * + * \sa reShape + */ +void DenseMatrix::reBuild(DataArrayDouble *array, int nbRows, int nbCols) +{ + int nbr(getNumberOfRowsExt(nbRows)),nbc(getNumberOfColsExt(nbCols)); + CheckArraySizes(array,nbr,nbc); + DataArrayDouble *data(_data); + if(data!=array) + { + _data=array; _data->incrRef(); + declareAsNew(); + } + if(nbr!=_nb_rows) + { + _nb_rows=nbr; + declareAsNew(); + } + if(nbc!=_nb_cols) + { + _nb_cols=nbc; + declareAsNew(); + } +} + +/*! + * 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). + * If the number of elements needs to be changed call reBuild method instead. + * + * \param [in] nbRows - The new number of rows (>0) + * \param [in] nbCols - The new number of columns (>0) + * \throw if the \c nbRows*nbCols is not equal to \c this->getNbOfElems() + * \sa reBuild + */ +void DenseMatrix::reShape(int nbRows, int nbCols) +{ + if(nbRows<0 || nbCols<0) + throw INTERP_KERNEL::Exception("DenseMatrix::reShape : number of rows and number of cols must be > 0 both !"); + if(nbRows*nbCols!=getNbOfElems()) + throw INTERP_KERNEL::Exception("DenseMatrix::reShape : This method is designed to change only the shape ! Number of elements must remain the same !"); + if(_nb_rows!=nbRows) + { + _nb_rows=nbRows; + declareAsNew(); + } + if(_nb_cols!=nbCols) + { + _nb_cols=nbCols; + declareAsNew(); + } +} + +void DenseMatrix::transpose() +{ + const MemArray& mem(getData()->accessToMemArray()); + double *pt(mem.toNoInterlace(getNumberOfCols())); + std::copy(pt,pt+getNbOfElems(),getData()->getPointer());//declareAsNew done here automatically by getPointer + free(pt); + std::swap(_nb_rows,_nb_cols); + updateTime(); +} + +bool DenseMatrix::isEqual(const DenseMatrix& other, double eps) const +{ + std::string tmp; + return isEqualIfNotWhy(other,eps,tmp); +} + +bool DenseMatrix::isEqualIfNotWhy(const DenseMatrix& other, double eps, std::string& reason) const +{ + if(_nb_rows!=other._nb_rows) + { + std::ostringstream oss; oss << "Number of rows differs (" << _nb_rows << "!=" << other._nb_rows << ") !"; + reason+=oss.str(); + return false; + } + if(_nb_cols!=other._nb_cols) + { + std::ostringstream oss; oss << "Number of cols differs (" << _nb_cols << "!=" << other._nb_cols << ") !"; + reason+=oss.str(); + return false; + } + std::string tmp1; + if(!_data->isEqualIfNotWhy(*other._data,eps,tmp1)) + { + reason+="Data differs : "+tmp1; + return false; + } + return true; +} + +DataArrayDouble *DenseMatrix::matVecMult(const DataArrayDouble *vec) const +{ + return MatVecMult(this,vec); +} + +DataArrayDouble *DenseMatrix::MatVecMult(const DenseMatrix *mat, const DataArrayDouble *vec) +{ + if(!mat || !vec) + throw INTERP_KERNEL::Exception("DenseMatrix::MatVecMult : input matrix or vec is NULL !"); + vec->checkAllocated(); + if(vec->getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("DenseMatrix::MatVecMult : input vector must have only one component !"); + if(vec->getNumberOfTuples()!=mat->getNumberOfCols()) + throw INTERP_KERNEL::Exception("DenseMatrix::MatVecMult : Number of columns of this must be equal to number of tuples of vec !"); + MEDCouplingAutoRefCountObjectPtr ret(DataArrayDouble::New()); ret->alloc(mat->getNumberOfRows(),1); + INTERP_KERNEL::matrixProduct(mat->getData()->begin(),mat->getNumberOfRows(),mat->getNumberOfCols(),vec->begin(),vec->getNumberOfTuples(),1,ret->getPointer()); + return ret.retn(); +} + +DenseMatrix *DenseMatrix::Add(const DenseMatrix *a1, const DenseMatrix *a2) +{ + if(!a1 || !a2) + throw INTERP_KERNEL::Exception("DenseMatrix::Add : input matrices must be not NULL !"); + CheckSameSize(a1,a2); + MEDCouplingAutoRefCountObjectPtr data(DataArrayDouble::Add(a1->getData(),a2->getData())); + MEDCouplingAutoRefCountObjectPtr ret(DenseMatrix::New(data,a1->getNumberOfRows(),a1->getNumberOfCols())); + return ret.retn(); +} + +void DenseMatrix::addEqual(const DenseMatrix *other) +{ + if(!other) + throw INTERP_KERNEL::Exception("DenseMatrix::addEqual : other must be not NULL !"); + CheckSameSize(this,other); + getData()->addEqual(other->getData()); +} + +DenseMatrix *DenseMatrix::Substract(const DenseMatrix *a1, const DenseMatrix *a2) +{ + if(!a1 || !a2) + throw INTERP_KERNEL::Exception("DenseMatrix::Substract : input matrices must be not NULL !"); + CheckSameSize(a1,a2); + MEDCouplingAutoRefCountObjectPtr data(DataArrayDouble::Substract(a1->getData(),a2->getData())); + MEDCouplingAutoRefCountObjectPtr ret(DenseMatrix::New(data,a1->getNumberOfRows(),a1->getNumberOfCols())); + return ret.retn(); +} + +void DenseMatrix::substractEqual(const DenseMatrix *other) +{ + if(!other) + throw INTERP_KERNEL::Exception("DenseMatrix::substractEqual : other must be not NULL !"); + CheckSameSize(this,other); + getData()->substractEqual(other->getData()); +} + +DenseMatrix *DenseMatrix::Multiply(const DenseMatrix *a1, const DenseMatrix *a2) +{ + if(!a1 || !a2) + throw INTERP_KERNEL::Exception("DenseMatrix::Multiply : input matrices must be not NULL !"); + CheckCompatibleSizeForMul(a1,a2); + int nbr(a1->getNumberOfRows()),nbc(a2->getNumberOfCols()); + MEDCouplingAutoRefCountObjectPtr data(DataArrayDouble::New()); data->alloc(nbr*nbc,1); + MEDCouplingAutoRefCountObjectPtr ret(DenseMatrix::New(data,a1->getNumberOfRows(),a2->getNumberOfCols())); + INTERP_KERNEL::matrixProduct(a1->getData()->begin(),a1->getNumberOfRows(),a1->getNumberOfCols(),a2->getData()->begin(),a2->getNumberOfRows(),a2->getNumberOfCols(),data->getPointer()); + return ret.retn(); +} + +DenseMatrix *DenseMatrix::Multiply(const DenseMatrix *a1, const DataArrayDouble *a2) +{ + if(!a1 || !a2 || !a2->isAllocated()) + throw INTERP_KERNEL::Exception("DenseMatrix::Multiply #2 : input matrices must be not NULL and a2 allocated !"); + if(a2->getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("DenseMatrix::Multiply #2 : The 2nd member must have exactly one component !"); + MEDCouplingAutoRefCountObjectPtr a2Bis(DenseMatrix::New(const_cast(a2),a2->getNumberOfTuples(),1)); + return DenseMatrix::Multiply(a1,a2Bis); +} + +DenseMatrix::~DenseMatrix() +{ +} + +DenseMatrix::DenseMatrix(int nbRows, int nbCols):_nb_rows(nbRows),_nb_cols(nbCols),_data(DataArrayDouble::New()) +{ + if(_nb_rows<0 || _nb_cols<0) + throw INTERP_KERNEL::Exception("constructor of DenseMatrix : number of rows and number of cols must be > 0 both !"); + int nbOfTuples(_nb_rows*_nb_cols); + _data->alloc(nbOfTuples,1); +} + +DenseMatrix::DenseMatrix(DataArrayDouble *array, int nbRows, int nbCols):_nb_rows(nbRows),_nb_cols(nbCols) +{ + CheckArraySizes(array,_nb_rows,_nb_cols); + _data=array; _data->incrRef(); +} + +int DenseMatrix::getNumberOfRowsExt(int nbRows) const +{ + if(nbRows<-1) + throw INTERP_KERNEL::Exception("DenseMatrix::getNumberOfRowsExt : invalid input must be >= -1 !"); + if(nbRows==-1) + return _nb_rows; + else + return nbRows; +} + +int DenseMatrix::getNumberOfColsExt(int nbCols) const +{ + if(nbCols<-1) + throw INTERP_KERNEL::Exception("DenseMatrix::getNumberOfColsExt : invalid input must be >= -1 !"); + if(nbCols==-1) + return _nb_cols; + else + return nbCols; +} + +void DenseMatrix::checkValidData() const +{ + if(!getData()) + throw INTERP_KERNEL::Exception("DenseMatrix::checkValidData : data is NULL !"); + if(!getData()->isAllocated()) + throw INTERP_KERNEL::Exception("DenseMatrix::checkValidData : data is not allocated !"); + if(getData()->getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("DenseMatrix::checkValidData : data has not 1 component !"); +} + +void DenseMatrix::CheckArraySizes(DataArrayDouble *array, int nbRows, int nbCols) +{ + if(nbRows<0 || nbCols<0) + throw INTERP_KERNEL::Exception("constructor #2 of DenseMatrix : number of rows and number of cols must be > 0 both !"); + if(!array || !array->isAllocated()) + throw INTERP_KERNEL::Exception("constructor #2 of DenseMatrix : input array is empty or not allocated !"); + if(array->getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("constructor #2 of DenseMatrix : input array must have exactly one component !"); + std::size_t nbr((std::size_t)nbRows),nbc((std::size_t)nbCols); + if(nbr*nbc!=array->getNbOfElems()) + 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 !"); +} + +void DenseMatrix::CheckSameSize(const DenseMatrix *a1, const DenseMatrix *a2) +{ + if(!a1 || !a2) + throw INTERP_KERNEL::Exception("DenseMatrix::CheckSameSize : a1 or a2 is NULL !"); + a1->checkValidData(); a2->checkValidData(); + if(a1->getNumberOfRows()!=a2->getNumberOfRows()) + throw INTERP_KERNEL::Exception("DenseMatrix::CheckSameSize : number of rows mismatches !"); + if(a1->getNumberOfCols()!=a2->getNumberOfCols()) + throw INTERP_KERNEL::Exception("DenseMatrix::CheckSameSize : number of columns mismatches !"); + +} + +void DenseMatrix::CheckCompatibleSizeForMul(const DenseMatrix *a1, const DenseMatrix *a2) +{ + if(!a1 || !a2) + throw INTERP_KERNEL::Exception("DenseMatrix::CheckCompatibleSizeForMul : a1 or a2 is NULL !"); + a1->checkValidData(); a2->checkValidData(); + if(a1->getNumberOfCols()!=a2->getNumberOfRows()) + throw INTERP_KERNEL::Exception("DenseMatrix::CheckCompatibleSizeForMul : number of cols of a1 must be equal to number of rows of a2 !"); +} + diff --git a/src/MEDCoupling/MEDCouplingMatrix.hxx b/src/MEDCoupling/MEDCouplingMatrix.hxx new file mode 100644 index 000000000..8f79784b6 --- /dev/null +++ b/src/MEDCoupling/MEDCouplingMatrix.hxx @@ -0,0 +1,86 @@ +// Copyright (C) 2007-2014 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay + +#ifndef __PARAMEDMEM_MEDCOUPLINGMATRIX_HXX__ +#define __PARAMEDMEM_MEDCOUPLINGMATRIX_HXX__ + +#include "MEDCoupling.hxx" +#include "MEDCouplingTimeLabel.hxx" +#include "MEDCouplingRefCountObject.hxx" +#include "MEDCouplingMemArray.hxx" +#include "MEDCouplingAutoRefCountObjectPtr.hxx" + +#include "InterpKernelException.hxx" + +namespace ParaMEDMEM +{ + /*! + * The aim of this class is \b NOT to reimplement all linear algebra but only to store a dense matrix. + * It only provides basic set/get and basic operations and bindings to linear algebra libraries (numpy/scipy) and a compatible format to Petsc. + */ + class DenseMatrix : public RefCountObject, public TimeLabel + { + public: + MEDCOUPLING_EXPORT static DenseMatrix *New(int nbRows, int nbCols); + MEDCOUPLING_EXPORT static DenseMatrix *New(DataArrayDouble *array, int nbRows, int nbCols); + MEDCOUPLING_EXPORT DenseMatrix *deepCpy() const; + MEDCOUPLING_EXPORT DenseMatrix *shallowCpy() const; + MEDCOUPLING_EXPORT std::size_t getHeapMemorySizeWithoutChildren() const; + MEDCOUPLING_EXPORT std::vector getDirectChildren() const; + MEDCOUPLING_EXPORT void updateTime() const; + // + MEDCOUPLING_EXPORT int getNumberOfRows() const { return _nb_rows; } + MEDCOUPLING_EXPORT int getNumberOfCols() const { return _nb_cols; } + MEDCOUPLING_EXPORT int getNbOfElems() const { return _nb_rows*_nb_cols; } + MEDCOUPLING_EXPORT void reBuild(DataArrayDouble *array, int nbRows=-1, int nbCols=-1); + MEDCOUPLING_EXPORT void reShape(int nbRows, int nbCols); + MEDCOUPLING_EXPORT void transpose(); + // + MEDCOUPLING_EXPORT bool isEqual(const DenseMatrix& other, double eps) const; + MEDCOUPLING_EXPORT bool isEqualIfNotWhy(const DenseMatrix& other, double eps, std::string& reason) const; + MEDCOUPLING_EXPORT DataArrayDouble *matVecMult(const DataArrayDouble *vec) const; + MEDCOUPLING_EXPORT static DataArrayDouble *MatVecMult(const DenseMatrix *mat, const DataArrayDouble *vec); + MEDCOUPLING_EXPORT static DenseMatrix *Add(const DenseMatrix *a1, const DenseMatrix *a2); + MEDCOUPLING_EXPORT void addEqual(const DenseMatrix *other); + MEDCOUPLING_EXPORT static DenseMatrix *Substract(const DenseMatrix *a1, const DenseMatrix *a2); + MEDCOUPLING_EXPORT void substractEqual(const DenseMatrix *other); + MEDCOUPLING_EXPORT static DenseMatrix *Multiply(const DenseMatrix *a1, const DenseMatrix *a2); + MEDCOUPLING_EXPORT static DenseMatrix *Multiply(const DenseMatrix *a1, const DataArrayDouble *a2); + // + MEDCOUPLING_EXPORT const DataArrayDouble *getData() const { return _data; } + MEDCOUPLING_EXPORT DataArrayDouble *getData() { return _data; } + private: + ~DenseMatrix(); + DenseMatrix(int nbRows, int nbCols); + DenseMatrix(DataArrayDouble *array, int nbRows, int nbCols); + int getNumberOfRowsExt(int nbRows) const; + int getNumberOfColsExt(int nbCols) const; + void checkValidData() const; + static void CheckArraySizes(DataArrayDouble *array, int nbRows, int nbCols); + static void CheckSameSize(const DenseMatrix *a1, const DenseMatrix *a2); + static void CheckCompatibleSizeForMul(const DenseMatrix *a1, const DenseMatrix *a2); + private: + int _nb_rows; + int _nb_cols; + MEDCouplingAutoRefCountObjectPtr _data; + }; +} + +#endif diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index c54ef2cec..e2ccfd8ea 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -4108,7 +4108,7 @@ DataArrayDouble *DataArrayDouble::computeAbs() const } /*! - * Apply a liner function to a given component of \a this array, so that + * Apply a linear function to a given component of \a this array, so that * an array element (x) becomes \f$ a * x + b \f$. * \param [in] a - the first coefficient of the function. * \param [in] b - the second coefficient of the function. @@ -4131,7 +4131,7 @@ void DataArrayDouble::applyLin(double a, double b, int compoId) } /*! - * Apply a liner function to all elements of \a this array, so that + * Apply a linear function to all elements of \a this array, so that * an element _x_ becomes \f$ a * x + b \f$. * \param [in] a - the first coefficient of the function. * \param [in] b - the second coefficient of the function. diff --git a/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx b/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx index 3f27e506a..718507f27 100644 --- a/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx @@ -697,6 +697,17 @@ void MEDCouplingTimeDiscretization::applyLin(double a, double b, int compoId) } } +void MEDCouplingTimeDiscretization::applyLin(double a, double b) +{ + std::vector arrays; + getArrays(arrays); + for(std::size_t j=0;japplyLin(a,b); + } +} + void MEDCouplingTimeDiscretization::applyFunc(int nbOfComp, FunctionToEvaluate func) { std::vector arrays; diff --git a/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx b/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx index cd028f491..4996c472d 100644 --- a/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx +++ b/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx @@ -138,6 +138,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT virtual void setUniformValue(int nbOfTuple, int nbOfCompo, double value); MEDCOUPLING_EXPORT virtual void setOrCreateUniformValueOnAllComponents(int nbOfTuple, double value); MEDCOUPLING_EXPORT virtual void applyLin(double a, double b, int compoId); + MEDCOUPLING_EXPORT virtual void applyLin(double a, double b); MEDCOUPLING_EXPORT virtual void applyFunc(int nbOfComp, FunctionToEvaluate func); MEDCOUPLING_EXPORT virtual void applyFunc(int nbOfComp, const std::string& func); MEDCOUPLING_EXPORT virtual void applyFunc2(int nbOfComp, const std::string& func); diff --git a/src/MEDCoupling_Swig/MEDCoupling.i b/src/MEDCoupling_Swig/MEDCoupling.i index a7080a02a..0cc0082dc 100644 --- a/src/MEDCoupling_Swig/MEDCoupling.i +++ b/src/MEDCoupling_Swig/MEDCoupling.i @@ -95,6 +95,12 @@ def ParaMEDMEMDataArrayIntTupleIdiv(self,*args): def ParaMEDMEMDataArrayIntTupleImod(self,*args): import _MEDCoupling return _MEDCoupling.DataArrayIntTuple____imod___(self, self, *args) +def ParaMEDMEMDenseMatrixIadd(self,*args): + import _MEDCoupling + return _MEDCoupling.DenseMatrix____iadd___(self, self, *args) +def ParaMEDMEMDenseMatrixIsub(self,*args): + import _MEDCoupling + return _MEDCoupling.DenseMatrix____isub___(self, self, *args) %} %include "MEDCouplingFinalize.i" diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index a7fbdca6b..aa6ee06da 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -14992,8 +14992,8 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertTrue(f.getArray().isEqual(DataArrayDouble([60.75,52.75,52.75,44.75,52.75,44.75,44.75,36.75]),1e-12)) f.applyLin(2.,0.,0)# here it is OK ! self.assertTrue(f.getArray().isEqual(DataArrayDouble([121.5,105.5,105.5,89.5,105.5,89.5,89.5,73.5]),1e-12)) - #f.applyLin(2.,0.) - #self.assertTrue(f.getArray().isEqual(DataArrayDouble([243.,211.,211.,179.,211.,179.,179.,127.]),1e-12)) + f.applyLin(2.,0.) + self.assertTrue(f.getArray().isEqual(DataArrayDouble([243.,211.,211.,179.,211.,179.,179.,147.]),1e-12)) pass def testSwig2StructurizeMe1(self): @@ -15015,6 +15015,67 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertTrue(f.isEqual(np)) pass + def testSwig2DenseMatrix1(self): + m0=DenseMatrix(DataArrayDouble([2,3,4,5,1,6]),2,3) + self.assertEqual(m0.getNumberOfRows(),2) + self.assertEqual(m0.getNumberOfCols(),3) + self.assertEqual(m0.getNbOfElems(),6) + ref=m0.getData().getHiddenCppPointer() + m00=m0.deepCpy() + self.assertTrue(m0.isEqual(m00,1e-12)) + m00.getData().setIJ(0,0,2.1) + self.assertTrue(not m0.isEqual(m00,1e-12)) + m00.getData().setIJ(0,0,2.) + self.assertTrue(m0.isEqual(m00,1e-12)) + self.assertTrue(m0.getData().isEqual(DataArrayDouble([2,3,4,5,1,6]),1e-12)) + # + m000=m0*DataArrayDouble([5,9,3]) + self.assertTrue(m000.getData().isEqual(DataArrayDouble([49.,52.]),1e-12)) + # + m0.reShape(3,2) + self.assertTrue(not m0.isEqual(m00,1e-12)) + self.assertEqual(m0.getNumberOfRows(),3) + self.assertEqual(m0.getNumberOfCols(),2) + self.assertEqual(ref,m0.getData().getHiddenCppPointer()) + self.assertTrue(m0.getData().isEqual(DataArrayDouble([2,3,4,5,1,6]),1e-12)) + m0.reShape(2,3) + self.assertTrue(m0.isEqual(m00,1e-12)) + self.assertEqual(ref,m0.getData().getHiddenCppPointer()) + self.assertEqual(m0.getNumberOfRows(),2) + self.assertEqual(m0.getNumberOfCols(),3) + self.assertTrue(m0.getData().isEqual(DataArrayDouble([2,3,4,5,1,6]),1e-12)) + #m0np=m0.getData().toNumPyArray() ; m0np=matrix(m0np.reshape(m0.getNumberOfRows(),m0.getNumberOfCols())) + m1=m0.deepCpy() + self.assertEqual(m1.getNumberOfRows(),2) + self.assertEqual(m1.getNumberOfCols(),3) + self.assertTrue(m1.getData().isEqual(DataArrayDouble([2,3,4,5,1,6]),1e-12)) + m11=m0.deepCpy() ; m11+=m1 + self.assertEqual(m11.getNumberOfRows(),2) + self.assertEqual(m11.getNumberOfCols(),3) + self.assertTrue(m11.getData().isEqual(DataArrayDouble([4,6,8,10,2,12]),1e-12)) + m11=m11+m1 + self.assertEqual(m11.getNumberOfRows(),2) + self.assertEqual(m11.getNumberOfCols(),3) + self.assertTrue(m11.getData().isEqual(DataArrayDouble([6,9,12,15,3,18]),1e-12)) + m11=m11-m1 + self.assertEqual(m11.getNumberOfRows(),2) + self.assertEqual(m11.getNumberOfCols(),3) + self.assertTrue(m11.getData().isEqual(DataArrayDouble([4,6,8,10,2,12]),1e-12)) + m11-=m1 + self.assertEqual(m1.getNumberOfRows(),2) + self.assertEqual(m1.getNumberOfCols(),3) + self.assertTrue(m1.getData().isEqual(DataArrayDouble([2,3,4,5,1,6]),1e-12)) + m1.transpose() + self.assertEqual(m1.getNumberOfRows(),3) + self.assertEqual(m1.getNumberOfCols(),2) + self.assertTrue(m1.getData().isEqual(DataArrayDouble([2,5,3,1,4,6]),1e-12)) + #m1np=m0np.transpose() + m2=m0*m1 + self.assertEqual(m2.getNumberOfRows(),2) + self.assertEqual(m2.getNumberOfCols(),2) + self.assertTrue(m2.getData().isEqual(DataArrayDouble([29,37,37,62]),1e-12)) + pass + def setUp(self): pass pass diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index e5add4eaa..8e476cce0 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -41,6 +41,7 @@ #include "MEDCouplingDefinitionTime.hxx" #include "MEDCouplingFieldDiscretization.hxx" #include "MEDCouplingCartesianAMRMesh.hxx" +#include "MEDCouplingMatrix.hxx" #include "MEDCouplingTypemaps.i" #include "InterpKernelAutoPtr.hxx" @@ -319,6 +320,15 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::MEDCouplingCartesianAMRMesh::getFather; %newobject ParaMEDMEM::MEDCouplingCartesianAMRMesh::getPatch; %newobject ParaMEDMEM::MEDCouplingCartesianAMRMesh::__getitem__; +%newobject ParaMEDMEM::DenseMatrix::New; +%newobject ParaMEDMEM::DenseMatrix::deepCpy; +%newobject ParaMEDMEM::DenseMatrix::shallowCpy; +%newobject ParaMEDMEM::DenseMatrix::getData; +%newobject ParaMEDMEM::DenseMatrix::matVecMult; +%newobject ParaMEDMEM::DenseMatrix::MatVecMult; +%newobject ParaMEDMEM::DenseMatrix::__add__; +%newobject ParaMEDMEM::DenseMatrix::__sub__; +%newobject ParaMEDMEM::DenseMatrix::__mul__; %feature("unref") MEDCouplingPointSet "$this->decrRef();" %feature("unref") MEDCouplingMesh "$this->decrRef();" @@ -342,6 +352,7 @@ using namespace INTERP_KERNEL; %feature("unref") MEDCouplingMultiFields "$this->decrRef();" %feature("unref") MEDCouplingCartesianAMRMesh "$this->decrRef();" %feature("unref") MEDCouplingCartesianAMRPatch "$this->decrRef();" +%feature("unref") DenseMatrix "$this->decrRef();" %rename(assign) *::operator=; %ignore ParaMEDMEM::MEDCouplingGaussLocalization::pushTinySerializationIntInfo; @@ -3291,6 +3302,7 @@ namespace ParaMEDMEM void setStartTime(double val, int iteration, int order) throw(INTERP_KERNEL::Exception); void setEndTime(double val, int iteration, int order) throw(INTERP_KERNEL::Exception); void applyLin(double a, double b, int compoId) throw(INTERP_KERNEL::Exception); + void applyLin(double a, double b) throw(INTERP_KERNEL::Exception); int getNumberOfComponents() const throw(INTERP_KERNEL::Exception); int getNumberOfTuples() const throw(INTERP_KERNEL::Exception); int getNumberOfValues() const throw(INTERP_KERNEL::Exception); @@ -4753,6 +4765,92 @@ namespace ParaMEDMEM } } }; + + class DenseMatrix : public RefCountObject, public TimeLabel + { + public: + static DenseMatrix *New(int nbRows, int nbCols) throw(INTERP_KERNEL::Exception); + static DenseMatrix *New(DataArrayDouble *array, int nbRows, int nbCols) throw(INTERP_KERNEL::Exception); + DenseMatrix *deepCpy() const throw(INTERP_KERNEL::Exception); + DenseMatrix *shallowCpy() const throw(INTERP_KERNEL::Exception); + // + int getNumberOfRows() const throw(INTERP_KERNEL::Exception); + int getNumberOfCols() const throw(INTERP_KERNEL::Exception); + int getNbOfElems() const throw(INTERP_KERNEL::Exception); + void reBuild(DataArrayDouble *array, int nbRows=-1, int nbCols=-1) throw(INTERP_KERNEL::Exception); + void reShape(int nbRows, int nbCols) throw(INTERP_KERNEL::Exception); + void transpose() throw(INTERP_KERNEL::Exception); + // + bool isEqual(const DenseMatrix& other, double eps) const throw(INTERP_KERNEL::Exception); + DataArrayDouble *matVecMult(const DataArrayDouble *vec) const throw(INTERP_KERNEL::Exception); + static DataArrayDouble *MatVecMult(const DenseMatrix *mat, const DataArrayDouble *vec) throw(INTERP_KERNEL::Exception); + %extend + { + DenseMatrix(int nbRows, int nbCols) throw(INTERP_KERNEL::Exception) + { + return DenseMatrix::New(nbRows,nbCols); + } + + DenseMatrix(DataArrayDouble *array, int nbRows, int nbCols) throw(INTERP_KERNEL::Exception) + { + return DenseMatrix::New(array,nbRows,nbCols); + } + + PyObject *isEqualIfNotWhy(const DenseMatrix& other, double eps) const throw(INTERP_KERNEL::Exception) + { + std::string ret1; + bool ret0=self->isEqualIfNotWhy(other,eps,ret1); + PyObject *ret=PyTuple_New(2); + PyObject *ret0Py=ret0?Py_True:Py_False; + Py_XINCREF(ret0Py); + PyTuple_SetItem(ret,0,ret0Py); + PyTuple_SetItem(ret,1,PyString_FromString(ret1.c_str())); + return ret; + } + + DataArrayDouble *getData() throw(INTERP_KERNEL::Exception) + { + DataArrayDouble *ret(self->getData()); + if(ret) + ret->incrRef(); + return ret; + } + + DenseMatrix *__add__(const DenseMatrix *other) throw(INTERP_KERNEL::Exception) + { + return ParaMEDMEM::DenseMatrix::Add(self,other); + } + + DenseMatrix *__sub__(const DenseMatrix *other) throw(INTERP_KERNEL::Exception) + { + return ParaMEDMEM::DenseMatrix::Substract(self,other); + } + + DenseMatrix *__mul__(const DenseMatrix *other) throw(INTERP_KERNEL::Exception) + { + return ParaMEDMEM::DenseMatrix::Multiply(self,other); + } + + DenseMatrix *__mul__(const DataArrayDouble *other) throw(INTERP_KERNEL::Exception) + { + return ParaMEDMEM::DenseMatrix::Multiply(self,other); + } + + PyObject *___iadd___(PyObject *trueSelf, const DenseMatrix *other) throw(INTERP_KERNEL::Exception) + { + self->addEqual(other); + Py_XINCREF(trueSelf); + return trueSelf; + } + + PyObject *___isub___(PyObject *trueSelf, const DenseMatrix *other) throw(INTERP_KERNEL::Exception) + { + self->substractEqual(other); + Py_XINCREF(trueSelf); + return trueSelf; + } + } + }; } %pythoncode %{ diff --git a/src/MEDCoupling_Swig/MEDCouplingFinalize.i b/src/MEDCoupling_Swig/MEDCouplingFinalize.i index e54cc904d..029ab3f97 100644 --- a/src/MEDCoupling_Swig/MEDCouplingFinalize.i +++ b/src/MEDCoupling_Swig/MEDCouplingFinalize.i @@ -48,6 +48,9 @@ DataArrayIntTuple.__imul__=ParaMEDMEMDataArrayIntTupleImul DataArrayIntTuple.__idiv__=ParaMEDMEMDataArrayIntTupleIdiv DataArrayIntTuple.__imod__=ParaMEDMEMDataArrayIntTupleImod +DenseMatrix.__iadd__=ParaMEDMEMDenseMatrixIadd +DenseMatrix.__isub__=ParaMEDMEMDenseMatrixIsub + del ParaMEDMEMDataArrayDoubleIadd del ParaMEDMEMDataArrayDoubleIsub del ParaMEDMEMDataArrayDoubleImul @@ -71,4 +74,6 @@ del ParaMEDMEMDataArrayIntTupleIsub del ParaMEDMEMDataArrayIntTupleImul del ParaMEDMEMDataArrayIntTupleIdiv del ParaMEDMEMDataArrayIntTupleImod +del ParaMEDMEMDenseMatrixIadd +del ParaMEDMEMDenseMatrixIsub %} diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapper.i b/src/MEDCoupling_Swig/MEDCouplingRemapper.i index 5154fc43e..10d9f9b38 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapper.i +++ b/src/MEDCoupling_Swig/MEDCouplingRemapper.i @@ -178,6 +178,12 @@ def ParaMEDMEMDataArrayIntTupleIdiv(self,*args): def ParaMEDMEMDataArrayIntTupleImod(self,*args): import _MEDCouplingRemapper return _MEDCouplingRemapper.DataArrayIntTuple____imod___(self, self, *args) +def ParaMEDMEMDenseMatrixIadd(self,*args): + import _MEDCouplingRemapper + return _MEDCouplingRemapper.DenseMatrix____iadd___(self, self, *args) +def ParaMEDMEMDenseMatrixIsub(self,*args): + import _MEDCouplingRemapper + return _MEDCouplingRemapper.DenseMatrix____isub___(self, self, *args) %} %include "MEDCouplingFinalize.i" diff --git a/src/MEDLoader/Swig/MEDLoader.i b/src/MEDLoader/Swig/MEDLoader.i index 24aa39736..7ca7af2ff 100644 --- a/src/MEDLoader/Swig/MEDLoader.i +++ b/src/MEDLoader/Swig/MEDLoader.i @@ -96,6 +96,12 @@ def ParaMEDMEMDataArrayIntTupleIdiv(self,*args): def ParaMEDMEMDataArrayIntTupleImod(self,*args): import _MEDLoader return _MEDLoader.DataArrayIntTuple____imod___(self, self, *args) +def ParaMEDMEMDenseMatrixIadd(self,*args): + import _MEDLoader + return _MEDLoader.DenseMatrix____iadd___(self, self, *args) +def ParaMEDMEMDenseMatrixIsub(self,*args): + import _MEDLoader + return _MEDLoader.DenseMatrix____isub___(self, self, *args) %} %include "MEDCouplingFinalize.i" -- 2.39.2