vtkdisplay/src/linear_algebra.cc

// nUtil - An utility Library for gnurbs
// Copyright (C) 2008-2019 Eric Bechet
//
// See the LICENSE file for contributions and license information.
// Please report all bugs and problems to <bechet@cadxfem.org>.
//

#include <iostream>
#include <iomanip>
#include "linear_algebra.hxx"


void Vector::Display(std::ostream& out) const
{
//  out << std::setiosflags(std::ios::fixed);
//  out << std::setprecision(8);
  out << std::endl;
  for (int j=0;j<Size();++j)
    out << coord[j] << " ";
  out << std::endl;
}

void Matrix::Display(std::ostream& out) const
{
//  out << std::setiosflags(std::ios::fixed);
//  out << std::setprecision(8);
  for (int i=0;i<this->SizeL();++i)
  {
    out << std::endl;
    for (int j=0;j<this->SizeC();++j)
      out << (*this)(i,j) << " ";
  }
  out << std::endl;
}

#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
 a[k][l]=h+s*(g-h*tau);

int Square_Matrix::Eigen_Vectors(Matrix &M,Vector &d)
{
  int j,iq,ip,i;
  double tresh,theta,tau,t,sm,s,h,g,c;

  std::vector<double> b(Size());
  std::vector<double> z(Size());
  std::vector<std::vector<double> > a(Size());
  for (i=0;i<Size();i++)
  {
    a[i].resize(Size());
  }

  int nrot=0;

  for (ip=0;ip<Size();ip++)
  {
    for (iq=0;iq<Size();iq++)
    {
      mat[ip][iq]=0.0;
      a[ip][iq]=M(ip,iq);
    }
    mat[ip][ip]=1.0;
  }
  for (ip=0;ip<Size();ip++)
  {
    b[ip]=d[ip]=a[ip][ip];
    z[ip]=0.0;
  }
  nrot=0;
  for (i=0;i<100;i++)
  {
    sm=0.0;
    for (ip=0;ip<Size()-1;ip++)
    {
      for (iq=ip+1;iq<Size();iq++)
        sm += fabs(a[ip][iq]);
    }
    if (sm == 0.0)
    {
      break;
    }
    if (i < 3)
      tresh=0.2*sm/ (Size() *Size());
    else
      tresh=0.0;
    for (ip=0;ip<Size()-1;ip++)
    {
      for (iq=ip+1;iq<Size();iq++)
      {
        g=100.0*fabs(a[ip][iq]);
        if ((i > 3) && ((fabs(d[ip]) +g) == fabs(d[ip])) && ((fabs(d[iq]) +g) == fabs(d[iq])))
          a[ip][iq]=0.0;
        else
        {
          if (fabs(a[ip][iq]) > tresh)
          {
            h=d[iq]-d[ip];
            if ((fabs(h) +g) == fabs(h))
              t= (a[ip][iq]) /h;
            else
            {
              theta=0.5*h/ (a[ip][iq]);
              t=1.0/ (fabs(theta) +sqrt(1.0+theta*theta));
              if (theta < 0.0) t = -t;
            }
            c=1.0/sqrt(1+t*t);
            s=t*c;
            tau=s/ (1.0+c);
            h=t*a[ip][iq];
            z[ip] -= h;
            z[iq] += h;
            d[ip] -= h;
            d[iq] += h;
            a[ip][iq]=0.0;
            for (j=0;j<=ip-1;j++)
            {
              ROTATE(a,j,ip,j,iq)
            }
            for (j=ip+1;j<=iq-1;j++)
            {
              ROTATE(a,ip,j,j,iq)
            }
            for (j=iq+1;j<Size();j++)
            {
              ROTATE(a,ip,j,iq,j)
            }
            for (j=0;j<Size();j++)
            {
              ROTATE(mat,j,ip,j,iq)
            }
            ++nrot;
          }
        }
      }
    }
    for (ip=0;ip<Size();ip++)
    {
      b[ip] += z[ip];
      d[ip]=b[ip];
      z[ip]=0.0;
    }
  }
  if (i==100)
  {
    return 1; // there was a problem.
  }
  return 0;
}

#undef ROTATE


