#include <base/subscriptor.h>
#include <base/smartpointer.h>
-#include <utility>
+#include <vector>
//forward declarations
template <typename number> class SparseMatrix;
class ostream;
-//TODO: Documentation on full rows
+
/**
* Structure representing the sparsity pattern of a sparse matrix.
* Initialize a rectangular matrix with
* #m# rows and #n# columns.
* The matrix may contain at most #max_per_row#
- * nonzero entries per row, with the exception
- * of the first #n_long_rows# rows. These may even
- * be full.
+ * nonzero entries per row.
*/
SparseMatrixStruct (const unsigned int m,
const unsigned int n,
- const unsigned int max_per_row,
- const unsigned int n_long_rows = 0);
+ const unsigned int max_per_row);
+
+ /**
+ * Initialize a rectangular
+ * matrix with #m# rows and #n#
+ * columns. The maximal number
+ * of nonzero entries for each
+ * row is given by the
+ * #row_lengths# array.
+ */
+ SparseMatrixStruct (const unsigned int m,
+ const unsigned int n,
+ const vector<unsigned int> &row_lengths);
/**
* Initialize a square matrix of dimension
- * #n# with at most #max_per_row#. There are
- * no full rows with this constructor
+ * #n# with at most #max_per_row#
+ * nonzero entries per row.
*/
SparseMatrixStruct (const unsigned int n,
const unsigned int max_per_row);
+ /**
+ * Initialize a square
+ * matrix with #m# rows and #m#
+ * columns. The maximal number
+ * of nonzero entries for each
+ * row is given by the
+ * #row_lengths# array.
+ */
+ SparseMatrixStruct (const unsigned int m,
+ const vector<unsigned int> &row_lengths);
+
/**
* Copy operator. For this the same holds
* as for the copy constructor: it is
* with at most #max_per_row#
* nonzero entries per row.
*
- * This matrix also allows for
- * the first rows to have up to
- * #n# entries, if #n_long_rows#
- * is different from zero. This
- * is useful in the case of
- * globally coupling degrees of
- * freedom.
+ * This function simply maps its
+ * operations to the other
+ * #reinit# function.
+ */
+ void reinit (const unsigned int m,
+ const unsigned int n,
+ const unsigned int max_per_row);
+
+ /**
+ * Reallocate memory for a matrix
+ * of size #m \times n#. The
+ * number of entries for each row
+ * is taken from the array
+ * #row_lengths# which has to
+ * give this number of each row
+ * #i=1...m#.
*
* If #m*n==0# all memory is freed,
* resulting in a total reinitialization
* memory is only allocated if the new
* size extends the old one. This is done
* to save time and to avoid fragmentation
- * of the heap. */
- void reinit (const unsigned int m,
- const unsigned int n,
- const unsigned int max_per_row,
- const unsigned int n_long_rows = 0);
-
+ * of the heap.
+ */
+ void reinit (const unsigned int m,
+ const unsigned int n,
+ const vector<unsigned int> &row_lengths);
+
/**
* This function compresses the sparsity
* structure that this object represents.
/**
* Maximum number of elements per
- * row. This is set to the value given
- * to the #reinit# function (or to the
- * constructor). Its value is more
- * or less meaningsless after #compress()#
- * has been called.
+ * row. This is set to the value
+ * given to the #reinit# function
+ * (or to the constructor), or to
+ * the maximum row length
+ * computed from the vectors in
+ * case the more flexible
+ * constructors or reinit
+ * versions are called. Its value
+ * is more or less meaningsless
+ * after #compress()# has been
+ * called.
*/
- unsigned int max_row_len;
+ unsigned int max_row_length;
- /**
- * Number of full rows. This
- * value only has a meaning
- * before compression and allows
- * the first rows to have more
- * elements than #max_per_row#
- */
- unsigned int n_long_rows;
-
/**
* Array which hold for each row which
* is the first element in #colnums#
* allocated memory may be larger than
* the region that is used. The actual
* number of elements that was allocated
- * is stored in #max_dim#. */
+ * is stored in #max_dim#.
+ */
unsigned int *rowstart;
/**
SparseMatrixStruct::SparseMatrixStruct (const unsigned int m,
const unsigned int n,
- const unsigned int max_per_row,
- const unsigned int long_rows)
+ const unsigned int max_per_row)
: max_dim(0),
max_vec_len(0),
rowstart(0),
colnums(0)
{
- reinit (m,n,max_per_row,long_rows);
+ reinit (m,n,max_per_row);
+};
+
+
+
+SparseMatrixStruct::SparseMatrixStruct (const unsigned int m,
+ const unsigned int n,
+ const vector<unsigned int> &row_lengths)
+ : max_dim(0),
+ max_vec_len(0),
+ rowstart(0),
+ colnums(0)
+{
+ reinit (m, n, row_lengths);
};
rowstart(0),
colnums(0)
{
- reinit (n,n,max_per_row,0);
+ reinit (n,n,max_per_row);
+};
+
+
+
+SparseMatrixStruct::SparseMatrixStruct (const unsigned int m,
+ const vector<unsigned int> &row_lengths)
+ : max_dim(0),
+ max_vec_len(0),
+ rowstart(0),
+ colnums(0)
+{
+ reinit (m, m, row_lengths);
};
void
SparseMatrixStruct::reinit (const unsigned int m,
const unsigned int n,
- const unsigned int max_per_row,
- const unsigned int long_rows)
+ const unsigned int max_per_row)
{
+ // simply map this function to the
+ // other #reinit# function
+ vector<unsigned int> row_lengths (m, max_per_row);
+ reinit (m, n, row_lengths);
+};
+
+
+
+void
+SparseMatrixStruct::reinit (const unsigned int m,
+ const unsigned int n,
+ const vector<unsigned int> &row_lengths)
+{
+ Assert (((m==0) && (n==0)) || (*max_element(row_lengths.begin(), row_lengths.end()) > 0),
+ ExcInvalidNumber(*max_element(row_lengths.begin(), row_lengths.end())));
+ Assert (row_lengths.size() == m, ExcInvalidNumber (m));
+
rows = m;
cols = n;
- vec_len = m * max_per_row + (n-max_per_row) * long_rows;
- max_row_len = max_per_row;
- n_long_rows = long_rows;
+ vec_len = accumulate (row_lengths.begin(), row_lengths.end(), 0);
+ max_row_length = (row_lengths.size() == 0 ?
+ 0 :
+ *max_element(row_lengths.begin(), row_lengths.end()));
// delete empty matrices
if ((m==0) || (n==0))
if (colnums) delete[] colnums;
rowstart = 0;
colnums = 0;
- max_vec_len = vec_len = max_dim = rows = cols = max_row_len = n_long_rows = 0;
+ max_vec_len = vec_len = max_dim = rows = cols = 0;
+ // if dimension is zero: ignore
+ // max_per_row
+ max_row_length = 0;
compressed = false;
return;
};
-
+
+ // allocate memory for the rowstart
+ // values, if necessary
if (rows > max_dim)
{
if (rowstart) delete[] rowstart;
max_dim = rows;
rowstart = new unsigned int[max_dim+1];
};
-
+
+ // allocate memory for the column
+ // numbers if necessary
if (vec_len > max_vec_len)
{
if (colnums) delete[] colnums;
max_vec_len = vec_len;
colnums = new int[max_vec_len];
};
-
- unsigned int start = 0;
- for (unsigned int i=0; i<=rows; i++)
- {
- rowstart[i] = start;
- start += (i<n_long_rows) ? n : max_per_row;
- }
-
+
+ // set the rowstart array
+ rowstart[0] = 0;
+ for (unsigned int i=1; i<=rows; i++)
+ rowstart[i] = rowstart[i-1]+row_lengths[i-1];
+ Assert (rowstart[rows]==vec_len, ExcInternalError());
+
+ // preset the column numbers by a
+ // value indicating it is not in
+ // use
fill_n (&colnums[0], vec_len, -1);
+ // if the matrix is square: let the
+ // first entry in each row be the
+ // diagonal value
if (rows == cols)
for (unsigned int i=0;i<rows;i++)
colnums[rowstart[i]] = i;
// reserve temporary storage to
// store the entries of one row
- int safe_row_len = (n_long_rows != 0) ? n_cols() : max_row_len;
- int *tmp_entries = new int[safe_row_len];
-
+ vector<int> tmp_entries (max_row_length);
+
// Traverse all rows
for (unsigned int line=0; line<rows; ++line)
{
vec_len = rowstart[rows] = next_row_start;
compressed = true;
-
- delete[] tmp_entries;
};
// number of elements per row using
// the stored value
if (!compressed)
- return max_row_len;
+ return max_row_length;
// if compress() was called, we
// use a better algorithm which