Index: include/storage/sparse.h
===================================================================
--- include/storage/sparse.h	(revision 404)
+++ include/storage/sparse.h	(revision 403)
@@ -13,7 +13,7 @@
 *
 *	Copyright (C) The University of Twente, 2004-2008.
 *	Copyright (C) RWTH Aachen, 2008-2009.
-*       Authors: Maneesh Khattri, Ivan Zapreev, David N. Jansen
+*       Author: David N. Jansen
 *
 *	This program is free software; you can redistribute it and/or
 *	modify it under the terms of the GNU General Public License
@@ -51,46 +51,53 @@
 #include "macro.h"
 #include "error.h"
 
+typedef int state_diff; /* difference between states */
+
+/*@abstract@*/ struct mtx_el {          /* a single matrix entry */
+        state_diff col_minus_row;       /* next-state minus previous-state */
+        double val;                     /* transition probability */
+};
+
+typedef /*@dependent@*/ const struct mtx_el * mtx_el_ptr;
+
 /*****************************************************************************
 			STRUCTURE
 name		: values
 purpose         : information about a single state in a transition probability
                   matrix.
-@member col	: column indices of non-zero elements.
-@member val	: values of non-zero elements.
+@member entries : column indices and values of non-zero elements.
+                  elements are sorted ascendingly with key col_minus_row.
 @member back_set: states from which this state is reachable.
+                  The elements are sorted *descendingly* with key col_minus_row.
+                  May be NULL if that information has not been collected.
 @member ncols   : number of next-states, i.e. entries in the row, not counting
                   the diagonal
 @member pcols   : number of previous-states, i.e. entries in the column, not
-                  counting the diagonal.
+                  counting the diagonal. Is == -1 if the information about
+                  previous-states is not up-to-date for this state.
 @member diag	: The value of the diagonal element in this row.
-remark		: There is a one-one relation between col and val. For instance
-		  col[0] indicates the column index of val[0]. This structure
-		  is similar to the compressed-row storage for sparse matrices.
 ******************************************************************************/
 /* The following splint comments only hold for normal and contiguous (ncolse)
 sparse matrices. In wrapper matrices, ``owned'' would have to be replaced by
 ``dependent''. However, this is not achievable without very deep changes of the
 sparse matrix implementation. */
-typedef struct values
+/*@abstract@*/ struct values            /* details of a single state */
 {
         /*@owned@*/ /*@relnull@*/
-	int *col;    /* the column ids */
+        struct mtx_el * entries;        /* list of next-states */
                                         /* is required to be the first field in
                                            the structure. See
                                            allocate_sparse_matrix_ncolse(). */
-        /*@owned@*/ /*@relnull@*/
-	double *val; /* the corresponding values */
         /*@only@*/ /*@null@*/
-	int *back_set;
+        mtx_el_ptr * back_set;
                                         /* list of pointers to previous-states*/
         state_count ncols;              /* number of next-states, i.e. size of
-                                           *col and *val */
+                                           *entries */
         state_count pcols;              /* number of previous-states, i.e. size
                                            of *back_set */
 	double diag;
                                         /* probability of self-loop */
-}values;
+};
 
 /*****************************************************************************
 			STRUCTURE
@@ -99,8 +106,7 @@
 @member cols            : number of columns in the matrix.
 @member valstruc        : this is a vector of rows, each element is a structure,
 			containing values of non-zero elements.
-remark			: There is a one-one relation between col and val. For instance
-			col[0] indicates the column index of val[0]. This structure
+remark                  : This structure
 			is similar to the compressed-row storage for sparse matrices.
 ******************************************************************************/
 typedef /*@abstract@*/ struct sparse
@@ -109,11 +115,88 @@
                                         /* number of rows */
 	int cols;
                                         /* numer of columns */
+        unsigned flags;                 /* flags */
         struct values valstruc[1];      /* struct hack: actually a rows-sized
                                            array of information about the matrix
                                            rows */
 }sparse;
 
+/* Flag values
+   Flags indicate which type of matrix the structure is; as not all functions
+   can handle all kinds of matrices, they have to check dynamically the type
+   correctness.
+   Additionally, flags store whether back_sets are up-to-date or not. In the
+   latter case, they may have to be regenerated. */
+#define mtx_flag_NORMAL ((unsigned) 1)
+#define mtx_flag_WRAPPER ((unsigned) 2)
+#define mtx_flag_DIAGONAL (mtx_flag_NORMAL | mtx_flag_WRAPPER)
+#define mtx_flag_NCOLSE ((unsigned) 4)
+#define mtx_flag_BACK_SET_CORRECT ((unsigned) 8)
+
+/* Flag tests */
+extern BOOL mtx_is_normal(/*@observer@*/ /*@temp@*/ const sparse * pM);
+#define mtx_is_normal(pM) ((void) (pM)->rows, (void) (pM)->cols, \
+                        (void) (pM)->valstruc[0].entries, \
+                        (void) (pM)->valstruc[0].back_set, \
+                        (void) ((pM)->valstruc[0].ncols + \
+                                (pM)->valstruc[0].pcols), \
+                        (void) (pM)->valstruc[0].diag, \
+        mtx_flag_NORMAL == \
+        test_flag((pM)->flags,mtx_flag_NORMAL|mtx_flag_WRAPPER|mtx_flag_NCOLSE))
+extern BOOL mtx_is_normal_or_diag(/*@observer@*/ /*@temp@*/ const sparse * pM);
+#define mtx_is_normal_or_diag(pM) ((void) (pM)->rows, (void) (pM)->cols, \
+                        (void) (pM)->valstruc[0].entries, \
+                        (void) (pM)->valstruc[0].back_set, \
+                        (void) ((pM)->valstruc[0].ncols + \
+                                (pM)->valstruc[0].pcols), \
+                        (void) (pM)->valstruc[0].diag, \
+        mtx_flag_NORMAL == \
+        test_flag((pM)->flags, mtx_flag_NORMAL | mtx_flag_NCOLSE))
+extern BOOL mtx_is_wrap(/*@observer@*/ /*@temp@*/ const sparse * pM);
+#define mtx_is_wrap(pM) ((void) (pM)->rows, (void) (pM)->cols, \
+                        (void) (pM)->valstruc[0].entries, \
+                        (void) (pM)->valstruc[0].back_set, \
+                        (void) ((pM)->valstruc[0].ncols + \
+                                (pM)->valstruc[0].pcols), \
+                        (void) (pM)->valstruc[0].diag, \
+        mtx_flag_WRAPPER == \
+        test_flag((pM)->flags,mtx_flag_NORMAL|mtx_flag_WRAPPER|mtx_flag_NCOLSE))
+extern BOOL mtx_is_wrap_or_diag(/*@observer@*/ /*@temp@*/ const sparse * pM);
+#define mtx_is_wrap_or_diag(pM) ((void) (pM)->rows, (void) (pM)->cols, \
+                        (void) (pM)->valstruc[0].entries, \
+                        (void) (pM)->valstruc[0].back_set, \
+                        (void) ((pM)->valstruc[0].ncols + \
+                                (pM)->valstruc[0].pcols), \
+                        (void) (pM)->valstruc[0].diag, \
+        mtx_flag_WRAPPER == \
+        test_flag((pM)->flags, mtx_flag_WRAPPER | mtx_flag_NCOLSE))
+extern BOOL mtx_is_ncolse(/*@observer@*/ /*@temp@*/ const sparse * pM);
+#define mtx_is_ncolse(pM) ((void) (pM)->rows, (void) (pM)->cols, \
+                        (void) (pM)->valstruc[0].entries, \
+                        (void) (pM)->valstruc[0].back_set, \
+                        (void) ((pM)->valstruc[0].ncols + \
+                                (pM)->valstruc[0].pcols), \
+                        (void) (pM)->valstruc[0].diag, \
+        mtx_flag_NCOLSE == \
+        test_flag((pM)->flags,mtx_flag_NORMAL|mtx_flag_WRAPPER|mtx_flag_NCOLSE))
+extern BOOL mtx_is_diag(/*@observer@*/ /*@temp@*/ const sparse * pM);
+#define mtx_is_diag(pM) ((void) (pM)->rows, (void) (pM)->cols, \
+                        (void) (pM)->valstruc[0].entries, \
+                        (void) (pM)->valstruc[0].back_set, \
+                        (void) ((pM)->valstruc[0].ncols + \
+                                (pM)->valstruc[0].pcols), \
+                        (void) (pM)->valstruc[0].diag, \
+        mtx_flag_DIAGONAL == \
+        test_flag((pM)->flags,mtx_flag_NORMAL|mtx_flag_WRAPPER|mtx_flag_NCOLSE))
+extern BOOL mtx_has_backsets(/*@observer@*/ /*@temp@*/ const sparse * pM);
+#define mtx_has_backsets(pM) ((void) (pM)->rows, (void) (pM)->cols, \
+                        (void) (pM)->valstruc[0].entries, \
+                        (void) (pM)->valstruc[0].back_set, \
+                        (void) ((pM)->valstruc[0].ncols + \
+                                (pM)->valstruc[0].pcols), \
+                        (void) (pM)->valstruc[0].diag, \
+        0 != test_flag((pM)->flags, mtx_flag_BACK_SET_CORRECT))
+
 /************************************************************************/
 /**************      SIMPLE MATRIX INFORMATION ACCESS      **************/
 /************************************************************************/
@@ -132,6 +215,7 @@
                         /*@observer@*/ /*@temp@*/ const sparse * p_mtx)
                         /*@modifies nothing@*/;
 #       define mtx_rows(p_mtx) ((void) (p_mtx)->cols, \
+                        (void) (p_mtx)->valstruc[0].entries, \
                         (void) (p_mtx)->valstruc[0].back_set, \
                         (void) ((p_mtx)->valstruc[0].ncols + \
                                 (p_mtx)->valstruc[0].pcols), \
@@ -142,6 +226,7 @@
                         /*@observer@*/ /*@temp@*/ const sparse * p_mtx)
                         /*@modifies nothing@*/;
 #       define mtx_cols(p_mtx) ((void) (p_mtx)->rows, \
+                        (void) (p_mtx)->valstruc[0].entries, \
                         (void) (p_mtx)->valstruc[0].back_set, \
                         (void) ((p_mtx)->valstruc[0].ncols + \
                                 (p_mtx)->valstruc[0].pcols), \
@@ -153,19 +238,17 @@
                         state_index row) /*@modifies nothing@*/;
 #       define mtx_next_num(p_mtx,row) ((void) (p_mtx)->rows, \
                         (void) (p_mtx)->cols, \
