Rodriguez -> g: some more speedup

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@584 4c865b51-4357-4376-afb4-474e03ccb993

4_grains/Rod2g.m

@@ -21,37 +21,34 @@ function invg = Rod2g(r)
 % AK - 13/11/2008 rewritten with the permutation and identity matrices hard
 % coded rather than derived each time (2x faster!)
 %
-% 2012-03-12, Modified by Nicola Vigano', nicola.vigano@esrf.fr
-%  - Made even faster, by saving partial computation: 6x faster than the latest :)
+% 2012, Modified by Nicola Vigano', nicola.vigano@esrf.fr
+%  - Made even faster, by saving partial computation: 15x faster than the latest :)
 
-    % Initialise the orientation matrix
-    invg = zeros(3, 3);
-    % initialise the permutation tensor
-    % e = calc_perm_tensor();
+    % initialise the permutation tensor (do the 3rd as 1st, to allocate it all)
+    e(:, :, 3) = [ 0     1     0; ...
+                  -1     0     0; ...
+                   0     0     0];
     e(:, :, 1) = [ 0     0     0; ...
                    0     0     1; ...
                    0    -1     0];
     e(:, :, 2) = [ 0     0    -1; ...
                    0     0     0; ...
                    1     0     0];
-    e(:, :, 3) = [ 0     1     0; ...
-                  -1     0     0; ...
-                   0     0     0];
 
-    delta = eye(3);
     % Saving partial computation
-    rSquareNorm = dot(r, r);
+    rSquareNorm = sum(r .^ 2);
 
+    % Initialise the orientation matrix
+    invg = zeros(3, 3);
     % Loop over the elements of the orientation matrix
     for ii = 1:3
         for jj = 1:3
-            invg(ii, jj) = (1 - rSquareNorm) * delta(ii, jj) + 2 * r(ii) * r(jj);
+            invg(ii, jj) = (1 - rSquareNorm) * (ii == jj) + 2 * r(ii) * r(jj);
             % Summing over the remaining component
             for k = 1:3
-                 invg(ii, jj) = invg(ii, jj) - 2 * e(ii, jj, k) * r(k);
+                invg(ii, jj) = invg(ii, jj) - 2 * r(k) * e(ii, jj, k);
             end
-            invg(ii, jj) = invg(ii, jj) / (1 + rSquareNorm);
         end
     end
-
+    invg = invg / (1 + rSquareNorm);
 end