In [1]:

    

%load_ext Cython


    

In [2]:

    

from math import sqrt
from array import array

def A(i, j):
    """i: row, j: column, ZERO-based."""
    return 1.0 / (((i + j) * (i + j + 1) >> 1) + i + 1)

def A_times_u(u, v):
    """v = Au, A: nxn Matrix, u: n vector."""
    u_len = len(u)
    
    for i in range(u_len):
        partial_sum = 0.0
        for j in range(u_len):
            partial_sum += A(i, j) * u[j]
            
        v[i] = partial_sum
        
def At_times_u(u, v):
    """v = A^T u, A: nxn Matrix, u: n vector."""
    u_len = len(u)
    
    for i in range(u_len):
        partial_sum = 0.0
        for j in range(u_len):
            partial_sum += A(j, i) * u[j]
            
        v[i] = partial_sum
        
def B_times_u(u, out, tmp):
    """Bu = A^T A u, tmp = A u, out = A^T tmp."""
    A_times_u(u, tmp)
    At_times_u(tmp, out)
    
def spectral_norm(n):
    u = array("d", [1.0] * n)
    v = array("d", [0.0] * n)
    tmp = array("d", [0.0] * n)
    
    # B^n u, set n=20 to adapt u to 
    # closely aligned with the principal eigenvector.
    for _ in range(10):
        B_times_u(u, v, tmp)
        B_times_u(v, u, tmp)
        
    vBv, vv = 0.0, 0.0
    
    for ue, ve in zip(u, v):
        vBv += ue * ve
        vv += ve * ve
    
    return sqrt(vBv / vv)

def main(n):
    print("%0.9f" % spectral_norm(n))


    

In [3]:

    

import cProfile as pf


    

In [4]:

    

pf.run('main(400)', sort='time')


    

    




1.274224081
         6400126 function calls in 3.533 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
  6400000    2.136    0.000    2.136    0.000 <ipython-input-2-f499672ca18f>:4(A)
       20    0.700    0.035    1.770    0.089 <ipython-input-2-f499672ca18f>:8(A_times_u)
       20    0.696    0.035    1.762    0.088 <ipython-input-2-f499672ca18f>:19(At_times_u)
        1    0.000    0.000    3.533    3.533 <ipython-input-2-f499672ca18f>:35(spectral_norm)
       20    0.000    0.000    3.533    0.177 <ipython-input-2-f499672ca18f>:30(B_times_u)
        1    0.000    0.000    3.533    3.533 {built-in method builtins.exec}
       40    0.000    0.000    0.000    0.000 {built-in method builtins.len}
        2    0.000    0.000    0.000    0.000 iostream.py:309(write)
        1    0.000    0.000    3.533    3.533 <ipython-input-2-f499672ca18f>:54(main)
        1    0.000    0.000    0.000    0.000 ioloop.py:932(add_callback)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.print}
        2    0.000    0.000    0.000    0.000 iostream.py:241(_schedule_flush)
        1    0.000    0.000    0.000    0.000 {method 'write' of '_io.FileIO' objects}
        1    0.000    0.000    0.000    0.000 stack_context.py:253(wrap)
        2    0.000    0.000    0.000    0.000 iostream.py:228(_is_master_process)
        1    0.000    0.000    0.000    0.000 posix.py:53(wake)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.hasattr}
        1    0.000    0.000    0.000    0.000 {built-in method math.sqrt}
        2    0.000    0.000    0.000    0.000 {method 'write' of '_io.StringIO' objects}
        2    0.000    0.000    0.000    0.000 {built-in method builtins.isinstance}
        2    0.000    0.000    0.000    0.000 {built-in method posix.getpid}
        1    0.000    0.000    0.000    0.000 {method 'append' of 'list' objects}
        1    0.000    0.000    0.000    0.000 {built-in method _thread.get_ident}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
        1    0.000    0.000    3.533    3.533 <string>:1(<module>)

In [5]:

    

%%cython
from math import sqrt
from array import array