+                        (void) (p_mtx)->valstruc[0].entries, \
                         (void) (p_mtx)->valstruc[0].back_set, \
                         (void) (p_mtx)->valstruc[0].pcols, \
                         (void) (p_mtx)->valstruc[0].diag, \
                         (const state_count) (p_mtx)->valstruc[(row)].ncols)
+
+        /* Because backsets are only generated upon explicit request, this
+           function cannot be implemented as a macro. */
         extern state_count mtx_prev_num(
                         /*@observer@*/ /*@temp@*/ const sparse * p_mtx,
                         state_index col);
-#       define mtx_prev_num(p_mtx,col) ((void) (p_mtx)->rows, \
-                        (void) (p_mtx)->cols, \
-                        (void) (p_mtx)->valstruc[0].back_set, \
-                        (void) (p_mtx)->valstruc[0].ncols, \
-                        (void) (p_mtx)->valstruc[0].diag, \
-                        (const state_count) (p_mtx)->valstruc[(col)].pcols)
 
 /*======================================================================*/
 /************************************************************************/
@@ -187,11 +270,6 @@
 /************************************************************************/
 /************************An alternative MATRIX ALLOCATION****************/
 /*======================================================================*/
-/**
-* WARNING: The way this sparse matrix is allocated in
-*       "allocate_sparse_matrix_ncolse" imposes some
-*	strict conditions on some operations on this matrix, such as adding a new element
-*/
 
 	/**
 	* Allocate a sparse matrix in a special way. Here we assume that we know the number of non-zero
@@ -220,6 +298,7 @@
 	*/
         extern err_state free_sparse_ncolse(/*@only@*/ sparse * pM)
                         /*@modifies pM@*/;
+#       define free_sparse_ncolse(pM) (free_mtx_sparse((pM)))
 
 	/**
 	* Sets a value to a pre-allocated element in the matrix.
@@ -233,6 +312,8 @@
 	*/
         extern err_state set_mtx_val_ncolse(sparse * pM, int row, int col,
                         double val) /*@modifies *pM@*/;
+#       define set_mtx_val_ncolse(pM,row,col,val) (set_mtx_val((pM), (row), \
+                        (col), (val)))
 
 	/**
 	* Searches for a pre-allocated matrix element and adds a value to it.
@@ -246,6 +327,8 @@
 	*/
         extern err_state add_mtx_val_ncolse(sparse * pM, int row, int col,
                         double val) /*@modifies *pM@*/;
+#       define add_mtx_val_ncolse(pM,row,col,val) (add_mtx_val((pM), (row), \
+                        (col), (val)))
 
 /*======================================================================*/
 /************************************************************************/
@@ -310,6 +393,7 @@
         @return         : err_ERROR: fail, err_OK: success
 	******************************************************************************/
         extern err_state free_mtx_wrap(/*@only@*/sparse * pM) /*@modifies pM@*/;
+#       define free_mtx_wrap(pM) (free_mtx_sparse((pM)))
 
 	/**
 	* Add a value to an existing one in the matrix or inserts a new value if the is none.
@@ -338,6 +422,7 @@
                         /*@temp@*/ /*@observer@*/ const sparse * pM,
                         int row) /*@modifies nothing@*/;
 #       define mtx_get_diag_val_nt(pM,row) ((void)(pM)->rows, (void)(pM)->cols,\
+                        (void) (pM)->valstruc[0].entries, \
                         (void) (pM)->valstruc[0].back_set, \
                         (void) ((pM)->valstruc[0].ncols + \
                                 (pM)->valstruc[0].pcols), \
@@ -367,6 +452,7 @@
                         double value) /*@modifies *pM@*/;
 #       define mtx_set_diag_val_nt(pM,row,value) ((void) (pM)->rows, \
                         (void) (pM)->cols, \
+                        (void) (pM)->valstruc[0].entries, \
                         (void) (pM)->valstruc[0].back_set, \
                         (void) ((pM)->valstruc[0].ncols + \
                                 (pM)->valstruc[0].pcols), \
@@ -881,8 +967,7 @@
                         state_index m_col = (row); \
                         double m_val = 0.0; \
                         state_count m_i__row; \
-                        const int * m_colrow; \
-                        const double * m_valrow; \
+                        const struct mtx_el * m_entries; \
                         if ( NULL == (p_mtx) || (unsigned) m_col >= \
                                                 (unsigned) mtx_rows((p_mtx)) ) \
                         { \
@@ -890,8 +975,7 @@
                                         "%d,col,val)", (const void *) (p_mtx), \
                                         m_col, EXIT_FAILURE)); \
                         } \
