Index: /issm/trunk/src/c/Bamgx/Metric.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Metric.h	(revision 2932)
+++ /issm/trunk/src/c/Bamgx/Metric.h	(revision 2933)
@@ -39,5 +39,5 @@
 		Real8 operator()(R2 x) const { return sqrt((x,x))/h;};
 		Real8 operator()(R2 x,R2 y) const { return ((x,y))/(h*h);};
-		int  IntersectWith(MetricIso M) {int r=0;if (M.h<h) r=1,h=M.h;return r;}
+		int  IntersectWith(MetricIso M) {int change=0;if (M.h<h) change=1,h=M.h;return change;}
 		MetricIso operator*(Real8 c) const {return  MetricIso(h/c);} 
 		MetricIso operator/(Real8 c) const {return  MetricIso(h*c);}
@@ -122,5 +122,5 @@
 	}
 
-	void  SimultaneousMatrixReduction( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V);
+	void  SimultaneousMatrixReduction( MetricAnIso M1,  MetricAnIso M2,D2xD2 &V);
 
 	inline MetricAnIso::MetricAnIso(const MatVVP2x2 M) {
Index: /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2932)
+++ /issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp	(revision 2933)
@@ -73,31 +73,49 @@
 		/*Get a new metric from an existing metric (M1=this)
 		 * and a new metric given in input M2 using a 
-		 * Simultaneous Matrix Reduction
-		 * */
-
-		int r=0;
-		MetricAnIso & M1 = *this;
-		D2xD2 M;
-		double l1,l2;
-
-		SimultaneousMatrixReduction(*this,M2,l1,l2,M);
-
-		R2 v1(M.x.x,M.y.x);
-		R2 v2(M.x.y,M.y.y);
-		double l11=M1(v1,v1);
-		double l12=M1(v2,v2);
-		double l21=M2(v1,v1);
-		double l22=M2(v2,v2);
-		if ( l11 < l21 )  r=1,l11=l21;
-		if ( l12 < l22 )  r=1,l12=l22; 
-		if (r) { // change
-			D2xD2 M_1(M.inv());
-			D2xD2 D(l11,0,0,l12); 
-			D2xD2 Mi(M_1.t()*D*M_1);
-			a11=Mi.x.x;
-			a21=0.5*(Mi.x.y+Mi.y.x);
-			a22=Mi.y.y;
-		}
-		return r;
+		 * Simultaneous Matrix Reduction:
+		 * If M1 and M2 are 2 metrics, we must build N=M1^-1 M2 (Alauzet2003 p16) 
+		 * the eigen vectors of N form a matrix P
+		 * The new metric M = M1 inter M2 is then given by:
+		 *
+		 *      -T [ max(lambda1, mu1)          0         ]  -1				 
+		 * M = P   [                                      ] P		 
+		 *         [        0            max(lambda2, mu2)] 
+		 *
+		 * where lambdai and mui can be computed using Rayleigh formula: 
+		 *    lambdai = vi' M1 vi
+		 * with vi eigen vectors of N (columns of P)
+		 */
+
+		int         change=0;
+		MetricAnIso &M1=*this;
+		D2xD2       P;
+
+		//Get P, eigen vectors of N=inv(M1) M2
+		SimultaneousMatrixReduction(*this,M2,P);
+
+		//extract the eigen vectors of P (columns)
+		R2 v1(P.x.x,P.y.x);
+		R2 v2(P.x.y,P.y.y);
+
+		//compute lambdai mui
+		double lambda1=M1(v1,v1);
+		double lambda2=M1(v2,v2);
+		double mu1=M2(v1,v1);
+		double mu2=M2(v2,v2);
+
+		//check where any change needs to be done on M1
+		if ( lambda1 < mu1 )  change=1,lambda1=mu1;
+		if ( lambda2 < mu2 )  change=1,lambda2=mu2; 
+
+		//update M1 if necessary
+		if (change) {
+			D2xD2 invP(P.inv());
+			D2xD2 D(lambda1,0,0,lambda2); 
+			D2xD2 M(invP.t()*D*invP);
+			a11=M.x.x;
+			a21=0.5*(M.x.y+M.y.x);
+			a22=M.y.y;
+		}
+		return change;
 	}
 	/*}}}1*/
@@ -171,53 +189,104 @@
 	/*}}}1*/
 	/*FUNCTION SimultaneousMatrixReduction{{{1*/
