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:	283
Committed:	Tue Sep 13 21:11:20 2011 UTC (13 years, 8 months ago) by francois
File size:	7719 byte(s)
Log Message:	structure de l'écriture
#	User	Rev	Content
1	francois	283	//------------------------------------------------------------
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) {
80			delete [] aaug;
81			return 0;
82			};
83			// j ai mon pivot max non nul
84			for (int ii=j+1;ii<n;ii++) aaug(ii,j)=aaug(ii,j)/aaug(j,j);
85
86			for (int iii=j+1;iii<n;iii++)
87			for (int k=j+1;k<n+1;k++)
88			aaug(iii,k)=aaug(iii,k)-aaug(iii,j)*aaug(j,k);
89			}
90
91			if (fabs(aaug(n-1,n-1))<0.00000000000001) {
92			delete [] aaug;
93			return 0;
94			};
95			// back substitution
96			x[n-1]=aaug(n-1,n)/aaug(n-1,n-1);
97			for (int j=n-2;j>-1;j--)
98			{
99			x[j]=aaug(j,n);
100			for (int k=n-1;k>j;k--)
101			x[j]=x[j]-aaug(j,k)*x[k];
102			x[j]=x[j]/aaug(j,j);
103			}
104
105
106			delete [] aaug;
107			return 1;
108			#undef a // enlever le (i,j) pour enlever le warning
109			#undef aaug
110			}
111
112			int OT_SYSTEME::resoud_QR(double m, double b, double *x, int dim1,int dim2)
113			{
114			double mat=new double[dim1(dim2+1)];
115			double *r=new double[dim2];
116			#define mat(i,j) ((mat+(i)(dim2+1)+(j)))
117			#define m(i,j) ((m+(i)dim2+(j)))
118			// initialisation avec mat augment�e
119			for (int i=0;i<dim1;i++)
120			{
121			for (int j=0;j<dim2;j++)
122			mat(i,j)=m(i,j);
123			mat(i,dim2)=b[i];
124			}
125			// QR
126			for (int k=0;k<dim2;k++)
127			{
128			double somme=0.;
129			for (int i=k;i<dim1;i++)
130			somme=somme+mat(i,k)*mat(i,k);
131			double alpha=sqrt(somme);
132			double pivot=mat(k,k);
133			if (pivot>0) alpha=-alpha;
134			r[k]=alpha;
135			if (OPERATEUR::egal(alpha,0.,0.000001))
136			{
137			delete [] mat;
138			delete [] r;
139			return 0;
140			}
141			double beta=somme-pivot*alpha;
142			pivot=pivot-alpha;
143			mat(k,k)=pivot;
144			for (int j=k+1;j<dim2+1;j++)
145			{
146			double somme=0.;
147			for (int i=k;i<dim1;i++)
148			somme=somme+mat(i,k)*mat(i,j);
149			double gamma=somme/beta;
150			for (int i=k;i<dim1;i++)
151			mat(i,j)=mat(i,j)-gamma*mat(i,k);
152			}
153			}
154			// resolution
155			x[dim2-1]=mat(dim2-1,dim2)/r[dim2-1];
156			for (int j=dim2-2;j>-1;j--)
157			{
158			x[j]=mat(j,dim2);
159			for (int k=dim2-1;k>j;k--)
160			x[j]=x[j]-mat(j,k)*x[k];
161			x[j]=x[j]/r[j];
162			}
163
164
165
166			#undef mat
167			delete [] mat;
168			delete [] r;
169			return 1;
170			}
171
172
173
174
175
176
177			int OT_SYSTEME::decompose_LU(double a,double res,int dim)
178			{
179			#define a(i,j) ((a+(i)(dim)+(j)))
180			#define res(i,j) ((res+(i)(dim)+(j)))
181
182			if (OPERATEUR::egal(a(0,0),0.,0.00001)==1) return 0;
183			for (int i=0;i<dim;i++)
184			res(i,0)=a(i,0);
185			for (int j=1;j<dim;j++)
186			res(0,j)=a(0,j)/res(0,0);
187
188
189			for (int j=1;j<dim;j++)
190			{
191			for (int i=j;i<dim;i++)
192			{
193			double sum=0.;
194			for (int k=0;k<j;k++)
195			sum=sum+res(i,k)*res(k,j);
196			res(i,j)=a(i,j)-sum;
197			}
198			for (int i=j+1;i<dim;i++)
199			{
200			double sum=0.;
201			if (OPERATEUR::egal(res(j,j),0.,0.00001)==1) return 0;
202			for (int k=0;k<j;k++)
203			sum=sum+res(j,k)*res(k,i);
204			res(j,i)=(a(j,i)-sum)/res(j,j);
205			}
206
207			}
208			return 1;
209
210			#undef a
211			#undef res
212			}
213
214
215			int OT_SYSTEME::resoud_LU(double matlu,double b,double *x,int dim)
216			{
217			#define matlu(i,j) ((matlu+(i)(dim)+(j)))
218
219			double* bb=new double[dim];
220			bb[0]=b[0]/matlu(0,0);
221			for (int i=1;i<dim;i++)
222			{
223			double sum=0.0;
224			for (int k=0;k<i;k++)
225			sum=sum+matlu(i,k)*bb[k];
226			bb[i]=(b[i]-sum)/matlu(i,i);
227			}
228
229			x[dim-1]=bb[dim-1];
230			for (int j=dim-2;j>=0;j--)
231			{
232			double sum=0.0;
233			for (int k=j+1;k<dim;k++)
234			sum=sum+matlu(j,k)*x[k];
235			x[j]=bb[j]-sum;
236
237			}
238
239			#undef lu
240
241			delete [] bb;
242			return 1;
243			}
244
245
246
247
248
249
250			void OT_SYSTEME::get_inv(double mat1,double mat2,int n,int *ierr)
251			{
252			#define mat(i,j) ((mataug+(j)n+i))
253			#define mat1(i,j) ((mat1+(j)n+i))
254			#define mat2(i,j) ((mat2+(j)n+i))
255
256			double mataug=(double)calloc(2nn,sizeof(double));
257
258			for (int i=0;i<n;i++)
259			{
260			for (int j=0;j<n;j++)
261			mat(i,j)=mat1(i,j);
262			mat(i,i+n)=1.;
263			}
264
265
266			*ierr=0;
267
268			for (int j=0;j<n;j++)
269			{
270			//recherche du pivot max total de la colonne j
271			double pivot=fabs(mat(j,j));
272			int ligne_pivot=j;
273			for (int i=j+1;i<n;i++)
274			{
275			if (fabs(mat(i,j))>pivot)
276			{
277			pivot=fabs(mat(i,j));
278			ligne_pivot=i;
279			}
280			}
281			//tester si le pivot max est nul
282			if (pivot<0.000000001) {
283			*ierr=1;
284			free(mataug);
285			return;
286			}
287			//tester si le pivot actuel est max
288			if (ligne_pivot!=j)
289			{
290			for (int i=j;i<2*n;i++)
291			{
292			double temp=mat(j,i);
293			mat(j,i)=mat(ligne_pivot,i);
294			mat(ligne_pivot,i)=temp;
295			}
296			}
297			//remplir la colonne j de 0
298			for (int i=j+1;i<n;i++)
299			mat(i,j)=-mat(i,j)/mat(j,j);
300			for (int i=j+1;i<n;i++)
301			for (int k=j+1;k<2*n;k++)
302			mat(i,k)=mat(i,k)+mat(i,j)*mat(j,k);
303			}
304
305			for (int i=0;i<n;i++)
306			{
307			mat2(n-1,i)=mat(n-1,n+i)/mat(n-1,n-1);
308			for (int j=n-2;j>-1;j--)
309			{
310			mat2(j,i)=mat(j,n+i);
311			for (int k=n-1;k>j;k--)
312			mat2(j,i)=mat2(j,i)-mat2(k,i)*mat(j,k);
313			mat2(j,i)=mat2(j,i)/mat(j,j);
314			}
315			}
316
317			free(mataug);
318			#undef mat
319			#undef mat1
320			#undef mat2
321			}
322