def A0(i, j):
    """i: row, j: column, ZERO-based."""
    return 1.0 / (((i + j) * (i + j + 1) >> 1) + i + 1)

def A_times_u0(u, v):
    """v = Au, A: nxn Matrix, u: n vector."""
    u_len = len(u)
    
    for i in range(u_len):
        partial_sum = 0.0
        for j in range(u_len):
            partial_sum += A0(i, j) * u[j]
            
        v[i] = partial_sum
        
def At_times_u0(u, v):
    """v = A^T u, A: nxn Matrix, u: n vector."""
    u_len = len(u)
    
    for i in range(u_len):
        partial_sum = 0.0
        for j in range(u_len):
            partial_sum += A0(j, i) * u[j]
            
        v[i] = partial_sum
        
def B_times_u0(u, out, tmp):
    """Bu = A^T A u, tmp = A u, out = A^T tmp."""
    A_times_u0(u, tmp)
    At_times_u0(tmp, out)
    
def spectral_norm0(n):
    u = array("d", [1.0] * n)
    v = array("d", [0.0] * n)
    tmp = array("d", [0.0] * n)
    
    # B^n u, set n=20 to adapt u to 
    # closely aligned with the principal eigenvector.
    for _ in range(10):
        B_times_u0(u, v, tmp)
        B_times_u0(v, u, tmp)
        
    vBv, vv = 0.0, 0.0
    
    for ue, ve in zip(u, v):
        vBv += ue * ve
        vv += ve * ve
    
    return sqrt(vBv / vv)

def main0(n):
    print("%0.9f" % spectral_norm0(n))


    

In [6]:

    

pf.run("main0(400)", sort='time')


    

    




1.274224081
         21 function calls in 1.483 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    1.483    1.483    1.483    1.483 {_cython_magic_f41508800d19f149a51a027ed7e2c3da.main0}
        2    0.000    0.000    0.000    0.000 iostream.py:309(write)
        1    0.000    0.000    1.483    1.483 {built-in method builtins.exec}
        2    0.000    0.000    0.000    0.000 iostream.py:228(_is_master_process)
        1    0.000    0.000    0.000    0.000 ioloop.py:932(add_callback)
        1    0.000    0.000    0.000    0.000 stack_context.py:253(wrap)
        2    0.000    0.000    0.000    0.000 iostream.py:241(_schedule_flush)
        1    0.000    0.000    1.483    1.483 <string>:1(<module>)
        2    0.000    0.000    0.000    0.000 {built-in method posix.getpid}
        2    0.000    0.000    0.000    0.000 {method 'write' of '_io.StringIO' objects}
        2    0.000    0.000    0.000    0.000 {built-in method builtins.isinstance}
        1    0.000    0.000    0.000    0.000 {built-in method _thread.get_ident}
        1    0.000    0.000    0.000    0.000 {built-in method builtins.hasattr}
        1    0.000    0.000    0.000    0.000 {method 'append' of 'list' objects}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

In [7]:

    

%%cython -a
# cython: cdivison=True, boundscheck=False, wrapround=False

from libc.math cimport sqrt
import numpy as np
#from cpython.array cimport array

cdef inline double A1(int i,int j):
    """i: row, j: column, ZERO-based."""
    return 1.0 / (((i + j) * (i + j + 1) >> 1) + i + 1)

cdef void A_times_u1(double[::1] u, double[::1] v):
    """v = Au, A: nxn Matrix, u: n vector."""
    cdef:
        int u_len = u.shape[0]
        int i
        double partial_sum
    
    for i in range(u_len):
        partial_sum = 0.0
        for j in range(u_len):
            partial_sum += A1(i, j) * u[j]
            
        v[i] = partial_sum
        
cdef void At_times_u1(double[::1] u, double[::1] v):
    """v = A^T u, A: nxn Matrix, u: n vector."""
    cdef:
        int u_len = u.shape[0]
        int i
        double partial_sum
    
    for i in range(u_len):
        partial_sum = 0.0
        for j in range(u_len):
            partial_sum += A1(j, i) * u[j]
            
        v[i] = partial_sum
        
