Changeset 2933

Timestamp:

01/29/10 11:33:45 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added some useful comments to Simultaneous reduction

Location:

issm/trunk/src/c/Bamgx

Files:

• Metric.h (modified) (2 diffs)
• objects/MetricAnIso.cpp (modified) (2 diffs)

issm/trunk/src/c/Bamgx/Metric.h

@@ -39 +39 @@
                 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 +122 @@
         }
 
-        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) {

issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp

@@ -73 +73 @@
                 /*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 +189 @@
         /*}}}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*/