Since it's been brought back to my attention, I'm going to try and fix
the Laspack post-init matrix addition bug. I'm needing to refactor
the sparsity pattern creation at the same time... and I noticed an
inconsistency in the way n_nz and n_oz were treated, particulary in
the case where user-added sparsity entries exist. So my first patch
just does the necessary factoring and fixes that inconsistency... but
since we don't actually *have* any user-added sparsity entries in any
of my apps or our examples, I don't trust my changes yet.
Would one of the INL guys mind testing the attached patch?
Thanks,
---
Roy
Index: include/base/dof_map.h
===================================================================
--- include/base/dof_map.h (revision 5966)
+++ include/base/dof_map.h (working copy)
@@ -22,6 +22,7 @@
// Local Includes -----------------------------------
#include "libmesh_common.h"
+#include "auto_ptr.h"
#include "enum_order.h"
#include "reference_counted_object.h"
#include "libmesh.h" // libMesh::invalid_uint
@@ -201,13 +202,20 @@
void distribute_dofs (MeshBase&);
/**
- * Computes the sparsity pattern for the matrix corresponding
- * to \p proc_id. Produces data that can be fed to Petsc for
+ * Computes the sparsity pattern for the matrices corresponding to
+ * \p proc_id and sends that data to Linear Algebra packages for
* preallocation of sparse matrices.
*/
void compute_sparsity (const MeshBase&);
/**
+ * Computes the sparsity pattern for the matrices corresponding to
+ * \p proc_id and sends that data to Linear Algebra packages for
+ * preallocation of a new sparse matrix.
+ */
+ void recompute_sparsity (const MeshBase&, SparseMatrix<Number>* m);
+
+ /**
* Attach an object to use to populate the
* sparsity pattern with extra entries.
*
@@ -828,6 +836,11 @@
private:
/**
+ * Builds a sparsity pattern
+ */
+ AutoPtr<SparsityPattern::Build> build_sparsity(const MeshBase& mesh) const;
+
+ /**
* Invalidates all active DofObject dofs for this system
*/
void invalidate_dofs(MeshBase& mesh) const;
Index: src/base/dof_map.C
===================================================================
--- src/base/dof_map.C (revision 5966)
+++ src/base/dof_map.C (working copy)
@@ -44,12 +44,93 @@
#include "threads.h"
+
namespace libMesh
{
// ------------------------------------------------------------
// DofMap member functions
+AutoPtr<SparsityPattern::Build> DofMap::build_sparsity
+ (const MeshBase& mesh) const
+{
+ libmesh_assert (mesh.is_prepared());
+ libmesh_assert (this->n_variables());
+ START_LOG("compute_sparsity()", "DofMap");
+
+ // Compute the sparsity structure of the global matrix. This can be
+ // fed into a PetscMatrix to allocate exacly the number of nonzeros
+ // necessary to store the matrix. This algorithm should be linear
+ // in the (# of elements)*(# nodes per element)
+
+ // We can be more efficient in the threaded sparsity pattern assembly
+ // if we don't need the exact pattern. For some sparse matrix formats
+ // a good upper bound will suffice.
+ bool need_full_sparsity_pattern=false;
+ std::vector<SparseMatrix<Number>* >::const_iterator
+ pos = this->_matrices.begin(),
+ end = this->_matrices.end();
+
+ for (; pos != end; ++pos)
+ if ((*pos)->need_full_sparsity_pattern())
+ need_full_sparsity_pattern = true;
+
+ // See if we need to include sparsity pattern entries for coupling
+ // between neighbor dofs
+ bool implicit_neighbor_dofs = this->use_coupled_neighbor_dofs(mesh);
+
+ // We can compute the sparsity pattern in parallel on multiple
+ // threads. The goal is for each thread to compute the full sparsity
+ // pattern for a subset of elements. These sparsity patterns can
+ // be efficiently merged in the SparsityPattern::Build::join()
+ // method, especially if there is not too much overlap between them.
+ // Even better, if the full sparsity pattern is not needed then
+ // the number of nonzeros per row can be estimated from the
+ // sparsity patterns created on each thread.
+ AutoPtr<SparsityPattern::Build> sp
+ (new SparsityPattern::Build (mesh,
+ *this,
+ this->_dof_coupling,
+ implicit_neighbor_dofs,
+ need_full_sparsity_pattern));
+
+ Threads::parallel_reduce (ConstElemRange (mesh.active_elements_begin(),
+ mesh.active_elements_end()), *sp);
+
+#ifndef NDEBUG
+ // Avoid declaring these variables unless asserts are enabled.
+ const unsigned int proc_id = mesh.processor_id();
+ const unsigned int n_dofs_on_proc = this->n_dofs_on_processor(proc_id);
+#endif
+ libmesh_assert (sp->sparsity_pattern.size() == n_dofs_on_proc);
+
+ STOP_LOG("compute_sparsity()", "DofMap");
+
+ // Check to see if we have any extra stuff to add to the sparsity_pattern
+ if (_extra_sparsity_function)
+ {
+ if (_augment_sparsity_pattern)
+ {
+ libmesh_here();
+ libMesh::out << "WARNING: You have specified both an extra sparsity function and object.\n"
+ << " Are you sure this is what you meant to do??"
+ << std::endl;
+ }
+
+ _extra_sparsity_function
+ (sp->sparsity_pattern, sp->n_nz,
+ sp->n_oz, _extra_sparsity_context);
+ }
+
+ if (_augment_sparsity_pattern)
+ _augment_sparsity_pattern->augment_sparsity_pattern
+ (sp->sparsity_pattern, sp->n_nz, sp->n_oz);
+
+ return sp;
+}
+
+
+
DofMap::DofMap(const unsigned int number) :
_dof_coupling(NULL),
_variables(),
@@ -1312,89 +1393,31 @@
void DofMap::compute_sparsity(const MeshBase& mesh)
{
- libmesh_assert (mesh.is_prepared());
- libmesh_assert (this->n_variables());
+ AutoPtr<SparsityPattern::Build> sp = this->build_sparsity(mesh);
- START_LOG("compute_sparsity()", "DofMap");
-
- // Compute the sparsity structure of the global matrix. This can be
- // fed into a PetscMatrix to allocate exacly the number of nonzeros
- // necessary to store the matrix. This algorithm should be linear
- // in the (# of elements)*(# nodes per element)
-
- // We can be more efficient in the threaded sparsity pattern assembly
- // if we don't need the exact pattern. For some sparse matrix formats
- // a good upper bound will suffice.
- bool need_full_sparsity_pattern=false;
+ // It is possible that some \p SparseMatrix implementations want to
+ // see it. Let them see it before we throw it away.
std::vector<SparseMatrix<Number>* >::const_iterator
pos = _matrices.begin(),
end = _matrices.end();
for (; pos != end; ++pos)
- if ((*pos)->need_full_sparsity_pattern())
- need_full_sparsity_pattern = true;
+ (*pos)->update_sparsity_pattern (sp->sparsity_pattern);
- // See if we need to include sparsity pattern entries for coupling
- // between neighbor dofs
- bool implicit_neighbor_dofs = this->use_coupled_neighbor_dofs(mesh);
-
- // We can compute the sparsity pattern in parallel on multiple
- // threads. The goal is for each thread to compute the full sparsity
- // pattern for a subset of elements. These sparsity patterns can
- // be efficiently merged in the SparsityPattern::Build::join()
- // method, especially if there is not too much overlap between them.
- // Even better, if the full sparsity pattern is not needed then
- // the number of nonzeros per row can be estimated from the
- // sparsity patterns created on each thread.
- SparsityPattern::Build sp (mesh,
- *this,
- _dof_coupling,
- implicit_neighbor_dofs,
- need_full_sparsity_pattern);
-
- Threads::parallel_reduce (ConstElemRange (mesh.active_elements_begin(),
- mesh.active_elements_end()), sp);
-
-#ifndef NDEBUG
- // Avoid declaring these variables unless asserts are enabled.
- const unsigned int proc_id = mesh.processor_id();
- const unsigned int n_dofs_on_proc = this->n_dofs_on_processor(proc_id);
-#endif
- libmesh_assert (sp.sparsity_pattern.size() == n_dofs_on_proc);
-
// steal the n_nz and n_oz arrays from sp -- it won't need them any more,
// and this is more efficient than copying them.
- _n_nz.swap(sp.n_nz);
- _n_oz.swap(sp.n_oz);
+ _n_nz.swap(sp->n_nz);
+ _n_oz.swap(sp->n_oz);
- STOP_LOG("compute_sparsity()", "DofMap");
+}
- // Check to see if we have any extra stuff to add to the sparsity_pattern
- if (_extra_sparsity_function)
- {
- if (_augment_sparsity_pattern)
- {
- libmesh_here();
- libMesh::out << "WARNING: You have specified both an extra sparsity function and object.\n"
- << " Are you sure this is what you meant to do??"
- << std::endl;
- }
- _extra_sparsity_function(sp.sparsity_pattern, _n_nz, _n_oz, _extra_sparsity_context);
- }
- if (_augment_sparsity_pattern)
- _augment_sparsity_pattern->augment_sparsity_pattern (sp.sparsity_pattern, _n_nz, _n_oz);
+void DofMap::recompute_sparsity(const MeshBase& mesh, SparseMatrix<Number>* m)
+{
+ AutoPtr<SparsityPattern::Build> sp = this->build_sparsity(mesh);
- // We are done with the sparsity_pattern. However, quite a
- // lot has gone into computing it. It is possible that some
- // \p SparseMatrix implementations want to see it. Let them
- // see it before we throw it away.
- pos = _matrices.begin();
- end = _matrices.end();
-
- for (; pos != end; ++pos)
- (*pos)->update_sparsity_pattern (sp.sparsity_pattern);
+ m->update_sparsity_pattern (sp->sparsity_pattern);
}
------------------------------------------------------------------------------
Live Security Virtual Conference
Exclusive live event will cover all the ways today's security and
threat landscape has changed and how IT managers can respond. Discussions
will include endpoint security, mobile security and the latest in malware
threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
