src/INTERP_KERNEL/LinearAlgebra/InterpKernelLUDecomp.cxx

   1 // Copyright (C) 2022-2024  CEA, EDF
   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
  20 // Implementation coming from Numerical Recipes in C of 1994 (version 2.04)
  21
  22 #include "InterpKernelLUDecomp.hxx"
  23 #include "InterpKernelException.hxx"
  24
  25 #include <cmath>
  26 #include <sstream>
  27
  28 using namespace INTERP_KERNEL;
  29
  30 LUDecomp::LUDecomp(const INTERP_KERNEL::DenseMatrix& a) : n(a.nrows()), lu(a), aref(a), indx(n)
  31 {
  32         const double TINY=1.0e-40;
  33         mcIdType i,imax,j,k;
  34         double big,temp;
  35         std::vector<double> vv(n);
  36         d=1.0;
  37         for (i=0;i<n;i++)
  38   {
  39                 big=0.0;
  40                 for (j=0;j<n;j++)
  41                         if ((temp=std::abs(lu[i][j])) > big) big=temp;
  42                 if (big == 0.0) THROW_IK_EXCEPTION("Singular matrix in LUDecomp");
  43                 vv[i]=1.0/big;
  44         }
  45         for (k=0;k<n;k++)
  46   {
  47                 big=0.0;
  48                 imax=k;
  49                 for (i=k;i<n;i++)
  50     {
  51                         temp=vv[i]*std::abs(lu[i][k]);
  52                         if (temp > big)
  53       {
  54                                 big=temp;
  55                                 imax=i;
  56                         }
  57                 }
  58                 if (k != imax) {
  59                         for (j=0;j<n;j++) {
  60                                 temp=lu[imax][j];
  61                                 lu[imax][j]=lu[k][j];
  62                                 lu[k][j]=temp;
  63                         }
  64                         d = -d;
  65                         vv[imax]=vv[k];
  66                 }
  67                 indx[k]=imax;
  68                 if (lu[k][k] == 0.0) lu[k][k]=TINY;
  69                 for (i=k+1;i<n;i++) {
  70                         temp=lu[i][k] /= lu[k][k];
  71                         for (j=k+1;j<n;j++)
  72                                 lu[i][j] -= temp*lu[k][j];
  73                 }
  74         }
  75 }
  76 void LUDecomp::solve(const std::vector<double>& b, std::vector<double> &x)
  77 {
  78         mcIdType i,ii=0,ip,j;
  79         double sum;
  80         if (b.size() != ToSizeT(n) || x.size() != ToSizeT(n))
  81                 THROW_IK_EXCEPTION("LUDecomp::solve bad sizes");
  82         for (i=0;i<n;i++) x[i] = b[i];
  83         for (i=0;i<n;i++) {
  84                 ip=indx[i];
  85                 sum=x[ip];
  86                 x[ip]=x[i];
  87                 if (ii != 0)
  88                         for (j=ii-1;j<i;j++) sum -= lu[i][j]*x[j];
  89                 else if (sum != 0.0)
  90                         ii=i+1;
  91                 x[i]=sum;
  92         }
  93         for (i=n-1;i>=0;i--) {
  94                 sum=x[i];
  95                 for (j=i+1;j<n;j++) sum -= lu[i][j]*x[j];
  96                 x[i]=sum/lu[i][i];
  97         }
  98 }
  99
 100 void LUDecomp::solve(const INTERP_KERNEL::DenseMatrix& b, INTERP_KERNEL::DenseMatrix &x)
 101 {
 102         mcIdType i,j,m(b.ncols());
 103         if (b.nrows() != n || x.nrows() != n || b.ncols() != x.ncols())
 104                 THROW_IK_EXCEPTION("LUDecomp::solve bad sizes");
 105         std::vector<double> xx(n);
 106         for (j=0;j<m;j++) {
 107                 for (i=0;i<n;i++) xx[i] = b[i][j];
 108                 solve(xx,xx);
 109                 for (i=0;i<n;i++) x[i][j] = xx[i];
 110         }
 111 }
 112
 113 void LUDecomp::inverse(INTERP_KERNEL::DenseMatrix &ainv)
 114 {
 115         mcIdType i,j;
 116         ainv.resize(n,n);
 117         for (i=0;i<n;i++) {
 118                 for (j=0;j<n;j++) ainv[i][j] = 0.;
 119                 ainv[i][i] = 1.;
 120         }
 121         solve(ainv,ainv);
 122 }
 123
 124 double LUDecomp::det()
 125 {
 126         double dd = d;
 127         for (mcIdType i=0;i<n;i++) dd *= lu[i][i];
 128         return dd;
 129 }
 130
 131 void LUDecomp::mprove(const std::vector<double>& b, std::vector<double> &x)
 132 {
 133         mcIdType i,j;
 134         std::vector<double> r(n);
 135         for (i=0;i<n;i++) {
 136                 long double sdp = -b[i];
 137                 for (j=0;j<n;j++)
 138                         sdp += (long double)aref[i][j] * (long double)x[j];
 139                 r[i]=double(sdp);
 140         }
 141         solve(r,r);
 142         for (i=0;i<n;i++) x[i] -= r[i];
 143 }