outil/src/ot_systeme.cpp

//------------------------------------------------------------
//------------------------------------------------------------
//  MAGiC
//  Jean Christophe Cuillière et Vincent FRANCOIS
//  Département de Génie Mécanique - UQTR
//------------------------------------------------------------
//  Le projet MAGIC est un projet de recherche du département
//  de génie mécanique de l'Université du Québec à
//  Trois Rivières
//  Les librairies ne peuvent être utilisées sans l'accord
//  des auteurs (contact : francois@uqtr.ca)
//------------------------------------------------------------
//------------------------------------------------------------
//
//       ot_systeme.cpp
//
//------------------------------------------------------------
//------------------------------------------------------------
//  COPYRIGHT 2000
//  Version du 02/03/2006 à 11H23
//------------------------------------------------------------
//------------------------------------------------------------


#include "gestionversion.h"

#pragma hdrstop
#include <stdlib.h>
#include <math.h>
#include "ot_systeme.h"
#include "ot_mathematique.h"

//---------------------------------------------------------------------------

OT_SYSTEME::OT_SYSTEME()
{
}


OT_SYSTEME::~OT_SYSTEME()
{
}


int OT_SYSTEME::resoud_gauss(double *a,double* b,double *x,int n)
{
double *aaug=new double[(n+1)*n];
#define a(i,j) (*(a+(i)*(n)+j))
#define aaug(i,j) (*(aaug+(i)*(n+1)+j))
for (int i=0;i<n;i++)
        for (int j=0;j<n;j++)
                aaug(i,j)=a(i,j);
for (int k=0;k<n;k++)
        aaug(k,n)=b[k];  
for (int j=0;j<n-1;j++)
        {
        // recherche le pivot
        double pvt=fabs(aaug(j,j));
        int ligne_pivot_debut=j;
        int ligne_pivot=j;

        for (int i=j+1;i<n;i++)
                if (fabs(aaug(i,j))>pvt)
                        {
                        pvt=fabs(aaug(i,j));
                        ligne_pivot=i;
                        }
        if (ligne_pivot!=ligne_pivot_debut) // trouver meilleur
                {
                for (int col=j;col<n+1;col++)
                        {
                        double temp=aaug(ligne_pivot_debut,col);
                        aaug(ligne_pivot_debut,col)=aaug(ligne_pivot,col);
                        aaug(ligne_pivot,col)=temp;
                        }
                }
        // j ai mon pivot max
        if (fabs(aaug(j,j))< 0.00000000000001) {delete [] aaug;return 0;};
        // j ai mon pivot max non nul
        for (int ii=j+1;ii<n;ii++) aaug(ii,j)=aaug(ii,j)/aaug(j,j);

        for (int iii=j+1;iii<n;iii++)
                for (int k=j+1;k<n+1;k++)
                        aaug(iii,k)=aaug(iii,k)-aaug(iii,j)*aaug(j,k);
        }

if (fabs(aaug(n-1,n-1))<0.00000000000001) {delete [] aaug;return 0;};
// back substitution
x[n-1]=aaug(n-1,n)/aaug(n-1,n-1);
for (int j=n-2;j>-1;j--)
        {
        x[j]=aaug(j,n);
        for (int k=n-1;k>j;k--)
                x[j]=x[j]-aaug(j,k)*x[k];
        x[j]=x[j]/aaug(j,j);
        }


delete [] aaug;
return 1;
#undef a // enlever le  (i,j) pour enlever le warning
#undef aaug
}