cdef inline void B_times_u1(double[::1] u, double[::1] out, double[::1] tmp):
    """Bu = A^T A u, tmp = A u, out = A^T tmp."""
    A_times_u1(u, tmp)
    At_times_u1(tmp, out)
    
cdef double spectral_norm1(int n):
    cdef:
        double[::1] u = np.ones((n,), dtype=np.double)
        double[::1] v = np.empty((n,), dtype=np.double)
        double[::1] tmp = np.empty((n,), dtype=np.double)
        #double[::1] u = array('d', [1.0]*n)
        #double[::1] v = array('d', [0.0]*n)
        #double[::1] tmp = array('d', [0.0]*n)
        int i, j
        double vBv=0.0, vv=0.0#, ue, ve
    
    # B^n u, set n=20 to adapt u to 
    # closely aligned with the principal eigenvector.
    for i in range(10):
        B_times_u1(u, v, tmp)
        B_times_u1(v, u, tmp)
        
    #for ue, ve in zip(u, v):
    #    vBv += ue * ve
    #    vv += ve * ve
    for j in range(n):
        vBv += u[j] * v[j]
        vv += v[j] * v[j]
    
    return sqrt(vBv / vv)

cpdef void main1(n):
    print("%0.9f" % spectral_norm1(n))


    

    
Out[7]:









    
    
    
    


Generated by Cython 0.25.2

    Yellow lines hint at Python interaction.

    Click on a line that starts with a "+" to see the C code that Cython generated for it.

+01: # cython: cdivison=True, boundscheck=False, wrapround=False
  __pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 02: 
 03: from libc.math cimport sqrt
+04: import numpy as np
  __pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 4, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) __PYX_ERR(0, 4, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 05: #from cpython.array cimport array
 06: 
+07: cdef inline double A1(int i,int j):
static CYTHON_INLINE double __pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_A1(int __pyx_v_i, int __pyx_v_j) {
  double __pyx_r;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("A1", 0);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_WriteUnraisable("_cython_magic_11ff13a900cc2fede1824e43f31ba3b4.A1", __pyx_clineno, __pyx_lineno, __pyx_filename, 0, 0);
  __pyx_r = 0;
  __pyx_L0:;
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
 08:     """i: row, j: column, ZERO-based."""
+09:     return 1.0 / (((i + j) * (i + j + 1) >> 1) + i + 1)
  __pyx_t_1 = (((((__pyx_v_i + __pyx_v_j) * ((__pyx_v_i + __pyx_v_j) + 1)) >> 1) + __pyx_v_i) + 1);
  if (unlikely(__pyx_t_1 == 0)) {
    PyErr_SetString(PyExc_ZeroDivisionError, "float division");
    __PYX_ERR(0, 9, __pyx_L1_error)
  }
  __pyx_r = (1.0 / ((double)__pyx_t_1));
  goto __pyx_L0;
 10: 
+11: cdef void A_times_u1(double[::1] u, double[::1] v):
static void __pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_A_times_u1(__Pyx_memviewslice __pyx_v_u, __Pyx_memviewslice __pyx_v_v) {
  int __pyx_v_u_len;
  int __pyx_v_i;
  double __pyx_v_partial_sum;
  int __pyx_v_j;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("A_times_u1", 0);
/* … */
  /* function exit code */
  __Pyx_RefNannyFinishContext();
}
 12:     """v = Au, A: nxn Matrix, u: n vector."""
 13:     cdef:
+14:         int u_len = u.shape[0]
  __pyx_v_u_len = (__pyx_v_u.shape[0]);
 15:         int i
 16:         double partial_sum
 17: 
+18:     for i in range(u_len):
  __pyx_t_1 = __pyx_v_u_len;
  for (__pyx_t_2 = 0; __pyx_t_2 < __pyx_t_1; __pyx_t_2+=1) {
    __pyx_v_i = __pyx_t_2;
+19:         partial_sum = 0.0
    __pyx_v_partial_sum = 0.0;
+20:         for j in range(u_len):
    __pyx_t_3 = __pyx_v_u_len;
    for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) {
      __pyx_v_j = __pyx_t_4;
+21:             partial_sum += A1(i, j) * u[j]
      __pyx_t_5 = __pyx_v_j;
      if (__pyx_t_5 < 0) __pyx_t_5 += __pyx_v_u.shape[0];
      __pyx_v_partial_sum = (__pyx_v_partial_sum + (__pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_A1(__pyx_v_i, __pyx_v_j) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_u.data) + __pyx_t_5)) )))));
    }
 22: 