void Square_Matrix::Complete_Base(int num)
{
  int i,j,k;
  std::vector<Vector> Tab(Size());
  for (i=0;i<Size();i++)
    Tab[i].Resize(Size());

  double pscal,norm;
  for (i=0;i<num;i++)
  {
    for (j=0;j<Size();j++)
      Tab[i][j]=mat[i][j];
    Tab[i].Normalize();
  }

  for (i=num;i<Size();i++)
  {
    for (j=0;j<num;j++)
    {
      for (k=0;k<Size();k++)
        Tab[i][k]=0;
      Tab[i][(i+j) %Size()]=1.0;
      for (k=0;k<i;k++)
      {
        pscal=Tab[i]*Tab[k];
        Tab[i].Add(Tab[k],-pscal);
      }
      norm=Tab[i].Norm();
      if (norm>1e-4)
        break;

    }
    Tab[i].Normalize();
  }

  for (i=0;i<Size();i++)
    for (j=0;j<Size();j++)
      mat[i][j]=Tab[i][j];
}


void Cholesky_Matrix::CholeskyDec_intern(void)
{
  int k,i,p;
  for (k=0;k<Size();++k)
  {
    double sumsqakp=0.;
    for (p=0;p<k;++p)
    {
      double akp= (*this)(k,p);
      sumsqakp+=akp*akp;
    }
    double& akk= (*this)(k,k);
    if (akk<sumsqakp)
    {
      std::cerr << "Cholesky : (akk<sumsqakp) " << akk << " < " << sumsqakp << std::endl;
      throw "Cholesky : (akk<sumsqakp)";
    }
    akk=sqrt(akk-sumsqakp);
    //    Display();
    //    cout << "a" << k << k << " " <<  akk << " sumsqakp " << sumsqakp << endl;
    for (i=k+1;i<Size();++i)
    {
      double sumaipakp=0.0;
      for (p=0;p<k;++p)
      {
        sumaipakp+= (*this)(i,p) * (*this)(k,p);
      }
      double& aik= (*this)(i,k);
      aik= (aik-sumaipakp) /akk;
      //      Display();
      //      cout << "a" << i << k << " " <<  aik <<" sumaipakp " << sumaipakp << endl;
      //      cout << aik << " " ;
    }
    //    cout << endl;
  }
}


void LU_Matrix::LUDec_intern(void)
{

  int i,imax,j,k;
  double big,dum,sum,temp;
  std::vector<double> vv(Size());


  d=1.;

  for (i=0;i<Size();++i)
  {
    big=0.0;
    for (j=0;j<Size();++j)
      if ((temp=fabs(mat[i][j])) > big) big=temp;

    if (big == 0.0)
    {
      d=0.0;
      return;
    }

    vv[i]=1.0/big;
  }


  for (j=0;j<Size();++j)
  {
    for (i=0;i<j;++i)
    {
      sum=mat[i][j];
      for (k=0;k<i;++k) sum -= mat[i][k]*mat[k][j];
      mat[i][j]=sum;
    }

    big=0.0;
    for (i=j;i<Size();++i)
    {
      sum=mat[i][j];
      for (k=0;k<j;++k)
        sum -= mat[i][k]*mat[k][j];
      mat[i][j]=sum;
      if ((dum=vv[i]*fabs(sum)) >= big)
      {
        big=dum;
        imax=i;
      }
    }
    if (j != imax)
    {
      for (k=0;k<Size();++k)
      {
        dum=mat[imax][k];
        mat[imax][k]=mat[j][k];
        mat[j][k]=dum;
      }
      d = - (d);
      vv[imax]=vv[j];
    }
    tab[j]=imax;

    if (mat[j][j] == 0.0)
    {
      mat[j][j] = 1e-50;
    }

    if (j != Size()-1)
    {
      dum=1.0/ (mat[j][j]);
      for (i=j+1;i<Size();++i) mat[i][j] *= dum;
    }
  }
  for (i=0;i<Size();++i)
    d*=mat[i][i];
}

void LU_Matrix::Solve_Linear_System(Vector &rhsv,Vector &sol)
{
  int i,im1,ii=-1,ip,j;
  double sum;

  if (d==0.0) return;

  sol=rhsv;

  for (i=0;i<Size();++i)
  {
    ip=tab[i];
    sum=sol[ip];
    sol[ip]=sol[i];
    im1=i-1;
    if (ii>=0)
      for (j=ii;j<=im1;++j) sum -= mat[i][j]*sol[j];
    else if (sum!=0.0) ii=i;
    sol[i]=sum;
  }

  for (i=Size()-1;i>=0;i--)
  {
    sum=sol[i];
    for (j=i+1;j<Size();++j) sum -= mat[i][j]*sol[j];
    sol[i]=sum/mat[i][i];
  }
}