int OT_SYSTEME::resoud_QR(double *m, double *b, double *x, int dim1,int dim2)
{
double *mat=new double[dim1*(dim2+1)];
double *r=new double[dim2];
#define mat(i,j) (*(mat+(i)*(dim2+1)+(j)))
#define m(i,j) (*(m+(i)*dim2+(j)))
// initialisation avec mat augmentée
for (int i=0;i<dim1;i++)
        {
        for (int j=0;j<dim2;j++)
                mat(i,j)=m(i,j);
        mat(i,dim2)=b[i];
        }
// QR
for (int k=0;k<dim2;k++)
        {
        double somme=0.;
        for (int i=k;i<dim1;i++)
                somme=somme+mat(i,k)*mat(i,k);
        double alpha=sqrt(somme);
        double pivot=mat(k,k);
        if (pivot>0) alpha=-alpha;
        r[k]=alpha;
        if (OPERATEUR::egal(alpha,0.,0.000001))
                {
                delete [] mat;
                delete [] r;
                return 0;
                }
        double beta=somme-pivot*alpha;
        pivot=pivot-alpha;
        mat(k,k)=pivot;
        for (int j=k+1;j<dim2+1;j++)
                {
                double somme=0.;
                for (int i=k;i<dim1;i++)
                        somme=somme+mat(i,k)*mat(i,j);
                double gamma=somme/beta;
                for (int i=k;i<dim1;i++)
                        mat(i,j)=mat(i,j)-gamma*mat(i,k);
                }
        }
// resolution
x[dim2-1]=mat(dim2-1,dim2)/r[dim2-1];
for (int j=dim2-2;j>-1;j--)
        {
        x[j]=mat(j,dim2);
        for (int k=dim2-1;k>j;k--)
                x[j]=x[j]-mat(j,k)*x[k];
        x[j]=x[j]/r[j];
        }


#undef mat
delete [] mat;
delete [] r;
return 1;
}


int OT_SYSTEME::decompose_LU(double *a,double* res,int dim)
{
#define a(i,j) (*(a+(i)*(dim)+(j)))
#define res(i,j) (*(res+(i)*(dim)+(j)))

if (OPERATEUR::egal(a(0,0),0.,0.00001)==1) return 0;
for (int i=0;i<dim;i++)
        res(i,0)=a(i,0);
for (int j=1;j<dim;j++)
        res(0,j)=a(0,j)/res(0,0);


for (int j=1;j<dim;j++)
        {
        for (int i=j;i<dim;i++)
                {
                double sum=0.;
                for (int k=0;k<j;k++)
                        sum=sum+res(i,k)*res(k,j);
                res(i,j)=a(i,j)-sum;
                }
        for (int i=j+1;i<dim;i++)
                {
                double sum=0.;
                if (OPERATEUR::egal(res(j,j),0.,0.00001)==1) return 0;
                for (int k=0;k<j;k++)
                        sum=sum+res(j,k)*res(k,i);
                res(j,i)=(a(j,i)-sum)/res(j,j);
                }

        }
return 1;

#undef a
#undef res
}


int OT_SYSTEME::resoud_LU(double *matlu,double *b,double *x,int dim)
{
#define matlu(i,j) (*(matlu+(i)*(dim)+(j)))

double* bb=new double[dim];
bb[0]=b[0]/matlu(0,0);
for (int i=1;i<dim;i++)
        {
        double sum=0.0;
        for (int k=0;k<i;k++)
                sum=sum+matlu(i,k)*bb[k];
        bb[i]=(b[i]-sum)/matlu(i,i);
        }

x[dim-1]=bb[dim-1];
for (int j=dim-2;j>=0;j--)
        {
        double sum=0.0;
        for (int k=j+1;k<dim;k++)
                sum=sum+matlu(j,k)*x[k];
        x[j]=bb[j]-sum;

        }

#undef lu

delete [] bb;
return 1;
}


