Changeset 236

Timestamp:

05/05/09 11:31:00 (16 years ago)

Author:

Mathieu Morlighem

Message:

Fixed gradient orthogonalisation

File:

• issm/trunk/src/c/Orthx/Orthx.cpp (modified) (1 diff)

issm/trunk/src/c/Orthx/Orthx.cpp

@@ -14 +14 @@
 
         /*intermediary:*/
-        Vec oldgradj2=NULL;
-        double norm,dot_product;;
+        double norm_new,norm_old,dot_product;;
 
+        /*Initialize output*/
+        VecDuplicate(gradj,&newgradj);
+        VecCopy(gradj,newgradj);
 
-        if (!oldgradj){
-                VecDuplicate(gradj,&oldgradj2);
-                VecCopy(gradj,oldgradj2);
-                /*scale to 1:*/
-                VecNorm(oldgradj2,NORM_INFINITY,&norm);
-                VecScale(oldgradj2,1.0/norm);
+        /*rough orthagonalization
+        gradient=gradient-(gradient'*oldgradient)*oldgradient /norm(oldgradient)^2; */
+        if(oldgradj){
+                VecNorm(oldgradj,NORM_2,&norm_old);
+                VecDot(newgradj,oldgradj,&dot_product);
+                VecAXPY(newgradj, -dot_product/pow(norm_old,2), oldgradj);
         }
 
-        /*% normalize gradient to 1:*/
-        VecDuplicate(gradj,&newgradj);
-        VecCopy(gradj,newgradj);
-        VecNorm(newgradj,NORM_INFINITY,&norm);
-        VecScale(newgradj,1.0/norm);
-
-        /*rough orthagonalization
-        gradient=gradient-(gradient'*oldgradient)*oldgradient; 
-        */
-
-        if(oldgradj){
-                VecDot(newgradj,oldgradj,&dot_product);
-                VecAXPY(newgradj, -dot_product,oldgradj);
-        }
-        else{
-                VecDot(newgradj,oldgradj2,&dot_product);
-                VecAXPY(newgradj, -dot_product,oldgradj2);
-        }
+        /*scale to 1: gradient=gradient/max(abs(gradient))*/
+        VecNorm(newgradj,NORM_INFINITY,&norm_new);
+        VecScale(newgradj,1.0/norm_new);
 
         /*Assign correct pointer*/
         *pnewgradj=newgradj;
-
-        /*Free ressources and return*/
-        VecFree(&oldgradj2);
 }