Bug fixed in LLSparseMatrix.put(i , j, val) and multiplication between CSR and CSC implemented

Author: counterclocker

examples/mul.py

@@ -0,0 +1,25 @@
+from sparse_lib.sparse.ll_mat import MakeLLSparseMatrix
+from sparse_lib.sparse.csc_mat import CSCSparseMatrix
+from sparse_lib.sparse.csr_mat import CSRSparseMatrix
+
+import sys
+
+A = MakeLLSparseMatrix(size=4)
+
+A.put_triplet([1.0, 2.0, 6.0, 3.0, 5.0, 4.0], [0, 1, 1, 2, 3, 3], [0, 1, 2, 2, 0, 3])
+
+A.print_to(sys.stdout)
+
+R = A.to_csr()
+
+R.print_to(sys.stdout)
+
+C = A.to_csc()
+
+C.debug_print()
+C.print_to(sys.stdout)
+
+print "*" * 80
+
+RC = R * C
+print RC

examples/umfpack.py

@@ -49,3 +49,13 @@
 print D
 print do_recip
 print R
+
+print "*" * 80
+print "L * U: "
+LU = L * U
+
+
+print "heheh"
+print LU
+import sys
+LU.print_to(sys.stdout)

sparse_lib/solvers/suitesparse/umfpack.c

@@ -209,19 +209,14 @@ cdef class UmfpackSolver:
 
         assert self.nrow == self.ncol, "Only square matrices are handled in UMFPACK"
 
-        # TODO: type csc in pxd file
         self.csc_mat  = self.A.to_csc()
 
-        #cdef double * val = self.val
-        #cdef int * col = <int *> self.col
-        #cdef int * ind = <int *> self.ind
-
         self.symbolic_computed = False
         self.numeric_computed = False
 
         # set default parameters for control
         umfpack_di_defaults(<double *>&self.control)
-        self.control[UMFPACK_PRL] = 3
+        self.set_verbosity(3)
 
     ####################################################################################################################
     # FREE MEMORY
@@ -422,6 +417,22 @@ cdef class UmfpackSolver:
         Returns:
             (L, U, P, Q, D, do_recip, R)
 
+            The original matrix A is factorized into
+
+                L U = P R A Q
+
+            where:
+             - L is unit lower triangular,
+             - U is upper triangular,
+             - P and Q are permutation matrices,
+             - R is a row-scaling diagonal matrix such that
+
+                  * the i-th row of A has been divided by R[i] if do_recip = True,
+                  * the i-th row of A has been multiplied by R[i] if do_recip = False.
+
+            L and U are returned as CSRSparseMatrix and CSCSparseMatrix sparse matrices respectively.
+            P, Q and R are returned as NumPy arrays.
+
         """
         # TODO: use properties?? we can only get matrices, not set them...
         self.create_numeric()

sparse_lib/solvers/suitesparse/umfpack.pyx

@@ -24,4 +24,7 @@ cdef class CSCSparseMatrix(ImmutableSparseMatrix):
         bint __col_indices_sorted  # are the column indices sorted in ascending order?
         int __first_row_not_ordered # first row that is not ordered
 
+    cdef at(self, int i, int j)
+    cdef safe_at(self, int i, int j)
+
 cdef MakeCSCSparseMatrix(int nrow, int ncol, int nnz, int * ind, int * row, double * val)

sparse_lib/sparse/csc_mat.c

@@ -10,6 +10,10 @@ from sparse_lib.sparse.sparse_mat cimport ImmutableSparseMatrix
 
 from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free
 
+cdef int CSC_MAT_PPRINT_ROW_THRESH = 500       # row threshold for choosing print format
+cdef int CSC_MAT_PPRINT_COL_THRESH = 20        # column threshold for choosing print format
+
+
 cdef class CSCSparseMatrix(ImmutableSparseMatrix):
     """
     Compressed Sparse Column Format matrix.