-                        m_colrow = (p_mtx)->valstruc[m_col].col; \
-                        m_valrow = (p_mtx)->valstruc[m_col].val; \
+                        m_entries = (p_mtx)->valstruc[m_col].entries; \
                         m_i__row = mtx_next_num((p_mtx), m_col) + 1; \
                         if ( m_col >= mtx_cols((p_mtx)) || \
                                         (m_val = mtx_get_diag_val_nt((p_mtx), \
@@ -902,12 +986,14 @@
                                 /* expression in the for statement; this */ \
                                 /* helps the compiler to optimize the code. */ \
                                 (void) (0 < --m_i__row && \
-                                        (m_col=*m_colrow++, m_val=*m_valrow++, \
+                                        (m_col=m_entries->col_minus_row+(row), \
+                                         m_val = m_entries++->val, \
                                          TRUE)); \
                         } \
                         for ( ; 0 < m_i__row ; \
                                 (void) (0 < --m_i__row && \
-                                        (m_col=*m_colrow++, m_val=*m_valrow++, \
+                                        (m_col=m_entries->col_minus_row+(row), \
+                                         m_val = m_entries++->val, \
                                          TRUE)) ) \
                         {
 
@@ -932,8 +1018,7 @@
                 { \
                         state_index m_col = (row); \
                         state_count m_i__row_nodiag; \
-                        const int * m_colrow; \
-                        const double * m_valrow; \
+                        const struct mtx_el * m_entries; \
                         if ( NULL == (p_mtx) || (unsigned) m_col >= \
                                                 (unsigned) mtx_rows((p_mtx)) ) \
                         { \
@@ -941,14 +1026,13 @@
                                         "nodiag(%p,%d,col,val)", (const void *)\
                                         (p_mtx), m_col, EXIT_FAILURE)); \
                         } \
-                        m_colrow = (p_mtx)->valstruc[m_col].col; \
-                        m_valrow = (p_mtx)->valstruc[m_col].val; \
+                        m_entries = (p_mtx)->valstruc[m_col].entries; \
                         for ( m_i__row_nodiag = mtx_next_num((p_mtx), m_col) ; \
                                         0 < m_i__row_nodiag-- ; \
-                                                m_colrow++, m_valrow++ ) \
+                                                m_entries++ ) \
                         { \
-                                const double m_val = *m_valrow; \
-                                m_col = *m_colrow;
+                                const double m_val = m_entries->val; \
+                                m_col = m_entries->col_minus_row + (row);
 
 #       define end_mtx_walk_row_nodiag \
                         } \
@@ -962,8 +1046,7 @@
                 { \
                         state_index m_col = (row); \
                         state_count m_i__row_sorted; \
-                        const int * m_colrow; \
-                        const double * m_valrow; \
+                        const struct mtx_el * m_entries; \
                         BOOL m_diag_passed = FALSE; \
                         if ( NULL == (p_mtx) || (unsigned) m_col >= \
                                                 (unsigned) mtx_rows((p_mtx)) ) \
@@ -972,8 +1055,7 @@
                                         "sorted(%p,%d,col,val)", (const void *)\
                                         (p_mtx), m_col, EXIT_FAILURE)); \
                         } \
-                        m_colrow = (p_mtx)->valstruc[m_col].col; \
-                        m_valrow = (p_mtx)->valstruc[m_col].val; \
+                        m_entries = (p_mtx)->valstruc[m_col].entries; \
                         for ( m_i__row_sorted = mtx_next_num((p_mtx), m_col) ; \
                                                 TRUE ; ) \
                         { \
@@ -985,7 +1067,7 @@
                                      /* surely not if there is another */ \
                                      /* element left of the diagonal */ \
                                         (0 < m_i__row_sorted \
-                                         && *m_colrow <= (row)) || \
+                                         && 0 >= m_entries->col_minus_row) || \
                                      /* yes, at least in principle */ \
                                         (m_diag_passed = TRUE, \
                                      /* except if the diagonal element */ \
@@ -1001,8 +1083,8 @@
                                         if ( 0 >= m_i__row_sorted ) \
                                                 break; \
                                         /* yes, there is: */ \
-                                        m_col = *m_colrow++; \
-                                        m_val = *m_valrow++; \
+                                        m_col = m_entries->col_minus_row+(row);\
+                                        m_val = m_entries++->val; \
                                         m_i__row_sorted--; \
                                 }
 
@@ -1018,35 +1100,34 @@
                 { \
                         state_index m_row = (clm); \
                         state_count m_i__column; \
-                        const int * m_rows; \
+                        const mtx_el_ptr * m_rows; \
                         double m_val = 0.0; \
+                        const char * m_error = err_PARAM; \
                         if ( NULL == (p_mtx) || (unsigned) m_row >= \
-                                                (unsigned) mtx_cols((p_mtx)) ) \
+                                                (unsigned) mtx_cols((p_mtx)) \
+                                        || (m_error = err_CALLBY, \
+                                            (m_i__column = mtx_prev_num( \
+                                                (p_mtx), m_row)) == \
+                                                state_index_ERROR) ) \
                         { \
-                                exit(err_macro_2(err_PARAM, "mtx_walk_column(" \
+                                exit(err_macro_2(m_error, "mtx_walk_column(" \
                                         "%p,row,%d,val)", (const void*)(p_mtx),\
                                         m_row, EXIT_FAILURE)); \
                         } \
-                        m_i__column = mtx_prev_num((p_mtx), m_row) + 1; \
+                        m_i__column++; \
                         m_rows = (p_mtx)->valstruc[m_row].back_set; \
                         if ( m_row >= mtx_rows((p_mtx)) || \
                                 (m_val = mtx_get_diag_val_nt((p_mtx), m_row)) \
                                         == 0.0 ) \
                         { \
                                 (void) (0 < --m_i__column \
-                                && (m_row = *m_rows++, err_state_iserror( \
-                                    get_mtx_val((p_mtx), m_row,(clm), &m_val)))\
-                                && (exit(err_macro_2(err_CALLBY, "mtx_walk_" \
-                                        "column(%p,row,%d,val)", (const void *)\
-                                        (p_mtx), (clm), EXIT_FAILURE)), TRUE));\
+                                && (m_row = (clm) - (*m_rows)->col_minus_row, \
+                                        m_val = (*m_rows++)->val, TRUE)); \
                         } \
                         for ( ; 0 < m_i__column ; \
                                 (void) (0 < --m_i__column \
-                                && (m_row = *m_rows++, err_state_iserror( \
-                                    get_mtx_val((p_mtx), m_row,(clm), &m_val)))\
-                                && (exit(err_macro_2(err_CALLBY, "mtx_walk_" \
-                                        "column(%p,row,%d,val)", (const void *)\
-                                        (p_mtx),(clm), EXIT_FAILURE)), TRUE)) )\
+                                && (m_row = (clm) - (*m_rows)->col_minus_row, \
+                                        m_val = (*m_rows++)->val, TRUE)) ) \
                         {
 
 #       define end_mtx_walk_column \
@@ -1061,63 +1142,64 @@
                 { \
                         state_index m_row = (clm); \
                         state_count m_i__column_nodiag; \
-                        const int * m_rows; \
+                        const mtx_el_ptr * m_rows; \
+                        const char * m_error = err_PARAM; \
                         if ( NULL == (p_mtx) || (unsigned) m_row >= \
-                                                (unsigned) mtx_cols((p_mtx)) ) \
+                                                (unsigned) mtx_cols((p_mtx)) \
+                                        || ((m_i__column_nodiag = \
+                                                mtx_prev_num((p_mtx),m_row)) ==\
+                                                state_index_ERROR \
+                                            && (m_error = err_CALLBY, TRUE)) ) \
                         { \
-                                exit(err_macro_2(err_PARAM, "mtx_walk_column_" \
+                                exit(err_macro_2(m_error, "mtx_walk_column_" \
                                         "nodiag(%p,row,%d,val)", (const void *)\
                                         (p_mtx), m_row, EXIT_FAILURE)); \
                         } \
                         m_rows = (p_mtx)->valstruc[m_row].back_set; \
-                        for ( m_i__column_nodiag=mtx_prev_num((p_mtx), m_row) ;\
+                        for ( ; \
                                 0 < m_i__column_nodiag ; m_rows++, \
                                                 m_i__column_nodiag-- ) \
                         { \
-                                double m_val; \
-                                m_row = *m_rows; \
-                                if ( err_state_iserror(get_mtx_val((p_mtx), \
-                                                m_row, (clm), &m_val)) ) \
-                                { \
-                                        exit(err_macro_2(err_CALLBY, "mtx_walk"\
-                                                "_column_nodiag(%p,row,%d,val)"\
-                                                , (const void *)(p_mtx), (clm),\
-                                                EXIT_FAILURE)); \
-                                }
+                                const double m_val = (*m_rows)->val; \
+                                m_row = (clm) - (*m_rows)->col_minus_row;
 
 #       define end_mtx_walk_column_nodiag \
                         } \
                         /*@-noeffect@*/(void)m_i__column_nodiag;/*@=noeffect@*/\
                 }
 
-        /* The _noval variant of mtx_walk_column is more efficient because the
-           call to get_mtx_val takes some time; if the value is not used, it is
-           not necessary to spend that. */
         /*@iter mtx_walk_column_noval(sef observer const sparse * p_mtx,
                         yield state_index m_row, sef unused state_index clm)@*/
 #       define mtx_walk_column_noval(p_mtx,m_row,clm) \
                 { \
                         state_index m_row = (clm); \
                         state_count m_i__column_noval; \
-                        const int * m_rows; \
+                        const mtx_el_ptr * m_rows; \
+                        const char * m_error = err_PARAM; \
                         if ( NULL == (p_mtx) || (unsigned) m_row >= \
-                                                (unsigned) mtx_cols((p_mtx)) ) \
+                                                (unsigned) mtx_cols((p_mtx)) \
+                                        || (m_error = err_CALLBY, \
+                                            (m_i__column_noval = mtx_prev_num( \
+                                                (p_mtx), m_row)) == \
+                                                state_index_ERROR) ) \
                         { \
-                                exit(err_macro_2(err_PARAM, "mtx_walk_column_" \
+                                exit(err_macro_2(m_error, "mtx_walk_column_" \
                                         "noval(%p,row,%d)", (const void *) \
                                         (p_mtx), m_row, EXIT_FAILURE)); \
                         } \
-                        m_i__column_noval = mtx_prev_num((p_mtx), m_row) + 1; \
+                        m_i__column_noval++; \
                         m_rows = (p_mtx)->valstruc[m_row].back_set; \
                         if ( m_row >= mtx_rows((p_mtx)) || \
                                 0.0 == mtx_get_diag_val_nt((p_mtx), m_row) ) \
                         { \
                                 (void) (0 < --m_i__column_noval \
-                                        && (m_row = *m_rows++, TRUE)); \
+                                && (m_row = (clm) - (*m_rows++)->col_minus_row,\
+                                        TRUE)); \
                         } \
                         for ( ; 0 < m_i__column_noval ; \
                                 (void) (0 < --m_i__column_noval \
-                                        && (m_row = *m_rows++, TRUE)) ) \
+                                && (m_row = (clm) - (*m_rows++)->col_minus_row,\
+                                        TRUE)) ) \
                         {
 
 #       define end_mtx_walk_column_noval \
@@ -1131,21 +1213,25 @@
                 { \
                         state_index m_row = (clm); \
                         state_count m_i__column_nodiag_noval; \
-                        const int * m_rows; \
+                        const mtx_el_ptr * m_rows; \
+                        const char * m_error = err_PARAM; \
                         if ( NULL == (p_mtx) || (unsigned) m_row >= \
-                                                (unsigned) mtx_cols((p_mtx)) ) \
+                                                (unsigned) mtx_cols((p_mtx)) \
+                                        || ((m_i__column_nodiag_noval = \
+                                                mtx_prev_num((p_mtx),m_row)) ==\
+                                                state_index_ERROR \
+                                            && (m_error = err_CALLBY, TRUE)) ) \
                         { \
-                                exit(err_macro_2(err_PARAM, "mtx_walk_column_" \
+                                exit(err_macro_2(m_error, "mtx_walk_column_" \
                                         "nodiag_noval(%p,row,%d)",(const void*)\
                                         (p_mtx), m_row, EXIT_FAILURE)); \
                         } \
                         m_rows = (p_mtx)->valstruc[m_row].back_set; \
-                        for ( m_i__column_nodiag_noval = mtx_prev_num((p_mtx), \
-                                        m_row) ; \
+                        for ( ; \
                                 0 < m_i__column_nodiag_noval ; m_rows++, \
                                                 m_i__column_nodiag_noval-- ) \
                         { \
-                                m_row = *m_rows;
+                                m_row = (clm) - (*m_rows)->col_minus_row;
 
 #       define end_mtx_walk_column_nodiag_noval \
                         } \
@@ -1162,8 +1248,7 @@
                         double m_val; \
                         state_count m_i__all; \
                         const struct values * m_values; \
-                        const int * m_colrow; \
-                        const double * m_valrow; \
+                        const struct mtx_el * m_entries; \
                         if ( NULL == (p_mtx) ) { \
                                 exit(err_macro_1(err_PARAM, "mtx_walk_all(%p," \
                                         "row,col,val)", (const void *) (p_mtx),\
@@ -1178,7 +1263,7 @@
                         /* of commas instead of semicola. */ \
                         (void) (m_col = --m_row, --m_values, \
                         m_i__all = m_values->ncols + 1, \
-                        m_colrow = m_values->col, m_valrow = m_values->val, \
+                        m_entries = m_values->entries, \
                         m_val = m_col<mtx_cols((p_mtx)) ? m_values->diag : 0.0,\
                         TRUE); \
                         for( ; 0 < m_i__all ; \
@@ -1186,13 +1271,13 @@
                                      ? (void) (0 < m_row && \
                                        (m_col = --m_row, --m_values, \
                                         m_i__all = m_values->ncols + 1, \
-                                        m_colrow = m_values->col, \
-                                        m_valrow = m_values->val, \
+                                        m_entries = m_values->entries, \
                                         m_val = m_col < mtx_cols((p_mtx)) \
                                                         ? m_values->diag : 0.0,\
                                         TRUE)) \
-                                     : (void) (m_col = *m_colrow++, \
-                                        m_val = *m_valrow++) ) \
+                                     : (void) (m_col = m_entries->col_minus_row\
+                                                        + m_row, \
+                                        m_val = m_entries++->val) ) \
                         { \
                                 /*@-realcompare@*/ \
                                 if ( 0.0 != m_val ) /*@=realcompare@*/ {
Index: include/algorithms/random_numbers/rand_num_generator.h
===================================================================
--- include/algorithms/random_numbers/rand_num_generator.h	(revision 404)
+++ include/algorithms/random_numbers/rand_num_generator.h	(revision 403)
@@ -50,6 +50,8 @@
 #ifndef RANDOM_NON_UNIF_H
 #define RANDOM_NON_UNIF_H
 
+#include "sparse.h"
+
 /**
 * The available random number generators for the simulation.
 * (1) Combined linear congruential generator taken from "Applied
@@ -122,8 +124,8 @@
 * @param prob the related probability distribution
 * @param size the size of the values and prob arrays
 */
-extern int generateRandNumberDiscrete(const int * values_local,
-                const double * distribution, int size);
+extern int generateRandNumberDiscrete(const struct mtx_el * distribution,
+                int size, state_index row);
 
 /**
 * This function frees all memory allocated for the random number generator
Index: include/modelchecking/simulation_common.h
===================================================================
--- include/modelchecking/simulation_common.h	(revision 404)
+++ include/modelchecking/simulation_common.h	(revision 403)
@@ -193,9 +193,9 @@
                 (0 != mtx_next_num((pM), (current_obs_state)) \
                         ? /* Compute to what state we will go */ \
                           generateRandNumberDiscrete( \
-                                (pM)->valstruc[(current_obs_state)].col, \
-                                (pM)->valstruc[(current_obs_state)].val, \
-                                mtx_next_num((pM), (current_obs_state))) \
+                                (pM)->valstruc[(current_obs_state)].entries, \
+                                mtx_next_num((pM), (current_obs_state)), \
+                                (current_obs_state)) \
                         : (current_obs_state))
 
 	/**
Index: src/storage/sparse.c
===================================================================
--- src/storage/sparse.c	(revision 404)
+++ src/storage/sparse.c	(revision 403)
@@ -13,7 +13,7 @@
 *
 *	Copyright (C) The University of Twente, 2004-2008.
 *	Copyright (C) RWTH Aachen, 2008-2009.
-*       Authors: Maneesh Khattri, Ivan Zapreev, David N. Jansen
+*       Author: David N. Jansen
 *
 *	This program is free software; you can redistribute it and/or
 *	modify it under the terms of the GNU General Public License
@@ -61,27 +61,28 @@
                         state_index m_col = (row); \
                         double * m_p_val; \
                         state_count m_i__c_row; \
-                        const int * m_colrow; \
-                        double * m_valrow; \
-                        if ( NULL == (p_mtx) || (unsigned) m_col >= \
+                        struct mtx_el * m_entries; \
+                        if ( NULL == (p_mtx) /* || mtx_is_wrap((p_mtx)) */ \
+                                        || (unsigned) m_col >= \
                                                 (unsigned) mtx_rows((p_mtx)) ) \
                         { \
                                 exit(err_macro_2(err_PARAM, "mtx_change_row(" \
                                         "%p,%d,col,val)", (void *) (p_mtx), \
                                         m_col, EXIT_FAILURE)); \
                         } \
-                        m_colrow = (p_mtx)->valstruc[m_col].col; \
-                        m_valrow = (p_mtx)->valstruc[m_col].val; \
+                        m_entries = (p_mtx)->valstruc[m_col].entries; \
                         m_i__c_row = mtx_next_num((p_mtx), m_col) + 1; \
                         m_p_val = &(p_mtx)->valstruc[m_col].diag; \
                         if ( m_col >= mtx_cols((p_mtx)) ) { \
                                 (void) (0 < --m_i__c_row && \
-                                        (m_col=*m_colrow++, m_p_val=m_valrow++,\
+                                        (m_col=m_entries->col_minus_row+(row), \
+                                         m_p_val = &m_entries++->val, \
                                          TRUE)); \
                         } \
                         for ( ; 0 < m_i__c_row ; \
                                 (void) (0 < --m_i__c_row && \
-                                        (m_col=*m_colrow++, m_p_val=m_valrow++,\
+                                        (m_col=m_entries->col_minus_row+(row), \
+                                         m_p_val = &m_entries++->val, \
                                          TRUE)) ) \
                         {
 
@@ -96,23 +97,22 @@
                 { \
                         state_index m_col = (row); \
                         state_count m_i__c_row_nodiag; \
-                        const int * m_colrow; \
-                        double * m_valrow; \
-                        if ( NULL == (p_mtx) || (unsigned) m_col >= \
+                        struct mtx_el * m_entries; \
+                        if ( NULL == (p_mtx) || mtx_is_wrap((p_mtx)) \
+                                             || (unsigned) m_col >= \
                                                 (unsigned) mtx_rows((p_mtx)) ) \
                         { \
                                 exit(err_macro_2(err_PARAM, "mtx_change_row_" \
                                         "nodiag(%p,%d,col,p_val)", \
                                         (void*)(p_mtx), m_col, EXIT_FAILURE)); \
                         } \
-                        m_colrow = (p_mtx)->valstruc[m_col].col; \
-                        m_valrow = (p_mtx)->valstruc[m_col].val; \
+                        m_entries = (p_mtx)->valstruc[m_col].entries; \
                         for ( m_i__c_row_nodiag = mtx_next_num((p_mtx),m_col) ;\
                                         0 < m_i__c_row_nodiag-- ; \
-                                                m_colrow++, m_valrow++ ) \
+                                                m_entries++ ) \
                         { \
-                                double * const m_p_val = m_valrow; \
-                                m_col = *m_colrow;
+                                double * const m_p_val = &m_entries->val; \
+                                m_col = m_entries->col_minus_row + (row);
 
 #define end_mtx_change_row_nodiag \
                         } \
@@ -127,9 +127,8 @@
                         double * m_p_val; \
                         state_count m_i__c_all; \
                         struct values * m_values; \
-                        const int * m_colrow; \
-                        double * m_valrow; \
-                        if ( NULL == (p_mtx) ) { \
+                        struct mtx_el * m_entries; \
+                        if ( NULL == (p_mtx) /* || mtx_is_wrap((p_mtx)) */ ) { \
                                 exit(err_macro_1(err_PARAM, "mtx_change_all(" \
                                         "%p,row,col,p_val)", (void *) (p_mtx), \
                                         EXIT_FAILURE)); \
@@ -138,7 +137,7 @@
                         m_values = &(p_mtx)->valstruc[m_row]; \
                         (void) (m_col = --m_row, --m_values, \
                         m_i__c_all = m_values->ncols + 1, \
-                        m_colrow = m_values->col, m_valrow = m_values->val, \
+                        m_entries = m_values->entries, \
                         m_p_val = m_col < mtx_cols((p_mtx)) ? &m_values->diag \
                                                             : NULL, \
                         TRUE); \
@@ -147,13 +146,13 @@
                                      ? (void) (0 < m_row && \
                                        (m_col = --m_row, --m_values, \
                                         m_i__c_all = m_values->ncols + 1, \
-                                        m_colrow = m_values->col, \
-                                        m_valrow = m_values->val, \
+                                        m_entries = m_values->entries, \
                                         m_p_val = m_col < mtx_cols((p_mtx)) \
                                                 ? &m_values->diag : NULL, \
                                         TRUE)) \
-                                     : (void) (m_col = *m_colrow++, \
-                                        m_p_val = m_valrow++) ) \
+                                     : (void) (m_col = \
+                                                m_entries->col_minus_row+m_row,\
+                                        m_p_val = &m_entries++->val) ) \
                         { \
                                 if ( NULL != m_p_val ) {
 
@@ -163,6 +162,99 @@
                 }
 
 
+/**
+* This function (re)generates the backset if needed and calculates the number
+* of predecessors of a state.
+* Parameter:    p_mtx = matrix to be used in the calculation
+*               req_col = column whose number of predecessors is needed
+*/
+state_count mtx_prev_num(/*@temp@*/ const sparse * pM,
+                        state_index req_col)
+{
+        /* Trick to avoid warning about cast to non-const sparse * */
+        sparse * p_mtx = * (sparse **) (void *) &pM;
+        state_index row, j;
+        struct values * valrow;
+        state_count pnum;
+
+        if ( NULL == p_mtx || mtx_rows(p_mtx) != mtx_cols(p_mtx)
+                        || mtx_cols(p_mtx) <= req_col || 0 > req_col )
+        {
+                err_msg_4(err_PARAM, "mtx_prev_num(%p[%dx%d],%d)", (void*)p_mtx,
+                        mtx_rows(p_mtx), mtx_cols(p_mtx), req_col,
+                        state_index_ERROR);
+        }
+
+        if ( mtx_has_backsets(p_mtx) && (pnum = p_mtx->valstruc[req_col].pcols)
+                                != state_index_NONE )
+        {
+                /* The required backset is up-to-date; do nothing. */
+                return pnum;
+        }
+
+        /* delete all old backsets */
+        row = mtx_rows(p_mtx);
+        valrow = &p_mtx->valstruc[0];
+        do {
+                free(valrow->back_set);
+                valrow->back_set = NULL;
+                valrow++->pcols = 0;
+        } while ( 0 < --row );
+
+        /* recreate new backsets */
+        row = 0;
+        valrow = &p_mtx->valstruc[0];
+        do {
+                /* we cannot use mtx_walk_... here because we need the pointer
+                   to the internal data structure. */
+                const struct mtx_el * entries = valrow->entries;
+                for ( j = valrow->ncols ; 0 < j ; j--, entries++ ) {
+                        const state_index col = entries->col_minus_row + row;
+                        struct values * valcol = &p_mtx->valstruc[col];
+                        /*@only@*/ /*@null@*/ mtx_el_ptr * temp_back_set;
+
+                        pnum = valcol->pcols;
+                        if ( state_index_NONE == pnum ) {
+                                /* there was not enough memory for this column's
+                                   backset in an earlier iteration. */
+                                continue;
+                        }
+
+                        temp_back_set = realloc(valcol->back_set,
+                                        (pnum + 1) * sizeof(mtx_el_ptr));
+                        if ( NULL == temp_back_set ) {
+                                if ( col != req_col ) {
+                                        /* Let's hope this backset is not
+                                           needed any more. */
+                                        /*@-usereleased@*/
+                                        free(valcol->back_set);
+                                        /*@=usereleased@*/
+                                        valcol->back_set = NULL;
+                                        valcol->pcols = state_index_NONE;
+                                        continue;
+                                }
+                                /*@-usereleased@*/ /*@-compdef@*/
+                                err_msg_4(err_MEMORY, "mtx_prev_num(%p[%dx%d],"
+                                        "%d)", (void *) p_mtx, mtx_rows(p_mtx),
+                                        mtx_cols(p_mtx), req_col,
+                                        state_index_ERROR);
+                                /*@=compdef@*/ /*@=usereleased@*/
+                        }
+                        temp_back_set[pnum] = entries;
+                        valcol->back_set = temp_back_set;
+                        valcol->pcols = pnum + 1;
+                  /*@-branchstate@*/
+                } /*@=branchstate@*/
+                valrow++;
+        } while ( ++row < /*@-compmempass@*/mtx_rows(p_mtx)/*@=compmempass@*/ );
+        set_flag(p_mtx->flags, mtx_flag_BACK_SET_CORRECT);
+
+        /*@-compdef@*/
+        return p_mtx->valstruc[req_col].pcols;
+        /*@=compdef@*/
+}
+
+
 /*======================================================================*/
 /************************************************************************/
 /************************A trivial MATRIX ALLOCATION*********************/
@@ -190,6 +282,8 @@
         }
 	pMatrix->rows = rows;
 	pMatrix->cols = cols;
+        pMatrix->flags = rows != cols ? mtx_flag_DIAGONAL
+                                : mtx_flag_DIAGONAL | mtx_flag_BACK_SET_CORRECT;
 	return pMatrix;
 }
 
@@ -219,8 +313,7 @@
 /*@only@*/ /*@null@*/
 sparse * allocate_sparse_matrix_ncolse(const int rows, const int cols, const int *ncolse){
 	int i, sum=0;
-	int *col_val;
-	double *val;
+        struct mtx_el * entries;
         sparse * pMatrix;
 
         if ( 0 >= rows || 0 >= cols || NULL == ncolse ) {
@@ -228,10 +321,10 @@
                                 rows, cols, (const void *) ncolse, NULL);
         }
 
-        /* In an ncolse matrix, I allocate one valstruc.col too much to
+        /* In an ncolse matrix, I allocate one valstruc.entries too much to
            allow for checking the space size even for the last row. */
         pMatrix = (sparse*) calloc(1, sizeof(sparse) - sizeof(pMatrix->valstruc)
-                                + sizeof(pMatrix->valstruc[rows].col)
+                                + sizeof(pMatrix->valstruc[rows].entries)
                                 + rows * sizeof(pMatrix->valstruc[0]));
         if ( NULL == pMatrix ) {
                 err_msg_3(err_MEMORY, "allocate_sparse_matrix_ncolse(%d,%d,%p)",
@@ -239,189 +332,246 @@
         }
 	pMatrix->rows = rows;
 	pMatrix->cols = cols;
+        pMatrix->flags = rows != cols ? mtx_flag_NCOLSE
+                                : mtx_flag_NCOLSE | mtx_flag_BACK_SET_CORRECT;
 	for(i=0;i<rows;i++){
 		sum+=ncolse[i];
 	}
 
         /* Allocate the continuous arrays for all non-zero column indices and
            values */
-        col_val = (int *) calloc((size_t) sum, sizeof(int));
-        val = (double *) calloc((size_t) sum, sizeof(double));
-        if ( NULL == col_val || NULL == val ) {
+        entries = (struct mtx_el *) calloc((size_t) sum, sizeof(struct mtx_el));
+        if ( NULL == entries ) {
                 err_msg_3(err_MEMORY, "allocate_sparse_matrix_ncolse(%d,%d,%p)",
                                 rows, cols, (const void *) ncolse,
-                                (free(val),free(col_val), free(pMatrix), NULL));
+                                (free(pMatrix), NULL));
         }
 
 	/* Initialize pointers for the matrix rows */
 	sum=0;
         for ( i = 0 ; ; i++ )
 	{
-		pMatrix->valstruc[i].col = &col_val[sum];
-                /* valstruc[rows].col is also allocated */
+		pMatrix->valstruc[i].entries = &entries[sum];
+                /* valstruc[rows].entries is also allocated */
                 if ( i >= rows )
                         break;
-                /* but valstruc[rows].val or ncolse[rows] may generate a bus
+                /* but ncolse[rows] may generate a bus
                    error */
-	 	pMatrix->valstruc[i].val = &val[sum];
 		sum+=ncolse[i];
 	}
 	return pMatrix;
 }
 
-/**
-* Frees the sparse matrix allocated with the allocate_sparse_matrix_ncolse(...) method.
-* WARNING: This method should be used with the matrices allocated with
-*	the allocate_sparse_matrix_ncolse(...) method.
-* @param pM the matrix to remove
+/***
+* The internal function looks for an entry in the array of entries provided.
+* entries[] - the array of entries to look in. May be NULL is size == 0.
+* size - the number of entries in the array
+* col_minus_row - the key to look for
+* found (out) - indicates whether the entry has been found.
+* Result: the array index of the entry; if not found, the array index where the
+*               entry should be inserted. state_index_ERROR if some error
+*               happened.
 */
-err_state free_sparse_ncolse(/*@only@*/ /*@i1@*/ /*@null@*/ sparse * pM) {
-	int i;
-	if( pM != NULL ) {
-		if(pM->rows==pM->cols){
-                        for ( i = 0 ; i< mtx_rows(pM) ; i++ ) {
-                                free(pM->valstruc[i].back_set);
-			}
-		}
-                free(pM->valstruc[0].col);
-                free(pM->valstruc[0].val);
-		free(pM);
-	}
-        else {
-                err_msg_3(err_PARAM, "free_sparse_ncolse(%p[%dx%d])",
-                                (void *) pM, NULL != pM ? mtx_rows(pM) : 0,
-                                NULL != pM ? mtx_cols(pM) : 0, err_ERROR);
+static state_index find_entry(const struct mtx_el entries[], state_count size,
+                state_diff col_minus_row, /*@out@*/ BOOL * found)
+{
+        state_index low;
+
+        if ( 0 > size || (0 < size && NULL == entries) || NULL == found ) {
+                /*@-mustdefine@*/
+                err_msg_4(err_PARAM, "find_entry(%p,%d,%d,%p)",
+                                (const void *) entries, size, col_minus_row,
+                                (void *) found, state_index_ERROR);
+                /*@=mustdefine@*/
         }
-        return err_OK;
+
+        /* binary search for the correct col_minus_row index */
+        low = 0; /* the lowest index that is not known to be too small */
+        /* size == the lowest index that is known to be too large */
+        while ( low < size ) {
+                state_index mid = (low + size) / 2;
+                if ( entries[mid].col_minus_row < col_minus_row ) {
+                        /* this check first because most insertions take place
+                           in the rightmost column */
+                        low = mid + 1;
+                } else if ( entries[mid].col_minus_row > col_minus_row ) {
+                        size = mid;
+                } else {
+                        *found = TRUE;
+                        return mid;
+                }
+        }
+        /* Now low >= size, so the element has not been found */
+        *found = FALSE;
+        return size;
 }
 
-/**
-* Sets a value to a pre-allocated element in the matrix.
-* WARNING: This method should be used with the matrices allocated with
-*	the allocate_sparse_matrix_ncolse(...) method.
-*       It assumes that the matrix be filled from left to right and from top to
-*       bottom: It appends new elements to the end of the respective arrays.
-* @param pM the matrix to add element to
-* @param row the row number in which the value is to be set.
-* @param col the column number in which the value is to be set.
-* @param val the value to be set.
+
+/***
+* The internal function looks for an entry in the matrix and inserts a new entry
+* if it hasn't been found.
+* pM - the matrix to look in
+* row, col - the row and column of the entry.
+* found (out) - indicates whether the entry has been found. May be NULL.
+* result (out) - pointer to the entry. If not found, result points to the next
+*               entry in the same row. If inserted, (*result)->col_minus_row is
+*               already set to the correct value and (*result)->val == 0.0.
 */
-err_state set_mtx_val_ncolse(/*@i1@*/ /*@null@*/ sparse * pM,
-                int row, int col, double val)
+static err_state insert_entry(sparse * pM, state_index row, state_index col,
+                /*@out@*/ struct mtx_el ** result)
 {
-        int idx, idxp = 0;
+        state_diff col_minus_row;
 
         if ( NULL == pM || row < 0 || row >= mtx_rows(pM) || col < 0
-                                || col >= mtx_cols(pM) )
+                                || col >= mtx_cols(pM) || mtx_is_wrap(pM)
+                                || row == col || NULL == result )
         {
-                err_msg_6(err_PARAM, "set_mtx_val_ncolse(%p[%dx%d],%d,%d,%g)",
+                err_msg_6(err_PARAM, "insert_entry(%p[%dx%d],%d,%d,%p)",
                                 (void *) pM, NULL != pM ? mtx_rows(pM) : 0,
                                 NULL != pM ? mtx_cols(pM) : 0, row, col,
-                                val, err_ERROR);
+                                (void *) result, err_ERROR);
         }
 
+        col_minus_row = col - row;
+        {
+                BOOL found;
+                struct mtx_el * entries = pM->valstruc[row].entries;
 		/* Here pM->ncols[row] is initially 0. */
-		if(val != 0 && row != col){
-                        idx = mtx_next_num(pM, row);
+                state_count idx = mtx_next_num(pM, row);
+                state_count found_idx = find_entry(entries, idx, col_minus_row,
+                                &found);
+                if ( state_index_ERROR == found_idx || 0 > found_idx
+                                                || idx < found_idx )
+                {
+                        /*@-mustdefine@*/
+                        err_msg_6(state_index_ERROR == found_idx ? err_CALLBY
+                                : err_INCONSISTENT, "insert_entry(%p[%dx%d],%d,"
+                                "%d,%p)",(void *) pM, mtx_rows(pM),mtx_cols(pM),
+                                row,col, (void *) result, err_ERROR);
+                        /*@=mustdefine@*/
+                }
+                if ( found ) {
+                        *result = &entries[found_idx];
+                        return err_OK;
+                }
+                if ( mtx_is_ncolse(pM) ) {
                         /* test if enough space for another entry has been
                            reserved */
-                        if ( &pM->valstruc[row].col[idx]
-                                                >= pM->valstruc[row+1].col )
+                        if ( &entries[idx] >= pM->valstruc[row + 1].entries )
                         {
                                 /* Too many elements in this row */
-                                err_msg_6(err_INCONSISTENT, "set_mtx_val_ncolse"
-                                        "(%p[%dx%d],%d,%d,%g)", (void *) pM,
+                                err_msg_6(err_INCONSISTENT, "insert_entry"
+                                        "(%p[%dx%d],%d,%d,%p)", (void *) pM,
                                         mtx_rows(pM), mtx_cols(pM), row, col,
-                                        val, err_ERROR);
+                                        (void *) result, err_ERROR);
                         }
-                        pM->valstruc[row].col[idx] = col; /* Set the new nonzero
-                                                        element column id */
-                        pM->valstruc[row].val[idx] = val; /* Set the new nonzero
-                                                        element value */
-                        /* Check if it is a square matrix */
-                        if ( mtx_rows(pM) == mtx_cols(pM) ) {
-                                int * temp_back_set;
+                } else {
+                        /* allocate space for the new entry */
+                        /*@null@*/ struct mtx_el *new_entries = realloc(entries,
+                                        (idx + 1) * sizeof(struct mtx_el));
+                        if ( NULL == new_entries ) {
+                                /*@-usereleased@*/ /*@-compdef@*/
+                                /*@-mustdefine@*/
+                                err_msg_6(err_MEMORY, "insert_entry(%p[%dx%d],"
+                                        "%d,%d,%p)", (void *) pM, mtx_rows(pM),
+                                        mtx_cols(pM), row, col, (void *) result,
+                                        err_ERROR);
+                                /*@=mustdefine@*/
+                                /*@=compdef@*/ /*@=usereleased@*/
+                        }
+                        /*@-usereleased@*/
+                        if ( entries != new_entries ) {
+                                if ( NULL != entries ) /*@=usereleased@*/ {
+                                        clear_flag(pM->flags,
+                                                mtx_flag_BACK_SET_CORRECT);
+                                }
+                                pM->valstruc[row].entries = entries=new_entries;
+                        }
+                        clear_flag(pM->flags, mtx_flag_WRAPPER);
+                }
+                pM->valstruc[row].ncols++; /* Increase the number of
+                                        nonzero entries in the current row */
+
+                /* Now let entries point to the new entry */
+                entries = *result = &entries[found_idx];
+                /* Move entries in further columns to the right */
+                if ( found_idx < idx ) {
+                        clear_flag(pM->flags, mtx_flag_BACK_SET_CORRECT);
+                        memmove(entries + 1, entries,
+                                (idx - found_idx) * sizeof(struct mtx_el));
+                }
+                /* Initialize the new entry */
+                entries->col_minus_row = col_minus_row;
+                entries->val = 0.0;
+        }
+
+        if ( mtx_has_backsets(pM) ) {
+                /* Try to keep the backset up-to-date */
+
+                /*@null@*/ mtx_el_ptr * temp_back_set = NULL;
+                state_index low, high, mid;
                                 /* Get the number of states from which we can
                                    reach state col */
-                                idxp = mtx_prev_num(pM, col);
+                state_count idxp = pM->valstruc[col].pcols;
+
+                if ( state_index_NONE == idxp ) {
+                        /* this column's backset is not up-to-date anyway, so
+                           don't spend too much effort */
+                        return err_OK;
+                }
                                 /* increase the size of the back_set array */
                                 temp_back_set = realloc(
                                                 pM->valstruc[col].back_set,
-                                                (idxp + 1) * sizeof(int));
+                                                (idxp+1) * sizeof(mtx_el_ptr));
                                 if ( NULL == temp_back_set ) {
-                                        /*@-usereleased@*/ /*@-compdef@*/
-                                        err_msg_6(err_MEMORY,
-                                                "set_mtx_val_ncolse(%p[%dx%d],"
-                                                "%d,%d,%g)", (void *) pM,
-                                                mtx_rows(pM), mtx_cols(pM), row,
-                                                col, val, err_ERROR);
-                                        /*@=usereleased@*/ /*@=compdef@*/
+                                        /* if there is not enough memory to keep
+                                           this backset, we just hope it's not
+                                           needed any more */
+                                        pM->valstruc[col].pcols
+                                                        = state_index_NONE;
+                                        free(pM->valstruc[col].back_set);
+                                        pM->valstruc[col].back_set = NULL;
+                                        return err_OK;
                                 }
                                 pM->valstruc[col].back_set = temp_back_set;
-                                /* Add the current row to this set */
-                                temp_back_set[idxp] = row;
                                 /* Increase the counter for the number of
                                    elements in the back_set */
 				pM->valstruc[col].pcols++;
-			}
-                        pM->valstruc[row].ncols++; /* Increase the number of
-                                        nonzero entries in the current row */
-		} else {
-			if(val != 0 && row == col){
-                                /* Handle the diagonal element */
-                                mtx_set_diag_val_nt(pM, row, val);
-			}
-		}
-        return err_OK;
-}
 
-/**
-* Searches for a pre-allocated matrix element and adds a value to it.
-* WARNING: This method should be used with the matrices allocated with
-*	the allocate_sparse_matrix_ncolse(...) method.
-* @param pM the matrix
-* @param row the row number in which the value is to be added.
-* @param col the column number in which the value is to be added.
-* @param val the value to be added.
-*/
-err_state add_mtx_val_ncolse(/*@i1@*/ /*@null@*/ sparse * pM,
-                int row, int col, double val)
-{
-        int ncols, i;
-	BOOL found = FALSE;
-        if ( NULL == pM || 0 > row || row >= mtx_rows(pM) || 0 > col
-                        || col >= mtx_cols(pM) )
-                err_msg_6(err_PARAM, "add_mtx_val_ncolse(%p[%dx%d],%d,%d,%g)",
-                                (void *) pM, NULL != pM ? mtx_rows(pM) : 0,
-                                NULL != pM ? mtx_cols(pM) : 0, row, col, val,
-                                err_ERROR);
-        if ( /*@-realcompare@*/ 0.0 != val /*@=realcompare@*/ ) {
-		if ( row != col ) {
-                        ncols = mtx_next_num(pM, row);
-			/* TODO: We could be smart and make at least a binary */
-			/* search here, as soon as all elements in pM->val[row].col */
-                        /* array are actually required to be ordered by value */
-			/* (increasing with the arrayindex) */
-			for( i=0; i < ncols; i++ ) {
-                                if ( pM->valstruc[row].col[i] == col ) {
-                                        pM->valstruc[row].val[i] += val;
-					found = TRUE;
-					break;
-				}
-			}
-		} else {
-                        mtx_set_diag_val_nt(pM, row,
-                                        mtx_get_diag_val_nt(pM, row) + val);
-			found = TRUE;
-		}
-		if( found == FALSE ) {
-                        if ( err_state_iserror(set_mtx_val_ncolse(pM, row, col,
-                                                                val)) )
-                                err_msg_6(err_CALLBY, "add_mtx_val_ncolse"
-                                        "(%p[%dx%d],%d,%d,%g)", (void *) pM,
+                /* another binary search to find the correct position
+                   in the backset. Note that back_set is ordered
+                   *descendingly*. */
+                low = 0;/* the lowest index that is not known to be too small */
+                high = idxp;/* the lowest index that is known to be too large */
+                /* Often, an element is inserted at the bottom end of the
+                   matrix; therefore, let's first compare with the bottommost
+                   element. */
+                mid = high - 1;
+                while ( low < high ) {
+                        if ( temp_back_set[mid]->col_minus_row > col_minus_row )
+                        {
+                                low = mid + 1;
+                        } else if ( temp_back_set[mid]->col_minus_row <
+                                                                col_minus_row )
+                        {
+                                high = mid;
+                        } else {
+                                /*@-compmempass@*/ /*@-mustdefine@*/
+                                err_msg_6(err_INCONSISTENT, "insert_entry(%p[%d"
+                                        "x%d],%d,%d,%p)", (void *) pM,
                                         mtx_rows(pM), mtx_cols(pM), row, col,
-                                        val, err_ERROR);
-		}
+                                        (void *) result, err_ERROR);
+                                /*@=mustdefine@*/ /*@=compmempass@*/
+                        }
+                        mid = (low + high) / 2;
+                }
+                /* Now let temp_back_set point to the new element of the
+                   back set */
+                temp_back_set = &temp_back_set[high];
+                /* move entries in further rows to the bottom */
+                memmove(temp_back_set + 1, temp_back_set,
+                                        (idxp - high) * sizeof(mtx_el_ptr));
+                *temp_back_set = *result;
         }
         return err_OK;
 }
@@ -430,9 +580,7 @@
 /************************************************************************/
 /************************General Sparse matrix methods*******************/
 /**********For matrices allocated with "allocate_sparse_matrix"**********/
-/*  NOTE: NOT ALL of these methods are applicable to matrices allocated
-	with the "allocate_sparse_matrix_ncolse" method! For example you can not
-	add a new row element with reallocating the row elements array!  */
+/*****************or with "allocate_sparse_matrix_ncolse"****************/
 /*======================================================================*/
 
 /**
@@ -450,7 +598,7 @@
 	int i;
 	int valid_state;
 
-        if ( NULL == pStateSpace || NULL == pQ
+        if ( NULL == pStateSpace || NULL == pQ || ! mtx_is_wrap_or_diag(pQ)
                         || NULL == pValidStates || 0 > pValidStates[0] )
         {
                 err_msg_8(err_PARAM, "initMatrix(%p[%dx%d],%p[%dx%d],%p[%d])",
@@ -494,7 +642,9 @@
 	/*Remove the non zero rows*/
         int i, length;
 
-        if ( NULL == pQ || NULL == pValidStates || 0 > pValidStates[0] ) {
+        if ( NULL == pQ || ! mtx_is_wrap_or_diag(pQ) || NULL == pValidStates
+                                || 0 > pValidStates[0] )
+        {
                 err_msg_5(err_PARAM, "cleanMatrix(%p[%dx%d],%p[%d])",
                                 (void *) pQ, NULL != pQ ? mtx_rows(pQ) : 0,
                                 NULL != pQ ? mtx_cols(pQ) : 0,
@@ -510,11 +660,8 @@
                 valrow->diag = 0.0;
                 if ( 0 != valrow->ncols ) {
                         valrow->ncols = 0;
-                        valrow->col = (int *) NULL;
-                        valrow->val = (double *) NULL;
-                        free(valrow->back_set);
-                        valrow->back_set = (int *) NULL;
-                        valrow->pcols = 0;
+                        valrow->entries = (struct mtx_el *) NULL;
+                        clear_flag(pQ->flags, mtx_flag_BACK_SET_CORRECT);
                 }
 	}
         /*Note: Because of the sparse matrix functions we do not need to
@@ -531,15 +678,24 @@
 err_state free_mtx_sparse(/*@only@*/ /*@i1@*/ /*@null@*/ sparse * pM)
 {
         state_index i, rows;
+        /*@dependent@*/ struct values * valrow;
+
 	if(pM)
 	{
                 rows = mtx_rows(pM);
 		for(i=0;i<rows;i++)
 		{
                         free(pM->valstruc[i].back_set);
-                        free(pM->valstruc[i].col);
-                        free(pM->valstruc[i].val);
 		}
+                if ( ! mtx_is_wrap_or_diag(pM) ) {
+                        valrow = &pM->valstruc[0];
+                        free(valrow->entries);
+                        if ( mtx_is_normal(pM) ) {
+                                for ( i = rows ; 0 < --i ; ) {
+                                        free((++valrow)->entries);
+                                }
+                        }
+                }
 		free(pM);
 	}
         else {
@@ -552,8 +708,6 @@
 
 /**
 * ADD a NEW value to the matrix.
-* WARNING: Should NOT be used with the matrices allocated with
-* allocate_sparse_matrix_ncolse(,,) method!
 * @param	: sparse * pM the matrix to add element to
 * @param	: int row: The row number in which the value is to be set.
 * @param	: int col: The column number in which the value is to be set.
@@ -577,75 +731,22 @@
 
 	if( val != 0 && row != col )
 	{
-                int idx = mtx_next_num(pM, row);
-                /*@null@*/ int *temp_col = (int*) realloc(pM->valstruc[row].col,
-                                (idx + 1) * sizeof(int));
-                /*@null@*/ double * temp_val = NULL;
-                if ( NULL != temp_col ) {
-                        pM->valstruc[row].col = temp_col;
-                        temp_val = (double *) realloc(pM->valstruc[row].val,
-                                (idx + 1) * sizeof(double));
-                        if ( NULL != temp_val ) {
-                                pM->valstruc[row].val = temp_val;
-                                if ( mtx_rows(pM) == mtx_cols(pM) ) {
-                                        int idxp = mtx_prev_num(pM, col);
-                                        /*@null@*/ int * temp_back_set =
-                                                (int *) realloc(
-                                                pM->valstruc[col].back_set,
-                                                (idxp + 1) * sizeof(int));
-                                        if ( NULL != temp_back_set ) {
-                                                pM->valstruc[col].back_set
-                                                        = temp_back_set;
-                                                temp_back_set[idxp] = row;
-                                                pM->valstruc[col].pcols++;
-                                        } else
-                                                temp_col = NULL;
-                                }
-                        }
-                }
-                if ( NULL == temp_col || NULL == temp_val ) {
-                        err_msg_6(err_MEMORY, "set_mtx_val(%p[%dx%d],%d,%d,%g)",
+                struct mtx_el * new_entry;
+                if ( err_state_iserror(insert_entry(pM, row, col, &new_entry)) )
+                {
+                        err_msg_6(err_CALLBY, "set_mtx_val(%p[%dx%d],%d,%d,%g)",
                                         (void *) pM, mtx_rows(pM), mtx_cols(pM),
                                         row, col, val, err_ERROR);
                 }
-                temp_col[idx] = col;
-                temp_val[idx] = val;
-                pM->valstruc[row].ncols++;
+                new_entry->val = val;
 	}
 	else if( val != 0 && row == col )
                 mtx_set_diag_val_nt(pM, row, val);
         return err_OK;
 }
 
-/*****************************************************************************
-name		: free_mtx_wrap
-role		: frees the sparse matrix. Does not clean rows.
-@param		: sparse * pM the matrix to remove
-******************************************************************************/
-err_state free_mtx_wrap(/*@only@*/ /*@i1@*/ /*@null@*/ sparse * pM)
-{
-	if(pM)
-	{
-                values * valrow = &pM->valstruc[0];
-                state_index i = mtx_rows(pM);
-                do {
-                        free(valrow++->back_set);
-                } while ( 0 < --i );
-		free(pM);
-	}
-        else {
-                err_msg_3(err_PARAM, "free_mtx_wrap(%p[%dx%d])",
-                                (void *) pM, NULL != pM ? mtx_rows(pM) : 0,
-                                NULL != pM ? mtx_cols(pM) : 0, err_ERROR);
-        }
-        return err_OK;
-}
-
 /**
 * Add a value to an existing one in the matrix or inserts a new value if the is none.
-* WARNING: Should NOT be used with the matrices allocated with
-* allocate_sparse_matrix_ncolse(,,) method! Unless you are 100% sure that
-* the pM[row,col] value already exists in the matrix!
 * @param	: sparse * pM the matrix
 * @param        : int row: The row number in which the value is to be added.
 * @param	: int col: The column number in which the value is to be added.
@@ -654,8 +755,6 @@
 err_state add_mtx_val(/*@i1@*/ /*@null@*/ sparse * pM,
                 int row, int col, double val)
 {
-        int ncols, i;
-        BOOL found = FALSE;
         /* Correction of an MRMC bug: It is an error if row == mtx_rows(pM) or
            col == mtx_cols(pM). David N. Jansen. */
         if ( NULL == pM || row < 0 || row >= mtx_rows(pM) || col < 0
@@ -672,47 +771,36 @@
         }
 	if ( row != col )
 	{
-                ncols = mtx_next_num(pM, row);
-		for( i=0; i < ncols; i++ )
-		{
-                        if ( pM->valstruc[row].col[i] == col )
-			{
-                                pM->valstruc[row].val[i] += val;
-				found=1;
-				break;
-			}
-		}
+                struct mtx_el * new_entry;
+                if ( err_state_iserror(insert_entry(pM, row, col, &new_entry)) )
+                {
+                        err_msg_6(err_CALLBY, "add_mtx_val(%p[%dx%d],%d,%d,%g)",
+                                        (void *) pM, mtx_rows(pM), mtx_cols(pM),
+                                        row, col, val, err_ERROR);
+                }
+                new_entry->val += val;
 	}
 	else
 	{
                 mtx_set_diag_val_nt(pM, row, mtx_get_diag_val_nt(pM,row) + val);
-		found = 1;
 	}
-	if( found == 0 )
-                if ( err_state_iserror(set_mtx_val(pM, row, col, val)) ) {
-                        err_msg_6(err_CALLBY, "add_mtx_val(%p[%dx%d],%d,%d,%g)",
-                                (void *) pM, mtx_rows(pM), mtx_cols(pM), row,
-                                col, val, err_ERROR);
-                }
         return err_OK;
 }
 
 /*****************************************************************************
 name		: get_mtx_val
 role		: get a value in the matrix. Perform binary search to get a
-		non-diagonal value in the matrix when the number of
-		non-diagonal elements in a row is >4.
+                  non-diagonal value in the matrix.
 @param		: sparse * pM the matrix to get element from
 @param		: int row: The row number in which the value is to be set.
 @param		: int col: The column number in which the value is to be set.
 @param		: double * val (by reference): The value.
-@return         : err_ERROR: fail, err_OK: success
+@return		: int: 0: fail, 1: success
+remark		: return 0: value was not found / row/col does not meet bounds.
 ******************************************************************************/
 err_state get_mtx_val(/*@observer@*/ /*@i1@*/ /*@null@*/ const sparse * pM,
                 int row, int col, /*@out@*/ /*@i1@*/ /*@null@*/ double * val)
 {
-        int i;
-
         /* Correction of an MRMC bug: It is an error if row == mtx_rows(pM) or
            col == mtx_cols(pM). David N. Jansen. */
         if ( NULL == pM || row < 0 || row >= mtx_rows(pM) || col < 0
@@ -728,38 +816,19 @@
 
 	if ( row != col )
 	{
-                int high = mtx_next_num(pM, row) - 1, low = 0;
-                                /* high and low are used for binary search */
-                *val = 0.0;
-		/* If there are <= than 4 elements in a row then we can do a regular search */
-		if(high > 4){
-			while(low < high){
-				i = (low+high)/2;
-                                if ( pM->valstruc[row].col[i] == col ) {
-                                        *val = pM->valstruc[row].val[i];
-					break;
-                                } else if ( col > pM->valstruc[row].col[i] ) {
-					low = i + 1;
-				}else{
-					high = i - 1;
-				}
-			}
-			/* low == high => there may be no required element in the (row,col) cell */
-                        if ( pM->valstruc[row].col[low] == col )
-			{
-                                *val = pM->valstruc[row].val[low];
-			}
-		}else{
-                        const int ncols = mtx_next_num(pM, row);
-			for( i=0; i < ncols; i++ )
-			{
-                                if ( pM->valstruc[row].col[i] == col )
-				{
-                                        *val = pM->valstruc[row].val[i];
-					break;
-				}
-			}
+                BOOL found;
+                const struct mtx_el * entries = pM->valstruc[row].entries;
+                state_index found_idx = find_entry(entries,mtx_next_num(pM,row),
+                                col - row, &found);
+                if ( state_index_ERROR == found_idx ) {
+                        /*@-mustdefine@*/
+                        err_msg_6(err_CALLBY, "get_mtx_val(%p[%dx%d],%d,%d,%p)",
+                                        (const void *) pM, mtx_rows(pM),
+                                        mtx_cols(pM), row, col, (void *) val,
+                                        err_ERROR);
+                        /*@=mustdefine@*/
 		}
+                *val = found ? entries[found_idx].val : 0.0;
 	}
         else
 	{
@@ -853,7 +922,7 @@
 err_state mult_mtx_const(/*@i1@*/ /*@null@*/ sparse * pM,
                 double constant)
 {
-        if ( NULL == pM ) {
+        if ( NULL == pM /* || mtx_is_wrap(pM) */ ) {
                 err_msg_4(err_PARAM, "mult_mtx_const(%p[%dx%d],%g)", (void *)pM,
                                 NULL != pM ? mtx_rows(pM) : 0,
                                 NULL != pM ? mtx_cols(pM) : 0, constant,
@@ -882,7 +951,9 @@
 {
         int i, id, length;
 
-        if ( NULL == pM || NULL == pValidStates || 0 > pValidStates[0] ) {
+        if ( NULL == pM /* || mtx_is_wrap(pM) */ || NULL == pValidStates
+                        || 0 > pValidStates[0] )
+        {
                 err_msg_6(err_PARAM, "mult_mtx_cer_const(%p[%dx%d],%g,%p[%d])",
                                 (void *) pM, NULL != pM ? mtx_rows(pM) : 0,
                                 NULL != pM ? mtx_cols(pM) : 0, constant,
@@ -919,7 +990,7 @@
         int i, id, length;
 	double cons;
 
-        if ( NULL == pM || NULL == constant
+        if ( NULL == pM /* || mtx_is_wrap(pM) */ || NULL == constant
                                 || NULL == pValidStates || 0 > pValidStates[0] )
         {
                 err_msg_7(err_PARAM,
@@ -970,7 +1041,7 @@
         int row_idx, length;
         double * pExitRates;
 
-        if ( NULL == pM || NULL == pRowSums
+        if ( NULL == pM || mtx_is_wrap(pM) || NULL == pRowSums
                                 || NULL == pValidStates || 0 > pValidStates[0] )
         {
                 err_msg_6(err_PARAM, "make_embedded_dtmc_out_of_rate_mtx_vs"
@@ -1054,7 +1125,7 @@
 {
         int row_idx, length;
 
-        if ( NULL == pM || NULL == pValidStates
+        if ( NULL == pM || mtx_is_wrap(pM) || NULL == pValidStates
                                 || 0 > pValidStates[0] || NULL == pExitRates )
         {
                 err_msg_6(err_PARAM, "restore_rate_mtx_out_of_embedded_dtmc_vs"
@@ -1113,7 +1184,7 @@
         int row_idx, length;
         double * pExitRates;
 
-        if ( NULL == pM || NULL == pRowSums ) {
+        if ( NULL == pM || mtx_is_wrap(pM) || NULL == pRowSums ) {
                 err_msg_4(err_PARAM, "make_embedded_dtmc_out_of_rate_mtx_all"
                                 "(%p[%dx%d],%p)", (void *) pM,
                                 NULL != pM ? mtx_rows(pM) : 0,
@@ -1182,7 +1253,7 @@
 {
         int row_idx, length;
 
-        if ( NULL == pM || NULL == pExitRates ) {
+        if ( NULL == pM || mtx_is_wrap(pM) || NULL == pExitRates ) {
                 err_msg_4(err_PARAM, "restore_rate_mtx_out_of_embedded_dtmc_all"
                                 "(%p[%dx%d],%p)", (void *) pM,
                                 NULL != pM ? mtx_rows(pM) : 0,
@@ -1421,11 +1492,11 @@
         /* Correction of an MRMC bug: It is an error if row == mtx_rows(pM).
            David N. Jansen. */
         if ( NULL == to || (unsigned) row >= (unsigned) mtx_rows(to)
+                                || ! mtx_is_wrap_or_diag(to)
                                 || NULL == from
                                 || mtx_rows(to) != mtx_rows(from)
                                 || mtx_cols(to) != mtx_cols(from)
-                                || NULL != to->valstruc[row].col
-                                || NULL != to->valstruc[row].val
+                                || NULL != to->valstruc[row].entries
                                 || 0 != mtx_next_num(to, row) )
         {
                 err_msg_7(err_PARAM,
@@ -1437,13 +1508,11 @@
                                 NULL != from ? mtx_cols(from) : 0, err_ERROR);
         }
 
-        free(to->valstruc[row].back_set);
-        to->valstruc[row].back_set = NULL;
-        to->valstruc[row].pcols = 0;
         mtx_set_diag_val_nt(to, row, mtx_get_diag_val_nt(from, row));
         if ( (to->valstruc[row].ncols = mtx_next_num(from, row)) != 0 ) {
-                to->valstruc[row].col = from->valstruc[row].col;
-                to->valstruc[row].val = from->valstruc[row].val;
+                to->valstruc[row].entries = from->valstruc[row].entries;
+                clear_flag(to->flags,
+                                mtx_flag_BACK_SET_CORRECT | mtx_flag_NORMAL);
         }
         /*@-compmempass@*/ /* pM is a wrapper matrix, but splint thinks it is
                               a normal one. */
@@ -1837,6 +1906,7 @@
 	int i;
 
         if ( NULL == pA || NULL == pDI || NULL == pLU
+                || ! mtx_is_diag(pDI) || ! mtx_is_diag(pLU)
                 || mtx_rows(pA)!=mtx_rows(pDI)||mtx_cols(pA)!=mtx_cols(pDI)
                 || mtx_rows(pA)!=mtx_rows(pLU)||mtx_cols(pA)!=mtx_cols(pLU))
         {
@@ -1891,6 +1961,7 @@
 	int i;
 
         if ( NULL==pA || NULL==pD || NULL==pLU || mtx_rows(pA)!=mtx_rows(pLU)
+                                        || ! mtx_is_diag(pLU)
                                         || mtx_cols(pA)!=mtx_cols(pLU) )
         {
                 /*@-mustdefine@*/
@@ -1940,10 +2011,10 @@
                 /*@i1@*/ /*@null@*/ sparse * pL,
                 /*@i1@*/ /*@null@*/ sparse * pU)
 {
-        BOOL isRFirst;
-	int i, j, theNColsL, theNColsU, col;
+        int i;
 
         if ( NULL == pA || NULL == pD || NULL == pL || NULL == pU
+                      || ! mtx_is_diag(pL) || ! mtx_is_diag(pU)
                       || mtx_rows(pA)!=mtx_rows(pL)||mtx_cols(pA)!=mtx_cols(pL)
                       || mtx_rows(pA)!=mtx_rows(pU)||mtx_cols(pA)!=mtx_cols(pU))
         {
@@ -1962,53 +2033,43 @@
         i = mtx_rows(pA);
         do
 	{
+                state_index split;
+                BOOL found;
+                const char * error_str;
+
                 --i;
 		/*Copy values to the pD vector*/
                 pD[i] = mtx_get_diag_val_nt(pA, i);
 
-		/*Initialize the counter for the number of row elements of the pL and pU matrixes*/
-		theNColsL = 0;
-		theNColsU = 0;
+                split = find_entry(pA->valstruc[i].entries, mtx_next_num(pA, i),
+                                                0, &found);
+                error_str = err_CALLBY;
+                if ( state_index_ERROR == split
+                                || (error_str = err_INCONSISTENT,
+                                    found) )
+                {
+                        err_msg_10(error_str, "split_A_into_D_L_U(%p[%dx%d],"
+                                "%p,%p[%dx%d],%p[%dx%d])", (const void *) pA,
+                                mtx_rows(pA),mtx_cols(pA), (void*)pD, (void*)pL,
+                                mtx_rows(pL), mtx_cols(pL), (void *) pU,
+                                mtx_rows(pU), mtx_cols(pU), err_ERROR);
+                }
 
-                isRFirst = TRUE;
-
                 /*Store possible Low diagonal values. Even if there are no
                 values, then the number of columns will be zero and thus this
                 will not
 		affect the entire process*/
-                pL->valstruc[i].col = pA->valstruc[i].col;
-                pL->valstruc[i].val = pA->valstruc[i].val;
+                pL->valstruc[i].ncols = split;
+		/*Copy values to the pL, pU matrixes*/
+                pL->valstruc[i].entries = pA->valstruc[i].entries;
+                pU->valstruc[i].entries = &pA->valstruc[i].entries[split];
+                pU->valstruc[i].ncols = /*@-compmempass@*/ mtx_next_num(pA, i) 
+                                /*@=compmempass@*/ - split;
 
-		/*Copy values to the pL, pU matrixes*/
-                for ( j = 0 ; j < /*@-compmempass@*/ mtx_next_num(pA, i)
-                                                /*@=compmempass@*/ ; j++ )
-		{
-			/*Get the column index of the next element in a row*/
-                        col = pA->valstruc[i].col[j];
-			/*Check and store element in pL or in pU*/
-			if( col < i )
-			{
-				theNColsL++;
-			}
-			else
-			{
-				/* Case of col > i*/
-				if( isRFirst )
-				{
-					/*Copy to the pU matrix*/
-                                        pU->valstruc[i].col =
-                                                        &pA->valstruc[i].col[j];
-                                        pU->valstruc[i].val =
-                                                        &pA->valstruc[i].val[j];
-					isRFirst = 0;
-				}
-				theNColsU++;
-			}
-		}
-		pL->valstruc[i].ncols = theNColsL;
-		pU->valstruc[i].ncols = theNColsU;
 	}
         while ( 0 < i );
+        clear_flag(pL->flags, mtx_flag_BACK_SET_CORRECT | mtx_flag_NORMAL);
+        clear_flag(pU->flags, mtx_flag_BACK_SET_CORRECT | mtx_flag_NORMAL);
         /*@-compmempass@*/
         return err_OK;
         /*@=compmempass@*/
Index: src/algorithms/random_numbers/rand_num_generator.c
===================================================================
--- src/algorithms/random_numbers/rand_num_generator.c	(revision 404)
+++ src/algorithms/random_numbers/rand_num_generator.c	(revision 403)
@@ -221,7 +221,8 @@
 /**
 * Helper function for choosing non-uniformly distributed random state.
 */
-static int returnDistribStateIndex(const double * distribution, int size) {
+static int returnDistribStateIndex(const struct mtx_el * distribution, int size)
+{
 	/* The cumulative probabilities form already considered states */
 	double cumulative_prob = 0;
 	/* Uniform distributed random number */
@@ -229,7 +230,7 @@
 	int i = 0;
 
 	do {
-		cumulative_prob += distribution[i];
+                cumulative_prob += distribution[i].val;
 		i++;
 	} while( i < size && cumulative_prob < unif_rand );
 
@@ -243,13 +244,13 @@
 * @param prob the related probability distribution
 * @param size the size of the values and prob arrays
 */
-int generateRandNumberDiscrete(const int * values, const double * distribution,
-                int size)
+int generateRandNumberDiscrete(const struct mtx_el * distribution,
+                int size, state_index row)
 {
-	IF_SAFETY( ( values != NULL ) && ( distribution != NULL )  && ( size != 0 ) )
+        IF_SAFETY( NULL != distribution  && 0 < size )
 		IF_SAFETY( pFGenRandNumDiscrete != NULL )
-                        return values[returnDistribStateIndex(distribution,
-                                                size)];
+                        return distribution[returnDistribStateIndex(
+                                distribution, size)].col_minus_row + row;
 		ELSE_SAFETY
 			printf("ERROR: The random-number generator for the discrete distribution is not set.\n");
                         exit(EXIT_FAILURE);
Index: src/lumping/lump.c
===================================================================
--- src/lumping/lump.c	(revision 404)
+++ src/lumping/lump.c	(revision 403)
@@ -65,8 +65,6 @@
                 /*@requires isnull P->first_Sp, P->first_PredCl@*/
                 /*@ensures notnull P->blocks@*/
                 /*@modifies *P->blocks@*/;
-static void quicksort(int * cols, double * vals, int left, int right)
-                /*@modifies *cols, *vals@*/;
 
 /**
 * This internal subroutine finds all predecessors of splitter and treats them
@@ -497,14 +495,6 @@
                         }
 		}
                 end_mtx_walk_row;
-		/* order row by column index */
-                if ( 1 < mtx_next_num(Q1, row) ) {
-                        /*@-type@*/
-                        /* sparse matrix internals used */
-                        quicksort(Q1->valstruc[row].col, Q1->valstruc[row].val,
-                                        0, mtx_next_num(Q1, row) - 1);
-                        /*@=type@*/
-		}
 	}
         end_part_walk_blocks;
         printf("Lumping: The number of partition blocks is %d\n", mtx_rows(Q1));
@@ -514,44 +504,6 @@
 }
 
 
-/*****************************************************************************
-name		: quicksort
-role		: This method sorts the row elements by column index.
-@param		: int *cols: column indices of non-zero elements.
-@param		: int *vals: values of non-zero elements.
-@param		: int left: first index.
-@param		: int right: last index.
-******************************************************************************/
-static void quicksort(int *cols, double* vals, int left, int right)
-{
-	int i, j, x, y;
-	double z;
-
-	i = left; j = right;
-	x = cols[(left+right)/2];
-
-	do {
-		while((cols[i] < x) && (i < right)) i++;
-		while((x < cols[j]) && (j > left)) j--;
-
-		if(i <= j) {
-			y = cols[i];
-			cols[i] = cols[j];
-			cols[j] = y;
-
-			z = vals[i];
-			vals[i] = vals[j];
-			vals[j] = z;
-
-			i++; j--;
-		}
-	} while(i <= j);
-
-	if(left < j) quicksort(cols, vals, left, j);
-	if(i < right) quicksort(cols, vals, i, right);
-}
-
-
 /**
 * This method changes the labelling structure such that it corresponds to the partition.
 * @param	: labelling *labellin: the labelling structure.
