def __init__(self,dim,dicto=None,shifts=None):
        """
        attributes: shifts, dims, data, is1D, length
        methods : dump
        """

        self.shifts=shifts
        
        if type(dim)==type(()) or type(dim)==type([]):
            self.data = spmatrix.ll_mat(reduce(operator.mul, dim), 1)
            self.dims = dim        
            if dicto:
                for k, v in dicto.iteritems():
                    shk = map(operator.__sub__, k, shifts)
                    self.data[self.comp(shk), 0]=v
        
        elif type(dim)==IntType:
            self.data = spmatrix.ll_mat(dim,1)
            self.dims = dim        
            if dicto:
                for k, v in dicto.iteritems():
                    shk = k - shifts
                    self.data[shk,0] = v
        
        self.is1D = type(self.dims)==IntType

    def __init__(self, portfolio):
        from pysparse.sparse import spmatrix
        # set up sparse matrices
        # these first two lists define all indices for asset and issuer arrays
        self.issuers = [i for i in portfolio.issuers()]
        self.assets = [a for a in portfolio.assets]
        self.asset_issuer_map = makeAssetIssuerIndexMap(self.issuers, self.assets)

        def ppfGen(assets):
            for ass in assets:
                yield norm.ppf(ass.dp)

        self.thresholds = np.fromiter(ppfGen(self.assets), np.double)

        self.n_issuers = len(self.issuers)
        self.n_assets = len(self.assets)
        factor_indices = portfolio.factor_indices() # on class for testing help
        self.n_factors = len(factor_indices.keys())
        
        #do the factor weights, also running sum for ideosyncratic weights
        wm = spmatrix.ll_mat(self.n_issuers, self.n_factors+self.n_issuers)
        for i, iss in enumerate(self.issuers):
            wsum = 0.0
            for f in iss.factors:
                j = factor_indices[f.name]
                w = np.sqrt(max(f.weight, 0.0))
                wm[i, j] = w
                wsum += w*w
            wm[i, self.n_factors+i] = np.sqrt(max(1.0 - wsum, 0.0))
            
        self.weights = wm.to_csr()

    def __init__(self, **kwargs):

        nrow = kwargs.get('nrow', 0)
        ncol = kwargs.get('ncol', 0)
        bandwidth = kwargs.get('bandwidth', 0)
        matrix = kwargs.get('matrix', None)
        sizeHint = kwargs.get('sizeHint', 0)
        symmetric = 'symmetric' in kwargs and kwargs['symmetric']
        size = kwargs.get('size',0)
        if size > 0:
            if nrow > 0 or ncol > 0:
                if size != nrow or size != ncol:
                    msg =  'size argument was given but does not match '
                    msg += 'nrow and ncol'
                raise ValueError, msg
            else:
                nrow = ncol = size

        if matrix is not None:
            self.matrix = matrix
        else:
            if symmetric and nrow==ncol:
                if sizeHint is None:
                    sizeHint = nrow
                    if bandwidth > 0:
                        sizeHint += 2*(bandwidth-1)*(2*nrow-bandwidth-2)
                self.matrix = spmatrix.ll_mat_sym(nrow, sizeHint)
            else:
                if sizeHint is None:
                    sizeHint = min(nrow,ncol)
                    if bandwidth > 0:
                        sizeHint = bandwidth * (2*sizeHint-bandwidth-1)/2
                self.matrix = spmatrix.ll_mat(nrow, ncol, sizeHint)

 def test_on_diagonal_matrix(self):
     a = spmatrix.ll_mat(3, 3)
     a.put([1, 2, 3])
     r = jdsym.jdsym(a, None, None, 3, 1, 1e-9, 100, krylov.qmrs)
     self.assertEqual(r[0], 3)
     self.assertTrue(numpy.allclose(r[1], [1, 2, 3]))
     self.assertTrue(numpy.allclose(numpy.abs(r[2]), numpy.identity(3)))

def construct_pysparse_matrix(n, nbr_elements):
    A = spmatrix.ll_mat(n, n, nbr_elements)

    for i in xrange(nbr_elements):
        A[i % n, (2 * i + 1) % n] = i / 3

    return A

    def LinearSystem(self):
        r"""
        Assembly linear system
        Depends on Velocity field and Gamma
        """
        # assembly matrix of linear system
        # using pysparse optimized matrix non zero elements 5*M         
        self.mUt = spmatrix.ll_mat(self.Nz*self.Nx, self.Nz*self.Nx, 5*self.Nz*self.Nx-2*self.Nz-2*self.Nx)

        for Ln in range(0, self.Nz*self.Nx, 1):
            # 1.0*u(x-1,z) + Gamma(x,z)*u(x,z) + 1.0*u(x+1,z) + 1.0*u(x,z-1) + 1.0*u(x,z+1)
            # turn the indices to the one of original matrix
            i = Ln%self.Nx
            k = Ln/self.Nx

            self.mUt[Ln,Ln] = self.Gamma(k, i)
            #is this right?
            if(i-1 >= 0): # u(x-1,z) inside grid in I
                self.mUt[Ln,Ln-1] = 1.0
            if(i+1 < self.Nx): # u(x+1,z) inside grid in I
                self.mUt[Ln,Ln+1] = 1.0
            if(k-1 >= 0): #u(x,z-1)
                self.mUt[Ln,Ln-self.Nx]= 1.0
            if(k+1 < self.Nz): #u(x,z+1)
                self.mUt[Ln,Ln+self.Nx]= 1.0

        
        return self.mUt

    def setUp(self):

        self.nbr_elements = 100000
        self.size = 1000000

        self.A_c = LLSparseMatrix(size=self.size, size_hint=self.nbr_elements, dtype=FLOAT64_T)
        self.A_p = spmatrix.ll_mat(self.size, self.size, self.nbr_elements)
        self.A_s = lil_matrix((self.size, self.size), dtype=np.float64) # how do we reserve space in advance?

