<h3>lac</h3>
<ol>
+ <li> <p>
+ New: There are now functions
+ <code class="member">SparsityPattern::symmetrize</code> and
+ <code class="member">SparseMatrix::symmetrize</code> that
+ generate a symmetric matrix from a non-symmetric one.
+ <br>
+ (WB 2000/12/02)
+ </p>
+
<li> <p>
New: almost all classes that store data now have a function
<code class="member">memory_consumption</code> that returns an
void add (const unsigned int i, const unsigned int j,
const number value);
+ /**
+ * Symmetrize the matrix by
+ * forming the mean value between
+ * the existing matrix and its
+ * transpose, $A = \frac 12(A+A^T)$.
+ *
+ * This operation assumes that
+ * the underlying sparsity
+ * pattern represents a symmetric
+ * object. If this is not the
+ * case, then the result of this
+ * operation will not be a
+ * symmetric matrix, since it
+ * only explicitely symmetrizes
+ * by looping over the lower left
+ * triangular part for efficiency
+ * reasons; if there are entries
+ * in the upper right triangle,
+ * then these elements are missed
+ * in the
+ * symmetrization. Symmetrization
+ * of the sparsity pattern can be
+ * obtain by the
+ * @ref{SparsityPattern}@p{::symmetrize}
+ * function.
+ */
+ void symmetrize ();
+
/**
* Copy the given matrix to this
* one. The operation throws an
{};
+
template <typename number>
SparseMatrix<number>::SparseMatrix (const SparseMatrix &m) :
Subscriptor (m),
};
+
template <typename number>
SparseMatrix<number>&
SparseMatrix<number>::operator = (const SparseMatrix<number> &m)
};
+
template <typename number>
SparseMatrix<number>::SparseMatrix (const SparsityPattern &c) :
cols(&c),
};
+
template <typename number>
SparseMatrix<number>::~SparseMatrix ()
{
};
+
template <typename number>
void
SparseMatrix<number>::reinit ()
}
+
template <typename number>
void
SparseMatrix<number>::reinit (const SparsityPattern &sparsity)
};
+
template <typename number>
void
SparseMatrix<number>::clear ()
};
+
template <typename number>
bool
SparseMatrix<number>::empty () const
};
+
template <typename number>
unsigned int
SparseMatrix<number>::n_nonzero_elements () const
};
+
template <typename number>
unsigned int
SparseMatrix<number>::n_actually_nonzero_elements () const
};
+
+template <typename number>
+void
+SparseMatrix<number>::symmetrize ()
+{
+ Assert (cols != 0, ExcMatrixNotInitialized());
+ Assert (cols->rows == cols->cols, ExcMatrixNotSquare());
+
+ const unsigned int n_rows = m();
+ for (unsigned int row=0; row<n_rows; ++row)
+ {
+ // first skip diagonal entry
+ number *val_ptr = &val[cols->rowstart[row]+1];
+ const unsigned int *colnum_ptr = &cols->colnums[cols->rowstart[row]+1];
+ const number *const val_end_of_row = &val[cols->rowstart[row+1]];
+
+ // treat lower left triangle
+ while ((val_ptr != val_end_of_row) && (*colnum_ptr<row))
+ {
+ // compute the mean of this
+ // and the transpose value
+ const number mean_value = (*val_ptr +
+ val[(*cols)(*colnum_ptr,row)]) / 2.0;
+ // set this value and the
+ // transpose one to the
+ // mean
+ *val_ptr = mean_value;
+ set (*colnum_ptr, row, mean_value);
+
+ // advance pointers
+ ++val_ptr;
+ ++colnum_ptr;
+ };
+ };
+};
+
+
+
template <typename number>
template <typename somenumber>
SparseMatrix<number> &
* because the column numbers are
* sorted.
*/
- unsigned 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.
* If the entry already exists, nothing
* bad happens.
*/
- void add (const unsigned int i, const unsigned int j);
+ void add (const unsigned int i,
+ const unsigned int j);
+
+ /**
+ * Make the sparsity pattern
+ * symmetric by adding the
+ * sparsity pattern of the
+ * transpose object.
+ *
+ * This function throws an
+ * exception if the sparsity
+ * pattern does not represent a
+ * square matrix.
+ */
+ void symmetrize ();
/**
* Print the sparsity of the matrix
// wrong: there was not enough space
// in this line
Assert (false, ExcNotEnoughSpace(i, rowstart[i+1]-rowstart[i]));
-}
+};
+
+
+
+void
+SparsityPattern::symmetrize ()
+{
+ Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
+ Assert (compressed==false, ExcMatrixIsCompressed());
+ Assert (rows==cols, ExcNotSquare());
+
+ // 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
+ //
+ // 2. that the @p{add} function can
+ // be called on elements that
+ // already exist without any harm
+ for (unsigned int row=0; row<rows; ++row)
+ for (unsigned int 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
+ if (colnums[k] == invalid_entry)
+ break;
+
+ // otherwise add the
+ // transpose entry
+ add (colnums[k], row);
+ };
+};
projects / dealii-svn.git / commitdiff