From: bangerth Date: Mon, 6 Aug 2007 21:23:32 +0000 (+0000) Subject: Comment on the CompressedSetSparsityPattern issue. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=31bb64ef51d0359d46106e7da9a512150603ce5a;p=dealii-svn.git Comment on the CompressedSetSparsityPattern issue. git-svn-id: https://svn.dealii.org/trunk@14911 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-27/doc/intro.dox b/deal.II/examples/step-27/doc/intro.dox index 0b6185294a..d0c2c8be48 100644 --- a/deal.II/examples/step-27/doc/intro.dox +++ b/deal.II/examples/step-27/doc/intro.dox @@ -187,7 +187,7 @@ the same time. -

Assembling matrices and vectors with $hp$ objects

+

Assembling matrices and vectors with hp objects

Following this, we have to set up matrices and vectors for the linear system of the correct size and assemble them. Setting them up works in exactly the @@ -268,7 +268,7 @@ usual fashion to assemble local contributions. -

A simple indicator for $hp$ refinement and estimating smoothness

+

A simple indicator for hp refinement and estimating smoothness

One of the central pieces of the adaptive finite element method is that we inspect the computed solution (a posteriori) with an indicator that tells us @@ -505,3 +505,66 @@ $\beta$, the formula above gives us a mean to calculate the value of the exponent $\mu$ that we can then use to determine that $\hat u(\hat x)$ is in $H^s(\hat K)$ with $s=\mu-\frac d2$. + + +

Complications with linear systems for hp discretizations

+ +One of the problems with $hp$ methods is that the high polynomial degree of +shape functions together with the large number of constrained degrees of +freedom leads to matrices with large numbers of nonzero entries in some +rows. At the same time, because there are areas where we use low polynomial +degree and consequently matrix rows with relatively few nonzero +entries. Consequently, allocating the sparsity pattern for these matrices is a +challenge. + +Most programs built on deal.II use the DoFTools::make_sparsity_pattern +function to allocate the sparsity pattern of a matrix, and later add a few +more entries necessary to handle constrained degrees of freedom using +ConstraintMatrix::condense. The sparsity pattern is then compressed using +SparsityPattern::compress. This method is explained in step-6 and used in +most tutorial programs. In order to work, it needs an initial upper estimate +for the maximal number of nonzero entries per row, something that can be had +from the DoFHandler::max_couplings_between_dofs function. This is necessary +due to the data structure used in the SparsityPattern class. + +Unfortunately, DoFHandler::max_couplings_between_dofs is unable to produce an +efficient upper estimate in 3d and for higher order elements. If used in these +situations, it therefore leads the SparsityPattern class to allocate much too +much memory, almost all of which will be released again when we call +SparsityPattern::compress. This deficiency, caused by the fact that +DoFHandler::max_couplings_between_dofs must produce a single number for the +maximal number of elements per row even though most rows will be significantly +shorter, can be so severe that the initial memory allocation for the +SparsityPattern exceeds the actual need by a factor of 10 or larger, and can +lead to a program running out of memory when in fact there would be plenty of +memory for all computations. + +A solution to the problem has already been discussed in @ref step_11 "step-11" +and @ref step_18 "step-18". It used an intermediate object of type +CompressedSparsityPattern. This class uses a different memory storage scheme +that is optimized to creating a sparsity pattern when maximal numbers +of entries per row are not accurately available, but is unsuitable for use as +the sparsity pattern actually underlying a sparse matrix. After building the +intermediate object, it is therefore copied into a true SparsityPattern +object, something that can be done very efficient and without having to +over-allocate memory. Typical code doing this is shown in the documentation of +the CompressedSparsityPattern class. This solution is slower than directly +building a SparsityPattern object, but only uses as much memory as is really +necessary. + +As it now turns out, the storage format used in the CompressedSparsityPattern +class is not very good for matrices with truly large numbers of entries per +row — where truly large numbers mean in the hundreds. This isn't +typically the case for lower order elements even in 3d, but happens for high +order elements in 3d; for example, a vertex degree of freedom of a $Q_5$ +element in 3d may couple to as many as 1700 other degrees of freedom. In such +a CompressedSparsityPattern will work, but by tuning the memory storage format +used internally in that class a bit will make it work several times +faster. This is what we did with the CompressedSetSparsityPattern class +— it has exactly the same interface as the CompressedSparsityPattern +class but internally stores things somewhat differently. For most cases, there +is not much of a difference in performance in the classes (though the old +class has a slight advantage for lower order elements in 3d), but for high +order and $hp$ elements in 3d, the CompressedSetSparsityPattern has a definite +edge. We will therefore use it later when we build the sparsity pattern in +this tutorial program.