void OT_SYSTEME::get_inv(double *mat1,double *mat2,int n,int *ierr)
{
#define mat(i,j) (*(mataug+(j)*n+i))
#define mat1(i,j) (*(mat1+(j)*n+i))
#define mat2(i,j) (*(mat2+(j)*n+i))

double *mataug=(double*)calloc(2*n*n,sizeof(double));
 
for (int i=0;i<n;i++)
 {
 for (int j=0;j<n;j++)
  mat(i,j)=mat1(i,j);
 mat(i,i+n)=1.;
 }
 

*ierr=0;

for (int j=0;j<n;j++)
 {
 //recherche du pivot max total de la colonne j
 double pivot=fabs(mat(j,j));
 int ligne_pivot=j;
 for (int i=j+1;i<n;i++)
  {
  if (fabs(mat(i,j))>pivot)
   {
   pivot=fabs(mat(i,j));
   ligne_pivot=i;
   }
  }
 //tester si le pivot max est nul
 if (pivot<0.000000001) {*ierr=1;free(mataug);return;}
 //tester si le pivot actuel est max
 if (ligne_pivot!=j)
  {
  for (int i=j;i<2*n;i++)
   {
   double temp=mat(j,i);
   mat(j,i)=mat(ligne_pivot,i);
   mat(ligne_pivot,i)=temp;
   }
  }
 //remplir la colonne j de 0
 for (int i=j+1;i<n;i++)
  mat(i,j)=-mat(i,j)/mat(j,j);
 for (int i=j+1;i<n;i++)
  for (int k=j+1;k<2*n;k++)
   mat(i,k)=mat(i,k)+mat(i,j)*mat(j,k);
 }

for (int i=0;i<n;i++)
 {
 mat2(n-1,i)=mat(n-1,n+i)/mat(n-1,n-1);
 for (int j=n-2;j>-1;j--)
  {
  mat2(j,i)=mat(j,n+i);
  for (int k=n-1;k>j;k--)
   mat2(j,i)=mat2(j,i)-mat2(k,i)*mat(j,k);
  mat2(j,i)=mat2(j,i)/mat(j,j);
  }
 }

free(mataug);
#undef mat
#undef mat1
#undef mat2
}

