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.resize(0);
}
OT_TENSEUR::OT_TENSEUR(int n)
{
 Dim1=n;
 Dim2=n;
 data.resize(n*n);
   for(int i=0;i<n*n;i++) data[i]=0.0;
}
OT_TENSEUR::OT_TENSEUR(int n,int m)
{
 Dim1=n;
 Dim2=m;
 data.resize(n*m);
   for(int i=0;i<n*m;i++) data[i]=0.0;
}
OT_TENSEUR:: OT_TENSEUR(vector<OT_VECTEUR_4DD>& v)
{
 Dim1=v.size();
 Dim2=Dim1;
 int Dim= Dim1*Dim2  ;
 data.resize(Dim);
 vector<OT_VECTEUR_4DD> ::iterator i1=v.begin();
 vector<OT_VECTEUR_4DD> ::iterator i2=v.begin();

  while( i1!=v.end())
  {
   OT_VECTEUR_4DD z1=*i1;
    while( i2!=v.end())
      {
      OT_VECTEUR_4DD z2=*i2;;
      data.insert(data.end(),z1*z2);

      i2++;
     }
     i1++;
  }


}

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.resize(Dim);
 vector<OT_VECTEUR_4DD> ::iterator i1=v1.begin();
 vector<OT_VECTEUR_4DD> ::iterator i2=v2.begin();

  while( i1!=v1.end())
  {
   OT_VECTEUR_4DD z1=*i1;
    while( i2!=v2.end())
      {
      OT_VECTEUR_4DD z2=*i2;;
      data.insert(data.end(),z1*z2);

      i2++;
     }
     i1++;
   }


}


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


OT_TENSEUR::OT_TENSEUR(const OT_TENSEUR& ot)
{
Dim1=ot.Dim1;
Dim2=ot.Dim2;
data=ot.data;
}


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

 if(this==&ot) return *this ;
  Dim1=ot.Dim1;
  Dim2=ot.Dim2;
  data=ot.data;
 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++)
       {
      if(t1(i,j)!=t2(i,j)) return 0;
       }

    }
     }
return 1;
}


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

vector<double2>::iterator k= data.begin();
k=k+i*Dim2+j;
return *k;
}


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

vector<const double2>::iterator k= data.begin();
k=k+i*Dim2+j;
return *k;
}
*/


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


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


int OT_TENSEUR:: is_equivalent_at( OT_TENSEUR& tns)
 {
vector<int> check;

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

//  bool verif=false;
  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;
    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::ostream& operator<<(std::ostream& os, OT_TENSEUR& tns){
  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;
            }
            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;
//#undef a
//#undef v
}


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(Dim1,Dim2);
 for(int i=0;i<Dim1;i++) {
    for(int j=0;j<Dim2;j++)
    {
TMP(i,j)=data[i*Dim2+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;
}


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