Source Code for Module lib.linear_algebra.kronecker_product

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2009-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 the calculation of the Kronecker product.""" 
 24   
 25  # Python imports. 
 26  from numpy import outer 
 27   
 28   
 29   
 30 -def kron_prod(A, B): 
 31      """Calculate the Kronecker product of the matrices A and B. 
 32   
 33      @param A:   ixj matrix. 
 34      @type A:    rank-2 numpy array 
 35      @param B:   kxl matrix. 
 36      @type B:    rank-2 numpy array 
 37      @return:    The Kronecker product. 
 38      @rtype:     ikxjl rank-2 numpy array 
 39      """ 
 40   
 41      # The outer product. 
 42      C = outer(A, B) 
 43   
 44      # Redefine the shape. 
 45      orig_shape = C.shape 
 46      C.shape = A.shape + B.shape 
 47   
 48      # Generate and return the Kronecker product matrix. 
 49      transpose_23(C) 
 50      C.shape = orig_shape 
 51      return C 
 52   
 53   
 54 -def transpose_12(matrix): 
 55      """Perform the 1,2 transpose of a rank-4, 3D tensor. 
 56   
 57      @param matrix:  The rank-4 tensor either in (9, 9) shape for a matrix or the (3, 3, 3, 3) shape 
 58                      for the tensor form. 
 59      @type matrix:   numpy array 
 60      """ 
 61   
 62      # Reshape if necessary. 
 63      reshape = False 
 64      if matrix.shape == (9, 9): 
 65          reshape = True 
 66          matrix.shape = (3, 3, 3, 3) 
 67   
 68      # Perform the transpose. 
 69      for i in range(3): 
 70          for j in range(i, 3): 
 71              for k in range(3): 
 72                  for l in range(3): 
 73                      # Store the element. 
 74                      element = matrix[i, j, k, l] 
 75   
 76                      # Overwrite. 
 77                      matrix[i, j, k, l] = matrix[j, i, k, l] 
 78                      matrix[j, i, k, l] = element 
 79   
 80      # Undo the reshape. 
 81      if reshape: 
 82          matrix.shape = (9, 9) 
 83   
 84   
 85 -def transpose_13(matrix): 
 86      """Perform the 1,3 transpose of a rank-4, 3D tensor. 
 87   
 88      @param matrix:  The rank-4 tensor either in (9, 9) shape for a matrix or the (3, 3, 3, 3) shape 
 89                      for the tensor form. 
 90      @type matrix:   numpy array 
 91      """ 
 92   
 93      # Reshape if necessary. 
 94      reshape = False 
 95      if matrix.shape == (9, 9): 
 96          reshape = True 
 97          matrix.shape = (3, 3, 3, 3) 
 98   
 99      # Perform the transpose. 
100      for i in range(3): 
101          for j in range(3): 
102              for k in range(i, 3): 
103                  for l in range(3): 
104                      # Store the element. 
105                      element = matrix[i, j, k, l] 
106   
107                      # Overwrite. 
108                      matrix[i, j, k, l] = matrix[k, j, i, l] 
109                      matrix[k, j, i, l] = element 
110   
111      # Undo the reshape. 
112      if reshape: 
113          matrix.shape = (9, 9) 
114   
115   
116 -def transpose_14(matrix): 
117      """Perform the 1,4 transpose of a rank-4, 3D tensor. 
118   
119      @param matrix:  The rank-4 tensor either in (9, 9) shape for a matrix or the (3, 3, 3, 3) shape 
120                      for the tensor form. 
121      @type matrix:   numpy array 
122      """ 
123   
124      # Reshape if necessary. 
125      reshape = False 
126      if matrix.shape == (9, 9): 
127          reshape = True 
128          matrix.shape = (3, 3, 3, 3) 
129   
130      # Perform the transpose. 
131      for i in range(3): 
132          for j in range(3): 
133              for k in range(3): 
134                  for l in range(i, 3): 
135                      # Store the element. 
136                      element = matrix[i, j, k, l] 
137   
138                      # Overwrite. 
139                      matrix[i, j, k, l] = matrix[l, j, k, i] 
140                      matrix[l, j, k, i] = element 
141   
142      # Undo the reshape. 
143      if reshape: 
144          matrix.shape = (9, 9) 
145   
146   
147 -def transpose_23(matrix): 
148      """Perform the 2,3 transpose of a rank-4, 3D tensor. 
149   
150      @param matrix:  The rank-4 tensor either in (9, 9) shape for a matrix or the (3, 3, 3, 3) shape 
151                      for the tensor form. 
152      @type matrix:   numpy array 
153      """ 
154   
155      # Reshape if necessary. 
156      reshape = False 
157      if matrix.shape == (9, 9): 
158          reshape = True 
159          matrix.shape = (3, 3, 3, 3) 
160   
161      # Perform the transpose. 
162      for i in range(3): 
163          for j in range(3): 
164              for k in range(j, 3): 
165                  for l in range(3): 
166                      # Store the element. 
167                      element = matrix[i, j, k, l] 
168   
169                      # Overwrite. 
170                      matrix[i, j, k, l] = matrix[i, k, j, l] 
171                      matrix[i, k, j, l] = element 
172   
173      # Undo the reshape. 
174      if reshape: 
175          matrix.shape = (9, 9) 
176   
177   
178 -def transpose_24(matrix): 
179      """Perform the 2,4 transpose of a rank-4, 3D tensor. 
180   
181      @param matrix:  The rank-4 tensor either in (9, 9) shape for a matrix or the (3, 3, 3, 3) shape 
182                      for the tensor form. 
183      @type matrix:   numpy array 
184      """ 
185   
186      # Reshape if necessary. 
187      reshape = False 
188      if matrix.shape == (9, 9): 
189          reshape = True 
190          matrix.shape = (3, 3, 3, 3) 
191   
192      # Perform the transpose. 
193      for i in range(3): 
194          for j in range(3): 
195              for k in range(3): 
196                  for l in range(j, 3): 
197                      # Store the element. 
198                      element = matrix[i, j, k, l] 
199   
200                      # Overwrite. 
201                      matrix[i, j, k, l] = matrix[i, l, k, j] 
202                      matrix[i, l, k, j] = element 
203   
204      # Undo the reshape. 
205      if reshape: 
206          matrix.shape = (9, 9) 
207   
208   
209 -def transpose_34(matrix): 
210      """Perform the 3,4 transpose of a rank-4, 3D tensor. 
211   
212      @param matrix:  The rank-4 tensor either in (9, 9) shape for a matrix or the (3, 3, 3, 3) shape 
213                      for the tensor form. 
214      @type matrix:   numpy array 
215      """ 
216   
217      # Reshape if necessary. 
218      reshape = False 
219      if matrix.shape == (9, 9): 
220          reshape = True 
221          matrix.shape = (3, 3, 3, 3) 
222   
223      # Perform the transpose. 
224      for i in range(3): 
225          for j in range(3): 
226              for k in range(3): 
227                  for l in range(k, 3): 
228                      # Store the element. 
229                      element = matrix[i, j, k, l] 
230   
231                      # Overwrite. 
232                      matrix[i, j, k, l] = matrix[i, j, l, k] 
233                      matrix[i, j, l, k] = element 
234   
235      # Undo the reshape. 
236      if reshape: 
237          matrix.shape = (9, 9) 
238