outil/src/ot_tenseur.cpp

//---------------------------------------------------------------------------
#include "gestionversion.h"

#pragma hdrstop

#include "ot_tenseur.h"

#include <iomanip>
//#include <conio.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
//---------------------------------------------------------------------------

#pragma package(smart_init)


OT_TENSEUR::OT_TENSEUR()
{
    Dim1=0;
    Dim2=0;
    data=NULL;
}

OT_TENSEUR::OT_TENSEUR(int n)
{
    Dim1=n;
    Dim2=n;
    double2 zero(0.);
    data =new double2[n*n];
    for (int i=0;i<Dim1*Dim2;i++) data[i]=zero;
}

OT_TENSEUR::OT_TENSEUR(int n,int m)
{
    Dim1=n;
    Dim2=m;
    double2 zero(0.);
    data =new double2[n*m];
    for (int i=0;i<Dim1*Dim2;i++) data[i]=zero;
}


OT_TENSEUR:: OT_TENSEUR(vector<OT_VECTEUR_4DD>& v)
{

    Dim1=v.size();
    Dim2=Dim1;
    double2 zero=0.;
    double2 prec=1e-6;
    prec.set_dx(1e-8);
    double2 c=0.5;
    OT_VECTEUR_4DD VCT_NUL(1e-7,1e-7,1e-7,1e-7), w;
    int Dim= Dim1*Dim2  ;
    data=new double2[Dim];
    for (int i=0;i<Dim1;i++)
    {
        OT_VECTEUR_4DD z1=v[i];
        for (int j=0;j<Dim2;j++)
        {
            OT_VECTEUR_4DD z2=v[j];
            double2 prscal=z1*z2;
            w=z1^z2;
            double2 wnorm=w.norme();
            if (wnorm<prec)
            {
                if (prscal<zero && prscal.get_fabs()>c) prscal=  -1.;
                if (prscal>zero && prscal.get_fabs()>c) prscal=   1.;
                if (prscal<zero && prscal.get_fabs()<c) prscal=   0.;
                if (prscal>zero && prscal.get_fabs()<c) prscal=   0.;
            }
            double2 val= prscal.get_fabs();
            if (val.get_x()<=val.get_dx())  {
                prscal=0.;
                prscal.set_dx(val.get_dx());
            }
            data[i*Dim2+j]=prscal;
        }
    }


}

OT_TENSEUR::OT_TENSEUR(vector<OT_VECTEUR_4DD>&v1,vector<OT_VECTEUR_4DD>&v2)
{

    Dim1=v1.size();
    Dim2=v2.size();
    int Dim= Dim1*Dim2  ;
    data=new double2[Dim];
    for (int i=0;i<Dim1;i++)
    {
        OT_VECTEUR_4DD z1=v1[i];
        for (int j=0;j<Dim2;j++)
        {
            OT_VECTEUR_4DD z2=v2[j];
            double2 prscal=z1*z2;
            double2 val= prscal.get_fabs();
            if (val.get_x()<=val.get_dx())  {
                prscal=0.;
                prscal.set_dx(val.get_dx());
            }
            data[i*Dim2+j]=prscal;
        }
    }


}


OT_TENSEUR::~OT_TENSEUR()
{
    Dim1=0;
    Dim2=0;

    if (data!=NULL) delete [] data;
}


OT_TENSEUR::OT_TENSEUR(const OT_TENSEUR& ot)
{
    Dim1=ot.Dim1;
    Dim2=ot.Dim2;
    data=new double2[Dim1*Dim2];

    for (int i=0;i<Dim1;i++)
    {
        for (int j=0;j<Dim2;j++)
        {
            data[i*Dim2+j]=ot(i,j);
        }
    }
}


OT_TENSEUR& OT_TENSEUR:: operator =(const OT_TENSEUR& ot)
{

    if (this==&ot) return *this ;

    if (data)  delete [] data;

    Dim1=ot.Dim1;
    Dim2=ot.Dim2;
    data=new double2[Dim1*Dim2];

    for (int i=0;i<Dim1;i++) {
        for (int j=0;j<Dim2;j++) {
            data[i*Dim2+j]=ot.data[i*Dim2+j];
        }
    }


    return *this;
}

int operator==( OT_TENSEUR& t1, OT_TENSEUR&t2)
{
    if (t1.Get_NbRows()!=t2.Get_NbRows()||t1.Get_NbCols()!=t2.Get_NbCols()) return 0;
    else {
        for (int i=0;i<t1.Get_NbRows();i++)
        {
            for (int j=0;j<t1.Get_NbCols();j++)
            {
                double2 val1=t1(i,j);
                double2 val2=t2(i,j);
                if (val1==val2) continue;
                else return 0;
            }

        }
    }
    return 1;
}

void OT_TENSEUR::identite(void)
{
    double2 un(1.0);
    for (int i=0;i<Dim1;i++)
        data[i*Dim2+i]=un;
}


