outil/src/ot_tenseur.cpp

//---------------------------------------------------------------------------

#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(std::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(std::vector<OT_VECTEUR_4DD>&v1,std::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;
}

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

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;
    std::vector<int> check;

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

        std::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))
            {
                std::vector<int>trouve;
                std::vector<double2>valeur;
                std::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 std::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;
}


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