// let c point to the constraint
// of this column
vector<ConstraintLine>::const_iterator c = lines.begin();
- while (c->line != (unsigned int)uncondensed_struct.get_column_numbers()[j])
+ while (c->line != uncondensed_struct.get_column_numbers()[j])
++c;
for (unsigned int q=0; q!=c->entries.size(); ++q)
// let c point to the constraint
// of this column
vector<ConstraintLine>::const_iterator c = lines.begin();
- while (c->line != (unsigned int)uncondensed_struct.get_column_numbers()[j])
+ while (c->line != uncondensed_struct.get_column_numbers()[j])
++c;
for (unsigned int p=0; p!=c->entries.size(); ++p)
for (unsigned int c=0; c<lines.size(); ++c)
distribute[lines[c].line] = static_cast<signed int>(c);
- int n_rows = sparsity.n_rows();
- for (int row=0; row<n_rows; ++row)
+ const unsigned int n_rows = sparsity.n_rows();
+ for (unsigned int row=0; row<n_rows; ++row)
{
if (distribute[row] == -1)
// regular line. loop over cols
for (unsigned int j=sparsity.get_rowstart_indices()[row];
j<sparsity.get_rowstart_indices()[row+1]; ++j)
// end of row reached?
- if (sparsity.get_column_numbers()[j] == -1)
+ if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
{
// this should not happen, since
// we only operate on compressed
for (unsigned int j=sparsity.get_rowstart_indices()[row];
j<sparsity.get_rowstart_indices()[row+1]; ++j)
// end of row reached?
- if (sparsity.get_column_numbers()[j] == -1)
+ if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
{
// this should not happen, since
// we only operate on compressed
// let c point to the constraint
// of this column
vector<ConstraintLine>::const_iterator c = lines.begin();
- while (c->line != (unsigned int)uncondensed.get_column_numbers()[j])
+ while (c->line != uncondensed.get_column_numbers()[j])
++c;
for (unsigned int q=0; q!=c->entries.size(); ++q)
// let c point to the constraint
// of this column
vector<ConstraintLine>::const_iterator c = lines.begin();
- while (c->line != (unsigned int)uncondensed.get_column_numbers()[j]) ++c;
+ while (c->line != uncondensed.get_column_numbers()[j]) ++c;
for (unsigned int p=0; p!=c->entries.size(); ++p)
for (unsigned int q=0; q!=next_constraint->entries.size(); ++q)
for (unsigned int j=sparsity.get_rowstart_indices()[row];
j<sparsity.get_rowstart_indices()[row+1]; ++j)
{
- const int column = sparsity.get_column_numbers()[j];
+ const unsigned int column = sparsity.get_column_numbers()[j];
// end of row reached?
- if (column == -1)
+ if (column == SparsityPattern::invalid_entry)
break;
else
if (distribute[column] != -1)
for (unsigned int j=sparsity.get_rowstart_indices()[row];
j<sparsity.get_rowstart_indices()[row+1]; ++j)
// end of row reached?
- if (sparsity.get_column_numbers()[j] == -1)
+ if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
break;
else
{
// of this row's entries
for (j=sparsity.get_rowstart_indices()[row];
j<sparsity.get_rowstart_indices()[row+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
+ if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
break;
// post-condition after loop:
// coordination, i.e. the number
for (unsigned int i=0; i<last_round_dofs.size(); ++i)
for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
+ if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
break;
else
next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
unsigned int coordination = 0;
for (unsigned int j=sparsity.get_rowstart_indices()[*s];
j<sparsity.get_rowstart_indices()[*s+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
+ if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
break;
else
++coordination;
unsigned int j;
for (j=sparsity.get_rowstart_indices()[row];
j<sparsity.get_rowstart_indices()[row+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
+ if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
break;
// post-condition after loop:
// coordination is now
for (unsigned int i=0; i<last_round_dofs.size(); ++i)
for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
+ if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
break;
else
next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
unsigned int coordination = 0;
for (unsigned int j=sparsity.get_rowstart_indices()[*s];
j<sparsity.get_rowstart_indices()[*s+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
+ if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
break;
else
++coordination;
const unsigned int n_dofs = matrix.m();
const SparsityPattern &sparsity = matrix.get_sparsity_pattern();
const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices();
- const int *sparsity_colnums = sparsity.get_column_numbers();
+ const unsigned int *sparsity_colnums = sparsity.get_column_numbers();
for (; dof != endd; ++dof)
{
- const int dof_number = dof->first;
+ Assert (dof->first >= 0, ExcInternalError());
+
+ const unsigned int dof_number = static_cast<unsigned int>(dof->first);
// for each boundary dof:
// set entries of this line
{
// we need not handle the
// row we have already cleared
- if ((signed int)row == dof_number)
+ if (row == dof_number)
continue;
// check whether the line has
//
// if this row contains an element
// for this dof, *p==dof_number
- const int * p = lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
- &sparsity_colnums[sparsity_rowstart[row+1]],
- dof_number);
+ const unsigned int * p = lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
+ &sparsity_colnums[sparsity_rowstart[row+1]],
+ dof_number);
// check whether this line has
// an entry in the regarding column
// (check for ==dof_number and
====== Things to be done for the lac library ================
=================================================================
-Use unsigned integers for the colnums array in dSMatrixStruct. This
- would enhance safety since colnum=-1 would no longer point to a
- valid address. Mark unused values by numeric_limits<uint>::max
- (or some such constant) instead, which probably has the same
- bit representation as (uint)-1, but is an unsigned, therefore
- eliminating the need for conversions.
-
Make ILU decomposition faster by following the comment in the
innermost loop (about the better use of the index j).
template <typename number, typename inverse_type>
template <typename number2>
-void PreconditionBlockSOR<number,inverse_type>::operator() (Vector<number2> &dst,
+void PreconditionBlockSOR<number,inverse_type>::operator() (Vector<number2> &dst,
const Vector<number2> &src) const
{
Assert(A!=0, ExcNoMatrixGivenToUse());
const SparsityPattern &spars = M.get_sparsity_pattern();
const unsigned int *rowstart = spars.get_rowstart_indices();
- const int *columns = spars.get_column_numbers();
+ const unsigned int *columns = spars.get_column_numbers();
Vector<number2> b_cell(blocksize), x_cell(blocksize);
{
b_cell_row=src(row);
for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
- if ((column=static_cast<unsigned int>(columns[j]))
- < begin_diag_block)
+ if ((column=columns[j]) < begin_diag_block)
b_cell_row -= M.global_entry(j) * dst(column);
b_cell(row_cell)=b_cell_row;
for (unsigned int column_cell=0, column=cell*blocksize;
{
b_cell_row=src(row);
for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
- if ((column=static_cast<unsigned int>(columns[j])) < begin_diag_block)
+ if ((column=columns[j]) < begin_diag_block)
{
b_cell_row -= M.global_entry(j) * dst(column);
}
// pattern of the old matrix!
const unsigned int * const rowstart_indices
= matrix.get_sparsity_pattern().get_rowstart_indices();
- const int * const column_numbers
+ const unsigned int * const column_numbers
= matrix.get_sparsity_pattern().get_column_numbers();
for (unsigned int row=0; row<m(); ++row)
- for (const int * col = &column_numbers[rowstart_indices[row]];
+ for (const unsigned int * col = &column_numbers[rowstart_indices[row]];
col != &column_numbers[rowstart_indices[row+1]]; ++col)
set (row, *col, matrix.global_entry(col-column_numbers));
};
// matrix
const SparsityPattern &sparsity = get_sparsity_pattern();
const unsigned int * const rowstart_indices = sparsity.get_rowstart_indices();
- const int * const column_numbers = sparsity.get_column_numbers();
+ const unsigned int * const column_numbers = sparsity.get_column_numbers();
/*
PSEUDO-ALGORITHM
// let k run over all lower-left
// elements of row i; skip
// diagonal element at start
- const int * first_of_row
+ const unsigned int * first_of_row
= &column_numbers[rowstart_indices[row]+1];
- const int * first_after_diagonal
+ const unsigned int * first_after_diagonal
= lower_bound (&column_numbers[rowstart_indices[row]+1],
&column_numbers[rowstart_indices[row+1]],
- static_cast<int>(row));
+ row);
// k := *col_ptr
- for (const int * col_ptr = first_of_row; col_ptr!=first_after_diagonal; ++col_ptr)
+ for (const unsigned int * col_ptr = first_of_row; col_ptr!=first_after_diagonal; ++col_ptr)
{
const unsigned int global_index_ik = col_ptr-column_numbers;
global_entry(global_index_ik) *= diag_element(*col_ptr);
// row correctly, but excluding
// the main diagonal entry
bool left_of_diagonal = true;
- for (const int * j = col_ptr+1;
+ for (const unsigned int * j = col_ptr+1;
j<&column_numbers[rowstart_indices[row+1]];
++j)
{
//
// note that j should never point
// to the diagonal itself!
- if (left_of_diagonal && (*j > static_cast<int>(row)))
+ if (left_of_diagonal && (*j > row))
{
- Assert (*j != static_cast<int>(row), ExcInternalError());
+ Assert (*j != row, ExcInternalError());
left_of_diagonal = false;
// a[i,i] -= a[i,k]*a[k,i]
Assert (dst.size() == m(), ExcSizeMismatch(dst.size(), m()));
const unsigned int N=dst.size();
- const unsigned int * const rowstart_indices
- = get_sparsity_pattern().get_rowstart_indices();
- const int * const column_numbers
- = get_sparsity_pattern().get_column_numbers();
+ const unsigned int * const rowstart_indices = get_sparsity_pattern().get_rowstart_indices();
+ const unsigned int * const column_numbers = get_sparsity_pattern().get_column_numbers();
// solve LUx=b in two steps:
// first Ly = b, then
// Ux = y
{
// get start of this row. skip the
// diagonal element
- const int * const rowstart = &column_numbers[rowstart_indices[row]+1];
+ const unsigned int * const rowstart = &column_numbers[rowstart_indices[row]+1];
// find the position where the part
// right of the diagonal starts
- const int * const first_after_diagonal
+ const unsigned int * const first_after_diagonal
= lower_bound (rowstart,
&column_numbers[rowstart_indices[row+1]],
- static_cast<int>(row));
+ row);
- for (const int * col=rowstart; col!=first_after_diagonal; ++col)
+ for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col)
dst(row) -= global_entry (col-column_numbers) * dst(*col);
};
for (int row=N-1; row>=0; --row)
{
// get end of this row
- const int * const rowend = &column_numbers[rowstart_indices[row+1]];
+ const unsigned int * const rowend = &column_numbers[rowstart_indices[row+1]];
// find the position where the part
// right of the diagonal starts
- const int * const first_after_diagonal
+ const unsigned int * const first_after_diagonal
= lower_bound (&column_numbers[rowstart_indices[row]+1],
&column_numbers[rowstart_indices[row+1]],
- row);
+ static_cast<unsigned int>(row));
- for (const int * col=first_after_diagonal; col!=rowend; ++col)
+ for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col)
dst(row) -= global_entry (col-column_numbers) * dst(*col);
// scale by the diagonal element.
template <typename number>
inline
-void SparseMatrix<number>::set (const unsigned int i, const unsigned int j,
- const number value) {
+void SparseMatrix<number>::set (const unsigned int i,
+ const unsigned int j,
+ const number value)
+{
Assert (cols != 0, ExcMatrixNotInitialized());
- Assert ((cols->operator()(i,j) != -1) || (value == 0.),
+ Assert ((cols->operator()(i,j) != SparsityPattern::invalid_entry) ||
+ (value == 0.),
ExcInvalidIndex(i,j));
- const int index = cols->operator()(i,j);
-
- if (index >= 0)
- val[index] = value;
+ if (value != 0.)
+ val[cols->operator()(i,j)] = value;
};
void SparseMatrix<number>::add (const unsigned int i, const unsigned int j,
const number value) {
Assert (cols != 0, ExcMatrixNotInitialized());
- Assert ((cols->operator()(i,j) != -1) || (value == 0.),
+ Assert ((cols->operator()(i,j) != SparsityPattern::invalid_entry) ||
+ (value == 0.),
ExcInvalidIndex(i,j));
- const int index = cols->operator()(i,j);
-
- if (index >= 0)
- val[index] += value;
+ if (value != 0.)
+ val[cols->operator()(i,j)] += value;
};
template <typename number>
inline
-number SparseMatrix<number>::operator () (const unsigned int i, const unsigned int j) const {
+number SparseMatrix<number>::operator () (const unsigned int i,
+ const unsigned int j) const
+{
Assert (cols != 0, ExcMatrixNotInitialized());
- Assert (cols->operator()(i,j) != -1,
+ Assert (cols->operator()(i,j) != SparsityPattern::invalid_entry,
ExcInvalidIndex(i,j));
return val[cols->operator()(i,j)];
};
template <typename number>
inline
-number SparseMatrix<number>::diag_element (const unsigned int i) const {
+number SparseMatrix<number>::diag_element (const unsigned int i) const
+{
Assert (cols != 0, ExcMatrixNotInitialized());
Assert (m() == n(), ExcMatrixNotSquare());
Assert (i<max_len, ExcInvalidIndex1(i));
template <typename number>
inline
number
-SparseMatrix<number>::raw_entry(unsigned int row, unsigned int index) const
+SparseMatrix<number>::raw_entry (const unsigned int row,
+ const unsigned int index) const
{
Assert(row<cols->rows, ExcIndexRange(row,0,cols->rows));
Assert(index<cols->row_length(row), ExcIndexRange(index,0,cols->row_length(row)));
template <typename number>
inline
-number SparseMatrix<number>::global_entry (const unsigned int j) const {
+number SparseMatrix<number>::global_entry (const unsigned int j) const
+{
Assert (cols != 0, ExcMatrixNotInitialized());
return val[j];
};
template <typename number>
inline
-number & SparseMatrix<number>::global_entry (const unsigned int j) {
+number & SparseMatrix<number>::global_entry (const unsigned int j)
+{
Assert (cols != 0, ExcMatrixNotInitialized());
return val[j];
};
// if not in MT mode or size<2000
// do it in an oldfashioned way
- const number *val_ptr = &val[cols->rowstart[0]];
- const int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
- somenumber *dst_ptr = &dst(0);
+ const number *val_ptr = &val[cols->rowstart[0]];
+ const unsigned int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
+ somenumber *dst_ptr = &dst(0);
for (unsigned int row=0; row<n_rows; ++row)
{
somenumber s = 0.;
const unsigned int end_row) const
{
#ifdef DEAL_II_USE_MT
- const number *val_ptr = &val[cols->rowstart[begin_row]];
- const int *colnum_ptr = &cols->colnums[cols->rowstart[begin_row]];
- somenumber *dst_ptr = &dst(begin_row);
+ const number *val_ptr = &val[cols->rowstart[begin_row]];
+ const unsigned int *colnum_ptr = &cols->colnums[cols->rowstart[begin_row]];
+ somenumber *dst_ptr = &dst(begin_row);
for (unsigned int row=begin_row; row<end_row; ++row)
{
somenumber s = 0.;
{
for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
- int p = cols->colnums[j];
+ const unsigned int p = cols->colnums[j];
dst(p) += val[j] * src(i);
}
}
Assert(m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size()));
Assert(n() == src.size(), ExcDimensionsDontMatch(n(),src.size()));
- const unsigned int n_rows = m();
- const number *val_ptr = &val[cols->rowstart[0]];
- const int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
- somenumber *dst_ptr = &dst(0);
+ const unsigned int n_rows = m();
+ const number *val_ptr = &val[cols->rowstart[0]];
+ const unsigned int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
+ somenumber *dst_ptr = &dst(0);
for (unsigned int row=0; row<n_rows; ++row)
{
somenumber s = 0.;
{
for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
- int p = cols->colnums[j];
+ const unsigned int p = cols->colnums[j];
dst(p) += val[j] * src(i);
}
}
#endif
// if not in MT mode or the matrix is
// too small: do it one-by-one
- somenumber sum = 0.;
- const number *val_ptr = &val[cols->rowstart[0]];
- const int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
+ somenumber sum = 0.;
+ const number *val_ptr = &val[cols->rowstart[0]];
+ const unsigned int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
for (unsigned int row=0; row<n_rows; ++row)
{
somenumber s = 0.;
{
#ifdef DEAL_II_USE_MT
somenumber sum = 0.;
- const number *val_ptr = &val[cols->rowstart[begin_row]];
- const int *colnum_ptr = &cols->colnums[cols->rowstart[begin_row]];
+ const number *val_ptr = &val[cols->rowstart[begin_row]];
+ const unsignedint *colnum_ptr = &cols->colnums[cols->rowstart[begin_row]];
for (unsigned int row=begin_row; row<end_row; ++row)
{
somenumber s = 0.;
somenumber s = b(i);
for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
- int p = cols->colnums[j];
+ const unsigned int p = cols->colnums[j];
s -= val[j] * u(p);
}
dst(i) = s;
somenumber s = b(i);
for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
- int p = cols->colnums[j];
+ const unsigned int p = cols->colnums[j];
s -= val[j] * u(p);
}
dst(i) = s;
const unsigned int first_right_of_diagonal_index
= (lower_bound (&cols->colnums[*rowstart_ptr+1],
&cols->colnums[*(rowstart_ptr+1)],
- static_cast<signed int>(row)) -
+ row) -
&cols->colnums[0]);
for (unsigned int j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
const unsigned int first_right_of_diagonal_index
= (lower_bound (&cols->colnums[*rowstart_ptr+1],
&cols->colnums[*(rowstart_ptr+1)],
- static_cast<signed int>(row)) -
+ static_cast<unsigned int>(row)) -
&cols->colnums[0]);
for (unsigned int j=first_right_of_diagonal_index; j<*(rowstart_ptr+1); ++j)
- if (cols->colnums[j] > row)
+ if (cols->colnums[j] > static_cast<unsigned int>(row))
*dst_ptr -= om* val[j] * dst(cols->colnums[j]);
*dst_ptr /= val[*rowstart_ptr];
{
somenumber s = dst(row);
for (unsigned int j=cols->rowstart[row]; j<cols->rowstart[row+1]; ++j)
- if (static_cast<unsigned int>(cols->colnums[j]) < row)
+ if (cols->colnums[j] < row)
s -= val[j] * dst(cols->colnums[j]);
dst(row) = s * om / val[cols->rowstart[row]];
Assert (m() == n(), ExcMatrixNotSquare());
Assert (m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size()));
- int p;
const unsigned int n = dst.size();
unsigned int j;
somenumber s;
s = 0.;
for (j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
- p = cols->colnums[j];
- if (p>=0)
+ const unsigned int p = cols->colnums[j];
+ if (p != SparsityPattern::invalid_entry)
{
if (i>j) s += val[j] * dst(p);
}
s = 0.;
for (j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
{
- p = cols->colnums[j];
- if (p>=0)
+ const unsigned int p = cols->colnums[j];
+ if (p != SparsityPattern::invalid_entry)
{
if (static_cast<unsigned int>(i)<j) s += val[j] * dst(p);
}
for (unsigned int i=0; i<m(); ++i)
{
for (unsigned int j=0; j<n(); ++j)
- if ((*cols)(i,j) != -1)
+ if ((*cols)(i,j) != SparsityPattern::invalid_entry)
out << setw(width)
<< val[cols->operator()(i,j)] << ' ';
else
* needed. Then #a_{24}=val[number of the found element] = 3#.
*
*
- * @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth; some documentation by Rasched Zamni
+ * @author Wolfgang Bangerth and others
*/
class SparsityPattern : public Subscriptor
{
public:
+ /**
+ * Define a value which is used
+ * to indicate that a certain
+ * value in the #colnums# array
+ * is unused, i.e. does not
+ * represent a certain column
+ * number index.
+ *
+ * Indices with this invalid
+ * value are used to insert new
+ * entries to the sparsity
+ * pattern using the #add# member
+ * function, and are removed when
+ * calling #compress#.
+ *
+ * You should not assume that the
+ * variable declared here has a
+ * certain value. The
+ * initialization is given here
+ * only to enable the compiler to
+ * perform some optimizations,
+ * but the actual value of the
+ * variable may change over time.
+ */
+ static const unsigned int invalid_entry = static_cast<unsigned int>(-1);
+
/**
* Initialize the matrix empty,
* that is with no memory
* element with row number #i# and
* column number #j#. If the matrix
* element is not a nonzero one,
- * return -1.
+ * return #SparsityPattern::invalid_entry#.
*
* This function is usually called
* by the #operator()# of the
* because the column numbers are
* sorted.
*/
- int operator() (const unsigned int i, const unsigned int j) const;
+ unsigned int operator() (const unsigned int i, const unsigned int j) const;
/**
* Add a nonzero entry to the matrix.
*/
void add (const unsigned int i, const unsigned int j);
- /**
- * This matrix adds a whole connectivity
- * list to the sparsity structure
- * respresented by this object. It assumes
- * the #rowcols# array to be a list of
- * indices which are all linked together,
- * i.e. all entries
- * #(rowcols[i], rowcols[j])# for all
- * #i,j=0...n# for this sparsity pattern
- * are created. #n# is assumed to be the
- * number of elements pointed to by
- * #rowcols#.
- */
- void add_matrix (const unsigned int n, const int* rowcols);
-
- //////////
- void add_matrix (const unsigned int m, const unsigned int n,
- const int* rows, const int* cols);
-
/**
* Print the sparsity of the matrix
* in a format that #gnuplot# understands
* avoid programs relying on outdated
* information!
*/
- const int * get_column_numbers () const;
+ const unsigned int * get_column_numbers () const;
/**
* element, i.e.
* #colnums[rowstart[r]]==r#.
*/
- int *colnums;
+ unsigned int *colnums;
/**
* Store whether the #compress# function
inline
-const int *
+const unsigned int *
SparsityPattern::get_column_numbers () const
{
return colnums;
};
+
inline
unsigned int
SparsityPattern::row_length (const unsigned int row) const
return rowstart[row+1]-rowstart[row];
}
+
+
inline
unsigned int
SparsityPattern::column_number (const unsigned int row,
- const unsigned int index) const
+ const unsigned int index) const
{
Assert(row<rows, ExcIndexRange(row,0,rows));
Assert(index<row_length(row), ExcIndexRange(index,0,row_length(row)));
}
+
inline
unsigned int
SparsityPattern::n_nonzero_elements () const
};
+
SparsityPattern::SparsityPattern (const SparsityPattern &original,
const unsigned int max_per_row,
const unsigned int extra_off_diagonals)
// the rest from that element onwards
// which is not a side-diagonal any
// more.
- const int * const
+ const unsigned int * const
original_row_start = &original.colnums[original.rowstart[row]] + 1;
- // thw following requires that
+ // the following requires that
// #original# be compressed since
- // otherwise there might be -1's
- const int * const
+ // otherwise there might be invalid_entry's
+ const unsigned int * const
original_row_end = &original.colnums[original.rowstart[row+1]];
- const int * const
+ const unsigned int * const
original_last_before_side_diagonals = lower_bound (original_row_start,
original_row_end,
- static_cast<int>
- (row
- -extra_off_diagonals));
- const int * const
+ row-extra_off_diagonals);
+ const unsigned int * const
original_first_after_side_diagonals = upper_bound (original_row_start,
original_row_end,
- static_cast<int>
- (row
- +extra_off_diagonals));
+ row+extra_off_diagonals);
- int * next_free_slot = &colnums[rowstart[row]] + 1;
+ unsigned int * next_free_slot = &colnums[rowstart[row]] + 1;
// copy elements before side-diagonals
next_free_slot = copy (original_row_start,
{
if (colnums) delete[] colnums;
max_vec_len = vec_len;
- colnums = new int[max_vec_len];
+ colnums = new unsigned int[max_vec_len];
};
// set the rowstart array
// preset the column numbers by a
// value indicating it is not in
// use
- fill_n (&colnums[0], vec_len, -1);
+ fill_n (&colnums[0], vec_len, invalid_entry);
// if the matrix is square: let the
// first entry in each row be the
SparsityPattern::compress ()
{
Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
+
+ // do nothing if already compressed
+ if (compressed)
+ return;
- if (compressed) return;
unsigned int next_free_entry = 0,
- next_row_start = 0,
- row_length = 0;
+ next_row_start = 0,
+ row_length = 0;
// first find out how many non-zero
// elements there are, in order to
const unsigned int
nonzero_elements = count_if (&colnums[rowstart[0]],
&colnums[rowstart[rows]],
- bind2nd(not_equal_to<int>(), -1));
+ bind2nd(not_equal_to<unsigned int>(), invalid_entry));
// now allocate the respective memory
- int *new_colnums = new int[nonzero_elements];
+ unsigned int *new_colnums = new unsigned int[nonzero_elements];
// reserve temporary storage to
// store the entries of one row
- vector<int> tmp_entries (max_row_length);
+ vector<unsigned int> tmp_entries (max_row_length);
// Traverse all rows
for (unsigned int line=0; line<rows; ++line)
// first unused entry is reached
row_length = 0;
for (unsigned int j=rowstart[line]; j<rowstart[line+1]; ++j,++row_length)
- if (colnums[j] != -1)
+ if (colnums[j] != invalid_entry)
tmp_entries[row_length] = colnums[j];
else
break;
// some internal checks
Assert ((rows!=cols) ||
- (new_colnums[rowstart[line]] == static_cast<signed int>(line)),
+ (new_colnums[rowstart[line]] == line),
ExcInternalError());
// assert that the first entry
// does not show up in
-int
-SparsityPattern::operator () (const unsigned int i, const unsigned int j) const
+unsigned int
+SparsityPattern::operator () (const unsigned int i,
+ const unsigned int j) const
{
Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
Assert (i<rows, ExcInvalidIndex(i,rows));
// if a first entry exists)
if (rowstart[i] != rowstart[i+1])
{
- if (static_cast<signed int>(j) == colnums[rowstart[i]])
+ if (j == colnums[rowstart[i]])
return rowstart[i];
}
else
// no first entry exists for this
// line
- return -1;
+ return invalid_entry;
// all other entries are sorted, so
// we can use a binary seach algorithm
// at the top of this function, so it
// may not be called for noncompressed
// structures.
- const int * const p = lower_bound (&colnums[rowstart[i]+1],
- &colnums[rowstart[i+1]],
- static_cast<signed int>(j));
- if ((*p == static_cast<signed int>(j)) &&
+ const unsigned int * const p = lower_bound (&colnums[rowstart[i]+1],
+ &colnums[rowstart[i+1]],
+ j);
+ if ((*p == j) &&
(p != &colnums[rowstart[i+1]]))
return (p - &colnums[0]);
else
- return -1;
+ return invalid_entry;
}
+
void
-SparsityPattern::add (const unsigned int i, const unsigned int j)
+SparsityPattern::add (const unsigned int i,
+ const unsigned int j)
{
Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
Assert (i<rows, ExcInvalidIndex(i,rows));
for (unsigned int k=rowstart[i]; k<rowstart[i+1]; k++)
{
// entry already exists
- if (colnums[k] == (signed int)j) return;
+ if (colnums[k] == j) return;
// empty entry found, put new
// entry here
- if (colnums[k] == -1)
+ if (colnums[k] == invalid_entry)
{
colnums[k] = j;
return;
-void
-SparsityPattern::add_matrix (const unsigned int n, const int* rowcols)
-{
- Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
- for (unsigned int i=0; i<n; ++i)
- for (unsigned int j=0; j<n; ++j)
- add(rowcols[i], rowcols[j]);
-}
-
-
-
-void
-SparsityPattern::add_matrix (const unsigned int m, const unsigned int n,
- const int* rows, const int* cols)
-{
- Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
- for (unsigned i=0; i<m; ++i)
- for (unsigned j=0; j<n; ++j)
- add(rows[i], cols[j]);
-}
-
-
-
void
SparsityPattern::print_gnuplot (ostream &out) const
{
Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
for (unsigned int i=0; i<rows; ++i)
for (unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j)
- if (colnums[j]>=0)
+ if (colnums[j] != invalid_entry)
out << i << " " << -colnums[j] << endl;
AssertThrow (out, ExcIO());
unsigned int b=0;
for (unsigned int i=0; i<rows; ++i)
for (unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j)
- if (colnums[j]>=0)
+ if (colnums[j] != invalid_entry)
{
if (static_cast<unsigned int>(abs(static_cast<int>(i-colnums[j]))) > b)
b = abs(static_cast<int>(i-colnums[j]));
projects / dealii-svn.git / commitdiff