TL: finished numpy example

examples/numpy/README.md

@@ -9,6 +9,13 @@ In particular, this script uses the [loadtxt](https://docs.scipy.org/doc/numpy-1
 and [savetxt](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.savetxt.html) numpy functions, to read the data stored
 in the [data.txt](data.txt) file, modify those and write the modified data into [data_double.txt](data_double.txt).
 
+## [callFortranRoutine.py](callFortranRoutine.py)
+
+Perform the [Rutishauser-Kahan-Pal-Walker algorithm](https://link.springer.com/content/pdf/10.1007/BF01405565.pdf) to transform an arrow matrix (with non-zero elements
+in the first row and column, and diagonal) into a symmetric trifdiagonal matrix.
+It uses the [f2py](https://www.numfys.net/howto/F2PY/) sub-module of numpy to compile a fortran subroutine that can be called in python. 
+
+
 ## [csv-read-interpol-plot.py](csv-read-interpol-plot.py)
 
 Read noisy data that are stored in a .csv file (comma separated values, like Excel files, ...), and build a polynomial that can fit the data.

examples/numpy/callFortranRoutine.py

+# -*- coding: utf-8 -*-
+"""
+Perform the Rutishauser-Kahan-Pal-Walker algorithm, using pure python
+and pure fortran.
+
+This algorithm allows to transform a matrix of the form
+
+[    1    sqrt(w1) sqrt(w2) ... sqrt(wn)]
+[sqrt(w1)   tau1                        ]
+[sqrt(w2)            tau2               ]
+[  :                        ...         ]
+[sqrt(wn)                         taun  ]
+
+into the tridiagonal matrix
+
+[    1    sqrt(b0)                                 ]
+[sqrt(b0)    a0    sqrt(b1)                        ]
+[         sqrt(b1)    a1     sqrt(b2)              ]
+[                    ...        ...        ...     ]
+[                            sqrt(b{n-1})   a{n-1} ]
+
+using givens rotations.
+"""
+import os
+import subprocess
+import numpy as np
+from time import time
+
+
+# Define Python function
+def tridiag_rpkw_python(nodes, weights):
+
+    n = nodes.size
+
+    alpha = np.array(nodes)
+    beta = np.zeros(n)
+    beta[0] = weights[0]
+
+    for i in range(n-1):
+
+        pi2 = weights[i+1]
+        lam = nodes[i+1]
+
+        gamma2 = 1.
+        sigma2 = 0.
+        tau = 0.
+
+        for k in range(i+2):
+
+            oldBeta = beta[k]
+            rho2 = oldBeta + pi2
+            beta[k] = gamma2*rho2
+            oldSigma2 = sigma2
+
+            if rho2 == 0:
+                gamma2 = 1.
+                sigma2 = 0.
+            else:
+                gamma2 = oldBeta/rho2
+                sigma2 = pi2/rho2
+
+            newTau = sigma2*(alpha[k]-lam) - gamma2*tau
+            alpha[k] = alpha[k] - (newTau-tau)
+            tau = newTau
+
+            if sigma2 == 0:
+                pi2 = oldSigma2*beta[k]
+            else:
+                pi2 = tau**2/sigma2
+
+    return alpha, beta
+
+
+# Compile fortran function
+print('Compiling fortran library')
+cmd = 'f2py -c -m libfor libfor.f95'
+with open(os.devnull, 'wb') as devnull:
+    subprocess.check_call(
+        cmd.split(), stdout=devnull, stderr=subprocess.STDOUT,
+        env=os.environ.copy())
+print(' -- done')
+
+
+# Define fortran function
+def tridiag_rpkw_fortran(nodes, weights):
+
+    # Import fortran library
+    import libfor
+
+    # Prepare input and output arrays
+    nodes = np.asarray(nodes)
+    weights = np.asarray(weights)
+    n = nodes.size
+    alpha = np.empty(n)
+    beta = np.empty(n)
+
+    # Call fortran library
+    libfor.tridiag_rpkw(nodes, weights, alpha, beta, n)
+
+    return alpha, beta
+
+
+# Performance test
+def runPerfo(nPoints):
+    print('Running performance test with a {0}x{0} matrix'.format(nPoints+1))
+
+    # Create nodes and weights
+    weights = np.full(nPoints, 1./nPoints)
+    nodes = np.linspace(-1, 1, num=nPoints, endpoint=False) + 1./nPoints/2.
+
+    # Python run
+    tBeg = time()
+    alpha, beta = tridiag_rpkw_python(nodes, weights)
+    tEnd = time()
+    tPython = tEnd - tBeg
+    print(' -- python routine: {:f}s'.format(tPython))
+
+    # Fortran run
+    tBeg = time()
+    alpha, beta = tridiag_rpkw_fortran(nodes, weights)
+    tEnd = time()
+    tFortran = tEnd - tBeg
+    print(' -- fortran function: {:f}s'.format(tFortran))
+    print(' -- speedUp factor: {:f}'.format(tPython/tFortran))
+
+
+runPerfo(100)
+runPerfo(500)
+runPerfo(1000)

examples/numpy/libfor.f95

+subroutine tridiag_rpkw(nodes, weights, alpha, beta, n)
+   implicit none
+
+   ! Argument variables
+   integer, intent(in) :: n
+   double precision, intent(in) :: nodes(n), weights(n)
+   double precision, intent(inout) :: alpha(n), beta(n)
+
+   ! Auxilliary variables
+   double precision :: pi2, lam, gamma2, sigma2, tau, rho2
+   double precision :: oldBeta, oldSigma2, newTau
+
+   ! Loop indexes
+   integer :: i, k
+
+   ! Intialize alpha, beta
+   do i=1,n
+      alpha(i) = nodes(i)
+      beta(i) = 0.
+   end do
+   beta(1) = weights(1)
+
+   ! Loop on each submatrices of the arrow matrix
+   do i=1,n-1
+
+      pi2 = weights(i+1)
+      lam = nodes(i+1)
+
+      gamma2 = 1.
+      sigma2 = 1.
+      tau = 0.
+
+      ! Apply rotations successively to tridiagonalize the submatrix
+      do k=1,i+1
+
+         oldBeta = beta(k)
+         rho2 = oldBeta + pi2
+         beta(k) = gamma2*rho2
+         oldSigma2 = sigma2
+
+         if (rho2 <= 0.) then
+            gamma2 = 1.
+            sigma2 = 0.
+         else
+            gamma2 = oldBeta/rho2
+            sigma2 = pi2/rho2
+         end if
+
+         newTau = sigma2*(alpha(k)-lam) - gamma2*tau
+         alpha(k) = alpha(k) - (newTau-tau)
+         tau = newTau
+
+         if (sigma2 <= 0.) then
+            pi2 = oldSigma2*beta(k)
+         else
+            pi2 = tau**2/sigma2
+         end if
+
+      end do
+
+   end do
+
+end subroutine
\ No newline at end of file