def poisson1d_vec(n):
    L = spmatrix.ll_mat(n, n, 3*n-2)
    e = numpy.ones(n)
    d = numpy.arange(n, dtype=numpy.int)
    L.put(2*e, d, d)
    L.put(-e[1:], d[1:], d[:-1])
    L.put(-e[1:], d[:-1], d[1:])
    return L

def poisson1d(n):
    L = spmatrix.ll_mat(n, n, 3*n-2)
    for i in range(n):
        L[i,i] = 2
        if i > 0:
            L[i,i-1] = -1
        if i < n-1:
            L[i,i+1] = -1
    return L

    def setUp(self):

        self.nbr_elements = 10000
        self.size = 100000

        self.A_c = LLSparseMatrix(size=self.size, size_hint=self.nbr_elements, dtype=FLOAT64_T)
        construct_sparse_matrix(self.A_c, self.size, self.nbr_elements)

        self.A_p = spmatrix.ll_mat(self.size, self.size, self.nbr_elements)
        construct_sparse_matrix(self.A_p, self.size, self.nbr_elements)

def lsq(lsq_ff):
    """ 

    :param lsq_ff:
    Convert the LSQP in the First Form(FF) ::
    
           minimize    c'x + 1/2|Qx-d|^2
           subject to  L <= Bx <= U,                       (LSQP-FF)
                       l <=  x <= u,
    to the Second Form (SF):: 
    
            minimize    c'x +1/2|r|^2
            subject to. [d] <= [Q  I][r] <= [d],
                        [L] <= [B  0][x] <= [U],            (LSQP-SF)
                        [l] <=       [x] <= [u],
                     -[inf] <=       [r] <= [inf].
    """
     
     
    p,n = lsq_ff.Q.shape
    m,n = lsq_ff.B.shape
     
    new_B = spmatrix.ll_mat(m+p, n+p, m+n+2*p+lsq_ff.B.nnz+lsq_ff.Q.nnz)
    new_B[:p,:n] = lsq_ff.Q
    new_B[p:,:n] = lsq_ff.B
    
    new_B.put(1, range(p), range(n,n+p))
   
    new_Lcon = np.zeros(p+m)    
    new_Lcon[:p] = lsq_ff.d   
    new_Lcon[p:] = lsq_ff.Lcon
    
    new_Ucon = np.zeros(p+m)    
    new_Ucon[:p] =  lsq_ff.d
    new_Ucon[p:] = lsq_ff.Ucon
    
    new_Lvar = -np.inf * np.ones(n+p)    
    new_Lvar[:n] = lsq_ff.Lvar
    
    new_Uvar = np.inf * np.ones(n+p)
    new_Uvar[:n] = lsq_ff.Uvar   
    
    new_Q = PysparseMatrix(nrow=n+p, ncol=n+p,\
                           sizeHint=p)
    new_Q.put(1, range(n,n+p), range(n,n+p))

    new_d = np.zeros(n+p)
    
    new_c = np.zeros(n+p)
    new_c[:n] = lsq_ff.c
    
    return LSQModel(Q=new_Q, B=new_B, d=new_d, c= new_c, Lcon=new_Lcon, \
                    Ucon=new_Ucon, Lvar=new_Lvar, Uvar=new_Uvar,
                    name= lsq_ff.name, dimQB=(p,n,m))

def poisson2d_vec(n):
    n2 = n*n
    L = spmatrix.ll_mat(n2, n2, 5*n2-4*n)
    d = np.arange(n2, dtype=np.int)
    L.put(4.0, d)
    L.put(-1.0, d[:-n], d[n:])
    L.put(-1.0, d[n:], d[:-n])
    for i in xrange(n):
        di = d[i*n:(i+1)*n]
        L.put(-1.0, di[1:], di[:-1])
        L.put(-1.0, di[:-1], di[1:])
    return L

 def testCompressStress(self):
     n = 20
     A = spmatrix.ll_mat(n, n)
     for k in range(20):
         for i in range(n * n / 2):
             i = random.randrange(n)
             j = random.randrange(n)
             A[i, j] = 1.0
         for i in range(n * n / 2):
             i = random.randrange(n)
             j = random.randrange(n)
             A[i, j] = 0.0

