outil/src/ot_systeme.cpp

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  MAGiC
//####//  Jean Christophe Cuilliere et Vincent FRANCOIS
//####//  Departement de Genie Mecanique - UQTR
//####//------------------------------------------------------------
//####//  MAGIC est un projet de recherche de l equipe ERICCA
//####//  du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
//####//  http://www.uqtr.ca/ericca
//####//  http://www.uqtr.ca/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       ot_systeme.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:53:59 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------


#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;
    };
    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)))
    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];
    }
    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);
        }
    }
    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++)
    {
        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;
            }
        }
        if (pivot<0.000000001) {
            *ierr=1;
            free(mataug);
            return;
        }
        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;
            }
        }
        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:	1156
Committed:	Thu Jun 13 22:02:48 2024 UTC (14 months, 2 weeks ago) by francois
File size:	7448 byte(s)
Log Message:	compatibilité Ubuntu 22.04 Suppression des refeences à Windows Ajout d'une banière
#	User	Rev	Content
1	francois	1156	//####//------------------------------------------------------------
2			//####//------------------------------------------------------------
3			//####// MAGiC
4			//####// Jean Christophe Cuilliere et Vincent FRANCOIS
5			//####// Departement de Genie Mecanique - UQTR
6			//####//------------------------------------------------------------
7			//####// MAGIC est un projet de recherche de l equipe ERICCA
8			//####// du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
9			//####// http://www.uqtr.ca/ericca
10			//####// http://www.uqtr.ca/
11			//####//------------------------------------------------------------
12			//####//------------------------------------------------------------
13			//####//
14			//####// ot_systeme.cpp
15			//####//
16			//####//------------------------------------------------------------
17			//####//------------------------------------------------------------
18			//####// COPYRIGHT 2000-2024
19			//####// jeu 13 jun 2024 11:53:59 EDT
20			//####//------------------------------------------------------------
21			//####//------------------------------------------------------------
22	francois	283
23
24
25			#pragma hdrstop
26			#include <stdlib.h>
27			#include <math.h>
28			#include "ot_systeme.h"
29			#include "ot_mathematique.h"
30
31
32			OT_SYSTEME::OT_SYSTEME()
33			{
34			}
35
36
37
38			OT_SYSTEME::~OT_SYSTEME()
39			{
40			}
41
42
43			int OT_SYSTEME::resoud_gauss(double a,double b,double *x,int n)
44			{
45			double aaug=new double[(n+1)n];
46			#define a(i,j) ((a+(i)(n)+j))
47			#define aaug(i,j) ((aaug+(i)(n+1)+j))
48			for (int i=0;i<n;i++)
49			for (int j=0;j<n;j++)
50			aaug(i,j)=a(i,j);
51			for (int k=0;k<n;k++)
52			aaug(k,n)=b[k];
53			for (int j=0;j<n-1;j++)
54			{
55			// recherche le pivot
56			double pvt=fabs(aaug(j,j));
57			int ligne_pivot_debut=j;
58			int ligne_pivot=j;
59
60			for (int i=j+1;i<n;i++)
61			if (fabs(aaug(i,j))>pvt)
62			{
63			pvt=fabs(aaug(i,j));
64			ligne_pivot=i;
65			}
66			if (ligne_pivot!=ligne_pivot_debut) // trouver meilleur
67			{
68			for (int col=j;col<n+1;col++)
69			{
70			double temp=aaug(ligne_pivot_debut,col);
71			aaug(ligne_pivot_debut,col)=aaug(ligne_pivot,col);
72			aaug(ligne_pivot,col)=temp;
73			}
74			}
75			// j ai mon pivot max
76			if (fabs(aaug(j,j))< 0.00000000000001) {
77			delete [] aaug;
78			return 0;
79			};
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) {
89			delete [] aaug;
90			return 0;
91			};
92			x[n-1]=aaug(n-1,n)/aaug(n-1,n-1);
93			for (int j=n-2;j>-1;j--)
94			{
95			x[j]=aaug(j,n);
96			for (int k=n-1;k>j;k--)
97			x[j]=x[j]-aaug(j,k)*x[k];
98			x[j]=x[j]/aaug(j,j);
99			}
100
101
102			delete [] aaug;
103			return 1;
104			#undef a // enlever le (i,j) pour enlever le warning
105			#undef aaug
106			}
107
108			int OT_SYSTEME::resoud_QR(double m, double b, double *x, int dim1,int dim2)
109			{
110			double mat=new double[dim1(dim2+1)];
111			double *r=new double[dim2];
112			#define mat(i,j) ((mat+(i)(dim2+1)+(j)))
113			#define m(i,j) ((m+(i)dim2+(j)))
114			for (int i=0;i<dim1;i++)
115			{
116			for (int j=0;j<dim2;j++)
117			mat(i,j)=m(i,j);
118			mat(i,dim2)=b[i];
119			}
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			x[dim2-1]=mat(dim2-1,dim2)/r[dim2-1];
149			for (int j=dim2-2;j>-1;j--)
150			{
151			x[j]=mat(j,dim2);
152			for (int k=dim2-1;k>j;k--)
153			x[j]=x[j]-mat(j,k)*x[k];
154			x[j]=x[j]/r[j];
155			}
156
157
158
159			#undef mat
160			delete [] mat;
161			delete [] r;
162			return 1;
163			}
164
165
166
167
168
169
170			int OT_SYSTEME::decompose_LU(double a,double res,int dim)
171			{
172			#define a(i,j) ((a+(i)(dim)+(j)))
173			#define res(i,j) ((res+(i)(dim)+(j)))
174
175			if (OPERATEUR::egal(a(0,0),0.,0.00001)==1) return 0;
176			for (int i=0;i<dim;i++)
177			res(i,0)=a(i,0);
178			for (int j=1;j<dim;j++)
179			res(0,j)=a(0,j)/res(0,0);
180
181
182			for (int j=1;j<dim;j++)
183			{
184			for (int i=j;i<dim;i++)
185			{
186			double sum=0.;
187			for (int k=0;k<j;k++)
188			sum=sum+res(i,k)*res(k,j);
189			res(i,j)=a(i,j)-sum;
190			}
191			for (int i=j+1;i<dim;i++)
192			{
193			double sum=0.;
194			if (OPERATEUR::egal(res(j,j),0.,0.00001)==1) return 0;
195			for (int k=0;k<j;k++)
196			sum=sum+res(j,k)*res(k,i);
197			res(j,i)=(a(j,i)-sum)/res(j,j);
198			}
199
200			}
201			return 1;
202
203			#undef a
204			#undef res
205			}
206
207
208			int OT_SYSTEME::resoud_LU(double matlu,double b,double *x,int dim)
209			{
210			#define matlu(i,j) ((matlu+(i)(dim)+(j)))
211
212			double* bb=new double[dim];
213			bb[0]=b[0]/matlu(0,0);
214			for (int i=1;i<dim;i++)
215			{
216			double sum=0.0;
217			for (int k=0;k<i;k++)
218			sum=sum+matlu(i,k)*bb[k];
219			bb[i]=(b[i]-sum)/matlu(i,i);
220			}
221
222			x[dim-1]=bb[dim-1];
223			for (int j=dim-2;j>=0;j--)
224			{
225			double sum=0.0;
226			for (int k=j+1;k<dim;k++)
227			sum=sum+matlu(j,k)*x[k];
228			x[j]=bb[j]-sum;
229
230			}
231
232			#undef lu
233
234			delete [] bb;
235			return 1;
236			}
237
238
239
240
241
242
243			void OT_SYSTEME::get_inv(double mat1,double mat2,int n,int *ierr)
244			{
245			#define mat(i,j) ((mataug+(j)n+i))
246			#define mat1(i,j) ((mat1+(j)n+i))
247			#define mat2(i,j) ((mat2+(j)n+i))
248
249			double mataug=(double)calloc(2nn,sizeof(double));
250
251			for (int i=0;i<n;i++)
252			{
253			for (int j=0;j<n;j++)
254			mat(i,j)=mat1(i,j);
255			mat(i,i+n)=1.;
256			}
257
258
259			*ierr=0;
260
261			for (int j=0;j<n;j++)
262			{
263			double pivot=fabs(mat(j,j));
264			int ligne_pivot=j;
265			for (int i=j+1;i<n;i++)
266			{
267			if (fabs(mat(i,j))>pivot)
268			{
269			pivot=fabs(mat(i,j));
270			ligne_pivot=i;
271			}
272			}
273			if (pivot<0.000000001) {
274			*ierr=1;
275			free(mataug);
276			return;
277			}
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			for (int i=j+1;i<n;i++)
288			mat(i,j)=-mat(i,j)/mat(j,j);
289			for (int i=j+1;i<n;i++)
290			for (int k=j+1;k<2*n;k++)
291			mat(i,k)=mat(i,k)+mat(i,j)*mat(j,k);
292			}
293
294			for (int i=0;i<n;i++)
295			{
296			mat2(n-1,i)=mat(n-1,n+i)/mat(n-1,n-1);
297			for (int j=n-2;j>-1;j--)
298			{
299			mat2(j,i)=mat(j,n+i);
300			for (int k=n-1;k>j;k--)
301			mat2(j,i)=mat2(j,i)-mat2(k,i)*mat(j,k);
302			mat2(j,i)=mat2(j,i)/mat(j,j);
303			}
304			}
305
306			free(mataug);
307			#undef mat
308			#undef mat1
309			#undef mat2
310			}
311