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

#pragma package(smart_init)


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

OT_TENSEUR::OT_TENSEUR(int n)
{
 Dim1=n;
 Dim2=n;

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

OT_TENSEUR::OT_TENSEUR(int n,int m)
{
 Dim1=n;
 Dim2=m;

 data =new double2[n*m];
      for(int i=0;i<Dim1*Dim2;i++) data[i]=0;
}


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

 Dim1=v.size();
 Dim2=Dim1;
 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;
       double2 prciseion=0.0;
       if(prscal.get_fabs()<=prciseion)  prscal=0.;
       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 prciseion=0.0;
       if(prscal.get_fabs()<=prciseion)  prscal=0.;
       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;
}


double2&OT_TENSEUR:: operator()(int i,int j)
{
assert(i<Dim1&&j<Dim2);

return data[i*Dim2+j];
}


double2& OT_TENSEUR:: operator()(int i,int j)const
{
assert(i<Dim1&&j<Dim2);

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)
 {
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();
            //val.set_dx(0.01);
            int indx1=-1;

            for(int s=0;s<Dim2;s++)
             {  int pos1;
                double2 vall=lig2[s];
                vall=vall.get_fabs();
                //vall.set_dx(0.01);
                if (val==vall&&!existe(trouve,s,indx1))//&&!existe(valeur,vall,pos1))
                  {
                  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;
 }


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:: 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;
}

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;

}


