Index: /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 19029)
+++ /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 19030)
@@ -333,4 +333,87 @@
 
 }/*}}}*/
+void Matrix2x2Eigen(IssmDouble* plambda1,IssmDouble* plambda2,IssmDouble* pvx, IssmDouble* pvy,IssmDouble a11, IssmDouble a21,IssmDouble a22){/*{{{*/
+	/*From symetric matrix (a11,a21;a21,a22), get eigen values lambda1 and lambda2 and one eigen vector v*/
+
+	/*Output*/
+	IssmDouble lambda1,lambda2;
+	IssmDouble vx,vy;
+
+	/*To get the eigen values, we must solve the following equation:
+	 *     | a11 - lambda    a21        |
+	 * det |                            | = 0
+	 *     | a21             a22-lambda |
+	 *
+	 * We have to solve the following polynom:
+	 *  lamda^2 + ( -a11 -a22)*lambda + (a11*a22-a21*a21) = 0*/
+
+	/*Compute polynom determinant*/
+	IssmDouble b=-a11-a22;
+	IssmDouble delta=b*b - 4*(a11*a22-a21*a21);
+
+	/*Compute norm of M to avoid round off errors*/
+	IssmDouble normM=a11*a11 + a22*a22 + a21*a21;
+
+	/*1: normM too small: eigen values = 0*/
+	if(normM<1.e-30){
+		lambda1=0.;
+		lambda2=0.;
+		vx=1.;
+		vy=0.;
+	}
+	/*2: delta is small -> double root*/
+	else if (delta < 1.e-5*normM){
+		lambda1=-b/2.;
+		lambda2=-b/2.;
+		vx=1.;
+		vy=0.;
+	}
+	/*3: general case -> two roots*/
+	else{
+		delta   = sqrt(delta);
+		lambda1 = (-b-delta)/2.;
+		lambda2 = (-b+delta)/2.;
+
+		/*Now, one must find the eigen vectors. For that we use the following property of the inner product
+		 *    <Ax,y> = <x,tAy>
+		 * Here, M'(M-lambda*Id) is symmetrical, which gives:
+		 *    ∀(x,y)∈R²xR² <M'x,y> = <M'y,x>
+		 * And we have the following:
+		 *    if y∈Ker(M'), ∀x∈R² <M'x,y> = <x,M'y> = 0
+		 * We have shown that
+		 *    Im(M') ⊥ Ker(M')
+		 *
+		 * To find the eigen vectors of M, we only have to find two vectors
+		 * of the image of M' and take their perpendicular as long as they are
+		 * not 0.
+		 * To do that, we take the images (1,0) and (0,1):
+		 *  x1 = (a11 - lambda)      x2 = a21
+		 *  y1 = a21                 y2 = (a22-lambda)
+		 *
+		 * We take the vector that has the larger norm and take its perpendicular.*/
+
+		IssmDouble norm1 = (a11-lambda1)*(a11-lambda1) + a21*a21; 
+		IssmDouble norm2 = a21*a21 + (a22-lambda1)*(a22-lambda1);
+
+		if(norm2<norm1){
+			norm1=sqrt(norm1);
+			vx = - a21/norm1;
+			vy = (a11-lambda1)/norm1;
+		}
+		else{
+			norm2=sqrt(norm2);
+			vx = - (a22-lambda1)/norm2;
+			vy = a21/norm2;
+		}
+	}
+
+	/*Assign output*/
+	*plambda1 = lambda1;
+	*plambda2 = lambda2;
+	*pvx      = vx;
+	*pvy      = vy;
+
+
+}/*}}}*/
 
 void Matrix3x3Determinant(IssmDouble* Adet,IssmDouble* A){/*{{{*/
Index: /issm/trunk-jpl/src/c/shared/Matrix/matrix.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Matrix/matrix.h	(revision 19029)
+++ /issm/trunk-jpl/src/c/shared/Matrix/matrix.h	(revision 19030)
@@ -14,4 +14,5 @@
 void Matrix2x2Invert(IssmDouble* Ainv, IssmDouble* A);
 void Matrix2x2Determinant(IssmDouble* Adet,IssmDouble* A);
+void Matrix2x2Eigen(IssmDouble* plambda1,IssmDouble* plambda2,IssmDouble* pvx, IssmDouble* pvy,IssmDouble a11, IssmDouble a21,IssmDouble a22);
 
 void Matrix3x3Invert(IssmDouble* Ainv, IssmDouble* A);