double det(double M[2][2])
{
  return M[0][0]*M[1][1]-M[1][0]*M[0][1];
}

double trace(double M[2][2])
{
  return M[0][0]+M[1][1];
}

double prodcc(double M1[2][2],double M2[2][2])
{
  double p=0.;
  for (int i=0;i<2;++i)
    for (int j=0;j<2;++j)
      p+=M1[i][j]*M2[i][j];
  return p;
}


void inverse(double M[2][2],double I[2][2])  // I : inverse de la matrice M
{
  double d = det(M);
  I[0][0] = M[1][1]/d;
  I[1][1] = M[0][0]/d;
  I[0][1] = -M[0][1]/d;
  I[1][0] = -M[1][0]/d;
}


void multi(double A[2][2],double B[2][2],double C[2][2]) // Effectue la multiplication de A et B et rentre le résultat dans la matrice C
{
  for (int i=0;i<2;i++)
  {
    for (int j=0;j<2;j++)
    {
      C[i][j] = 0;
      for (int k=0;k<2;k++)
      {
        C[i][j]+=A[i][k]*B[k][j];
      }
    }
  }
}

double multi(double V1[2],double M[2][2],double V2[2]) // Effectue la multiplication Ai Mij Cj 
{
  double res=0.;
  for (int i=0;i<2;i++)
  {
    for (int j=0;j<2;j++)
    {
      res+=V1[i]*M[i][j]*V2[j];
    }
  }
  return res;
}

