code: def tdot_blas(mat, out=None):
      
      if (mat.dtype != 'float64') or (len(mat.shape) != 2):
          return np.dot(mat, mat.T)
      nn = mat.shape[0]
      if out is None:
          out = np.zeros((nn, nn))
      else:
          assert(out.dtype == 'float64')
          assert(out.shape == (nn, nn))
          # FIXME: should allow non-contiguous out, and copy output into it:
          assert(8 in out.strides)
          # zeroing needed because of dumb way I copy across triangular answer
          out[:] = 0.0
  
      # # Call to DSYRK from BLAS
      mat = np.asfortranarray(mat)
      out = blas.dsyrk(alpha=1.0, a=mat, beta=0.0, c=out, overwrite_c=1,
                       trans=0, lower=0)
  
      symmetrify(out, upper=True)
      return np.ascontiguousarray(out)
docs: returns np.dot(mat, mat.T), but faster for large 2D arrays of doubles.