r24202 - /branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py -- June 20, 2014 - 15:22

Author: tlinnet
Date: Fri Jun 20 15:22:16 2014
New Revision: 24202

URL: http://svn.gna.org/viewcvs/relax?rev=24202&view=rev
Log:
Added the "dtype" argument to function matrix_exponential_rankN.

This is to force the conversion of dtype, if they are of other type.
This can be conversion from complex128 to complex64.

Task #7807 (https://gna.org/task/index.php?7807): Speed-up of dispersion 
models for Clustered analysis.

Modified:
    branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py

Modified: branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py?rev=24202&r1=24201&r2=24202&view=diff
==============================================================================
--- branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py   
(original)
+++ branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py   Fri 
Jun 20 15:22:16 2014
@@ -57,15 +57,17 @@
         return array(eA.real)
 
 
-def matrix_exponential_rankN(A):
+def matrix_exponential_rankN(A, dtype=None):
     """Calculate the exact matrix exponential using the eigenvalue 
decomposition approach, for higher dimensional data.
 
     Here X is the Row and Column length, of the outer square matrix.
 
-    @param A:   The square matrix to calculate the matrix exponential of.
-    @type A:    numpy float array of rank [NE][NS][NM][NO][ND][X][X]
-    @return:    The matrix exponential.  This will have the same 
dimensionality as the A matrix.
-    @rtype:     numpy float array of rank [NE][NS][NM][NO][ND][X][X]
+    @param A:       The square matrix to calculate the matrix exponential of.
+    @type A:        numpy float array of rank [NE][NS][NM][NO][ND][X][X]
+    @param dtype:   If provided, forces the calculation to use the data type 
specified.
+    @type type:     data-type, optional
+    @return:        The matrix exponential.  This will have the same 
dimensionality as the A matrix.
+    @rtype:         numpy float array of rank [NE][NS][NM][NO][ND][X][X]
     """
 
     # Set initial to None.
@@ -75,6 +77,15 @@
         NE, NS, NM, NO, ND, Row, Col = A.shape
     elif len(A.shape) == 6:
         NS, NM, NO, ND, Row, Col = A.shape
+
+    # Convert dtype, if specified.
+    if dtype != None:
+        dtype_mat = A.dtype
+
+        # If the dtype is different from the input.
+        if dtype_mat != dtype:
+            # This needs to be made as a copy.
+            A = A.astype(dtype)
 
     # Is the original matrix real?
     complex_flag = any(iscomplex(A))
@@ -102,7 +113,11 @@
         eye_mat = tile(eye(Row)[newaxis, newaxis, newaxis, newaxis, newaxis, 
...], (NE, NS, NM, NO, ND, 1, 1) )
 
     # Transform it to a diagonal matrix, with elements from vector down the 
diagonal.
-    W_exp_diag = multiply(W_exp, eye_mat )
+    # Use the dtype, if specified.
+    if dtype != None:
+        W_exp_diag = multiply(W_exp, eye_mat, dtype=dtype )
+    else:
+        W_exp_diag = multiply(W_exp, eye_mat)
 
     # Make dot products for higher dimension.
     # "...", the Ellipsis notation, is designed to mean to insert as many 
full slices (:)