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
#	Content
1	//####//------------------------------------------------------------
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
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