+23:         v[i] = partial_sum
    __pyx_t_6 = __pyx_v_i;
    if (__pyx_t_6 < 0) __pyx_t_6 += __pyx_v_v.shape[0];
    *((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_6)) )) = __pyx_v_partial_sum;
  }
 24: 
+25: cdef void At_times_u1(double[::1] u, double[::1] v):
static void __pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_At_times_u1(__Pyx_memviewslice __pyx_v_u, __Pyx_memviewslice __pyx_v_v) {
  int __pyx_v_u_len;
  int __pyx_v_i;
  double __pyx_v_partial_sum;
  int __pyx_v_j;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("At_times_u1", 0);
/* … */
  /* function exit code */
  __Pyx_RefNannyFinishContext();
}
 26:     """v = A^T u, A: nxn Matrix, u: n vector."""
 27:     cdef:
+28:         int u_len = u.shape[0]
  __pyx_v_u_len = (__pyx_v_u.shape[0]);
 29:         int i
 30:         double partial_sum
 31: 
+32:     for i in range(u_len):
  __pyx_t_1 = __pyx_v_u_len;
  for (__pyx_t_2 = 0; __pyx_t_2 < __pyx_t_1; __pyx_t_2+=1) {
    __pyx_v_i = __pyx_t_2;
+33:         partial_sum = 0.0
    __pyx_v_partial_sum = 0.0;
+34:         for j in range(u_len):
    __pyx_t_3 = __pyx_v_u_len;
    for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) {
      __pyx_v_j = __pyx_t_4;
+35:             partial_sum += A1(j, i) * u[j]
      __pyx_t_5 = __pyx_v_j;
      if (__pyx_t_5 < 0) __pyx_t_5 += __pyx_v_u.shape[0];
      __pyx_v_partial_sum = (__pyx_v_partial_sum + (__pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_A1(__pyx_v_j, __pyx_v_i) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_u.data) + __pyx_t_5)) )))));
    }
 36: 
+37:         v[i] = partial_sum
    __pyx_t_6 = __pyx_v_i;
    if (__pyx_t_6 < 0) __pyx_t_6 += __pyx_v_v.shape[0];
    *((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_6)) )) = __pyx_v_partial_sum;
  }
 38: 
+39: cdef inline void B_times_u1(double[::1] u, double[::1] out, double[::1] tmp):
static CYTHON_INLINE void __pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_B_times_u1(__Pyx_memviewslice __pyx_v_u, __Pyx_memviewslice __pyx_v_out, __Pyx_memviewslice __pyx_v_tmp) {
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("B_times_u1", 0);
/* … */
  /* function exit code */
  __Pyx_RefNannyFinishContext();
}
 40:     """Bu = A^T A u, tmp = A u, out = A^T tmp."""
+41:     A_times_u1(u, tmp)
  __pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_A_times_u1(__pyx_v_u, __pyx_v_tmp);
+42:     At_times_u1(tmp, out)
  __pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_At_times_u1(__pyx_v_tmp, __pyx_v_out);
 43: 
