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
//------------------------------------------------------------
//------------------------------------------------------------


#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:	1019
Committed:	Tue Jun 4 21:16:50 2019 UTC (6 years ago) by francois
File size:	7691 byte(s)
Log Message:	restructuration de magic outil est sorti de lib pour pouvoir etre utiliser en dehors de lib template est merge avec outil poly_occ et un sous projet de magic qui utilise le nouveau outil
#	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
26			#pragma hdrstop
27			#include <stdlib.h>
28			#include <math.h>
29			#include "ot_systeme.h"
30			#include "ot_mathematique.h"
31
32			//---------------------------------------------------------------------------
33
34			OT_SYSTEME::OT_SYSTEME()
35			{
36			}
37
38
39
40			OT_SYSTEME::~OT_SYSTEME()
41			{
42			}
43
44
45			int OT_SYSTEME::resoud_gauss(double a,double b,double *x,int n)
46			{
47			double aaug=new double[(n+1)n];
48			#define a(i,j) ((a+(i)(n)+j))
49			#define aaug(i,j) ((aaug+(i)(n+1)+j))
50			for (int i=0;i<n;i++)
51			for (int j=0;j<n;j++)
52			aaug(i,j)=a(i,j);
53			for (int k=0;k<n;k++)
54			aaug(k,n)=b[k];
55			for (int j=0;j<n-1;j++)
56			{
57			// recherche le pivot
58			double pvt=fabs(aaug(j,j));
59			int ligne_pivot_debut=j;
60			int ligne_pivot=j;
61
62			for (int i=j+1;i<n;i++)
63			if (fabs(aaug(i,j))>pvt)
64			{
65			pvt=fabs(aaug(i,j));
66			ligne_pivot=i;
67			}
68			if (ligne_pivot!=ligne_pivot_debut) // trouver meilleur
69			{
70			for (int col=j;col<n+1;col++)
71			{
72			double temp=aaug(ligne_pivot_debut,col);
73			aaug(ligne_pivot_debut,col)=aaug(ligne_pivot,col);
74			aaug(ligne_pivot,col)=temp;
75			}
76			}
77			// j ai mon pivot max
78			if (fabs(aaug(j,j))< 0.00000000000001) {
79			delete [] aaug;
80			return 0;
81			};
82			// j ai mon pivot max non nul
83			for (int ii=j+1;ii<n;ii++) aaug(ii,j)=aaug(ii,j)/aaug(j,j);
84
85			for (int iii=j+1;iii<n;iii++)
86			for (int k=j+1;k<n+1;k++)
87			aaug(iii,k)=aaug(iii,k)-aaug(iii,j)*aaug(j,k);
88			}
89
90			if (fabs(aaug(n-1,n-1))<0.00000000000001) {
91			delete [] aaug;
92			return 0;
93			};
94			// back substitution
95			x[n-1]=aaug(n-1,n)/aaug(n-1,n-1);
96			for (int j=n-2;j>-1;j--)
97			{
98			x[j]=aaug(j,n);
99			for (int k=n-1;k>j;k--)
100			x[j]=x[j]-aaug(j,k)*x[k];
101			x[j]=x[j]/aaug(j,j);
102			}
103
104
105			delete [] aaug;
106			return 1;
107			#undef a // enlever le (i,j) pour enlever le warning
108			#undef aaug
109			}
110
111			int OT_SYSTEME::resoud_QR(double m, double b, double *x, int dim1,int dim2)
112			{
113			double mat=new double[dim1(dim2+1)];
114			double *r=new double[dim2];
115			#define mat(i,j) ((mat+(i)(dim2+1)+(j)))
116			#define m(i,j) ((m+(i)dim2+(j)))
117			// initialisation avec mat augment�e
118			for (int i=0;i<dim1;i++)
119			{
120			for (int j=0;j<dim2;j++)
121			mat(i,j)=m(i,j);
122			mat(i,dim2)=b[i];
123			}
124			// QR
125			for (int k=0;k<dim2;k++)
126			{
127			double somme=0.;
128			for (int i=k;i<dim1;i++)
129			somme=somme+mat(i,k)*mat(i,k);
130			double alpha=sqrt(somme);
131			double pivot=mat(k,k);
132			if (pivot>0) alpha=-alpha;
133			r[k]=alpha;
134			if (OPERATEUR::egal(alpha,0.,0.000001))
135			{
136			delete [] mat;
137			delete [] r;
138			return 0;
139			}
140			double beta=somme-pivot*alpha;
141			pivot=pivot-alpha;
142			mat(k,k)=pivot;
143			for (int j=k+1;j<dim2+1;j++)
144			{
145			double somme=0.;
146			for (int i=k;i<dim1;i++)
147			somme=somme+mat(i,k)*mat(i,j);
148			double gamma=somme/beta;
149			for (int i=k;i<dim1;i++)
150			mat(i,j)=mat(i,j)-gamma*mat(i,k);
151			}
152			}
153			// resolution
154			x[dim2-1]=mat(dim2-1,dim2)/r[dim2-1];
155			for (int j=dim2-2;j>-1;j--)
156			{
157			x[j]=mat(j,dim2);
158			for (int k=dim2-1;k>j;k--)
159			x[j]=x[j]-mat(j,k)*x[k];
160			x[j]=x[j]/r[j];
161			}
162
163
164
165			#undef mat
166			delete [] mat;
167			delete [] r;
168			return 1;
169			}
170
171
172
173
174
175
176			int OT_SYSTEME::decompose_LU(double a,double res,int dim)
177			{
178			#define a(i,j) ((a+(i)(dim)+(j)))
179			#define res(i,j) ((res+(i)(dim)+(j)))
180
181			if (OPERATEUR::egal(a(0,0),0.,0.00001)==1) return 0;
182			for (int i=0;i<dim;i++)
183			res(i,0)=a(i,0);
184			for (int j=1;j<dim;j++)
185			res(0,j)=a(0,j)/res(0,0);
186
187
188			for (int j=1;j<dim;j++)
189			{
190			for (int i=j;i<dim;i++)
191			{
192			double sum=0.;
193			for (int k=0;k<j;k++)
194			sum=sum+res(i,k)*res(k,j);
195			res(i,j)=a(i,j)-sum;
196			}
197			for (int i=j+1;i<dim;i++)
198			{
199			double sum=0.;
200			if (OPERATEUR::egal(res(j,j),0.,0.00001)==1) return 0;
201			for (int k=0;k<j;k++)
202			sum=sum+res(j,k)*res(k,i);
203			res(j,i)=(a(j,i)-sum)/res(j,j);
204			}
205
206			}
207			return 1;
208
209			#undef a
210			#undef res
211			}
212
213
214			int OT_SYSTEME::resoud_LU(double matlu,double b,double *x,int dim)
215			{
216			#define matlu(i,j) ((matlu+(i)(dim)+(j)))
217
218			double* bb=new double[dim];
219			bb[0]=b[0]/matlu(0,0);
220			for (int i=1;i<dim;i++)
221			{
222			double sum=0.0;
223			for (int k=0;k<i;k++)
224			sum=sum+matlu(i,k)*bb[k];
225			bb[i]=(b[i]-sum)/matlu(i,i);
226			}
227
228			x[dim-1]=bb[dim-1];
229			for (int j=dim-2;j>=0;j--)
230			{
231			double sum=0.0;
232			for (int k=j+1;k<dim;k++)
233			sum=sum+matlu(j,k)*x[k];
234			x[j]=bb[j]-sum;
235
236			}
237
238			#undef lu
239
240			delete [] bb;
241			return 1;
242			}
243
244
245
246
247
248
249			void OT_SYSTEME::get_inv(double mat1,double mat2,int n,int *ierr)
250			{
251			#define mat(i,j) ((mataug+(j)n+i))
252			#define mat1(i,j) ((mat1+(j)n+i))
253			#define mat2(i,j) ((mat2+(j)n+i))
254
255			double mataug=(double)calloc(2nn,sizeof(double));
256
257			for (int i=0;i<n;i++)
258			{
259			for (int j=0;j<n;j++)
260			mat(i,j)=mat1(i,j);
261			mat(i,i+n)=1.;
262			}
263
264
265			*ierr=0;
266
267			for (int j=0;j<n;j++)
268			{
269			//recherche du pivot max total de la colonne j
270			double pivot=fabs(mat(j,j));
271			int ligne_pivot=j;
272			for (int i=j+1;i<n;i++)
273			{
274			if (fabs(mat(i,j))>pivot)
275			{
276			pivot=fabs(mat(i,j));
277			ligne_pivot=i;
278			}
279			}
280			//tester si le pivot max est nul
281			if (pivot<0.000000001) {
282			*ierr=1;
283			free(mataug);
284			return;
285			}
286			//tester si le pivot actuel est max
287			if (ligne_pivot!=j)
288			{
289			for (int i=j;i<2*n;i++)
290			{
291			double temp=mat(j,i);
292			mat(j,i)=mat(ligne_pivot,i);
293			mat(ligne_pivot,i)=temp;
294			}
295			}
296			//remplir la colonne j de 0
297			for (int i=j+1;i<n;i++)
298			mat(i,j)=-mat(i,j)/mat(j,j);
299			for (int i=j+1;i<n;i++)
300			for (int k=j+1;k<2*n;k++)
301			mat(i,k)=mat(i,k)+mat(i,j)*mat(j,k);
302			}
303
304			for (int i=0;i<n;i++)
305			{
306			mat2(n-1,i)=mat(n-1,n+i)/mat(n-1,n-1);
307			for (int j=n-2;j>-1;j--)
308			{
309			mat2(j,i)=mat(j,n+i);
310			for (int k=n-1;k>j;k--)
311			mat2(j,i)=mat2(j,i)-mat2(k,i)*mat(j,k);
312			mat2(j,i)=mat2(j,i)/mat(j,j);
313			}
314			}
315
316			free(mataug);
317			#undef mat
318			#undef mat1
319			#undef mat2
320			}
321