Revision:	1156
Committed:	Thu Jun 13 22:02:48 2024 UTC (14 months ago) by francois
File size:	7515 byte(s)
Log Message:	compatibilité Ubuntu 22.04 Suppression des refeences à Windows Ajout d'une banière
#	User	Rev	Content
1	francois	1061	// nUtil - An utility Library for gnurbs
2			// Copyright (C) 2008-2019 Eric Bechet
3			//
4			// See the LICENSE file for contributions and license information.
5			// Please report all bugs and problems to <bechet@cadxfem.org>.
6			//
7
8			#include <iostream>
9			#include <iomanip>
10	francois	1156	#include "linear_algebra.hxx"
11	francois	1061
12
13			void Vector::Display(std::ostream& out) const
14			{
15			// out << std::setiosflags(std::ios::fixed);
16			// out << std::setprecision(8);
17			out << std::endl;
18			for (int j=0;j<Size();++j)
19			out << coord[j] << " ";
20			out << std::endl;
21			}
22
23			void Matrix::Display(std::ostream& out) const
24			{
25			// out << std::setiosflags(std::ios::fixed);
26			// out << std::setprecision(8);
27			for (int i=0;i<this->SizeL();++i)
28			{
29			out << std::endl;
30			for (int j=0;j<this->SizeC();++j)
31			out << (*this)(i,j) << " ";
32			}
33			out << std::endl;
34			}
35
36			#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s(h+gtau);\
37			a[k][l]=h+s(g-htau);
38
39			int Square_Matrix::Eigen_Vectors(Matrix &M,Vector &d)
40			{
41			int j,iq,ip,i;
42			double tresh,theta,tau,t,sm,s,h,g,c;
43
44			std::vector<double> b(Size());
45			std::vector<double> z(Size());
46			std::vector<std::vector<double> > a(Size());
47			for (i=0;i<Size();i++)
48			{
49			a[i].resize(Size());
50			}
51
52			int nrot=0;
53
54			for (ip=0;ip<Size();ip++)
55			{
56			for (iq=0;iq<Size();iq++)
57			{
58			mat[ip][iq]=0.0;
59			a[ip][iq]=M(ip,iq);
60			}
61			mat[ip][ip]=1.0;
62			}
63			for (ip=0;ip<Size();ip++)
64			{
65			b[ip]=d[ip]=a[ip][ip];
66			z[ip]=0.0;
67			}
68			nrot=0;
69			for (i=0;i<100;i++)
70			{
71			sm=0.0;
72			for (ip=0;ip<Size()-1;ip++)
73			{
74			for (iq=ip+1;iq<Size();iq++)
75			sm += fabs(a[ip][iq]);
76			}
77			if (sm == 0.0)
78			{
79			break;
80			}
81			if (i < 3)
82			tresh=0.2sm/ (Size() Size());
83			else
84			tresh=0.0;
85			for (ip=0;ip<Size()-1;ip++)
86			{
87			for (iq=ip+1;iq<Size();iq++)
88			{
89			g=100.0*fabs(a[ip][iq]);
90			if ((i > 3) && ((fabs(d[ip]) +g) == fabs(d[ip])) && ((fabs(d[iq]) +g) == fabs(d[iq])))
91			a[ip][iq]=0.0;
92			else
93			{
94			if (fabs(a[ip][iq]) > tresh)
95			{
96			h=d[iq]-d[ip];
97			if ((fabs(h) +g) == fabs(h))
98			t= (a[ip][iq]) /h;
99			else
100			{
101			theta=0.5*h/ (a[ip][iq]);
102			t=1.0/ (fabs(theta) +sqrt(1.0+theta*theta));
103			if (theta < 0.0) t = -t;
104			}
105			c=1.0/sqrt(1+t*t);
106			s=t*c;
107			tau=s/ (1.0+c);
108			h=t*a[ip][iq];
109			z[ip] -= h;
110			z[iq] += h;
111			d[ip] -= h;
112			d[iq] += h;
113			a[ip][iq]=0.0;
114			for (j=0;j<=ip-1;j++)
115			{
116			ROTATE(a,j,ip,j,iq)
117			}
118			for (j=ip+1;j<=iq-1;j++)
119			{
120			ROTATE(a,ip,j,j,iq)
121			}
122			for (j=iq+1;j<Size();j++)
123			{
124			ROTATE(a,ip,j,iq,j)
125			}
126			for (j=0;j<Size();j++)
127			{
128			ROTATE(mat,j,ip,j,iq)
129			}
130			++nrot;
131			}
132			}
133			}
134			}
135			for (ip=0;ip<Size();ip++)
136			{
137			b[ip] += z[ip];
138			d[ip]=b[ip];
139			z[ip]=0.0;
140			}
141			}
142			if (i==100)
143			{
144			return 1; // there was a problem.
145			}
146			return 0;
147			}
148
149			#undef ROTATE
150
151
152			void Square_Matrix::Complete_Base(int num)
153			{
154			int i,j,k;
155			std::vector<Vector> Tab(Size());
156			for (i=0;i<Size();i++)
157			Tab[i].Resize(Size());
158
159			double pscal,norm;
160			for (i=0;i<num;i++)
161			{
162			for (j=0;j<Size();j++)
163			Tab[i][j]=mat[i][j];
164			Tab[i].Normalize();
165			}
166
167			for (i=num;i<Size();i++)
168			{
169			for (j=0;j<num;j++)
170			{
171			for (k=0;k<Size();k++)
172			Tab[i][k]=0;
173			Tab[i][(i+j) %Size()]=1.0;
174			for (k=0;k<i;k++)
175			{
176			pscal=Tab[i]*Tab[k];
177			Tab[i].Add(Tab[k],-pscal);
178			}
179			norm=Tab[i].Norm();
180			if (norm>1e-4)
181			break;
182
183			}
184			Tab[i].Normalize();
185			}
186
187			for (i=0;i<Size();i++)
188			for (j=0;j<Size();j++)
189			mat[i][j]=Tab[i][j];
190			}
191
192
193
194
195
196			void Cholesky_Matrix::CholeskyDec_intern(void)
197			{
198			int k,i,p;
199			for (k=0;k<Size();++k)
200			{
201			double sumsqakp=0.;
202			for (p=0;p<k;++p)
203			{
204			double akp= (*this)(k,p);
205			sumsqakp+=akp*akp;
206			}
207			double& akk= (*this)(k,k);
208			if (akk<sumsqakp)
209			{
210			std::cerr << "Cholesky : (akk<sumsqakp) " << akk << " < " << sumsqakp << std::endl;
211			throw "Cholesky : (akk<sumsqakp)";
212			}
213			akk=sqrt(akk-sumsqakp);
214			// Display();
215			// cout << "a" << k << k << " " << akk << " sumsqakp " << sumsqakp << endl;
216			for (i=k+1;i<Size();++i)
217			{
218			double sumaipakp=0.0;
219			for (p=0;p<k;++p)
220			{
221			sumaipakp+= (this)(i,p) (*this)(k,p);
222			}
223			double& aik= (*this)(i,k);
224			aik= (aik-sumaipakp) /akk;
225			// Display();
226			// cout << "a" << i << k << " " << aik <<" sumaipakp " << sumaipakp << endl;
227			// cout << aik << " " ;
228			}
229			// cout << endl;
230			}
231			}
232
233
234			void LU_Matrix::LUDec_intern(void)
235			{
236
237			int i,imax,j,k;
238			double big,dum,sum,temp;
239			std::vector<double> vv(Size());
240
241
242			d=1.;
243
244			for (i=0;i<Size();++i)
245			{
246			big=0.0;
247			for (j=0;j<Size();++j)
248			if ((temp=fabs(mat[i][j])) > big) big=temp;
249
250			if (big == 0.0)
251			{
252			d=0.0;
253			return;
254			}
255
256			vv[i]=1.0/big;
257			}
258
259
260			for (j=0;j<Size();++j)
261			{
262			for (i=0;i<j;++i)
263			{
264			sum=mat[i][j];
265			for (k=0;k<i;++k) sum -= mat[i][k]*mat[k][j];
266			mat[i][j]=sum;
267			}
268
269			big=0.0;
270			for (i=j;i<Size();++i)
271			{
272			sum=mat[i][j];
273			for (k=0;k<j;++k)
274			sum -= mat[i][k]*mat[k][j];
275			mat[i][j]=sum;
276			if ((dum=vv[i]*fabs(sum)) >= big)
277			{
278			big=dum;
279			imax=i;
280			}
281			}
282			if (j != imax)
283			{
284			for (k=0;k<Size();++k)
285			{
286			dum=mat[imax][k];
287			mat[imax][k]=mat[j][k];
288			mat[j][k]=dum;
289			}
290			d = - (d);
291			vv[imax]=vv[j];
292			}
293			tab[j]=imax;
294
295			if (mat[j][j] == 0.0)
296			{
297			mat[j][j] = 1e-50;
298			}
299
300			if (j != Size()-1)
301			{
302			dum=1.0/ (mat[j][j]);
303			for (i=j+1;i<Size();++i) mat[i][j] *= dum;
304			}
305			}
306			for (i=0;i<Size();++i)
307			d*=mat[i][i];
308			}
309
310			void LU_Matrix::Solve_Linear_System(Vector &rhsv,Vector &sol)
311			{
312			int i,im1,ii=-1,ip,j;
313			double sum;
314
315			if (d==0.0) return;
316
317			sol=rhsv;
318
319			for (i=0;i<Size();++i)
320			{
321			ip=tab[i];
322			sum=sol[ip];
323			sol[ip]=sol[i];
324			im1=i-1;
325			if (ii>=0)
326			for (j=ii;j<=im1;++j) sum -= mat[i][j]*sol[j];
327			else if (sum!=0.0) ii=i;
328			sol[i]=sum;
329			}
330
331			for (i=Size()-1;i>=0;i--)
332			{
333			sum=sol[i];
334			for (j=i+1;j<Size();++j) sum -= mat[i][j]*sol[j];
335			sol[i]=sum/mat[i][i];
336			}
337			}
338
339
340
341
342			double det(double M[2][2])
343			{
344			return M[0][0]M[1][1]-M[1][0]M[0][1];
345			}
346
347			double trace(double M[2][2])
348			{
349			return M[0][0]+M[1][1];
350			}
351
352			double prodcc(double M1[2][2],double M2[2][2])
353			{
354			double p=0.;
355			for (int i=0;i<2;++i)
356			for (int j=0;j<2;++j)
357			p+=M1[i][j]*M2[i][j];
358			return p;
359			}
360
361
362			void inverse(double M[2][2],double I[2][2]) // I : inverse de la matrice M
363			{
364			double d = det(M);
365			I[0][0] = M[1][1]/d;
366			I[1][1] = M[0][0]/d;
367			I[0][1] = -M[0][1]/d;
368			I[1][0] = -M[1][0]/d;
369			}
370
371
372			void multi(double A[2][2],double B[2][2],double C[2][2]) // Effectue la multiplication de A et B et rentre le résultat dans la matrice C
373			{
374			for (int i=0;i<2;i++)
375			{
376			for (int j=0;j<2;j++)
377			{
378			C[i][j] = 0;
379			for (int k=0;k<2;k++)
380			{
381			C[i][j]+=A[i][k]*B[k][j];
382			}
383			}
384			}
385			}
386
387			double multi(double V1[2],double M[2][2],double V2[2]) // Effectue la multiplication Ai Mij Cj
388			{
389			double res=0.;
390			for (int i=0;i<2;i++)
391			{
392			for (int j=0;j<2;j++)
393			{
394			res+=V1[i]M[i][j]V2[j];
395			}
396			}
397			return res;
398			}
399