Explanation of issue with multiplying matrices element-wise and numpy ndarray views. This is what caused to the difference between the numpy and fortran pulse propagator functions in my blog post (http://themodernscientist.com/posts/2013/2013-06-09-simulation_of_nmr_shaped_pulses/)

numpy_ndarray_views.py

import numpy as np
a = np.eye(3)
b = np.eye(3)
c = np.eye(3)

for i in range(10):
    d = np.random.random((3,3))
    
    # copy each matrix element individually
    d00 = d[0,0].copy()
    d01 = d[0,1].copy()
    d02 = d[0,2].copy()
    d10 = d[1,0].copy()
    d11 = d[1,1].copy()
    d12 = d[1,2].copy()
    d20 = d[2,0].copy()
    d21 = d[2,1].copy()
    d22 = d[2,2].copy()
    c00 = c[0,0].copy()
    c01 = c[0,1].copy()
    c02 = c[0,2].copy()
    c10 = c[1,0].copy()
    c11 = c[1,1].copy()
    c12 = c[1,2].copy()
    c20 = c[2,0].copy()
    c21 = c[2,1].copy()
    c22 = c[2,2].copy()
    
    # a = d * a
    a = np.dot(d,a)
    
    # b = d * b --> this operates on variables that are views
    # which causes their values to change during the multiplication
    b[0,0] = d[0,0]*b[0,0] + d[0,1]*b[1,0] + d[0,2]*b[2,0]
    b[0,1] = d[0,0]*b[0,1] + d[0,1]*b[1,1] + d[0,2]*b[2,1]
    b[0,2] = d[0,0]*b[0,2] + d[0,1]*b[1,2] + d[0,2]*b[2,2]
    b[1,0] = d[1,0]*b[0,0] + d[1,1]*b[1,0] + d[1,2]*b[2,0]
    b[1,1] = d[1,0]*b[0,1] + d[1,1]*b[1,1] + d[1,2]*b[2,1]
    b[1,2] = d[1,0]*b[0,2] + d[1,1]*b[1,2] + d[1,2]*b[2,2]
    b[2,0] = d[2,0]*b[0,0] + d[2,1]*b[1,0] + d[2,2]*b[2,0]
    b[2,1] = d[2,0]*b[0,1] + d[2,1]*b[1,1] + d[2,2]*b[2,1]
    b[2,2] = d[2,0]*b[0,2] + d[2,1]*b[1,2] + d[2,2]*b[2,2]
    
    # c = d * c
    c[0,0] = d00*c00 + d01*c10 + d02*c20
    c[0,1] = d00*c01 + d01*c11 + d02*c21
    c[0,2] = d00*c02 + d01*c12 + d02*c22
    c[1,0] = d10*c00 + d11*c10 + d12*c20
    c[1,1] = d10*c01 + d11*c11 + d12*c21
    c[1,2] = d10*c02 + d11*c12 + d12*c22
    c[2,0] = d20*c00 + d21*c10 + d22*c20
    c[2,1] = d20*c01 + d21*c11 + d22*c21
    c[2,2] = d20*c02 + d21*c12 + d22*c22

print 'max(abs(a-b)) ', np.max(np.abs(a-b))
print 'max(abs(a-c)) ', np.max(np.abs(a-c))
print 'a', a
print 'b', b #oops
print 'c', c

Here are the results:

max(abs(a-b))  1489.53863842
max(abs(a-c))  0.0
a [[ 14.29088931  17.46476968  27.7205024 ]
 [ 28.49985681  34.82941377  55.28208071]
 [ 42.15901178  51.52210923  81.77711907]]
b [[  5.46872177e-01   2.01405996e+02   4.85068475e+02]
 [  1.04308000e+00   3.84152960e+02   9.25198329e+02]
 [  1.77152075e+00   6.52428328e+02   1.57131576e+03]]
c [[ 14.29088931  17.46476968  27.7205024 ]
 [ 28.49985681  34.82941377  55.28208071]
 [ 42.15901178  51.52210923  81.77711907]]