Revision:	169
Committed:	Fri Feb 13 23:04:14 2009 UTC (16 years, 3 months ago) by francois
Original Path:	*magic/lib/outil/outil/src/ot_systeme.cpp*
File size:	7306 byte(s)
Log Message:	Resolution de bug avec la version de gcc 4.3.2 --> mise la norme de c++
#	Rev	Content
1	5	//------------------------------------------------------------
2		//------------------------------------------------------------
3		// MAGiC
4		// Jean Christophe Cuillière et Vincent FRANCOIS
5		// Département de Génie Mécanique - UQTR
6		//------------------------------------------------------------
7		// Le projet MAGIC est un projet de recherche du département
8		// de génie mécanique de l'Université du Québec à
9		// Trois Rivières
10		// Les librairies ne peuvent être utilisées sans l'accord
11		// des auteurs (contact : francois@uqtr.ca)
12		//------------------------------------------------------------
13		//------------------------------------------------------------
14		//
15		// ot_systeme.cpp
16		//
17		//------------------------------------------------------------
18		//------------------------------------------------------------
19		// COPYRIGHT 2000
20		// Version du 02/03/2006 à 11H23
21		//------------------------------------------------------------
22		//------------------------------------------------------------
23
24
25		#include "gestionversion.h"
26
27		#pragma hdrstop
28		#include <stdlib.h>
29		#include <math.h>
30		#include "ot_systeme.h"
31		#include "ot_mathematique.h"
32
33		//---------------------------------------------------------------------------
34
35		OT_SYSTEME::OT_SYSTEME()
36		{
37		}
38
39
40
41		OT_SYSTEME::~OT_SYSTEME()
42		{
43		}
44
45
46		int OT_SYSTEME::resoud_gauss(double a,double b,double *x,int n)
47		{
48		double aaug=new double[(n+1)n];
49		#define a(i,j) ((a+(i)(n)+j))
50		#define aaug(i,j) ((aaug+(i)(n+1)+j))
51		for (int i=0;i<n;i++)
52		for (int j=0;j<n;j++)
53		aaug(i,j)=a(i,j);
54		for (int k=0;k<n;k++)
55		aaug(k,n)=b[k];
56		for (int j=0;j<n-1;j++)
57		{
58		// recherche le pivot
59		double pvt=fabs(aaug(j,j));
60		int ligne_pivot_debut=j;
61		int ligne_pivot=j;
62
63		for (int i=j+1;i<n;i++)
64		if (fabs(aaug(i,j))>pvt)
65		{
66		pvt=fabs(aaug(i,j));
67		ligne_pivot=i;
68		}
69		if (ligne_pivot!=ligne_pivot_debut) // trouver meilleur
70		{
71		for (int col=j;col<n+1;col++)
72		{
73		double temp=aaug(ligne_pivot_debut,col);
74		aaug(ligne_pivot_debut,col)=aaug(ligne_pivot,col);
75		aaug(ligne_pivot,col)=temp;
76		}
77		}
78		// j ai mon pivot max
79		if (fabs(aaug(j,j))< 0.00000000000001) {delete [] aaug;return 0;};
80		// j ai mon pivot max non nul
81		for (int ii=j+1;ii<n;ii++) aaug(ii,j)=aaug(ii,j)/aaug(j,j);
82
83		for (int iii=j+1;iii<n;iii++)
84		for (int k=j+1;k<n+1;k++)
85		aaug(iii,k)=aaug(iii,k)-aaug(iii,j)*aaug(j,k);
86		}
87
88		if (fabs(aaug(n-1,n-1))<0.00000000000001) {delete [] aaug;return 0;};
89		// back substitution
90		x[n-1]=aaug(n-1,n)/aaug(n-1,n-1);
91		for (int j=n-2;j>-1;j--)
92		{
93		x[j]=aaug(j,n);
94		for (int k=n-1;k>j;k--)
95		x[j]=x[j]-aaug(j,k)*x[k];
96		x[j]=x[j]/aaug(j,j);
97		}
98
99
100		delete [] aaug;
101		return 1;
102		#undef a // enlever le (i,j) pour enlever le warning
103		#undef aaug
104		}
105
106		int OT_SYSTEME::resoud_QR(double m, double b, double *x, int dim1,int dim2)
107		{
108		double mat=new double[dim1(dim2+1)];
109		double *r=new double[dim2];
110		#define mat(i,j) ((mat+(i)(dim2+1)+(j)))
111		#define m(i,j) ((m+(i)dim2+(j)))
112		// initialisation avec mat augmentée
113		for (int i=0;i<dim1;i++)
114		{
115		for (int j=0;j<dim2;j++)
116		mat(i,j)=m(i,j);
117		mat(i,dim2)=b[i];
118		}
119		// QR
120		for (int k=0;k<dim2;k++)
121		{
122		double somme=0.;
123		for (int i=k;i<dim1;i++)
124		somme=somme+mat(i,k)*mat(i,k);
125		double alpha=sqrt(somme);
126		double pivot=mat(k,k);
127		if (pivot>0) alpha=-alpha;
128		r[k]=alpha;
129		if (OPERATEUR::egal(alpha,0.,0.000001))
130		{
131		delete [] mat;
132		delete [] r;
133		return 0;
134		}
135		double beta=somme-pivot*alpha;
136		pivot=pivot-alpha;
137		mat(k,k)=pivot;
138		for (int j=k+1;j<dim2+1;j++)
139		{
140		double somme=0.;
141		for (int i=k;i<dim1;i++)
142		somme=somme+mat(i,k)*mat(i,j);
143		double gamma=somme/beta;
144		for (int i=k;i<dim1;i++)
145		mat(i,j)=mat(i,j)-gamma*mat(i,k);
146		}
147		}
148		// resolution
149		x[dim2-1]=mat(dim2-1,dim2)/r[dim2-1];
150		for (int j=dim2-2;j>-1;j--)
151		{
152		x[j]=mat(j,dim2);
153		for (int k=dim2-1;k>j;k--)
154		x[j]=x[j]-mat(j,k)*x[k];
155		x[j]=x[j]/r[j];
156		}
157
158
159
160		#undef mat
161		delete [] mat;
162		delete [] r;
163		return 1;
164		}
165
166
167
168
169
170
171		int OT_SYSTEME::decompose_LU(double a,double res,int dim)
172		{
173		#define a(i,j) ((a+(i)(dim)+(j)))
174		#define res(i,j) ((res+(i)(dim)+(j)))
175
176		if (OPERATEUR::egal(a(0,0),0.,0.00001)==1) return 0;
177		for (int i=0;i<dim;i++)
178		res(i,0)=a(i,0);
179		for (int j=1;j<dim;j++)
180		res(0,j)=a(0,j)/res(0,0);
181
182
183		for (int j=1;j<dim;j++)
184		{
185		for (int i=j;i<dim;i++)
186		{
187		double sum=0.;
188		for (int k=0;k<j;k++)
189		sum=sum+res(i,k)*res(k,j);
190		res(i,j)=a(i,j)-sum;
191		}
192		for (int i=j+1;i<dim;i++)
193		{
194		double sum=0.;
195		if (OPERATEUR::egal(res(j,j),0.,0.00001)==1) return 0;
196		for (int k=0;k<j;k++)
197		sum=sum+res(j,k)*res(k,i);
198		res(j,i)=(a(j,i)-sum)/res(j,j);
199		}
200
201		}
202		return 1;
203
204		#undef a
205		#undef res
206		}
207
208
209		int OT_SYSTEME::resoud_LU(double matlu,double b,double *x,int dim)
210		{
211		#define matlu(i,j) ((matlu+(i)(dim)+(j)))
212
213		double* bb=new double[dim];
214		bb[0]=b[0]/matlu(0,0);
215		for (int i=1;i<dim;i++)
216		{
217		double sum=0.0;
218		for (int k=0;k<i;k++)
219		sum=sum+matlu(i,k)*bb[k];
220		bb[i]=(b[i]-sum)/matlu(i,i);
221		}
222
223		x[dim-1]=bb[dim-1];
224		for (int j=dim-2;j>=0;j--)
225		{
226		double sum=0.0;
227		for (int k=j+1;k<dim;k++)
228		sum=sum+matlu(j,k)*x[k];
229		x[j]=bb[j]-sum;
230
231		}
232
233		#undef lu
234
235		delete [] bb;
236		return 1;
237		}
238
239
240
241
242
243
244		void OT_SYSTEME::get_inv(double mat1,double mat2,int n,int *ierr)
245		{
246		#define mat(i,j) ((mataug+(j)n+i))
247		#define mat1(i,j) ((mat1+(j)n+i))
248		#define mat2(i,j) ((mat2+(j)n+i))
249
250		double mataug=(double)calloc(2nn,sizeof(double));
251
252		for (int i=0;i<n;i++)
253		{
254		for (int j=0;j<n;j++)
255		mat(i,j)=mat1(i,j);
256		mat(i,i+n)=1.;
257		}
258
259
260		*ierr=0;
261
262		for (int j=0;j<n;j++)
263		{
264		//recherche du pivot max total de la colonne j
265		double pivot=fabs(mat(j,j));
266		int ligne_pivot=j;
267		for (int i=j+1;i<n;i++)
268		{
269		if (fabs(mat(i,j))>pivot)
270		{
271		pivot=fabs(mat(i,j));
272		ligne_pivot=i;
273		}
274		}
275		//tester si le pivot max est nul
276		if (pivot<0.000000001) {*ierr=1;free(mataug);return;}
277		//tester si le pivot actuel est max
278		if (ligne_pivot!=j)
279		{
280		for (int i=j;i<2*n;i++)
281		{
282		double temp=mat(j,i);
283		mat(j,i)=mat(ligne_pivot,i);
284		mat(ligne_pivot,i)=temp;
285		}
286		}
287		//remplir la colonne j de 0
288		for (int i=j+1;i<n;i++)
289		mat(i,j)=-mat(i,j)/mat(j,j);
290		for (int i=j+1;i<n;i++)
291		for (int k=j+1;k<2*n;k++)
292		mat(i,k)=mat(i,k)+mat(i,j)*mat(j,k);
293		}
294
295		for (int i=0;i<n;i++)
296		{
297		mat2(n-1,i)=mat(n-1,n+i)/mat(n-1,n-1);
298		for (int j=n-2;j>-1;j--)
299		{
300		mat2(j,i)=mat(j,n+i);
301		for (int k=n-1;k>j;k--)
302		mat2(j,i)=mat2(j,i)-mat2(k,i)*mat(j,k);
303		mat2(j,i)=mat2(j,i)/mat(j,j);
304		}
305		}
306
307		free(mataug);
308		#undef mat
309		#undef mat1
310		#undef mat2
311		}
312