+44: cdef double spectral_norm1(int n):
static double __pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_spectral_norm1(int __pyx_v_n) {
  __Pyx_memviewslice __pyx_v_u = { 0, 0, { 0 }, { 0 }, { 0 } };
  __Pyx_memviewslice __pyx_v_v = { 0, 0, { 0 }, { 0 }, { 0 } };
  __Pyx_memviewslice __pyx_v_tmp = { 0, 0, { 0 }, { 0 }, { 0 } };
  CYTHON_UNUSED int __pyx_v_i;
  int __pyx_v_j;
  double __pyx_v_vBv;
  double __pyx_v_vv;
  double __pyx_r;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("spectral_norm1", 0);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4);
  __Pyx_XDECREF(__pyx_t_5);
  __PYX_XDEC_MEMVIEW(&__pyx_t_6, 1);
  __Pyx_WriteUnraisable("_cython_magic_11ff13a900cc2fede1824e43f31ba3b4.spectral_norm1", __pyx_clineno, __pyx_lineno, __pyx_filename, 0, 0);
  __pyx_r = 0;
  __pyx_L0:;
  __PYX_XDEC_MEMVIEW(&__pyx_v_u, 1);
  __PYX_XDEC_MEMVIEW(&__pyx_v_v, 1);
  __PYX_XDEC_MEMVIEW(&__pyx_v_tmp, 1);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
 45:     cdef:
+46:         double[::1] u = np.ones((n,), dtype=np.double)
  __pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_ones); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_3 = PyTuple_New(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_GIVEREF(__pyx_t_1);
  PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_1);
  __pyx_t_1 = 0;
  __pyx_t_1 = PyTuple_New(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_GIVEREF(__pyx_t_3);
  PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_t_3);
  __pyx_t_3 = 0;
  __pyx_t_3 = PyDict_New(); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_double); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  if (PyDict_SetItem(__pyx_t_3, __pyx_n_s_dtype, __pyx_t_5) < 0) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_1, __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_5);
  if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_v_u = __pyx_t_6;
  __pyx_t_6.memview = NULL;
  __pyx_t_6.data = NULL;
+47:         double[::1] v = np.empty((n,), dtype=np.double)
  __pyx_t_5 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 47, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_empty); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 47, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 47, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_t_1 = PyTuple_New(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 47, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_GIVEREF(__pyx_t_5);
  PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_t_5);
  __pyx_t_5 = 0;
  __pyx_t_5 = PyTuple_New(1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 47, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_GIVEREF(__pyx_t_1);
  PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_1);
  __pyx_t_1 = 0;
  __pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 47, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 47, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_double); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 47, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  if (PyDict_SetItem(__pyx_t_1, __pyx_n_s_dtype, __pyx_t_4) < 0) __PYX_ERR(0, 47, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_5, __pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 47, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_4);
  if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 47, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_v_v = __pyx_t_6;
  __pyx_t_6.memview = NULL;
  __pyx_t_6.data = NULL;
+48:         double[::1] tmp = np.empty((n,), dtype=np.double)
  __pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_empty); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_5 = PyTuple_New(1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_GIVEREF(__pyx_t_4);
  PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_4);
  __pyx_t_4 = 0;
  __pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_GIVEREF(__pyx_t_5);
  PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_5);
  __pyx_t_5 = 0;
  __pyx_t_5 = PyDict_New(); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_double); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  if (PyDict_SetItem(__pyx_t_5, __pyx_n_s_dtype, __pyx_t_2) < 0) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = __Pyx_PyObject_Call(__pyx_t_1, __pyx_t_4, __pyx_t_5); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_t_2);
  if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_v_tmp = __pyx_t_6;
  __pyx_t_6.memview = NULL;
  __pyx_t_6.data = NULL;
 49:         #double[::1] u = array('d', [1.0]*n)
 50:         #double[::1] v = array('d', [0.0]*n)
 51:         #double[::1] tmp = array('d', [0.0]*n)
 52:         int i, j
+53:         double vBv=0.0, vv=0.0#, ue, ve
  __pyx_v_vBv = 0.0;
  __pyx_v_vv = 0.0;
 54: 
 55:     # B^n u, set n=20 to adapt u to 
 56:     # closely aligned with the principal eigenvector.
