#include <string>
#include <exception>
+#include <sstream>
namespace INTERP_KERNEL
{
ExprEval/InterpKernelValue.cxx
ExprEval/InterpKernelAsmX86.cxx
GaussPoints/InterpKernelGaussCoords.cxx
+ LinearAlgebra/InterpKernelDenseMatrix.cxx
+ LinearAlgebra/InterpKernelLUDecomp.cxx
+ LinearAlgebra/InterpKernelQRDecomp.cxx
)
INCLUDE_DIRECTORIES(
${CMAKE_CURRENT_SOURCE_DIR}/Geometric2D
${CMAKE_CURRENT_SOURCE_DIR}/ExprEval
${CMAKE_CURRENT_SOURCE_DIR}/GaussPoints
+ ${CMAKE_CURRENT_SOURCE_DIR}/LinearAlgebra
)
IF (NOT DEFINED MSVC)
--- /dev/null
+// Copyright (C) 2022 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
+//
+
+#pragma once
+
+#include "INTERPKERNELDefines.hxx"
+#include "MCIdType.hxx"
+
+namespace INTERP_KERNEL
+{
+ template <class T>
+ class INTERPKERNEL_EXPORT DenseMatrixT
+ {
+ private:
+ mcIdType nn;
+ mcIdType mm;
+ T **v;
+ public:
+ DenseMatrixT();
+ DenseMatrixT(mcIdType n, mcIdType m);
+ DenseMatrixT(mcIdType n, mcIdType m, const T &a);
+ DenseMatrixT(mcIdType n, mcIdType m, const T *a);
+ DenseMatrixT(const DenseMatrixT &rhs);
+ DenseMatrixT & operator=(const DenseMatrixT &rhs);
+ using value_type = T;
+ //! subscripting: pointer to row i
+ T* operator[](const mcIdType i) { return v[i]; }
+ const T* operator[](const mcIdType i) const { return v[i]; }
+ mcIdType nrows() const { return nn; }
+ mcIdType ncols() const { return mm; }
+ void resize(mcIdType newn, mcIdType newm);
+ void assign(mcIdType newn, mcIdType newm, const T &a);
+ ~DenseMatrixT();
+ };
+
+ using DenseMatrix = DenseMatrixT<double>;
+}
--- /dev/null
+// Copyright (C) 2022 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
+//
+
+#pragma once
+
+#include "InterpKernelDenseMatrix.hxx"
+
+namespace INTERP_KERNEL
+{
+ template <class T>
+ DenseMatrixT<T>::DenseMatrixT() : nn(0), mm(0), v(nullptr) {}
+
+ template <class T>
+ DenseMatrixT<T>::DenseMatrixT(mcIdType n, mcIdType m) : nn(n), mm(m), v(n>0 ? new T*[n] : nullptr)
+ {
+ mcIdType i,nel=m*n;
+ if (v) v[0] = nel>0 ? new T[nel] : nullptr;
+ for (i=1;i<n;i++) v[i] = v[i-1] + m;
+ }
+
+ template <class T>
+ DenseMatrixT<T>::DenseMatrixT(mcIdType n, mcIdType m, const T &a) : nn(n), mm(m), v(n>0 ? new T*[n] : nullptr)
+ {
+ mcIdType i,j,nel=m*n;
+ if (v) v[0] = nel>0 ? new T[nel] : nullptr;
+ for (i=1; i< n; i++) v[i] = v[i-1] + m;
+ for (i=0; i< n; i++) for (j=0; j<m; j++) v[i][j] = a;
+ }
+
+ template <class T>
+ DenseMatrixT<T>::DenseMatrixT(mcIdType n, mcIdType m, const T *a) : nn(n), mm(m), v(n>0 ? new T*[n] : nullptr)
+ {
+ mcIdType i,j,nel=m*n;
+ if (v) v[0] = nel>0 ? new T[nel] : nullptr;
+ for (i=1; i< n; i++) v[i] = v[i-1] + m;
+ for (i=0; i< n; i++) for (j=0; j<m; j++) v[i][j] = *a++;
+ }
+
+ template <class T>
+ DenseMatrixT<T>::DenseMatrixT(const DenseMatrixT &rhs) : nn(rhs.nn), mm(rhs.mm), v(nn>0 ? new T*[nn] : nullptr)
+ {
+ mcIdType i,j,nel=mm*nn;
+ if (v) v[0] = nel>0 ? new T[nel] : nullptr;
+ for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
+ for (i=0; i< nn; i++) for (j=0; j<mm; j++) v[i][j] = rhs[i][j];
+ }
+
+ /*!
+ * postcondition: normal assignment via copying has been performed
+ * if matrix and rhs were different sizes, matrix
+ * has been resized to match the size of rhs
+ */
+ template <class T>
+ DenseMatrixT<T> & DenseMatrixT<T>::operator=(const DenseMatrixT<T> &rhs)
+ {
+ if (this != &rhs) {
+ mcIdType i,j,nel;
+ if (nn != rhs.nn || mm != rhs.mm) {
+ if ( v ) {
+ delete[] (v[0]);
+ delete[] (v);
+ }
+ nn=rhs.nn;
+ mm=rhs.mm;
+ v = nn>0 ? new T*[nn] : nullptr;
+ nel = mm*nn;
+ if (v) v[0] = nel>0 ? new T[nel] : nullptr;
+ for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
+ }
+ for (i=0; i< nn; i++) for (j=0; j<mm; j++) v[i][j] = rhs[i][j];
+ }
+ return *this;
+ }
+
+ template <class T>
+ void DenseMatrixT<T>::resize(mcIdType newn, mcIdType newm)
+ {
+ mcIdType i,nel;
+ if (newn != nn || newm != mm)
+ {
+ if ( v )
+ {
+ delete[] (v[0]);
+ delete[] (v);
+ }
+ nn = newn;
+ mm = newm;
+ v = nn>0 ? new T*[nn] : nullptr;
+ nel = mm*nn;
+ if (v) v[0] = nel>0 ? new T[nel] : nullptr;
+ for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
+ }
+ }
+
+ template <class T>
+ void DenseMatrixT<T>::assign(mcIdType newn, mcIdType newm, const T& a)
+ {
+ mcIdType i,j,nel;
+ if (newn != nn || newm != mm)
+ {
+ if ( v )
+ {
+ delete[] (v[0]);
+ delete[] (v);
+ }
+ nn = newn;
+ mm = newm;
+ v = nn>0 ? new T*[nn] : nullptr;
+ nel = mm*nn;
+ if (v) v[0] = nel>0 ? new T[nel] : nullptr;
+ for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
+ }
+ for (i=0; i< nn; i++) for (j=0; j<mm; j++) v[i][j] = a;
+ }
+
+ template <class T>
+ DenseMatrixT<T>::~DenseMatrixT()
+ {
+ if ( v )
+ {
+ delete[] (v[0]);
+ delete[] (v);
+ }
+ }
+}
--- /dev/null
+// Copyright (C) 2022 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
+//
+
+// Implementation coming from Numerical Recipes in C of 1994 (version 2.04) p 384
+// Root finding and non linear sets of equation using line searches and backtracking
+
+#include "InterpKernelDenseMatrix.hxx"
+#include "InterpKernelLUDecomp.hxx"
+#include "InterpKernelQRDecomp.hxx"
+#include "InterpKernelException.hxx"
+#include "MCIdType.hxx"
+
+#include <vector>
+#include <limits>
+#include <cmath>
+
+template<class T>
+inline T sqr(const T a) { return a*a; }
+
+namespace INTERP_KERNEL
+{
+ template <class T>
+ void LineSearches(const std::vector<double>& xold, const double fold, const std::vector<double>& g, std::vector<double> &p,
+ std::vector<double> &x, double &f, const double stpmax, bool &check, T &func)
+ {
+ const double ALF=1.0e-4, TOLX=std::numeric_limits<double>::epsilon();
+ double a,alam,alam2=0.0,alamin,b,disc,f2=0.0;
+ double rhs1,rhs2,slope=0.0,sum=0.0,temp,test,tmplam;
+ mcIdType i,n=xold.size();
+ check=false;
+ for (i=0;i<n;i++) sum += p[i]*p[i];
+ sum=std::sqrt(sum);
+ if (sum > stpmax)
+ for (i=0;i<n;i++)
+ p[i] *= stpmax/sum;
+ for (i=0;i<n;i++)
+ slope += g[i]*p[i];
+ if (slope >= 0.0) THROW_IK_EXCEPTION("Roundoff problem in LineSearches.");
+ test=0.0;
+ for (i=0;i<n;i++)
+ {
+ temp=std::abs(p[i])/std::max(std::abs(xold[i]),1.0);
+ if (temp > test) test=temp;
+ }
+ alamin=TOLX/test;
+ alam=1.0;
+ for (;;)
+ {
+ for (i=0;i<n;i++) x[i]=xold[i]+alam*p[i];
+ f=func(x);
+ if (alam < alamin)
+ {
+ for (i=0;i<n;i++) x[i]=xold[i];
+ check=true;
+ return;
+ } else if (f <= fold+ALF*alam*slope) return;
+ else
+ {
+ if (alam == 1.0)
+ tmplam = -slope/(2.0*(f-fold-slope));
+ else
+ {
+ rhs1=f-fold-alam*slope;
+ rhs2=f2-fold-alam2*slope;
+ a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
+ b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
+ if (a == 0.0) tmplam = -slope/(2.0*b);
+ else {
+ disc=b*b-3.0*a*slope;
+ if (disc < 0.0) tmplam=0.5*alam;
+ else if (b <= 0.0) tmplam=(-b+std::sqrt(disc))/(3.0*a);
+ else tmplam=-slope/(b+std::sqrt(disc));
+ }
+ if (tmplam>0.5*alam)
+ tmplam=0.5*alam;
+ }
+ }
+ alam2=alam;
+ f2 = f;
+ alam=std::max(tmplam,0.1*alam);
+ }
+ }
+ template <class T>
+ class JacobianCalculator
+ {
+ private:
+ const double EPS;
+ T &func;
+ public:
+ JacobianCalculator(T &funcc) : EPS(1.0e-8),func(funcc) {}
+ INTERP_KERNEL::DenseMatrix operator() (const std::vector<double>& x, const std::vector<double>& fvec)
+ {
+ mcIdType n=x.size();
+ INTERP_KERNEL::DenseMatrix df(n,n);
+ std::vector<double> xh=x;
+ for (mcIdType j=0;j<n;j++)
+ {
+ double temp=xh[j];
+ double h=EPS*std::abs(temp);
+ if (h == 0.0) h=EPS;
+ xh[j]=temp+h;
+ h=xh[j]-temp;
+ std::vector<double> f=func(xh);
+ xh[j]=temp;
+ for (mcIdType i=0;i<n;i++)
+ df[i][j]=(f[i]-fvec[i])/h;
+ }
+ return df;
+ }
+ };
+
+ template <class T>
+ class FMin
+ {
+ private:
+ std::vector<double> fvec;
+ T &func;
+ mcIdType n;
+ public:
+ FMin(T &funcc) : func(funcc){}
+ double operator() (const std::vector<double>& x)
+ {
+ n=x.size();
+ double sum=0;
+ fvec=func(x);
+ for (mcIdType i=0;i<n;i++) sum += sqr(fvec[i]);
+ return 0.5*sum;
+ }
+ std::vector<double>& getVector() { return fvec; }
+ };
+
+ template <class T>
+ void SolveWithNewton(std::vector<double> &x, bool &check, T &vecfunc)
+ {
+ const mcIdType MAXITS=200;
+ const double TOLF=1.0e-8,TOLMIN=1.0e-12,STPMX=100.0;
+ const double TOLX=std::numeric_limits<double>::epsilon();
+ mcIdType i,j,its,n=x.size();
+ double den,f,fold,stpmax,sum,temp,test;
+ std::vector<double> g(n),p(n),xold(n);
+ INTERP_KERNEL::DenseMatrix fjac(n,n);
+ FMin<T> fmin(vecfunc);
+ JacobianCalculator<T> fdjac(vecfunc);
+ std::vector<double> &fvec=fmin.getVector();
+ f=fmin(x);
+ test=0.0;
+ for (i=0;i<n;i++)
+ if (std::abs(fvec[i]) > test) test=std::abs(fvec[i]);
+ if (test < 0.01*TOLF)
+ {
+ check=false;
+ return;
+ }
+ sum=0.0;
+ for (i=0;i<n;i++) sum += sqr(x[i]);
+ stpmax=STPMX*std::max(std::sqrt(sum),double(n));
+ for (its=0;its<MAXITS;its++)
+ {
+ fjac=fdjac(x,fvec);
+ for (i=0;i<n;i++)
+ {
+ sum=0.0;
+ for (j=0;j<n;j++) sum += fjac[j][i]*fvec[j];
+ g[i]=sum;
+ }
+ for (i=0;i<n;i++) xold[i]=x[i];
+ fold=f;
+ for (i=0;i<n;i++) p[i] = -fvec[i];
+ INTERP_KERNEL::LUDecomp alu(fjac);
+ alu.solve(p,p);
+ LineSearches(xold,fold,g,p,x,f,stpmax,check,fmin);
+ test=0.0;
+ for (i=0;i<n;i++)
+ if (std::abs(fvec[i]) > test) test=std::abs(fvec[i]);
+ if (test < TOLF)
+ {
+ check=false;
+ return;
+ }
+ if (check)
+ {
+ test=0.0;
+ den=std::max(f,0.5*n);
+ for (i=0;i<n;i++) {
+ temp=std::abs(g[i])*std::max(std::abs(x[i]),1.0)/den;
+ if (temp > test) test=temp;
+ }
+ check=(test < TOLMIN);
+ return;
+ }
+ test=0.0;
+ for (i=0;i<n;i++)
+ {
+ temp=(std::abs(x[i]-xold[i]))/std::max(std::abs(x[i]),1.0);
+ if (temp > test) test=temp;
+ }
+ if (test < TOLX)
+ return;
+ }
+ THROW_IK_EXCEPTION("MAXITS exceeded in newt");
+ }
+}
--- /dev/null
+// Copyright (C) 2022 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
+//
+
+#include "InterpKernelDenseMatrix.txx"
+
+template class INTERPKERNEL_EXPORT INTERP_KERNEL::DenseMatrixT<double>;
--- /dev/null
+// Copyright (C) 2022 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
+//
+
+// Implementation coming from Numerical Recipes in C of 1994 (version 2.04)
+
+#include "InterpKernelLUDecomp.hxx"
+#include "InterpKernelException.hxx"
+
+#include <sstream>
+
+using namespace INTERP_KERNEL;
+
+LUDecomp::LUDecomp(const INTERP_KERNEL::DenseMatrix& a) : n(a.nrows()), lu(a), aref(a), indx(n)
+{
+ const double TINY=1.0e-40;
+ mcIdType i,imax,j,k;
+ double big,temp;
+ std::vector<double> vv(n);
+ d=1.0;
+ for (i=0;i<n;i++)
+ {
+ big=0.0;
+ for (j=0;j<n;j++)
+ if ((temp=std::abs(lu[i][j])) > big) big=temp;
+ if (big == 0.0) THROW_IK_EXCEPTION("Singular matrix in LUDecomp");
+ vv[i]=1.0/big;
+ }
+ for (k=0;k<n;k++)
+ {
+ big=0.0;
+ imax=k;
+ for (i=k;i<n;i++)
+ {
+ temp=vv[i]*std::abs(lu[i][k]);
+ if (temp > big)
+ {
+ big=temp;
+ imax=i;
+ }
+ }
+ if (k != imax) {
+ for (j=0;j<n;j++) {
+ temp=lu[imax][j];
+ lu[imax][j]=lu[k][j];
+ lu[k][j]=temp;
+ }
+ d = -d;
+ vv[imax]=vv[k];
+ }
+ indx[k]=imax;
+ if (lu[k][k] == 0.0) lu[k][k]=TINY;
+ for (i=k+1;i<n;i++) {
+ temp=lu[i][k] /= lu[k][k];
+ for (j=k+1;j<n;j++)
+ lu[i][j] -= temp*lu[k][j];
+ }
+ }
+}
+void LUDecomp::solve(const std::vector<double>& b, std::vector<double> &x)
+{
+ mcIdType i,ii=0,ip,j;
+ double sum;
+ if (b.size() != ToSizeT(n) || x.size() != ToSizeT(n))
+ THROW_IK_EXCEPTION("LUDecomp::solve bad sizes");
+ for (i=0;i<n;i++) x[i] = b[i];
+ for (i=0;i<n;i++) {
+ ip=indx[i];
+ sum=x[ip];
+ x[ip]=x[i];
+ if (ii != 0)
+ for (j=ii-1;j<i;j++) sum -= lu[i][j]*x[j];
+ else if (sum != 0.0)
+ ii=i+1;
+ x[i]=sum;
+ }
+ for (i=n-1;i>=0;i--) {
+ sum=x[i];
+ for (j=i+1;j<n;j++) sum -= lu[i][j]*x[j];
+ x[i]=sum/lu[i][i];
+ }
+}
+
+void LUDecomp::solve(const INTERP_KERNEL::DenseMatrix& b, INTERP_KERNEL::DenseMatrix &x)
+{
+ mcIdType i,j,m(b.ncols());
+ if (b.nrows() != n || x.nrows() != n || b.ncols() != x.ncols())
+ THROW_IK_EXCEPTION("LUDecomp::solve bad sizes");
+ std::vector<double> xx(n);
+ for (j=0;j<m;j++) {
+ for (i=0;i<n;i++) xx[i] = b[i][j];
+ solve(xx,xx);
+ for (i=0;i<n;i++) x[i][j] = xx[i];
+ }
+}
+
+void LUDecomp::inverse(INTERP_KERNEL::DenseMatrix &ainv)
+{
+ mcIdType i,j;
+ ainv.resize(n,n);
+ for (i=0;i<n;i++) {
+ for (j=0;j<n;j++) ainv[i][j] = 0.;
+ ainv[i][i] = 1.;
+ }
+ solve(ainv,ainv);
+}
+
+double LUDecomp::det()
+{
+ double dd = d;
+ for (mcIdType i=0;i<n;i++) dd *= lu[i][i];
+ return dd;
+}
+
+void LUDecomp::mprove(const std::vector<double>& b, std::vector<double> &x)
+{
+ mcIdType i,j;
+ std::vector<double> r(n);
+ for (i=0;i<n;i++) {
+ long double sdp = -b[i];
+ for (j=0;j<n;j++)
+ sdp += (long double)aref[i][j] * (long double)x[j];
+ r[i]=double(sdp);
+ }
+ solve(r,r);
+ for (i=0;i<n;i++) x[i] -= r[i];
+}
--- /dev/null
+// Copyright (C) 2022 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
+//
+
+#pragma once
+
+#include "InterpKernelDenseMatrix.hxx"
+#include <vector>
+
+namespace INTERP_KERNEL
+{
+ class LUDecomp
+ {
+ private:
+ mcIdType n;
+ INTERP_KERNEL::DenseMatrix lu;
+ std::vector<mcIdType> indx;
+ double d;
+ const INTERP_KERNEL::DenseMatrix& aref;
+ public:
+ LUDecomp (const INTERP_KERNEL::DenseMatrix& a);
+ void solve(const std::vector<double>& b, std::vector<double> &x);
+ void solve(const INTERP_KERNEL::DenseMatrix& b, INTERP_KERNEL::DenseMatrix &x);
+ void inverse(INTERP_KERNEL::DenseMatrix &ainv);
+ double det();
+ void mprove(const std::vector<double>& b, std::vector<double> &x);
+ };
+}
--- /dev/null
+// Copyright (C) 2022 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
+//
+
+// Implementation coming from Numerical Recipes in C of 1994 (version 2.04)
+
+#include "InterpKernelQRDecomp.hxx"
+#include "InterpKernelException.hxx"
+
+#include <cmath>
+
+using namespace INTERP_KERNEL;
+
+template<class T>
+inline T sqr(const T a) { return a*a; }
+
+template<class T>
+inline T sign(const T &a, const T &b) { return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a); }
+
+QRDecomp::QRDecomp(const INTERP_KERNEL::DenseMatrix& a): n(a.nrows()), qt(n,n), r(a), sing(false)
+{
+ mcIdType i,j,k;
+ std::vector<double> c(n), d(n);
+ double scale,sigma,sum,tau;
+ for (k=0;k<n-1;k++) {
+ scale=0.0;
+ for (i=k;i<n;i++) scale=std::max(scale,std::abs(r[i][k]));
+ if (scale == 0.0) {
+ sing=true;
+ c[k]=d[k]=0.0;
+ } else {
+ for (i=k;i<n;i++) r[i][k] /= scale;
+ for (sum=0.0,i=k;i<n;i++) sum += sqr(r[i][k]);
+ sigma=sign(std::sqrt(sum),r[k][k]);
+ r[k][k] += sigma;
+ c[k]=sigma*r[k][k];
+ d[k] = -scale*sigma;
+ for (j=k+1;j<n;j++) {
+ for (sum=0.0,i=k;i<n;i++) sum += r[i][k]*r[i][j];
+ tau=sum/c[k];
+ for (i=k;i<n;i++) r[i][j] -= tau*r[i][k];
+ }
+ }
+ }
+ d[n-1]=r[n-1][n-1];
+ if (d[n-1] == 0.0) sing=true;
+ for (i=0;i<n;i++) {
+ for (j=0;j<n;j++) qt[i][j]=0.0;
+ qt[i][i]=1.0;
+ }
+ for (k=0;k<n-1;k++) {
+ if (c[k] != 0.0) {
+ for (j=0;j<n;j++) {
+ sum=0.0;
+ for (i=k;i<n;i++)
+ sum += r[i][k]*qt[i][j];
+ sum /= c[k];
+ for (i=k;i<n;i++)
+ qt[i][j] -= sum*r[i][k];
+ }
+ }
+ }
+ for (i=0;i<n;i++) {
+ r[i][i]=d[i];
+ for (j=0;j<i;j++) r[i][j]=0.0;
+ }
+}
+
+void QRDecomp::solve(const std::vector<double>& b, std::vector<double> &x)
+{
+ qtmult(b,x);
+ rsolve(x,x);
+}
+
+void QRDecomp::qtmult(const std::vector<double>& b, std::vector<double> &x)
+{
+ mcIdType i,j;
+ double sum;
+ for (i=0;i<n;i++)
+ {
+ sum = 0.;
+ for (j=0;j<n;j++) sum += qt[i][j]*b[j];
+ x[i] = sum;
+ }
+}
+
+void QRDecomp::rsolve(const std::vector<double>& b, std::vector<double> &x)
+{
+ mcIdType i,j;
+ double sum;
+ if (sing) THROW_IK_EXCEPTION("attempting solve in a singular QR");
+ for (i=n-1;i>=0;i--) {
+ sum=b[i];
+ for (j=i+1;j<n;j++) sum -= r[i][j]*x[j];
+ x[i]=sum/r[i][i];
+ }
+}
+void QRDecomp::update(const std::vector<double>& u, const std::vector<double>& v)
+{
+ mcIdType i,k;
+ std::vector<double> w(u);
+ for (k=n-1;k>=0;k--)
+ if (w[k] != 0.0) break;
+ if (k < 0) k=0;
+ for (i=k-1;i>=0;i--) {
+ rotate(i,w[i],-w[i+1]);
+ if (w[i] == 0.0)
+ w[i]=std::abs(w[i+1]);
+ else if (std::abs(w[i]) > std::abs(w[i+1]))
+ w[i]=std::abs(w[i])*std::sqrt(1.0+sqr(w[i+1]/w[i]));
+ else w[i]=std::abs(w[i+1])*std::sqrt(1.0+sqr(w[i]/w[i+1]));
+ }
+ for (i=0;i<n;i++) r[0][i] += w[0]*v[i];
+ for (i=0;i<k;i++)
+ rotate(i,r[i][i],-r[i+1][i]);
+ for (i=0;i<n;i++)
+ if (r[i][i] == 0.0) sing=true;
+}
+
+void QRDecomp::rotate(const mcIdType i, const double a, const double b)
+{
+ mcIdType j;
+ double c,fact,s,w,y;
+ if (a == 0.0)
+ {
+ c=0.0;
+ s=(b >= 0.0 ? 1.0 : -1.0);
+ }
+ else if (std::abs(a) > std::abs(b))
+ {
+ fact=b/a;
+ c=sign(1.0/std::sqrt(1.0+(fact*fact)),a);
+ s=fact*c;
+ }
+ else
+ {
+ fact=a/b;
+ s=sign(1.0/std::sqrt(1.0+(fact*fact)),b);
+ c=fact*s;
+ }
+ for (j=i;j<n;j++)
+ {
+ y=r[i][j];
+ w=r[i+1][j];
+ r[i][j]=c*y-s*w;
+ r[i+1][j]=s*y+c*w;
+ }
+ for (j=0;j<n;j++)
+ {
+ y=qt[i][j];
+ w=qt[i+1][j];
+ qt[i][j]=c*y-s*w;
+ qt[i+1][j]=s*y+c*w;
+ }
+}
--- /dev/null
+// Copyright (C) 2022 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
+//
+
+// Implementation coming from Numerical Recipes in C of 1994 (version 2.04)
+
+#pragma once
+
+#include "InterpKernelDenseMatrix.hxx"
+#include <vector>
+
+namespace INTERP_KERNEL
+{
+ struct QRDecomp
+ {
+ private:
+ mcIdType n;
+ INTERP_KERNEL::DenseMatrix qt;
+ INTERP_KERNEL::DenseMatrix r;
+ bool sing;
+ public:
+ QRDecomp(const INTERP_KERNEL::DenseMatrix& a);
+ void solve(const std::vector<double>& b, std::vector<double> &x);
+ void qtmult(const std::vector<double>& b, std::vector<double> &x);
+ void rsolve(const std::vector<double>& b, std::vector<double> &x);
+ void update(const std::vector<double>& u, const std::vector<double>& v);
+ void rotate(const mcIdType i, const double a, const double b);
+ };
+}
