From 2b1e9d5e6a7c31dd6b94aa53a50aa64641c85b39 Mon Sep 17 00:00:00 2001
From: Anthony Geay <anthony.geay@edf.fr>
Date: Tue, 23 Aug 2022 11:33:45 +0200
Subject: [PATCH] Integration of newton solver coming from numerical recipes to
 solve

---
 .../Bases/InterpKernelException.hxx           |   1 +
 src/INTERP_KERNEL/CMakeLists.txt              |   4 +
 src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx |  53 +++++
 src/INTERP_KERNEL/InterpKernelDenseMatrix.txx | 141 +++++++++++
 .../InterpKernelRootsMultiDim.hxx             | 218 ++++++++++++++++++
 .../LinearAlgebra/InterpKernelDenseMatrix.cxx |  22 ++
 .../LinearAlgebra/InterpKernelLUDecomp.cxx    | 142 ++++++++++++
 .../LinearAlgebra/InterpKernelLUDecomp.hxx    |  43 ++++
 .../LinearAlgebra/InterpKernelQRDecomp.cxx    | 170 ++++++++++++++
 .../LinearAlgebra/InterpKernelQRDecomp.hxx    |  44 ++++
 10 files changed, 838 insertions(+)
 create mode 100644 src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx
 create mode 100644 src/INTERP_KERNEL/InterpKernelDenseMatrix.txx
 create mode 100644 src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx
 create mode 100644 src/INTERP_KERNEL/LinearAlgebra/InterpKernelDenseMatrix.cxx
 create mode 100644 src/INTERP_KERNEL/LinearAlgebra/InterpKernelLUDecomp.cxx
 create mode 100644 src/INTERP_KERNEL/LinearAlgebra/InterpKernelLUDecomp.hxx
 create mode 100644 src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.cxx
 create mode 100644 src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.hxx

diff --git a/src/INTERP_KERNEL/Bases/InterpKernelException.hxx b/src/INTERP_KERNEL/Bases/InterpKernelException.hxx
index 7703b7b46..f4f05bd62 100644
--- a/src/INTERP_KERNEL/Bases/InterpKernelException.hxx
+++ b/src/INTERP_KERNEL/Bases/InterpKernelException.hxx
@@ -25,6 +25,7 @@
 
 #include <string>
 #include <exception>
+#include <sstream>
 
 namespace INTERP_KERNEL
 {
diff --git a/src/INTERP_KERNEL/CMakeLists.txt b/src/INTERP_KERNEL/CMakeLists.txt
index 5de054ed0..bc07afc1c 100644
--- a/src/INTERP_KERNEL/CMakeLists.txt
+++ b/src/INTERP_KERNEL/CMakeLists.txt
@@ -62,6 +62,9 @@ SET(interpkernel_SOURCES
   ExprEval/InterpKernelValue.cxx
   ExprEval/InterpKernelAsmX86.cxx
   GaussPoints/InterpKernelGaussCoords.cxx
+  LinearAlgebra/InterpKernelDenseMatrix.cxx
+  LinearAlgebra/InterpKernelLUDecomp.cxx
+  LinearAlgebra/InterpKernelQRDecomp.cxx
   )
 
 INCLUDE_DIRECTORIES(
@@ -70,6 +73,7 @@ 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)
diff --git a/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx b/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx
new file mode 100644
index 000000000..38e4a6128
--- /dev/null
+++ b/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx
@@ -0,0 +1,53 @@
+// 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>;
+}
diff --git a/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx b/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx
new file mode 100644
index 000000000..9705f4fdb
--- /dev/null
+++ b/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx
@@ -0,0 +1,141 @@
+// 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);
+    }
+  }
+}
diff --git a/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx b/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx
new file mode 100644
index 000000000..3fba59fdb
--- /dev/null
+++ b/src/INTERP_KERNEL/InterpKernelRootsMultiDim.hxx
@@ -0,0 +1,218 @@
+// 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");
+  }
+}
diff --git a/src/INTERP_KERNEL/LinearAlgebra/InterpKernelDenseMatrix.cxx b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelDenseMatrix.cxx
new file mode 100644
index 000000000..141b3b935
--- /dev/null
+++ b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelDenseMatrix.cxx
@@ -0,0 +1,22 @@
+// 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>;
diff --git a/src/INTERP_KERNEL/LinearAlgebra/InterpKernelLUDecomp.cxx b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelLUDecomp.cxx
new file mode 100644
index 000000000..3416e5ad2
--- /dev/null
+++ b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelLUDecomp.cxx
@@ -0,0 +1,142 @@
+// 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];
+}
diff --git a/src/INTERP_KERNEL/LinearAlgebra/InterpKernelLUDecomp.hxx b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelLUDecomp.hxx
new file mode 100644
index 000000000..4c980a06f
--- /dev/null
+++ b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelLUDecomp.hxx
@@ -0,0 +1,43 @@
+// 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);
+  };
+}
diff --git a/src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.cxx b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.cxx
new file mode 100644
index 000000000..78afd9829
--- /dev/null
+++ b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.cxx
@@ -0,0 +1,170 @@
+// 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;
+	}
+}
diff --git a/src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.hxx b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.hxx
new file mode 100644
index 000000000..9f05a0ef8
--- /dev/null
+++ b/src/INTERP_KERNEL/LinearAlgebra/InterpKernelQRDecomp.hxx
@@ -0,0 +1,44 @@
+// 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);
+  };
+}
-- 
2.39.2