#!/usr/bin/env python

# Matrix Singular Value Decomposition on CPU using numpy and scipy

import numpy as np
from scipy import dot as sp_dot
from scipy import linalg as sp_linalg
from six import print_      # for python 2 and 3 compatibility
import time

col = 2048
row = col + 1
thres = 25
demo_types = [np.float32, np.float64, np.complex64, np.complex128]

for t in demo_types:
    a = np.arange(row*col).reshape(row, col).astype(t)
    
    # Matrix SVD with numpy
    print_('\n\n\nNumpy svd for type ' + str(np.dtype(t)))
    t_0 = time.time()
    # decomposition
    u, s, vh = np.linalg.svd(a)
    t_1 = time.time()
    # reconstruction
    S = np.zeros((row, col)).astype(t)
    S[:min(row, col),:min(row, col)] = np.diag(s)
    a_rec = np.dot(u[:,:thres], np.dot(S[:thres,:thres], vh[:thres,:]))
    t_2 = time.time()
    print_('Decomposition time:\t\t{0:8.2f}s'.format(t_1 - t_0))
    print_('Resconstruction time:\t\t{0:8.2f}s'.format(t_2 - t_1))
    print_('Maximum relative error:\t\t  {0:10.2e}'\
        .format(np.max(np.abs(a-a_rec))/np.max(np.abs(a))))
    
    
    # Matrix SVD with scipy
    print_('\nScipy svd for type ' + str(np.dtype(t)))
    t_0 = time.time()
    # decomposition
    u, s, vh = sp_linalg.svd(a)
    t_1 = time.time()
    # reconstruction
    S = sp_linalg.diagsvd(s, row, col)
    a_rec = sp_dot(u[:,:thres], sp_dot(S[:thres,:thres], vh[:thres,:]))
    t_2 = time.time()
    print_('Decomposition time:\t\t{0:8.2f}s'.format(t_1 - t_0))
    print_('Resconstruction time:\t\t{0:8.2f}s'.format(t_2 - t_1))
    print_('Maximum relative error:\t\t  {0:10.2e}'\
        .format(np.max(np.abs(a-a_rec))/np.max(np.abs(a))))