Revision:	75
Committed:	Wed Apr 2 12:52:33 2008 UTC (17 years, 1 month ago) by souaissa
Original Path:	*magic/lib/outil/outil/src/ot_tenseur.cpp*
File size:	9563 byte(s)
Log Message:	Mise à jour de ot_tenseur et ot_mathematique et ot_doubleprecision
#	User	Rev	Content
1	souaissa	58	//---------------------------------------------------------------------------
2	francois	17	#include "gestionversion.h"
3
4			#pragma hdrstop
5	souaissa	58
6	francois	17	#include "ot_tenseur.h"
7	souaissa	58
8			#include <iomanip>
9			#include <conio.h>
10			#include <stdio.h>
11	souaissa	69	#include <assert.h>
12	francois	17	//---------------------------------------------------------------------------
13	souaissa	58
14	francois	17	#pragma package(smart_init)
15
16
17
18	souaissa	69
19			OT_TENSEUR::OT_TENSEUR()
20	souaissa	58	{
21	souaissa	69	Dim1=0;
22			Dim2=0;
23	souaissa	70	data=0;
24	souaissa	58	}
25	souaissa	70
26	souaissa	69	OT_TENSEUR::OT_TENSEUR(int n)
27	souaissa	58	{
28	souaissa	69	Dim1=n;
29			Dim2=n;
30	souaissa	70
31			data =new double2[n*n];
32			for(int i=0;i<Dim1*Dim2;i++) data[i]=0;
33			}
34
35	souaissa	69	OT_TENSEUR::OT_TENSEUR(int n,int m)
36	francois	17	{
37	souaissa	69	Dim1=n;
38			Dim2=m;
39	souaissa	70
40			data =new double2[n*m];
41			for(int i=0;i<Dim1*Dim2;i++) data[i]=0;
42	souaissa	58	}
43	souaissa	70
44
45
46	souaissa	69	OT_TENSEUR:: OT_TENSEUR(vector<OT_VECTEUR_4DD>& v)
47	francois	17	{
48	souaissa	70
49	souaissa	69	Dim1=v.size();
50			Dim2=Dim1;
51			int Dim= Dim1*Dim2 ;
52	souaissa	70	data=new double2[Dim];
53			for(int i=0;i<Dim1;i++)
54			{
55			OT_VECTEUR_4DD z1=v[i];
56			for (int j=0;j<Dim2;j++)
57			{
58			OT_VECTEUR_4DD z2=v[j];
59			double2 prscal=z1*z2;
60	souaissa	75	double2 prciseion=0.0;
61			if(prscal.get_fabs()<=prciseion) prscal=0.;
62	souaissa	70	data[i*Dim2+j]=prscal;
63			}
64	souaissa	69	}
65	francois	17
66
67			}
68
69	souaissa	69	OT_TENSEUR::OT_TENSEUR(vector<OT_VECTEUR_4DD>&v1,vector<OT_VECTEUR_4DD>&v2)
70	souaissa	58	{
71	francois	17
72	souaissa	69	Dim1=v1.size();
73			Dim2=v2.size();
74			int Dim= Dim1*Dim2 ;
75	souaissa	70	data=new double2[Dim];
76			for(int i=0;i<Dim1;i++)
77			{
78			OT_VECTEUR_4DD z1=v1[i];
79			for (int j=0;j<Dim2;j++)
80			{
81			OT_VECTEUR_4DD z2=v2[j];
82			double2 prscal=z1*z2;
83	souaissa	75	double2 prciseion=0.0;
84			if(prscal.get_fabs()<=prciseion) prscal=0.;
85	souaissa	70	data[i*Dim2+j]=prscal;
86			}
87	souaissa	69	}
88	francois	17
89	souaissa	69
90	souaissa	58	}
91	francois	17
92	souaissa	69
93			OT_TENSEUR::~OT_TENSEUR()
94	souaissa	58	{
95	souaissa	69	Dim1=0;
96			Dim2=0;
97	souaissa	70
98			if (data!=NULL) delete [] data;
99	souaissa	58	}
100	francois	17
101
102	souaissa	69	OT_TENSEUR::OT_TENSEUR(const OT_TENSEUR& ot)
103	francois	17	{
104	souaissa	69	Dim1=ot.Dim1;
105			Dim2=ot.Dim2;
106	souaissa	70	data=new double2[Dim1*Dim2];
107
108			for(int i=0;i<Dim1;i++)
109			{
110			for (int j=0;j<Dim2;j++)
111			{
112			data[i*Dim2+j]=ot(i,j);
113			}
114			}
115	francois	17	}
116
117	souaissa	58
118	souaissa	69	OT_TENSEUR& OT_TENSEUR:: operator =(const OT_TENSEUR& ot)
119	francois	17	{
120	souaissa	58
121	souaissa	69	if(this==&ot) return *this ;
122	souaissa	70
123			if(data) delete [] data;
124
125	souaissa	69	Dim1=ot.Dim1;
126			Dim2=ot.Dim2;
127	souaissa	70	data=new double2[Dim1*Dim2];
128
129			for(int i=0;i<Dim1;i++) {
130			for (int j=0;j<Dim2;j++) {
131			data[iDim2+j]=ot.data[iDim2+j];
132			}
133			}
134
135
136	souaissa	69	return *this;
137	francois	17	}
138
139	souaissa	69	int operator==( OT_TENSEUR& t1, OT_TENSEUR&t2)
140	francois	17	{
141	souaissa	69	if(t1.Get_NbRows()!=t2.Get_NbRows()\|\|t1.Get_NbCols()!=t2.Get_NbCols()) return 0;
142	souaissa	58	else {
143	souaissa	69	for(int i=0;i<t1.Get_NbRows();i++)
144	souaissa	58	{
145	souaissa	69	for(int j=0;j<t1.Get_NbCols();j++)
146	souaissa	58	{
147	souaissa	75	double2 val1=t1(i,j);
148			double2 val2=t2(i,j);
149			if(val1==val2) continue;
150			else return 0;
151	souaissa	58	}
152	francois	17
153	souaissa	58	}
154			}
155			return 1;
156			}
157	francois	17
158
159	souaissa	69
160			double2&OT_TENSEUR:: operator()(int i,int j)
161	souaissa	58	{
162	souaissa	69	assert(i<Dim1&&j<Dim2);
163
164	souaissa	70	return data[i*Dim2+j];
165	souaissa	58	}
166	francois	17
167
168	souaissa	70
169			double2& OT_TENSEUR:: operator()(int i,int j)const
170	souaissa	58	{
171	souaissa	69	assert(i<Dim1&&j<Dim2);
172
173	souaissa	70	return data[i*Dim2+j];
174	souaissa	58	}
175	francois	17
176	souaissa	58
177	souaissa	70
178	souaissa	69	int OT_TENSEUR:: Get_NbRows()const
179	francois	17	{
180	souaissa	69	return Dim1;
181	francois	17	}
182
183	souaissa	69
184
185			int OT_TENSEUR:: Get_NbCols()const
186	souaissa	58	{
187	souaissa	69	return Dim2;
188	souaissa	58	}
189	francois	17
190
191	souaissa	58
192	souaissa	70	int OT_TENSEUR:: est_til_equivalent( OT_TENSEUR& tns)
193	souaissa	58	{
194			vector<int> check;
195
196	souaissa	69	for(int i=0;i<tns.Get_NbRows();i++)
197	souaissa	58	{
198
199	souaissa	69	vector<double2> lig1;
200			for(int ix1=0;ix1<Dim2;ix1++) lig1.insert(lig1.end(),data[i*Dim2+ix1]);
201
202			for(int j=0;j<Dim2;j++)
203	francois	17	{
204	souaissa	69	int p=-1;
205	souaissa	58	if(!existe(check,j,p))
206			{
207			vector<int>trouve;
208			vector<double2>valeur;
209			vector<double2> lig2;
210	souaissa	69	for(int ix2=0;ix2<Dim2;ix2++) lig2.insert(lig2.end(),tns(j,ix2));
211			for(int r=0;r<Dim2;r++)
212	souaissa	58	{
213	souaissa	69	double2 val=lig1[r];
214			val=val.get_fabs();
215	souaissa	70	//val.set_dx(0.01);
216	souaissa	58	int indx1=-1;
217	francois	17
218	souaissa	69	for(int s=0;s<Dim2;s++)
219	souaissa	58	{ int pos1;
220	souaissa	69	double2 vall=lig2[s];
221			vall=vall.get_fabs();
222	souaissa	70	//vall.set_dx(0.01);
223	souaissa	58	if (val==vall&&!existe(trouve,s,indx1))//&&!existe(valeur,vall,pos1))
224			{
225			trouve.insert(trouve.end(),s);
226			valeur.insert(valeur.end(),val);
227			break;
228			}
229
230			}
231
232			}
233			int t_trouve=trouve.size();
234	souaissa	69	if (t_trouve==Dim2)
235	souaissa	58	{
236			check.insert(check.end(),j);
237			break;
238
239	francois	17	}
240	souaissa	58
241			}
242			}
243			continue;
244			}
245			int t_chek=check.size();
246	souaissa	69	if (t_chek==Dim1)return 1;
247	souaissa	58	else return 0;
248	francois	17
249	souaissa	58	}
250
251
252
253	souaissa	69	template <class T> int OT_TENSEUR:: existe(const vector<T>& vct,T val,int& idx)
254	souaissa	58	{
255	francois	17
256	souaissa	58	if(val<0) val=-val;
257	souaissa	69	for(unsigned int i=0;i<vct.size();i++)
258	souaissa	58	{
259			double x=vct[i];
260			if(x<0) val=-val;
261			if(x==val){
262			idx=i; return 1; }
263			}
264			idx=-1;
265			return 0;
266			}
267	souaissa	69
268
269			OT_TENSEUR operator+( OT_TENSEUR& ot1, OT_TENSEUR& ot2)
270			{
271			assert (ot1.Get_NbRows()==ot2.Get_NbRows()&&ot1.Get_NbCols()==ot2.Get_NbCols());
272	souaissa	70	OT_TENSEUR res(ot1.Get_NbRows(),ot1.Get_NbCols());
273	souaissa	69	for(int i=0;i<ot1.Get_NbRows();i++)
274			{
275			for(int j=0;j<ot2.Get_NbCols();j++)
276			{
277			res(i,j)=ot1(i,j)+ot2(i,j);
278			}
279			}
280
281			return res;
282			}
283
284
285	souaissa	70	OT_TENSEUR operator*( double2& coef, OT_TENSEUR& ot)
286			{
287			OT_TENSEUR T=ot;
288			for(int i=0;i<ot.Get_NbRows();i++)
289			for(int j=0;j<ot.Get_NbCols();j++)
290	souaissa	69
291	souaissa	70	T(i,j)=coef*T(i,j);
292
293			return T;
294
295			}
296
297			OT_TENSEUR operator*(OT_TENSEUR& ot , double2& coef)
298			{
299			return coef*ot;
300			}
301
302			OT_TENSEUR OT_TENSEUR:: transpose()
303			{
304			OT_TENSEUR TMP(Dim2,Dim1);
305			for(int i=0;i<Dim1;i++) {
306			for(int j=0;j<Dim2;j++)
307			{
308			TMP(j,i)=(*this)(i,j); // ot(i,j);
309			}
310			}
311			return TMP;
312			}
313
314			OT_TENSEUR operator*( OT_TENSEUR& ot1, OT_TENSEUR& ot2)
315			{
316			int dim1=ot1.Get_NbRows();
317			int dim3= ot1.Get_NbCols();
318			int dim2=ot2.Get_NbCols();
319			OT_TENSEUR TMP (dim1,dim2);
320
321			for(int i=0;i<dim1;i++)
322			{
323			for(int j=0;j<dim2;j++)
324			{
325			for(int k=0;k<dim3;k++)
326			TMP(i,j)=TMP(i,j)+ot1(i,k)*ot2(k,j);
327			}
328			}
329			return TMP;
330			}
331
332	souaissa	69	std::ostream& operator<<(std::ostream& os, OT_TENSEUR& tns){
333	souaissa	70	os<<"------------------------------------------------------------------------------------"<<endl;
334			os<<"# TENSEUR: DIM1= "<< tns.Get_NbRows()<<" , DIM2: "<<tns.Get_NbCols()<<" #"<<endl;
335			os<<"------------------------------------------------------------------------------------"<<endl;
336	souaissa	69	for(int i=0;i<tns.Get_NbRows();i++) {
337			for(int j=0;j<tns.Get_NbCols();j++) {
338			os<<setw(18)<<tns(i,j);
339			}
340			os<<endl;
341			}
342	souaissa	70	os<<"------------------------------------------------------------------------------------"<<endl;
343	souaissa	69	return os;
344			}
345
346
347
348			void OT_TENSEUR::get_orthogonalisation(OT_TENSEUR& a,OT_VECTEUR_4DD& d,OT_TENSEUR& v,int n,int&nrot)
349			{
350
351			static const int ITMAX=50;
352			double2 zro=0.0;
353			int n1=n-1,n2=1/(n*n),ip1,iq1;// ,ip,j,i,iq// n=3,
354			double2 tresh,theta,tau,t,sm,s,h,g,c;
355			double2 *b=new double2[n]; // OT_VECTEUR_3D b,z;
356			double2 *z=new double2[n];
357
358			for(int ip=0;ip<n;ip++)
359			{
360			for(int iq=0;iq<n;iq++)
361			//v[ip*n+iq]=0.0;
362			v(ip,iq)=0.0;
363			v(ip,ip)=1.0;
364			b[(unsigned int)ip]=d[(unsigned int)ip]=a(ip,ip);
365			z[(unsigned int)ip]=0.0;
366			}
367			nrot=0;
368			for(int i=1;i<=ITMAX;i++)
369			{
370			sm=0.0;
371			for(int ip=0;ip<n1;ip++)
372			for(int iq=ip+1;iq<n;iq++)
373			sm =sm+(a(ip,iq)).get_fabs();
374			if(sm==zro)
375			return;
376			if(i<4)
377			tresh=0.2smn2;
378			else
379			tresh=0.0;
380			for(int ip=0;ip<n1;ip++)
381			for(int iq=ip+1;iq<n;iq++)
382			{
383			g=100*(a(ip,iq)).get_fabs();
384			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())
385			a(ip,iq)=0.0;
386			else if((a(ip,iq)).get_fabs()>tresh)
387			{
388			h=d[(unsigned int)iq]-d[(unsigned int)ip];
389			if(((h).get_fabs()+g)==(h).get_fabs())
390			t=a(ip,iq)/h;
391			else
392			{
393			theta=0.5*h/a(ip,iq);
394			t=1.0/((theta).get_fabs()+sqrt(1.0+theta*theta));
395			if(theta<zro)t=zro-t;
396			}
397			c=1.0/sqrt(1+t*t);
398			s=t*c;
399			tau=s/(1.0+c);
400			h=t*a(ip,iq);
401			z[(unsigned int)ip]=z[(unsigned int)ip]-h;
402			z[(unsigned int)iq]=z[(unsigned int)iq]+h;
403			d[(unsigned int)ip]=d[(unsigned int)ip]-h;
404			d[(unsigned int)iq]=d[(unsigned int)iq]+h;
405			a(ip,iq)=0.0;
406			ip1=ip-1;
407			iq1=iq-1;
408			for(int j=0;j<=ip1;j++)
409			{
410			g=a(j,ip);
411			h=a(j,iq);
412			a(j,ip)=g-s(h+gtau);
413			a(j,iq)=h+s(g-htau);
414			}
415			for(int j=ip+1;j<=iq1;j++)
416			{
417			g=a(ip,j);
418			h=a(j,iq);
419			a(ip,j)=g-s(h+gtau);
420			a(j,iq)=h+s(g-htau);
421			}
422			for(int j=iq+1;j<n;j++)
423			{
424			g=a(ip,j);
425			h=a(iq,j);
426			a(ip,j)=g-s(h+gtau);
427			a(iq,j)=h+s(g-htau);
428			}
429			for(int j=0;j<n;j++)
430			{
431			g=v(j,ip);
432			h=v(j,iq);
433			v(j,ip)=g-s(h+gtau);
434			v(j,iq)=h+s(g-htau);
435			}
436			++nrot;
437			}
438			}
439			for(int ip=0;ip<n;ip++)
440			{
441			b[(unsigned int)ip]=b[(unsigned int)ip]+z[(unsigned int)ip];
442			d[(unsigned int)ip]=b[(unsigned int)ip];
443			z[(unsigned int)ip]=0.0;
444			}
445			}
446
447
448			delete [] b;
449			delete [] z;
450	souaissa	70
451	souaissa	58	}
452	francois	17
453	souaissa	69
454	souaissa	58
455
456
457
458
459
460
461
462