+57:     for i in range(10):
  for (__pyx_t_7 = 0; __pyx_t_7 < 10; __pyx_t_7+=1) {
    __pyx_v_i = __pyx_t_7;
+58:         B_times_u1(u, v, tmp)
    __pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_B_times_u1(__pyx_v_u, __pyx_v_v, __pyx_v_tmp);
+59:         B_times_u1(v, u, tmp)
    __pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_B_times_u1(__pyx_v_v, __pyx_v_u, __pyx_v_tmp);
  }
 60: 
 61:     #for ue, ve in zip(u, v):
 62:     #    vBv += ue * ve
 63:     #    vv += ve * ve
+64:     for j in range(n):
  __pyx_t_7 = __pyx_v_n;
  for (__pyx_t_8 = 0; __pyx_t_8 < __pyx_t_7; __pyx_t_8+=1) {
    __pyx_v_j = __pyx_t_8;
+65:         vBv += u[j] * v[j]
    __pyx_t_9 = __pyx_v_j;
    if (__pyx_t_9 < 0) __pyx_t_9 += __pyx_v_u.shape[0];
    __pyx_t_10 = __pyx_v_j;
    if (__pyx_t_10 < 0) __pyx_t_10 += __pyx_v_v.shape[0];
    __pyx_v_vBv = (__pyx_v_vBv + ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_u.data) + __pyx_t_9)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_10)) )))));
+66:         vv += v[j] * v[j]
    __pyx_t_11 = __pyx_v_j;
    if (__pyx_t_11 < 0) __pyx_t_11 += __pyx_v_v.shape[0];
    __pyx_t_12 = __pyx_v_j;
    if (__pyx_t_12 < 0) __pyx_t_12 += __pyx_v_v.shape[0];
    __pyx_v_vv = (__pyx_v_vv + ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_11)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_v.data) + __pyx_t_12)) )))));
  }
 67: 
+68:     return sqrt(vBv / vv)
  if (unlikely(__pyx_v_vv == 0)) {
    PyErr_SetString(PyExc_ZeroDivisionError, "float division");
    __PYX_ERR(0, 68, __pyx_L1_error)
  }
  __pyx_r = sqrt((__pyx_v_vBv / __pyx_v_vv));
  goto __pyx_L0;
 69: 
+70: cpdef void main1(n):
static PyObject *__pyx_pw_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_1main1(PyObject *__pyx_self, PyObject *__pyx_v_n); /*proto*/
static void __pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_main1(PyObject *__pyx_v_n, CYTHON_UNUSED int __pyx_skip_dispatch) {
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("main1", 0);
/* … */
  /* function exit code */
  goto __pyx_L0;
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_WriteUnraisable("_cython_magic_11ff13a900cc2fede1824e43f31ba3b4.main1", __pyx_clineno, __pyx_lineno, __pyx_filename, 0, 0);
  __pyx_L0:;
  __Pyx_RefNannyFinishContext();
}

