"""
well described here http://www.netlib.org/lapack/lug/node123.html
"""
from __future__ import division, print_function
import numpy as np
import numba as nb
[docs]@nb.jit(nopython=True)
def pack2den_u(pack, dtype=np.float64):
"""
Unpacks a packed format to dense format
"""
dim = size2dim(len(pack))
den = np.zeros((dim,dim), dtype=dtype)
for i in range(dim):
for j in range(i+1):
den[i,j] = pack[ij2pack_u(i,j)]
den[j,i] = den[i,j]
return den
[docs]@nb.jit(nopython=True)
def pack2den_l(pack, dtype=np.float64):
"""
Unpacks a packed format to dense format
"""
dim = size2dim(len(pack))
den = np.zeros((dim,dim), dtype=dtype)
for j in range(dim):
for i in range(j,dim):
den[j,i] = pack[ij2pack_l(i,j,dim)]
den[i,j] = den[j,i]
return den
[docs]@nb.jit(nopython=True)
def size2dim(pack_size):
"""
Computes dimension of the matrix by it's packed size. Packed size is n*(n+1)//2
"""
rdim = (np.sqrt(8.0*pack_size+1.0) - 1.0)/2.0
ndim = int(rdim)
if ndim != rdim : raise SystemError('!ndim/=rdim')
if ndim*(ndim+1)//2 != pack_size: SystemError('!pack_size')
return ndim
[docs]@nb.jit(nopython=True)
def ij2pack_u(i,j):
ma = max(i,j)
return ma*(ma+1)//2+min(i,j)
[docs]@nb.jit(nopython=True)
def ij2pack_l(i,j,dim):
mi = min(i,j)
return max(i,j)+((2*dim-mi-1)*mi)//2
[docs]@nb.jit(nopython=True)
def cp_block_pack_u(bmat, s1, f1, s2, f2, gmat, add=False):
if add:
for i in range(s1,f1):
for j in range(s2,f2):
if j>i:
continue
ma = max(i,j)
gmat[ma*(ma+1)//2+min(i,j)] += bmat[i-s1, j-s2]
else:
for i in range(s1,f1):
for j in range(s2,f2):
ma = max(i,j)
gmat[ma*(ma+1)//2+min(i,j)] = bmat[i-s1, j-s2]