@@ -30,6 +34,126 @@ cdef class CSCSparseMatrix(ImmutableSparseMatrix):
         PyMem_Free(self.row)
         PyMem_Free(self.ind)
 
+    ####################################################################################################################
+    # Set/Get items
+    ####################################################################################################################
+    ####################################################################################################################
+    #                                            *** SET ***
+    def __setitem__(self, tuple key, value):
+        raise SyntaxError("Assign individual elements is not allowed")
+
+    #                                            *** GET ***
+    cdef at(self, int i, int j):
+        """
+        Direct access to element ``(i, j)``.
+
+        Warning:
+            There is not out of bounds test.
+
+        See:
+            :meth:`safe_at`.
+
+        """
+        cdef int k
+
+        if self.is_symmetric:
+            raise NotImplemented("Access to csr_mat(i, j) not (yet) implemented")
+
+        # TODO: column indices are NOT necessarily sorted... what do we do about it?
+        for k from self.ind[j] <= k < self.ind[j+1]:
+            if i == self.row[k]:
+                return self.val[k]
+
+        return 0.0
+
+    cdef safe_at(self, int i, int j):
+        """
+        Return element ``(i, j)`` but with check for out of bounds indices.
+
+        Raises:
+            IndexError: when index out of bound.
+
+        """
+        if not 0 <= i < self.nrow or not 0 <= j < self.ncol:
+            raise IndexError("Index out of bounds")
+
+        return self.at(i, j)
+
+    def __getitem__(self, tuple key):
+        """
+        Return ``csr_mat[i, j]``.
+
+        Args:
+          key = (i,j): Must be a couple of integers.
+
+        Raises:
+            IndexError: when index out of bound.
+
+        Returns:
+            ``csr_mat[i, j]``.
+        """
+        if len(key) != 2:
+            raise IndexError('Index tuple must be of length 2 (not %d)' % len(key))
+
+        cdef int i = key[0]
+        cdef int j = key[1]
+
+        return self.safe_at(i, j)
+
+    ####################################################################################################################
+    # String representations
+    ####################################################################################################################
+    def __repr__(self):
+        s = "CSCSparseMatrix of size %d by %d with %d non zero values" % (self.nrow, self.ncol, self.nnz)
+        return s
+
+    def print_to(self, OUT):
+        """
+        Print content of matrix to output stream.
+
+        Args:
+            OUT: Output stream that print (Python3) can print to.
+        """
+        # TODO: adapt to any numbers... and allow for additional parameters to control the output
+        cdef int i, k, first = 1;
+
+        cdef double *mat
+        cdef int j
+        cdef double val
+
+        print('CSCSparseMatrix ([%d,%d]):' % (self.nrow, self.ncol), file=OUT)
+
+        if not self.nnz or not self.__status_ok:
+            return
+
+        if self.nrow <= CSC_MAT_PPRINT_COL_THRESH and self.ncol <= CSC_MAT_PPRINT_ROW_THRESH:
+            # create linear vector presentation
+            # TODO: Skip this creation
+            mat = <double *> PyMem_Malloc(self.nrow * self.ncol * sizeof(double))
+
+            if not mat:
+                raise MemoryError()
+
+            for j from 0 <= j < self.ncol:
+                for i from 0 <= i < self.nrow:
+                    mat[j* self.nrow + i] = 0.0
+
+                k = self.ind[j]
+                while k < self.ind[j+1]:
+                    mat[(j*self.nrow)+self.row[k]] = self.val[k]
+                    k += 1
+
+            for i from 0 <= i < self.nrow:
+                for j from 0 <= j < self.ncol:
+                    val = mat[(j*self.nrow)+i]
+                    #print('%9.*f ' % (6, val), file=OUT, end='')
+                    print('{0:9.6f} '.format(val), end='')
+                print()
+
+            PyMem_Free(mat)
+
+        else:
+            print('Matrix too big to print out', file=OUT)
 
     ####################################################################################################################
     # DEBUG