double2& OT_TENSEUR:: operator()(int i,int j)
{
    assert(i<Dim1&&j<Dim2);
    double2 prciseion=0.0;
    if (data[i*Dim2+j].get_fabs()<=prciseion)
    {
        data[i*Dim2+j]=0.;

    }
    return data[i*Dim2+j];
}


double2& OT_TENSEUR:: operator()(int i,int j)const
{
    assert(i<Dim1&&j<Dim2);
    double2 prciseion=0.0;
    if (data[i*Dim2+j].get_fabs()<=prciseion)
    {
        data[i*Dim2+j]=0.;

    }
    return data[i*Dim2+j];
}


int OT_TENSEUR:: Get_NbRows()const
{
    return Dim1;
}


int OT_TENSEUR:: Get_NbCols()const
{
    return Dim2;
}


int OT_TENSEUR:: est_til_equivalent( OT_TENSEUR& tns)
{

    if (Dim1!=tns.Get_NbRows()||Dim2!=tns.Get_NbCols()) return 0;
    vector<int> check;

    for (int i=0;i<tns.Get_NbRows();i++)
    {

        vector<double2> lig1;
        for (int ix1=0;ix1<Dim2;ix1++) lig1.insert(lig1.end(),data[i*Dim2+ix1]);

        for (int j=0;j<Dim2;j++)
        {
            int p=-1;
            if (!existe(check,j,p))
            {
                vector<int>trouve;
                vector<double2>valeur;
                vector<double2> lig2;
                for (int ix2=0;ix2<Dim2;ix2++) lig2.insert(lig2.end(),tns(j,ix2));
                for (int r=0;r<Dim2;r++)
                {
                    double2 val=lig1[r];
                    val=val.get_fabs();

                    int indx1=-1;

                    for (int s=0;s<Dim2;s++)
                    {
                        int pos1;
                        double2 vall=lig2[s];
                        vall=vall.get_fabs();

                        if (val==vall&&!existe(trouve,s,indx1))
                        {
                            trouve.insert(trouve.end(),s);
                            valeur.insert(valeur.end(),val);
                            break;
                        }

                    }

                }
                int t_trouve=trouve.size();
                if (t_trouve==Dim2)
                {
                    check.insert(check.end(),j);
                    break;

                }

            }
        }
        continue;
    }
    int t_chek=check.size();
    if (t_chek==Dim1)return 1;
    else return 0;

}


template <class T> int OT_TENSEUR:: existe(const vector<T>& vct,T val,int& idx)
{

    if (val<0)   val=-val;
    for (unsigned int i=0;i<vct.size();i++)
    {
        double x=vct[i];
        if (x<0)   val=-val;
        if (x==val) {
            idx=i;
            return 1;
        }
    }
    idx=-1;
    return 0;
}


OT_TENSEUR operator+( OT_TENSEUR& ot1, OT_TENSEUR& ot2)
{
    assert (ot1.Get_NbRows()==ot2.Get_NbRows()&&ot1.Get_NbCols()==ot2.Get_NbCols());
    OT_TENSEUR res(ot1.Get_NbRows(),ot1.Get_NbCols());
    for (int i=0;i<ot1.Get_NbRows();i++)
    {
        for (int j=0;j<ot2.Get_NbCols();j++)
        {
            res(i,j)=ot1(i,j)+ot2(i,j);
        }
    }

    return res;
}


vector<double2> operator*( OT_TENSEUR& ot1, vector<double2>& v)
{
    assert (ot1.Get_NbCols()==v.size());
    vector<double2> res(ot1.Get_NbCols());
    for (int i=0;i<ot1.Get_NbRows();i++)
    {
        res[i]=double2(0.0);
        for (int j=0;j<ot1.Get_NbCols();j++)
        {
            res[i]=res[i]+ot1(i,j)+v[j];
        }
    }

    return res;


}

OT_TENSEUR operator*( double2& coef, OT_TENSEUR& ot)
{
    OT_TENSEUR T=ot;
    for (int i=0;i<ot.Get_NbRows();i++)
        for (int j=0;j<ot.Get_NbCols();j++)

            T(i,j)=coef*T(i,j);

    return T;

}

OT_TENSEUR operator*(OT_TENSEUR& ot , double2& coef)
{
    return coef*ot;
}

OT_TENSEUR OT_TENSEUR::inverse_homogene(void)
{
    OT_TENSEUR TMP(4,4);
    if (Dim1!=4) return TMP;
    if (Dim2!=4) return TMP;
    double2 zero(0.);
    double2 un(1.);
    OT_VECTEUR_4DD x1((*this)(0,0),(*this)(1,0),(*this)(2,0),zero);
    OT_VECTEUR_4DD x2((*this)(0,1),(*this)(1,1),(*this)(2,1),zero);
    OT_VECTEUR_4DD x3((*this)(0,2),(*this)(1,2),(*this)(2,2),zero);
    OT_VECTEUR_4DD x4((*this)(0,3)/(*this)(3,3),(*this)(1,3)/(*this)(3,3),(*this)(2,3)/(*this)(3,3),un);
    TMP(0,0)=(*this)(0,0);
    TMP(1,0)=(*this)(0,1);
    TMP(2,0)=(*this)(0,2);
    TMP(3,0)=zero;
    TMP(0,1)=(*this)(1,0);
    TMP(1,1)=(*this)(1,1);
    TMP(2,1)=(*this)(1,2);
    TMP(3,1)=zero;
    TMP(0,2)=(*this)(2,0);
    TMP(1,2)=(*this)(2,1);
    TMP(2,2)=(*this)(2,2);
    TMP(3,2)=zero;
    TMP(0,3)=zero-x4*x1;
    TMP(1,3)=zero-x4*x2;
    TMP(2,3)=zero-x4*x3;
    TMP(3,3)=un;
    return TMP;
}


