reinit (original.rows, original.cols, max_per_row);
- // now copy the entries from
- // the other object
+ // now copy the entries from the other object
for (size_type row=0; row<original.rows; ++row)
{
- // copy the elements of this row
- // of the other object
+ // copy the elements of this row of the other object
//
- // note that the first object actually
- // is the main-diagonal element,
+ // note that the first object actually is the main-diagonal element,
// which we need not copy
//
- // we do the copying in two steps:
- // first we note that the elements in
- // @p{original} are sorted, so we may
- // first copy all the elements up to
- // the first side-diagonal one which
- // is to be filled in. then we insert
- // the side-diagonals, finally copy
- // the rest from that element onwards
- // which is not a side-diagonal any
- // more.
+ // we do the copying in two steps: first we note that the elements in
+ // @p{original} are sorted, so we may first copy all the elements up to
+ // the first side-diagonal one which is to be filled in. then we insert
+ // the side-diagonals, finally copy the rest from that element onwards
+ // which is not a side-diagonal any more.
const size_type *const
original_row_start = &original.colnums[original.rowstart[row]] + 1;
- // the following requires that
- // @p{original} be compressed since
+ // the following requires that @p{original} be compressed since
// otherwise there might be invalid_entry's
const size_type *const
original_row_end = &original.colnums[original.rowstart[row+1]];
- // find pointers before and
- // after extra
- // off-diagonals. if at top or
- // bottom of matrix, then set
- // these pointers such that no
- // copying is necessary (see
- // the @p{copy} commands)
+ // find pointers before and after extra off-diagonals. if at top or
+ // bottom of matrix, then set these pointers such that no copying is
+ // necessary (see the @p{copy} commands)
const size_type *const
original_last_before_side_diagonals
= (row > extra_off_diagonals ?
row+extra_off_diagonals) :
original_row_end);
- // find first free slot. the
- // first slot in each row is
- // the diagonal element
+ // find first free slot. the first slot in each row is the diagonal
+ // element
size_type *next_free_slot = &colnums[rowstart[row]] + 1;
// copy elements before side-diagonals
original_row_end,
next_free_slot);
- // this error may happen if the
- // sum of previous elements per row
- // and those of the new diagonals
- // exceeds the maximum number of
- // elements per row given to this
- // constructor
+ // this error may happen if the sum of previous elements per row and
+ // those of the new diagonals exceeds the maximum number of elements per
+ // row given to this constructor
Assert (next_free_slot <= &colnums[rowstart[row+1]],
ExcNotEnoughSpace (0,rowstart[row+1]-rowstart[row]));
};
const unsigned int max_per_row,
const bool)
{
- // simply map this function to the
- // other @p{reinit} function
+ // simply map this function to the other @p{reinit} function
const std::vector<unsigned int> row_lengths (m, max_per_row);
reinit (m, n, row_lengths);
}
const size_type n,
const unsigned int max_per_row)
{
- // simply map this function to the
- // other @p{reinit} function
+ // simply map this function to the other @p{reinit} function
const std::vector<unsigned int> row_lengths (m, max_per_row);
reinit (m, n, row_lengths);
}
rowstart = 0;
colnums = 0;
max_vec_len = max_dim = rows = cols = 0;
- // if dimension is zero: ignore
- // max_per_row
+ // if dimension is zero: ignore max_per_row
max_row_length = 0;
compressed = false;
return;
}
- // first, if the matrix is
- // quadratic, we
- // will have to make sure that each
- // row has at least one entry for
- // the diagonal element. make this
- // more obvious by having a
- // variable which we can query
+ // first, if the matrix is quadratic, we will have to make sure that each
+ // row has at least one entry for the diagonal element. make this more
+ // obvious by having a variable which we can query
store_diagonal_first_in_row = (m == n);
- // find out how many entries we
- // need in the @p{colnums} array. if
- // this number is larger than
- // @p{max_vec_len}, then we will need
- // to reallocate memory
+ // find out how many entries we need in the @p{colnums} array. if this
+ // number is larger than @p{max_vec_len}, then we will need to reallocate
+ // memory
//
- // note that the number of elements
- // per row is bounded by the number
- // of columns
+ // note that the number of elements per row is bounded by the number of
+ // columns
//
std::size_t vec_len = 0;
for (size_type i=0; i<m; ++i)
row_lengths[i]),
n);
- // sometimes, no entries are
- // requested in the matrix (this
- // most often happens when blocks
- // in a block matrix are simply
- // zero). in that case, allocate
- // exactly one element, to have a
- // valid pointer to some memory
+ // sometimes, no entries are requested in the matrix (this most often
+ // happens when blocks in a block matrix are simply zero). in that case,
+ // allocate exactly one element, to have a valid pointer to some memory
if (vec_len == 0)
{
vec_len = 1;
if (store_diagonal_first_in_row && (max_row_length==0) && (m!=0))
max_row_length = 1;
- // allocate memory for the rowstart
- // values, if necessary. even
- // though we re-set the pointers
- // again immediately after deleting
- // their old content, set them to
- // zero in between because the
- // allocation might fail, in which
- // case we get an exception and the
- // destructor of this object will
- // be called -- where we look at
- // the non-nullness of the (now
- // invalid) pointer again and try
- // to delete the memory a second
- // time.
+ // allocate memory for the rowstart values, if necessary. even though we
+ // re-set the pointers again immediately after deleting their old content,
+ // set them to zero in between because the allocation might fail, in which
+ // case we get an exception and the destructor of this object will be called
+ // -- where we look at the non-nullness of the (now invalid) pointer again
+ // and try to delete the memory a second time.
if (rows > max_dim)
{
if (rowstart)
rowstart = new std::size_t[max_dim+1];
}
- // allocate memory for the column
- // numbers if necessary
+ // allocate memory for the column numbers if necessary
if (vec_len > max_vec_len)
{
if (colnums)
((vec_len == 1) && (rowstart[rows] == 0)),
ExcInternalError());
- // preset the column numbers by a
- // value indicating it is not in
- // use
+ // preset the column numbers by a value indicating it is not in use
std::fill_n (&colnums[0], vec_len, invalid_entry);
- // if diagonal elements are
- // special: let the first entry in
- // each row be the diagonal value
+ // if diagonal elements are special: let the first entry in each row be the
+ // diagonal value
if (store_diagonal_first_in_row)
for (size_type i=0; i<rows; i++)
colnums[rowstart[i]] = i;
next_row_start = 0,
row_length = 0;
- // first find out how many non-zero
- // elements there are, in order to
- // allocate the right amount of
- // memory
+ // first find out how many non-zero elements there are, in order to allocate
+ // the right amount of memory
const std::size_t nonzero_elements
= std::count_if (&colnums[rowstart[0]],
&colnums[rowstart[rows]],
size_type *new_colnums = new size_type[nonzero_elements];
- // reserve temporary storage to
- // store the entries of one row
+ // reserve temporary storage to store the entries of one row
std::vector<size_type> tmp_entries (max_row_length);
// Traverse all rows
for (size_type line=0; line<rows; ++line)
{
- // copy used entries, break if
- // first unused entry is reached
+ // copy used entries, break if first unused entry is reached
row_length = 0;
for (size_type j=rowstart[line]; j<rowstart[line+1]; ++j,++row_length)
if (colnums[j] != invalid_entry)
tmp_entries[row_length] = colnums[j];
else
break;
- // now @p{rowstart} is
- // the number of entries in
- // this line
+ // now @p{rowstart} is the number of entries in this line
- // Sort only beginning at the
- // second entry, if optimized
- // storage of diagonal entries
- // is on.
+ // Sort only beginning at the second entry, if optimized storage of
+ // diagonal entries is on.
- // if this line is empty or has
- // only one entry, don't sort
+ // if this line is empty or has only one entry, don't sort
if (row_length > 1)
std::sort ((store_diagonal_first_in_row)
? tmp_entries.begin()+1
: tmp_entries.begin(),
tmp_entries.begin()+row_length);
- // insert column numbers
- // into the new field
+ // insert column numbers into the new field
for (size_type j=0; j<row_length; ++j)
new_colnums[next_free_entry++] = tmp_entries[j];
- // note new start of this and
- // the next row
+ // note new start of this and the next row
rowstart[line] = next_row_start;
next_row_start = next_free_entry;
- // some internal checks: either
- // the matrix is not quadratic,
- // or if it is, then the first
- // element of this row must be
- // the diagonal element
- // (i.e. with column
- // index==line number)
+ // some internal checks: either the matrix is not quadratic, or if it
+ // is, then the first element of this row must be the diagonal element
+ // (i.e. with column index==line number)
Assert ((!store_diagonal_first_in_row) ||
(new_colnums[rowstart[line]] == line),
ExcInternalError());
- // assert that the first entry
- // does not show up in
- // the remaining ones and that
- // the remaining ones are unique
- // among themselves (this handles
- // both cases, quadratic and
- // rectangular matrices)
+ // assert that the first entry does not show up in the remaining ones
+ // and that the remaining ones are unique among themselves (this handles
+ // both cases, quadratic and rectangular matrices)
//
- // the only exception here is
- // if the row contains no
- // entries at all
+ // the only exception here is if the row contains no entries at all
Assert ((rowstart[line] == next_row_start)
||
(std::find (&new_colnums[rowstart[line]+1],
ExcInternalError());
};
- // assert that we have used all
- // allocated space, no more and no
- // less
+ // assert that we have used all allocated space, no more and no less
Assert (next_free_entry == nonzero_elements,
ExcInternalError());
// set iterator-past-the-end
rowstart[rows] = next_row_start;
- // set colnums to the newly
- // allocated array and delete the
- // old one
+ // set colnums to the newly allocated array and delete the old one
delete[] colnums;
colnums = new_colnums;
void
SparsityPattern::copy_from (const CSP &csp)
{
- // first determine row lengths for
- // each row. if the matrix is
- // quadratic, then we might have to
- // add an additional entry for the
- // diagonal, if that is not yet
- // present. as we have to call
- // compress anyway later on, don't
- // bother to check whether that
- // diagonal entry is in a certain
- // row or not
+ // first determine row lengths for each row. if the matrix is quadratic,
+ // then we might have to add an additional entry for the diagonal, if that
+ // is not yet present. as we have to call compress anyway later on, don't
+ // bother to check whether that diagonal entry is in a certain row or not
const bool do_diag_optimize = (csp.n_rows() == csp.n_cols());
std::vector<unsigned int> row_lengths (csp.n_rows());
for (size_type i=0; i<csp.n_rows(); ++i)
}
reinit (csp.n_rows(), csp.n_cols(), row_lengths);
- // now enter all the elements into
- // the matrix, if there are
- // any. note that if the matrix is
- // quadratic, then we already have
- // the diagonal element
+ // now enter all the elements into the matrix, if there are any. note that
+ // if the matrix is quadratic, then we already have the diagonal element
// preallocated
if (n_rows() != 0 && n_cols() != 0)
for (size_type row = 0; row<csp.n_rows(); ++row)
}
}
- // do not need to compress the sparsity
- // pattern since we already have
- // allocated the right amount of data,
- // and the CSP data is sorted, too.
+ // do not need to compress the sparsity pattern since we already have
+ // allocated the right amount of data, and the CSP data is sorted, too.
compressed = true;
}
template <typename number>
void SparsityPattern::copy_from (const FullMatrix<number> &matrix)
{
- // first init with the number of entries
- // per row. if this matrix is square then we
- // also have to allocate memory for the
- // diagonal entry, unless we have already
- // counted it
+ // first init with the number of entries per row. if this matrix is square
+ // then we also have to allocate memory for the diagonal entry, unless we
+ // have already counted it
std::vector<unsigned int> entries_per_row (matrix.m(), 0);
for (size_type row=0; row<matrix.m(); ++row)
{
bool
SparsityPattern::empty () const
{
- // let's try to be on the safe side of
- // life by using multiple possibilities in
- // the check for emptiness... (sorry for
- // this kludge -- emptying matrices and
- // freeing memory was not present in the
- // original implementation and I don't
- // know at how many places I missed
- // something in adding it, so I try to
- // be cautious. wb)
+ // let's try to be on the safe side of life by using multiple possibilities
+ // in the check for emptiness... (sorry for this kludge -- emptying matrices
+ // and freeing memory was not present in the original implementation and I
+ // don't know at how many places I missed something in adding it, so I try
+ // to be cautious. wb)
if ((rowstart==0) || (rows==0) || (cols==0))
{
Assert (rowstart==0, ExcInternalError());
SparsityPattern::size_type
SparsityPattern::max_entries_per_row () const
{
- // if compress() has not yet been
- // called, we can get the maximum
- // number of elements per row using
- // the stored value
+ // if compress() has not yet been called, we can get the maximum number of
+ // elements per row using the stored value
if (!compressed)
return max_row_length;
- // if compress() was called, we
- // use a better algorithm which
- // gives us a sharp bound
+ // if compress() was called, we use a better algorithm which gives us a
+ // sharp bound
size_type m = 0;
- for (size_type i=1; i<rows; ++i)
+ for (size_type i=1; i<=rows; ++i)
m = std::max (m, static_cast<size_type>(rowstart[i]-rowstart[i-1]));
return m;
Assert (j<cols, ExcIndexRange(j,0,cols));
Assert (compressed, ExcNotCompressed());
- // let's see whether there is
- // something in this line
+ // let's see whether there is something in this line
if (rowstart[i] == rowstart[i+1])
return invalid_entry;
- // If special storage of diagonals
- // was requested, we can get the
- // diagonal element faster by this
- // query.
+ // If special storage of diagonals was requested, we can get the diagonal
+ // element faster by this query.
if (store_diagonal_first_in_row && (i==j))
return rowstart[i];
- // all other entries are sorted, so
- // we can use a binary search
- // algorithm
+ // all other entries are sorted, so we can use a binary search algorithm
//
- // note that the entries are only
- // sorted upon compression, so this
- // would fail for non-compressed
- // sparsity patterns; however, that
- // is why the Assertion is at the
- // top of this function, so it may
- // not be called for noncompressed
- // structures.
+ // note that the entries are only sorted upon compression, so this would
+ // fail for non-compressed sparsity patterns; however, that is why the
+ // Assertion is at the top of this function, so it may not be called for
+ // noncompressed structures.
const size_type *sorted_region_start = (store_diagonal_first_in_row ?
&colnums[rowstart[i]+1] :
&colnums[rowstart[i]]);
Assert (global_index < n_nonzero_elements(),
ExcIndexRange (global_index, 0, n_nonzero_elements()));
- // first find the row in which the
- // entry is located. for this note
- // that the rowstart array indexes
- // the global indices at which each
- // row starts. since it is sorted,
- // and since there is an element
- // for the one-past-last row, we
- // can simply use a bisection
- // search on it
+ // first find the row in which the entry is located. for this note that the
+ // rowstart array indexes the global indices at which each row starts. since
+ // it is sorted, and since there is an element for the one-past-last row, we
+ // can simply use a bisection search on it
const size_type row
= (std::upper_bound (&rowstart[0], &rowstart[rows], global_index)
- &rowstart[0] - 1);
- // now, the column index is simple
- // since that is what the colnums
- // array stores:
+ // now, the column index is simple since that is what the colnums array
+ // stores:
const size_type col = colnums[global_index];
// so return the respective pair
{
Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
Assert (compressed==false, ExcMatrixIsCompressed());
- // Note that we only require a
- // quadratic matrix here, no
- // special treatment of diagonals
+ // Note that we only require a quadratic matrix here, no special treatment
+ // of diagonals
Assert (rows==cols, ExcNotQuadratic());
- // loop over all elements presently
- // in the sparsity pattern and add
- // the transpose element. note:
+ // loop over all elements presently in the sparsity pattern and add the
+ // transpose element. note:
//
- // 1. that the sparsity pattern
- // changes which we work on, but
- // not the present row
+ // 1. that the sparsity pattern changes which we work on, but not the
+ // present row
//
- // 2. that the @p{add} function can
- // be called on elements that
- // already exist without any harm
+ // 2. that the @p{add} function can be called on elements that already exist
+ // without any harm
for (size_type row=0; row<rows; ++row)
for (size_type k=rowstart[row]; k<rowstart[row+1]; k++)
{
- // check whether we are at
- // the end of the entries of
- // this row. if so, go to
- // next row
+ // check whether we are at the end of the entries of this row. if so,
+ // go to next row
if (colnums[k] == invalid_entry)
break;
- // otherwise add the
- // transpose entry if this is
- // not the diagonal (that
- // would not harm, only take
- // time to check up)
+ // otherwise add the transpose entry if this is not the diagonal (that
+ // would not harm, only take time to check up)
if (colnums[k] != row)
add (colnums[k], row);
};
for (size_type i=0; i<rows; ++i)
for (size_type j=rowstart[i]; j<rowstart[i+1]; ++j)
if (colnums[j] != invalid_entry)
- // while matrix entries are
- // usually written (i,j),
- // with i vertical and j
- // horizontal, gnuplot output
- // is x-y, that is we have to
- // exchange the order of
- // output
+ // while matrix entries are usually written (i,j), with i vertical and
+ // j horizontal, gnuplot output is x-y, that is we have to exchange
+ // the order of output
out << colnums[j] << " " << -static_cast<signed int>(i) << std::endl;
AssertThrow (out, ExcIO());
b = std::abs(static_cast<signed int>(i-colnums[j]));
}
else
- // leave if at the end of
- // the entries of this line
+ // leave if at the end of the entries of this line
break;
return b;
}
{
AssertThrow (out, ExcIO());
- // first the simple objects,
- // bracketed in [...]
+ // first the simple objects, bracketed in [...]
out << '[' << max_dim << ' '
<< rows << ' '
<< cols << ' '
projects / dealii-svn.git / commitdiff