def ll_mat_rand(n, m, density):
    """return a ll_mat object representing a general n-by-m sparse matrix filled with random non-zero values

    The number of non-zero entries is less or equal than
    n*m*density. The values of the non-zero entries are in the range
    [0.0, 1.0)."""
    nnz = int(density * n * m)
    A = spmatrix.ll_mat(n, m, nnz)
    for k in xrange(nnz):
        i = random.randrange(n)
        j = random.randrange(m)
        A[i, j] = random.random()
    return A

    def setUp(self):

        self.nbr_elements = 1000
        self.size = 10000

        self.A_c = LLSparseMatrix(size=self.size, size_hint=self.nbr_elements, dtype=FLOAT64_T)
        construct_sparse_matrix(self.A_c, self.size, self.nbr_elements)

        self.A_p = spmatrix.ll_mat(self.size, self.size, self.nbr_elements)
        construct_sparse_matrix(self.A_p, self.size, self.nbr_elements)

        self.stride = 10

        self.v = np.arange(0, self.size * self.stride, dtype=np.float64)

def poisson2d(n):
    L = spmatrix.ll_mat(n * n, n * n)
    for i in range(n):
        for j in range(n):
            k = i + n * j
            L[k, k] = 4
            if i > 0:
                L[k, k - 1] = -1
            if i < n - 1:
                L[k, k + 1] = -1
            if j > 0:
                L[k, k - n] = -1
            if j < n - 1:
                L[k, k + n] = -1
    return L

def poisson2d_sym_blk_vec(n):
    n2 = n*n
    L = spmatrix.ll_mat_sym(n2, 3*n2-2*n)
    D = spmatrix.ll_mat_sym(n, 2*n-1)
    e = numpy.ones(n)
    d = numpy.arange(n, dtype=numpy.int)
    D.put(4*e, d, d)
    D.put(-e[1:], d[1:], d[:-1])
    P = spmatrix.ll_mat(n, n, n-1)
    P.put(-e,d,d)
    for i in xrange(n-1):
        L[i*n:(i+1)*n, i*n:(i+1)*n] = D
        L[(i+1)*n:(i+2)*n, i*n:(i+1)*n] = P
    # Last diagonal block
    L[n2-n:n2, n2-n:n2] = D
    return L

    def setUp(self):

        self.nbr_elements = 1000
        self.size = 10000

        self.A_c = LLSparseMatrix(size=self.size, size_hint=self.nbr_elements, itype=INT32_T, dtype=FLOAT64_T)
        self.A_p = spmatrix.ll_mat(self.size, self.size, self.nbr_elements)
        self.A_s = lil_matrix((self.size, self.size), dtype=np.float64)

        self.list_of_matrices = []
        self.list_of_matrices.append(self.A_c)
        self.list_of_matrices.append(self.A_p)
        self.list_of_matrices.append(self.A_s)

        construct_random_matrices(self.list_of_matrices, self.size, self.nbr_elements)

        self.v = np.arange(0, self.size, dtype=np.float64)

    def setUp(self):

        self.nbr_elements = 1000
        self.size = 10000
        self.nbr_elements_to_add = 1000

        assert self.nbr_elements_to_add < self.size

        self.A_c = LLSparseMatrix(size=self.size, size_hint=self.nbr_elements, itype=INT32_T, dtype=FLOAT64_T)
        construct_sparse_matrix(self.A_c, self.size, self.nbr_elements)

        self.A_p = spmatrix.ll_mat(self.size, self.size, self.nbr_elements)
        construct_sparse_matrix(self.A_p, self.size, self.nbr_elements)

        self.id1 = np.ones(self.nbr_elements_to_add, dtype=np.int32)
        self.id2 = self.v = np.arange(0, self.nbr_elements_to_add, dtype=np.int32)

        self.val = np.empty(self.nbr_elements_to_add, dtype=np.float64)

    def setUp(self):

        self.nbr_elements = 1000
        self.size = 10000
        self.put_size = 10000

        assert self.put_size <= self.size

        self.A_c = LLSparseMatrix(size=self.size, size_hint=self.nbr_elements, itype=INT32_T, dtype=FLOAT64_T)
        construct_sparse_matrix(self.A_c, self.size, self.nbr_elements)

        self.A_p = spmatrix.ll_mat(self.size, self.size, self.nbr_elements)
        construct_sparse_matrix(self.A_p, self.size, self.nbr_elements)

        self.id1 = np.arange(0, self.put_size, dtype=np.int32)
        self.id2 = np.full(self.put_size, 37, dtype=np.int32)

        self.b = np.arange(0, self.put_size,dtype=np.float64)