-	void SimultaneousMatrixReduction( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) {
+	void SimultaneousMatrixReduction( MetricAnIso M1,  MetricAnIso M2, D2xD2 &V) {
+		/*In this routine we must return a matrix V that is composed of the 
+		 * eigen vectors of N=inv(M1) M2.
+		 * Instead of looking at N directly, we are going to use the fact that
+		 * M1 and M2 are symmetrical, positive definite. 
+		 * The eigen values of N are given by solving
+		 *    inv(M1) M2 V = lambda V
+		 * which is equivalent to
+		 *    M2 V = lambda M1 V
+		 * and we will hence solve
+		 *    (M2 - lambda M1) V = 0
+		 */
+
+		int i;
+
+		//M1 and M2 components
 		double a11=M1.a11,a21=M1.a21,a22=M1.a22;
 		double b11=M2.a11,b21=M2.a21,b22=M2.a22;
-		//  M1 v = l M2 v
-		// (M1 - l M2) v =0
-		// det (M1 - l M2) =0
-		// det (M1 - l M2) = a l^2 + b l + c;
-		// = (a11 - l * b11) * (a22 - l * b22) - (a21 - l * b21 ) ^2
-		const double /*c11 = a11*a11,*/ c21= a21*a21;
-		const double /*d11 = b11*b11,*/ d21= b21*b21;
-		const double a=b11*b22 - d21;
-		const double b=-a11*b22-a22*b11+2*a21*b21;
-		const double c=-c21+a11*a22;
-		const double bb = b*b,ac= a*c;
-		const double b2a = b/(2*a);
-		const double delta = bb - 4 * ac;
+
+		/*To get the eigen values, we solve the following problem:
+		 *    det(M2-lambda M1) = 0
+		 *    (b11 - lambda a11)(b22-lambda a22) - (b21-lambda a21)^2
+		 * and we have the following trinome:
+		 *    a lambda^2 + b lambda + c =0
+		 * with:
+		 *    a = a11 a22 - a21 a21 (=det(M1))
+		 *    b = -a11 b22 -b11 a22 + 2 b21 a21
+		 *    c = b11 b22 - b21 b21 (=det(M2))
+		 *    */
+		const double a= a11*a22  - a21*a21;
+		const double b=-a11*b22 - b11*a22+2*b21*a21;
+		const double c=-b21*b21 + b11*b22;
+		const double bb=b*b,ac=a*c;
+		const double delta= bb-4*ac;
+
+		// first, case of a double root if:
+		//  - all the terms are very small (a??)
+		//  - or : delta is very small
 		if ( (bb + Abs(ac) < 1.0e-34 ) ||  (delta < 1.0e-6*bb) ){
-			// racine double;
-			if (Abs(a) < 1.e-30 )
-			 l1 = l2 = 0;
-			else 
-			 l1=l2=-b/(2*a); 
+			//all vectors are eigen vectors -> choose 1,0 and 0,1
 			V= D2xD2(1,0,0,1);
 		}
+
+		//general case: two distinct roots: lambda1 and lambda2
 		else {
+
+			/*Compute eigen values*/
 			const double delta2 = sqrt(delta);
-			l1= (-b - delta2)/(2*a);
-			l2= (-b + delta2)/(2*a);
-			// M1 v = l M2 v
-			//  ( (M1 - I M2) x,y)  = (x,(M1 - I M2) y) \forall y
-			// so Ker((M1 - I M2)) = Im((M1 - I M2))^\perp
-			double v0 = a11-l1*b11, v1 = a21-l1*b21,v2 = a22 - l1*b22;
-			double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
-			double vp1x,vp1y,vp2x,vp2y;
-
-			if(s1 < s0)
-			 s0=sqrt(s0),vp1x=v1/s0,vp1y=-v0/s0;
-			else
-			 s1=sqrt(s1),vp1x=v2/s1,vp1y=-v1/s1;
-
-			v0 = a11-l2*b11, v1 = a21-l2*b21,v2 = a22 - l2*b22;
-			s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
-			if(s1 < s0)
-			 s0=sqrt(s0),vp2x=v1/s0,vp2y=-v0/s0;
-			else
-			 s1=sqrt(s1),vp2x=v2/s1,vp2y=-v1/s1;
-			V=D2xD2(vp1x,vp2x,vp1y,vp2y);
-		}
-		return;
+			double lambda[2];
+			lambda[0]= (-b - delta2)/(2*a);
+			lambda[1]= (-b + delta2)/(2*a);
+
+			/*compute eigen vectors*/
+			double vp[2][2];
+			double v0,v1,v2;
+			double s0,s1;
+
+			for(i=0;i<2;i++){
+				/*Now, one must find the eigen vectors. For that we use the 
+				 * following property of the scalare product
+				 *    (Ax,b) = transp(b) Ax = transp(x) transp(A) b
+				 *           = (transp(A) b ,x)
+				 * Here we are dealing with A= M2 - lambda M1 which is symmetrical:
+				 *    for all (x,y) in R2 
+				 *       ((M2 - lambda M1)x,y)=((M2 - lambda M1)y,x)
+				 * If y is in Ker(M2 - lambda M1):
+				 *    for all x in R2
+				 *       ((M2 - lambda M1)y,x)=0
+				 * This shows that:
+				 *    Ker(M2 - lambda M1) is orthogonal to Im(M2 - lambda M1)
+				 * To find the eigen vectors, we only have to find two vectors
+				 * of the image and take their perpendicular as long as they are
+				 * not 0.
+				 * To do that, we take (1,0) and (0,1) and take the larger norm*/
+
+				//compute V = M2 - lambdai M1
+				v0 = b11 - lambda[i]*a11;
+				v1 = b21 - lambda[i]*a21;
+				v2 = b22 - lambda[i]*a22;
+
+				// compute s1=norm(V(1,0)) and s0=norm(V(0,1))
+				s0 = v0*v0 + v1*v1;
+				s1 = v1*v1 + v2*v2;
+
+				//compute vp1 = (vp1x,vp1y)
+				if(s1 < s0){
+					s0=sqrt(s0);
+					vp[0][i]=   v1/s0;
+					vp[1][i]= - v0/s0;
+				}
+				else{
+					s1=sqrt(s1);
+					vp[0][i]=   v2/s0;
+					vp[1][i]= - v1/s0;
+				}
+			}
+
+			//compute V from vp
+			V=D2xD2(vp[0][0],vp[0][1],vp[1][0],vp[1][1]);
+		}
 	}
 	/*}}}1*/