OT_TENSEUR OT_TENSEUR:: transpose()
{
    OT_TENSEUR TMP(Dim2,Dim1);
    for (int i=0;i<Dim1;i++) {
        for (int j=0;j<Dim2;j++)
        {
            TMP(j,i)=(*this)(i,j);  // ot(i,j);
        }
    }
    return TMP;
}


OT_TENSEUR operator*( OT_TENSEUR& ot1, OT_TENSEUR& ot2)
{
    int dim1=ot1.Get_NbRows();
    int dim3= ot1.Get_NbCols();
    int dim2=ot2.Get_NbCols();
    OT_TENSEUR TMP (dim1,dim2);

    for (int i=0;i<dim1;i++)
    {
        for (int j=0;j<dim2;j++)
        {
            for (int k=0;k<dim3;k++)
                TMP(i,j)=TMP(i,j)+ot1(i,k)*ot2(k,j);
        }
    }
    return TMP;
}

OT_VECTEUR_4DD operator*( OT_TENSEUR& ot1, OT_VECTEUR_4DD& ot2)
{
    OT_VECTEUR_4DD TMP(0.,0.,0.,0.);
    if (ot1.Get_NbRows()!=4) return TMP;
    if (ot1.Get_NbCols()!=4) return TMP;
    for (int i=0;i<4;i++)
    {
        TMP[i]=TMP[i]+ot1(i,0)*ot2[0]+ot1(i,1)*ot2[1]+ot1(i,2)*ot2[2]+ot1(i,3)*ot2[3];
    }
    return TMP;
}

OT_VECTEUR_3DD operator*( OT_TENSEUR& ot1, OT_VECTEUR_3DD& ot2)
{
    OT_VECTEUR_3DD TMP(0.,0.,0.);
    if (ot1.Get_NbRows()!=3) return TMP;
    if (ot1.Get_NbCols()!=3) return TMP;
    for (int i=0;i<3;i++)
    {
        TMP[i]=TMP[i]+ot1(i,0)*ot2[0]+ot1(i,1)*ot2[1]+ot1(i,2)*ot2[2];
    }
    return TMP;
}

std::ostream& operator<<(std::ostream& os, OT_TENSEUR& tns) {
    os<<"------------------------------------------------------------------------------------"<<endl;
    os<<"#  TENSEUR:   DIM1=  "<< tns.Get_NbRows()<<"    , DIM2:  "<<tns.Get_NbCols()<<"  #"<<endl;
    os<<"------------------------------------------------------------------------------------"<<endl;
    for (int i=0;i<tns.Get_NbRows();i++) {
        for (int j=0;j<tns.Get_NbCols();j++) {
            os<<setw(18)<<tns(i,j);
        }
        os<<endl;
    }
    os<<"------------------------------------------------------------------------------------"<<endl;
    return os;
}


