forward-mode eigh derivatives return NaN's
First off, thanks for developing JAX, it's awesome!

I recently noticed that when one computes higher order derivatives of expressions involving eigenvectors obtained from `jax.numpy.eigh` or `jax.scipy.eigh`, NaNs are returned. However, it only appears to be a problem with forward-mode autodiff. This is unfortunate, since forward-mode tends to be more efficient for higher order derivatives.

A minimal example:
```python
import jax
import jax.numpy as np

# Create a random symmetric matrix
seed = jax.random.PRNGKey(0)
tmp = jax.random.uniform(seed, (2,2))
a = np.dot(tmp, tmp.T)

def test(inp):
    val, vec = np.linalg.eigh(inp)
    return np.dot(np.dot(vec, inp), vec.T)

def test_deriv(func, inp):
    grad_func = jax.jacfwd(func)
    hess_func = jax.jacfwd(grad_func)
    cube_func = jax.jacfwd(hess_func) # This derivative returns NaN, but jax.jacrev works!
    print(grad_func(inp))
    print(hess_func(inp))
    print(cube_func(inp))

test_deriv(test, a)
```

Additional notes:
- Expressions involving eigenvalues do not appear to have this issue
- Any array contractions involving eigenvectors and other arrays (`dot`, `einsum`, etc) also exhibit similar behavior: `jacfwd` fails, while `jacrev` works.
- I'm using CPU
- JAX version 0.1.46

Thanks for your time!
diff --git a/jax/lax/lax.py b/jax/lax/lax.py
index 713f94f55..a0318eed9 100644
--- a/jax/lax/lax.py
+++ b/jax/lax/lax.py
@@ -1720,7 +1720,7 @@ def _div_transpose_rule(cotangent, x, y):
 div_p = standard_binop([_num, _num], 'div')
 ad.defjvp(div_p,
           lambda g, x, y: div(_brcast(g, y), y),
-          lambda g, x, y: div(mul(neg(_brcast(g, x)), x), pow(y, _two(y))))
+          lambda g, x, y: div(mul(neg(_brcast(g, x)), x), square(y)))
 ad.primitive_transposes[div_p] = _div_transpose_rule
 
 rem_p = standard_binop([_num, _num], 'rem')
diff --git a/tests/linalg_test.py b/tests/linalg_test.py
index 18f3688df..5762f4120 100644
--- a/tests/linalg_test.py
+++ b/tests/linalg_test.py
@@ -538,6 +538,20 @@ class NumpyLinalgTest(jtu.JaxTestCase):
     jac0 = jax.jacobian(np.linalg.solve, argnums=0)(A[0], b[0])
     jac1 = jax.jacobian(np.linalg.solve, argnums=1)(A[0], b[0])
 
+  def testIssue1383(self):
+    seed = jax.random.PRNGKey(0)
+    tmp = jax.random.uniform(seed, (2,2))
+    a = np.dot(tmp, tmp.T)
+
+    def f(inp):
+      val, vec = np.linalg.eigh(inp)
+      return np.dot(np.dot(vec, inp), vec.T)
+
+    grad_func = jax.jacfwd(f)
+    hess_func = jax.jacfwd(grad_func)
+    cube_func = jax.jacfwd(hess_func)
+    self.assertFalse(onp.any(onp.isnan(cube_func(a))))
+
 
 class ScipyLinalgTest(jtu.JaxTestCase):
 
['tests/linalg_test.py', 'jax/lax/lax.py']