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_coldiag : 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;ivalstruc[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;ivalstruc[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.