/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_1main1(PyObject *__pyx_self, PyObject *__pyx_v_n); /*proto*/
static PyObject *__pyx_pw_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_1main1(PyObject *__pyx_self, PyObject *__pyx_v_n) {
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("main1 (wrapper)", 0);
  __pyx_r = __pyx_pf_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_main1(__pyx_self, ((PyObject *)__pyx_v_n));

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_main1(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_n) {
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("main1", 0);
  __Pyx_XDECREF(__pyx_r);
  __pyx_t_1 = __Pyx_void_to_None(__pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_main1(__pyx_v_n, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 70, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_r = __pyx_t_1;
  __pyx_t_1 = 0;
  goto __pyx_L0;

  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_AddTraceback("_cython_magic_11ff13a900cc2fede1824e43f31ba3b4.main1", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
+71:     print("%0.9f" % spectral_norm1(n))
  __pyx_t_1 = __Pyx_PyInt_As_int(__pyx_v_n); if (unlikely((__pyx_t_1 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 71, __pyx_L1_error)
  __pyx_t_2 = PyFloat_FromDouble(__pyx_f_46_cython_magic_11ff13a900cc2fede1824e43f31ba3b4_spectral_norm1(__pyx_t_1)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 71, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = PyUnicode_Format(__pyx_kp_u_0_9f, __pyx_t_2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 71, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = PyTuple_New(1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 71, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_GIVEREF(__pyx_t_3);
  PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_t_3);
  __pyx_t_3 = 0;
  __pyx_t_3 = __Pyx_PyObject_Call(__pyx_builtin_print, __pyx_t_2, NULL); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 71, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;

In [9]:

    

pf.run("main1(400)", sort='time')


    

    




1.274224081
         24 function calls in 0.011 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.011    0.011    0.011    0.011 {built-in method _cython_magic_11ff13a900cc2fede1824e43f31ba3b4.main1}
        2    0.000    0.000    0.000    0.000 iostream.py:309(write)
        1    0.000    0.000    0.011    0.011 {built-in method builtins.exec}
        1    0.000    0.000    0.000    0.000 ioloop.py:932(add_callback)
        1    0.000    0.000    0.000    0.000 {built-in method numpy.core.multiarray.empty}
        1    0.000    0.000    0.000    0.000 {built-in method numpy.core.multiarray.copyto}
        1    0.000    0.000    0.000    0.000 stack_context.py:253(wrap)
        2    0.000    0.000    0.000    0.000 iostream.py:241(_schedule_flush)
        1    0.000    0.000    0.000    0.000 numeric.py:150(ones)
        2    0.000    0.000    0.000    0.000 iostream.py:228(_is_master_process)
        1    0.000    0.000    0.011    0.011 <string>:1(<module>)
        2    0.000    0.000    0.000    0.000 {method 'write' of '_io.StringIO' objects}
        1    0.000    0.000    0.000    0.000 {built-in method builtins.hasattr}
        2    0.000    0.000    0.000    0.000 {built-in method posix.getpid}
        1    0.000    0.000    0.000    0.000 {method 'append' of 'list' objects}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
        2    0.000    0.000    0.000    0.000 {built-in method builtins.isinstance}
        1    0.000    0.000    0.000    0.000 {built-in method _thread.get_ident}

In [11]:

    

pf.run("main1(5500)", sort='time')


    

    




1.274224153
         26 function calls in 1.917 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    1.917    1.917    1.917    1.917 {built-in method _cython_magic_11ff13a900cc2fede1824e43f31ba3b4.main1}
        2    0.000    0.000    0.000    0.000 iostream.py:309(write)
        1    0.000    0.000    1.917    1.917 {built-in method builtins.exec}
        1    0.000    0.000    0.000    0.000 {built-in method numpy.core.multiarray.copyto}
        1    0.000    0.000    0.000    0.000 ioloop.py:932(add_callback)
        1    0.000    0.000    0.000    0.000 {built-in method numpy.core.multiarray.empty}
        1    0.000    0.000    0.000    0.000 {method 'write' of '_io.FileIO' objects}
        1    0.000    0.000    0.000    0.000 numeric.py:150(ones)
        2    0.000    0.000    0.000    0.000 iostream.py:241(_schedule_flush)
        2    0.000    0.000    0.000    0.000 iostream.py:228(_is_master_process)
        1    0.000    0.000    0.000    0.000 stack_context.py:253(wrap)
        1    0.000    0.000    0.000    0.000 posix.py:53(wake)
        1    0.000    0.000    1.917    1.917 <string>:1(<module>)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.hasattr}
        2    0.000    0.000    0.000    0.000 {method 'write' of '_io.StringIO' objects}
        1    0.000    0.000    0.000    0.000 {method 'append' of 'list' objects}
        2    0.000    0.000    0.000    0.000 {built-in method posix.getpid}
        1    0.000    0.000    0.000    0.000 {built-in method _thread.get_ident}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
        2    0.000    0.000    0.000    0.000 {built-in method builtins.isinstance}

In [ ]: