Source Code for Module lib.linear_algebra.matrix_power

 1  ############################################################################### 
 2  #                                                                             # 
 3  # Copyright (C) 2013 Edward d'Auvergne                                        # 
 4  #                                                                             # 
 5  # This file is part of the program relax (http://www.nmr-relax.com).          # 
 6  #                                                                             # 
 7  # This program is free software: you can redistribute it and/or modify        # 
 8  # it under the terms of the GNU General Public License as published by        # 
 9  # the Free Software Foundation, either version 3 of the License, or           # 
10  # (at your option) any later version.                                         # 
11  #                                                                             # 
12  # This program is distributed in the hope that it will be useful,             # 
13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
15  # GNU General Public License for more details.                                # 
16  #                                                                             # 
17  # You should have received a copy of the GNU General Public License           # 
18  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
19  #                                                                             # 
20  ############################################################################### 
21   
22  # Module docstring. 
23  """Module for matrix power operations.""" 
24   
25  # Python module imports. 
26  from numpy import diag, dot, eye 
27  from numpy.linalg import eig, inv 
28   
29  # relax module imports. 
30  from lib.errors import RelaxError 
31   
32   
33 -def square_matrix_power(x, y): 
34      """Compute x raised to the power y when x is a square matrix and y is a scalar. 
35   
36      @param x:   The square matrix. 
37      @type x:    numpy rank-2 array 
38      @param y:   The power. 
39      @type y:    float 
40      @return:    The matrix power of x. 
41      @rtype:     numpy rank-2 array 
42      """ 
43   
44      # Sanity check. 
45      s = x.shape 
46      if len(s) != 2 or s[0] != s[1]: 
47          raise RelaxError("The matrix '%s' must be square." % x) 
48   
49      # Catch the zeroth power. 
50      if y == 0: 
51          return eye(s[0]) 
52   
53      # The eigensystem of x. 
54      e, v = eig(x) 
55      d = diag(e) 
56   
57      # Return the matrix power. 
58      return dot(dot(v, d**y), inv(v)) 
59