void OT_TENSEUR::get_orthogonalisation(OT_TENSEUR& a,OT_VECTEUR_4DD& d,OT_TENSEUR& v,int n,int&nrot)
{

    static const int ITMAX=50;
    double2 zro=0.0;
    int n1=n-1,n2=1/(n*n),ip1,iq1;// ,ip,j,i,iq// n=3,
    double2 tresh,theta,tau,t,sm,s,h,g,c;
    double2 *b=new double2[n];                 // OT_VECTEUR_3D b,z;
    double2 *z=new double2[n];

    for (int ip=0;ip<n;ip++)
    {
        for (int iq=0;iq<n;iq++)
            //v[ip*n+iq]=0.0;
            v(ip,iq)=0.0;
        v(ip,ip)=1.0;
        b[(unsigned int)ip]=d[(unsigned int)ip]=a(ip,ip);
        z[(unsigned int)ip]=0.0;
    }
    nrot=0;
    for (int i=1;i<=ITMAX;i++)
    {
        sm=0.0;
        for (int ip=0;ip<n1;ip++)
            for (int iq=ip+1;iq<n;iq++)
                sm =sm+(a(ip,iq)).get_fabs();
        if (sm==zro)
            return;
        if (i<4)
            tresh=0.2*sm*n2;
        else
            tresh=0.0;
        for (int ip=0;ip<n1;ip++)
            for (int iq=ip+1;iq<n;iq++)
            {
                g=100*(a(ip,iq)).get_fabs();
                if (i>4&&((d[(unsigned int)ip]+g).get_fabs()==d[(unsigned int)ip].get_fabs())&&((d[(unsigned int)iq])+g).get_fabs()==d[(unsigned int)iq].get_fabs())
                    a(ip,iq)=0.0;
                else if ((a(ip,iq)).get_fabs()>tresh)
                {
                    h=d[(unsigned int)iq]-d[(unsigned int)ip];
                    if (((h).get_fabs()+g)==(h).get_fabs())
                        t=a(ip,iq)/h;
                    else
                    {
                        theta=0.5*h/a(ip,iq);
                        t=1.0/((theta).get_fabs()+sqrt(1.0+theta*theta));
                        if (theta<zro)t=zro-t;
                    }
                    c=1.0/sqrt(1+t*t);
                    s=t*c;
                    tau=s/(1.0+c);
                    h=t*a(ip,iq);
                    z[(unsigned int)ip]=z[(unsigned int)ip]-h;
                    z[(unsigned int)iq]=z[(unsigned int)iq]+h;
                    d[(unsigned int)ip]=d[(unsigned int)ip]-h;
                    d[(unsigned int)iq]=d[(unsigned int)iq]+h;
                    a(ip,iq)=0.0;
                    ip1=ip-1;
                    iq1=iq-1;
                    for (int j=0;j<=ip1;j++)
                    {
                        g=a(j,ip);
                        h=a(j,iq);
                        a(j,ip)=g-s*(h+g*tau);
                        a(j,iq)=h+s*(g-h*tau);
                    }
                    for (int j=ip+1;j<=iq1;j++)
                    {
                        g=a(ip,j);
                        h=a(j,iq);
                        a(ip,j)=g-s*(h+g*tau);
                        a(j,iq)=h+s*(g-h*tau);
                    }
                    for (int j=iq+1;j<n;j++)
                    {
                        g=a(ip,j);
                        h=a(iq,j);
                        a(ip,j)=g-s*(h+g*tau);
                        a(iq,j)=h+s*(g-h*tau);
                    }
                    for (int j=0;j<n;j++)
                    {
                        g=v(j,ip);
                        h=v(j,iq);
                        v(j,ip)=g-s*(h+g*tau);
                        v(j,iq)=h+s*(g-h*tau);
                    }
                    ++nrot;
                }
            }
        for (int ip=0;ip<n;ip++)
        {
            b[(unsigned int)ip]=b[(unsigned int)ip]+z[(unsigned int)ip];
            d[(unsigned int)ip]=b[(unsigned int)ip];
            z[(unsigned int)ip]=0.0;
        }
    }


    delete [] b;
    delete [] z;

}


int  OT_TENSEUR::listes_equivalentes(std::vector<double2>& list1,std::vector<double2>& list2,std::vector<unsigned int>& indx)

{
    int compt=0;
    int pos=-1;
    int exist=0;
    for (int i=0;i<list1.size();i++)
    {
        double2 val1=list1[i];
        double prc=1e-3;
        double prcc;
        if (pos==-1)
        {
            for (int k=0;k< list2.size();k++)
            {
                double2 val2=list2[k];
                prcc =val2.get_x()-val1.get_x();
                if (fabs(val2.get_x())>val2.get_dx())
                    prcc =prcc/val2.get_x();
                if ( fabs(prcc)<prc) {
                    indx.insert(indx.end(),k);
                    pos=k;
                    compt++;
                    break;
                }
            }
        }


        else  {
            for ( int r=0;r<list2.size();r++)
            {
                double2 vval2=list2[r];
                for (int j=0;j<indx.size();j++)
                {
                    exist=0;
                    if (r==indx[j])
                    {
                        exist=1;
                        break;
                    }
                }

                prcc =vval2.get_x()-val1.get_x();
                if (fabs(vval2.get_x())>vval2.get_dx())
                    prcc =prcc/vval2.get_x();

                if (fabs(prcc)<prc&&!exist)  // vval2==val1
                {
                    indx.insert(indx.end(),r);
                    pos=r;
                    compt++;
                    break;
                }

            }
        }
    }

    if (compt==list1.size())return 1;
    else return 0;


}


Revision:	283
Committed:	Tue Sep 13 21:11:20 2011 UTC (13 years, 8 months ago) by francois
File size:	15193 byte(s)
Log Message:	structure de l'écriture
#	User	Rev	Content
1	francois	283	//---------------------------------------------------------------------------
2			#include "gestionversion.h"
3
4			#pragma hdrstop
5
6			#include "ot_tenseur.h"
7
8			#include <iomanip>
9			//#include <conio.h>
10			#include <stdio.h>
11			#include <assert.h>
12			#include <math.h>
13			//---------------------------------------------------------------------------
14
15			#pragma package(smart_init)
16
17
18
19
20			OT_TENSEUR::OT_TENSEUR()
21			{
22			Dim1=0;
23			Dim2=0;
24			data=NULL;
25			}
26
27			OT_TENSEUR::OT_TENSEUR(int n)
28			{
29			Dim1=n;
30			Dim2=n;
31			double2 zero(0.);
32			data =new double2[n*n];
33			for (int i=0;i<Dim1*Dim2;i++) data[i]=zero;
34			}
35
36			OT_TENSEUR::OT_TENSEUR(int n,int m)
37			{
38			Dim1=n;
39			Dim2=m;
40			double2 zero(0.);
41			data =new double2[n*m];
42			for (int i=0;i<Dim1*Dim2;i++) data[i]=zero;
43			}
44
45
46
47			OT_TENSEUR:: OT_TENSEUR(vector<OT_VECTEUR_4DD>& v)
48			{
49
50			Dim1=v.size();
51			Dim2=Dim1;
52			double2 zero=0.;
53			double2 prec=1e-6;
54			prec.set_dx(1e-8);
55			double2 c=0.5;
56			OT_VECTEUR_4DD VCT_NUL(1e-7,1e-7,1e-7,1e-7), w;
57			int Dim= Dim1*Dim2 ;
58			data=new double2[Dim];
59			for (int i=0;i<Dim1;i++)
60			{
61			OT_VECTEUR_4DD z1=v[i];
62			for (int j=0;j<Dim2;j++)
63			{
64			OT_VECTEUR_4DD z2=v[j];
65			double2 prscal=z1*z2;
66			w=z1^z2;
67			double2 wnorm=w.norme();
68			if (wnorm<prec)
69			{
70			if (prscal<zero && prscal.get_fabs()>c) prscal= -1.;
71			if (prscal>zero && prscal.get_fabs()>c) prscal= 1.;
72			if (prscal<zero && prscal.get_fabs()<c) prscal= 0.;
73			if (prscal>zero && prscal.get_fabs()<c) prscal= 0.;
74			}
75			double2 val= prscal.get_fabs();
76			if (val.get_x()<=val.get_dx()) {
77			prscal=0.;
78			prscal.set_dx(val.get_dx());
79			}
80			data[i*Dim2+j]=prscal;
81			}
82			}
83
84
85			}
86
87			OT_TENSEUR::OT_TENSEUR(vector<OT_VECTEUR_4DD>&v1,vector<OT_VECTEUR_4DD>&v2)
88			{
89
90			Dim1=v1.size();
91			Dim2=v2.size();
92			int Dim= Dim1*Dim2 ;
93			data=new double2[Dim];
94			for (int i=0;i<Dim1;i++)
95			{
96			OT_VECTEUR_4DD z1=v1[i];
97			for (int j=0;j<Dim2;j++)
98			{
99			OT_VECTEUR_4DD z2=v2[j];
100			double2 prscal=z1*z2;
101			double2 val= prscal.get_fabs();
102			if (val.get_x()<=val.get_dx()) {
103			prscal=0.;
104			prscal.set_dx(val.get_dx());
105			}
106			data[i*Dim2+j]=prscal;
107			}
108			}
109
110
111			}
112
113
114			OT_TENSEUR::~OT_TENSEUR()
115			{
116			Dim1=0;
117			Dim2=0;
118
119			if (data!=NULL) delete [] data;
120			}
121
122
123			OT_TENSEUR::OT_TENSEUR(const OT_TENSEUR& ot)
124			{
125			Dim1=ot.Dim1;
126			Dim2=ot.Dim2;
127			data=new double2[Dim1*Dim2];
128
129			for (int i=0;i<Dim1;i++)
130			{
131			for (int j=0;j<Dim2;j++)
132			{
133			data[i*Dim2+j]=ot(i,j);
134			}
135			}
136			}
137
138
139			OT_TENSEUR& OT_TENSEUR:: operator =(const OT_TENSEUR& ot)
140			{
141
142			if (this==&ot) return *this ;
143
144			if (data) delete [] data;
145
146			Dim1=ot.Dim1;
147			Dim2=ot.Dim2;
148			data=new double2[Dim1*Dim2];
149
150			for (int i=0;i<Dim1;i++) {
151			for (int j=0;j<Dim2;j++) {
152			data[iDim2+j]=ot.data[iDim2+j];
153			}
154			}
155
156
157			return *this;
158			}
159
160			int operator==( OT_TENSEUR& t1, OT_TENSEUR&t2)
161			{
162			if (t1.Get_NbRows()!=t2.Get_NbRows()\|\|t1.Get_NbCols()!=t2.Get_NbCols()) return 0;
163			else {
164			for (int i=0;i<t1.Get_NbRows();i++)
165			{
166			for (int j=0;j<t1.Get_NbCols();j++)
167			{
168			double2 val1=t1(i,j);
169			double2 val2=t2(i,j);
170			if (val1==val2) continue;
171			else return 0;
172			}
173
174			}
175			}
176			return 1;
177			}
178
179			void OT_TENSEUR::identite(void)
180			{
181			double2 un(1.0);
182			for (int i=0;i<Dim1;i++)
183			data[i*Dim2+i]=un;
184			}
185
186
187
188			double2& OT_TENSEUR:: operator()(int i,int j)
189			{
190			assert(i<Dim1&&j<Dim2);
191			double2 prciseion=0.0;
192			if (data[i*Dim2+j].get_fabs()<=prciseion)
193			{
194			data[i*Dim2+j]=0.;
195
196			}
197			return data[i*Dim2+j];
198			}
199
200
201
202			double2& OT_TENSEUR:: operator()(int i,int j)const
203			{
204			assert(i<Dim1&&j<Dim2);
205			double2 prciseion=0.0;
206			if (data[i*Dim2+j].get_fabs()<=prciseion)
207			{
208			data[i*Dim2+j]=0.;
209
210			}
211			return data[i*Dim2+j];
212			}
213
214
215
216			int OT_TENSEUR:: Get_NbRows()const
217			{
218			return Dim1;
219			}
220
221
222
223			int OT_TENSEUR:: Get_NbCols()const
224			{
225			return Dim2;
226			}
227
228
229
230			int OT_TENSEUR:: est_til_equivalent( OT_TENSEUR& tns)
231			{
232
233			if (Dim1!=tns.Get_NbRows()\|\|Dim2!=tns.Get_NbCols()) return 0;
234			vector<int> check;
235
236			for (int i=0;i<tns.Get_NbRows();i++)
237			{
238
239			vector<double2> lig1;
240			for (int ix1=0;ix1<Dim2;ix1++) lig1.insert(lig1.end(),data[i*Dim2+ix1]);
241
242			for (int j=0;j<Dim2;j++)
243			{
244			int p=-1;
245			if (!existe(check,j,p))
246			{
247			vector<int>trouve;
248			vector<double2>valeur;
249			vector<double2> lig2;
250			for (int ix2=0;ix2<Dim2;ix2++) lig2.insert(lig2.end(),tns(j,ix2));
251			for (int r=0;r<Dim2;r++)
252			{
253			double2 val=lig1[r];
254			val=val.get_fabs();
255
256			int indx1=-1;
257
258			for (int s=0;s<Dim2;s++)
259			{
260			int pos1;
261			double2 vall=lig2[s];
262			vall=vall.get_fabs();
263
264			if (val==vall&&!existe(trouve,s,indx1))
265			{
266			trouve.insert(trouve.end(),s);
267			valeur.insert(valeur.end(),val);
268			break;
269			}
270
271			}
272
273			}
274			int t_trouve=trouve.size();
275			if (t_trouve==Dim2)
276			{
277			check.insert(check.end(),j);
278			break;
279
280			}
281
282			}
283			}
284			continue;
285			}
286			int t_chek=check.size();
287			if (t_chek==Dim1)return 1;
288			else return 0;
289
290			}
291
292
293
294			template <class T> int OT_TENSEUR:: existe(const vector<T>& vct,T val,int& idx)
295			{
296
297			if (val<0) val=-val;
298			for (unsigned int i=0;i<vct.size();i++)
299			{
300			double x=vct[i];
301			if (x<0) val=-val;
302			if (x==val) {
303			idx=i;
304			return 1;
305			}
306			}
307			idx=-1;
308			return 0;
309			}
310
311
312			OT_TENSEUR operator+( OT_TENSEUR& ot1, OT_TENSEUR& ot2)
313			{
314			assert (ot1.Get_NbRows()==ot2.Get_NbRows()&&ot1.Get_NbCols()==ot2.Get_NbCols());
315			OT_TENSEUR res(ot1.Get_NbRows(),ot1.Get_NbCols());
316			for (int i=0;i<ot1.Get_NbRows();i++)
317			{
318			for (int j=0;j<ot2.Get_NbCols();j++)
319			{
320			res(i,j)=ot1(i,j)+ot2(i,j);
321			}
322			}
323
324			return res;
325			}
326
327
328			vector<double2> operator*( OT_TENSEUR& ot1, vector<double2>& v)
329			{
330			assert (ot1.Get_NbCols()==v.size());
331			vector<double2> res(ot1.Get_NbCols());
332			for (int i=0;i<ot1.Get_NbRows();i++)
333			{
334			res[i]=double2(0.0);
335			for (int j=0;j<ot1.Get_NbCols();j++)
336			{
337			res[i]=res[i]+ot1(i,j)+v[j];
338			}
339			}
340
341			return res;
342
343
344			}
345
346			OT_TENSEUR operator*( double2& coef, OT_TENSEUR& ot)
347			{
348			OT_TENSEUR T=ot;
349			for (int i=0;i<ot.Get_NbRows();i++)
350			for (int j=0;j<ot.Get_NbCols();j++)
351
352			T(i,j)=coef*T(i,j);
353
354			return T;
355
356			}
357
358			OT_TENSEUR operator*(OT_TENSEUR& ot , double2& coef)
359			{
360			return coef*ot;
361			}
362
363			OT_TENSEUR OT_TENSEUR::inverse_homogene(void)
364			{
365			OT_TENSEUR TMP(4,4);
366			if (Dim1!=4) return TMP;
367			if (Dim2!=4) return TMP;
368			double2 zero(0.);
369			double2 un(1.);
370			OT_VECTEUR_4DD x1((this)(0,0),(this)(1,0),(*this)(2,0),zero);
371			OT_VECTEUR_4DD x2((this)(0,1),(this)(1,1),(*this)(2,1),zero);
372			OT_VECTEUR_4DD x3((this)(0,2),(this)(1,2),(*this)(2,2),zero);
373			OT_VECTEUR_4DD x4((this)(0,3)/(this)(3,3),(this)(1,3)/(this)(3,3),(this)(2,3)/(this)(3,3),un);
374			TMP(0,0)=(*this)(0,0);
375			TMP(1,0)=(*this)(0,1);
376			TMP(2,0)=(*this)(0,2);
377			TMP(3,0)=zero;
378			TMP(0,1)=(*this)(1,0);
379			TMP(1,1)=(*this)(1,1);
380			TMP(2,1)=(*this)(1,2);
381			TMP(3,1)=zero;
382			TMP(0,2)=(*this)(2,0);
383			TMP(1,2)=(*this)(2,1);
384			TMP(2,2)=(*this)(2,2);
385			TMP(3,2)=zero;
386			TMP(0,3)=zero-x4*x1;
387			TMP(1,3)=zero-x4*x2;
388			TMP(2,3)=zero-x4*x3;
389			TMP(3,3)=un;
390			return TMP;
391			}
392
393
394
395			OT_TENSEUR OT_TENSEUR:: transpose()
396			{
397			OT_TENSEUR TMP(Dim2,Dim1);
398			for (int i=0;i<Dim1;i++) {
399			for (int j=0;j<Dim2;j++)
400			{
401			TMP(j,i)=(*this)(i,j); // ot(i,j);
402			}
403			}
404			return TMP;
405			}
406
407
408			OT_TENSEUR operator*( OT_TENSEUR& ot1, OT_TENSEUR& ot2)
409			{
410			int dim1=ot1.Get_NbRows();
411			int dim3= ot1.Get_NbCols();
412			int dim2=ot2.Get_NbCols();
413			OT_TENSEUR TMP (dim1,dim2);
414
415			for (int i=0;i<dim1;i++)
416			{
417			for (int j=0;j<dim2;j++)
418			{
419			for (int k=0;k<dim3;k++)
420			TMP(i,j)=TMP(i,j)+ot1(i,k)*ot2(k,j);
421			}
422			}
423			return TMP;
424			}
425
426			OT_VECTEUR_4DD operator*( OT_TENSEUR& ot1, OT_VECTEUR_4DD& ot2)
427			{
428			OT_VECTEUR_4DD TMP(0.,0.,0.,0.);
429			if (ot1.Get_NbRows()!=4) return TMP;
430			if (ot1.Get_NbCols()!=4) return TMP;
431			for (int i=0;i<4;i++)
432			{
433			TMP[i]=TMP[i]+ot1(i,0)ot2[0]+ot1(i,1)ot2[1]+ot1(i,2)ot2[2]+ot1(i,3)ot2[3];
434			}
435			return TMP;
436			}
437
438			OT_VECTEUR_3DD operator*( OT_TENSEUR& ot1, OT_VECTEUR_3DD& ot2)
439			{
440			OT_VECTEUR_3DD TMP(0.,0.,0.);
441			if (ot1.Get_NbRows()!=3) return TMP;
442			if (ot1.Get_NbCols()!=3) return TMP;
443			for (int i=0;i<3;i++)
444			{
445			TMP[i]=TMP[i]+ot1(i,0)ot2[0]+ot1(i,1)ot2[1]+ot1(i,2)*ot2[2];
446			}
447			return TMP;
448			}
449
450			std::ostream& operator<<(std::ostream& os, OT_TENSEUR& tns) {
451			os<<"------------------------------------------------------------------------------------"<<endl;
452			os<<"# TENSEUR: DIM1= "<< tns.Get_NbRows()<<" , DIM2: "<<tns.Get_NbCols()<<" #"<<endl;
453			os<<"------------------------------------------------------------------------------------"<<endl;
454			for (int i=0;i<tns.Get_NbRows();i++) {
455			for (int j=0;j<tns.Get_NbCols();j++) {
456			os<<setw(18)<<tns(i,j);
457			}
458			os<<endl;
459			}
460			os<<"------------------------------------------------------------------------------------"<<endl;
461			return os;
462			}
463
464
465
466			void OT_TENSEUR::get_orthogonalisation(OT_TENSEUR& a,OT_VECTEUR_4DD& d,OT_TENSEUR& v,int n,int&nrot)
467			{
468
469			static const int ITMAX=50;
470			double2 zro=0.0;
471			int n1=n-1,n2=1/(n*n),ip1,iq1;// ,ip,j,i,iq// n=3,
472			double2 tresh,theta,tau,t,sm,s,h,g,c;
473			double2 *b=new double2[n]; // OT_VECTEUR_3D b,z;
474			double2 *z=new double2[n];
475
476			for (int ip=0;ip<n;ip++)
477			{
478			for (int iq=0;iq<n;iq++)
479			//v[ip*n+iq]=0.0;
480			v(ip,iq)=0.0;
481			v(ip,ip)=1.0;
482			b[(unsigned int)ip]=d[(unsigned int)ip]=a(ip,ip);
483			z[(unsigned int)ip]=0.0;
484			}
485			nrot=0;
486			for (int i=1;i<=ITMAX;i++)
487			{
488			sm=0.0;
489			for (int ip=0;ip<n1;ip++)
490			for (int iq=ip+1;iq<n;iq++)
491			sm =sm+(a(ip,iq)).get_fabs();
492			if (sm==zro)
493			return;
494			if (i<4)
495			tresh=0.2smn2;
496			else
497			tresh=0.0;
498			for (int ip=0;ip<n1;ip++)
499			for (int iq=ip+1;iq<n;iq++)
500			{
501			g=100*(a(ip,iq)).get_fabs();
502			if (i>4&&((d[(unsigned int)ip]+g).get_fabs()==d[(unsigned int)ip].get_fabs())&&((d[(unsigned int)iq])+g).get_fabs()==d[(unsigned int)iq].get_fabs())
503			a(ip,iq)=0.0;
504			else if ((a(ip,iq)).get_fabs()>tresh)
505			{
506			h=d[(unsigned int)iq]-d[(unsigned int)ip];
507			if (((h).get_fabs()+g)==(h).get_fabs())
508			t=a(ip,iq)/h;
509			else
510			{
511			theta=0.5*h/a(ip,iq);
512			t=1.0/((theta).get_fabs()+sqrt(1.0+theta*theta));
513			if (theta<zro)t=zro-t;
514			}
515			c=1.0/sqrt(1+t*t);
516			s=t*c;
517			tau=s/(1.0+c);
518			h=t*a(ip,iq);
519			z[(unsigned int)ip]=z[(unsigned int)ip]-h;
520			z[(unsigned int)iq]=z[(unsigned int)iq]+h;
521			d[(unsigned int)ip]=d[(unsigned int)ip]-h;
522			d[(unsigned int)iq]=d[(unsigned int)iq]+h;
523			a(ip,iq)=0.0;
524			ip1=ip-1;
525			iq1=iq-1;
526			for (int j=0;j<=ip1;j++)
527			{
528			g=a(j,ip);
529			h=a(j,iq);
530			a(j,ip)=g-s(h+gtau);
531			a(j,iq)=h+s(g-htau);
532			}
533			for (int j=ip+1;j<=iq1;j++)
534			{
535			g=a(ip,j);
536			h=a(j,iq);
537			a(ip,j)=g-s(h+gtau);
538			a(j,iq)=h+s(g-htau);
539			}
540			for (int j=iq+1;j<n;j++)
541			{
542			g=a(ip,j);
543			h=a(iq,j);
544			a(ip,j)=g-s(h+gtau);
545			a(iq,j)=h+s(g-htau);
546			}
547			for (int j=0;j<n;j++)
548			{
549			g=v(j,ip);
550			h=v(j,iq);
551			v(j,ip)=g-s(h+gtau);
552			v(j,iq)=h+s(g-htau);
553			}
554			++nrot;
555			}
556			}
557			for (int ip=0;ip<n;ip++)
558			{
559			b[(unsigned int)ip]=b[(unsigned int)ip]+z[(unsigned int)ip];
560			d[(unsigned int)ip]=b[(unsigned int)ip];
561			z[(unsigned int)ip]=0.0;
562			}
563			}
564
565
566			delete [] b;
567			delete [] z;
568
569			}
570
571
572
573
574			int OT_TENSEUR::listes_equivalentes(std::vector<double2>& list1,std::vector<double2>& list2,std::vector<unsigned int>& indx)
575
576			{
577			int compt=0;
578			int pos=-1;
579			int exist=0;
580			for (int i=0;i<list1.size();i++)
581			{
582			double2 val1=list1[i];
583			double prc=1e-3;
584			double prcc;
585			if (pos==-1)
586			{
587			for (int k=0;k< list2.size();k++)
588			{
589			double2 val2=list2[k];
590			prcc =val2.get_x()-val1.get_x();
591			if (fabs(val2.get_x())>val2.get_dx())
592			prcc =prcc/val2.get_x();
593			if ( fabs(prcc)<prc) {
594			indx.insert(indx.end(),k);
595			pos=k;
596			compt++;
597			break;
598			}
599			}
600			}
601
602
603			else {
604			for ( int r=0;r<list2.size();r++)
605			{
606			double2 vval2=list2[r];
607			for (int j=0;j<indx.size();j++)
608			{
609			exist=0;
610			if (r==indx[j])
611			{
612			exist=1;
613			break;
614			}
615			}
616
617			prcc =vval2.get_x()-val1.get_x();
618			if (fabs(vval2.get_x())>vval2.get_dx())
619			prcc =prcc/vval2.get_x();
620
621			if (fabs(prcc)<prc&&!exist) // vval2==val1
622			{
623			indx.insert(indx.end(),r);
624			pos=r;
625			compt++;
626			break;
627			}
628
629			}
630			}
631			}
632
633			if (compt==list1.size())return 1;
634			else return 0;
635
636
637			}
638
639
640
641
642
643