From: Wolfgang Bangerth Date: Fri, 18 Apr 2025 20:00:20 +0000 (-0600) Subject: Make it possible to create sparse outputs from MatrixOut. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5abd7a6ca79be59bcfb13e9f817dddff9465eda6;p=dealii.git Make it possible to create sparse outputs from MatrixOut. --- diff --git a/include/deal.II/lac/matrix_out.h b/include/deal.II/lac/matrix_out.h index ee9b1d5f05..f82ac2f1b0 100644 --- a/include/deal.II/lac/matrix_out.h +++ b/include/deal.II/lac/matrix_out.h @@ -106,17 +106,50 @@ public: unsigned int block_size; /** - * If true, plot discontinuous patches, one for each entry. + * By default, this class shows each matrix entry as one value, + * and then produces a bilinear display of all of these + * values. This results in a continuous plot on a mesh with + * $N\times N$ vertices for a matrix of size $N$. If, on the other + * hand, the current variable is set to `true`, then this class + * instead creates a plot in which each matrix entry corresponds + * to a single patch on which the value of the matrix entry is + * shown as a constant function. This creates a discontinuous plot + * with $N\times N$ cells on a mesh of size $(N+1)\times(N+1)$. + * + * The default value is false. */ bool discontinuous; + /** + * By default, the MatrixOut class creates a plot in which each matrix + * entry is actually shown, even if it is zero. For large + * matrices, this results in very large output files, or indeed + * exhausts the available memory. On the other hand, if the + * current flags is set to `true`, then the class only outputs + * patches whenever there are non-zero matrix entries to be + * shown. For sparse matrices, this leads to an output size that + * is proportional to the number of nonzero entries, rather than + * proportional to $N^2$. + * + * @note Internally, the current implementation continues to loop + * over all matrix entries, whether they are zero or not. As a + * consequence, the *run time* of outputting sparse matrices + * continues to be proportional to $N^2$, rather than proportional + * to the number of nonzero entries. For large matrices, this means + * that you still have to have patience -- but at least it is possible + * to output information about matrices of size $10,000\times 10,000$ + * or $50,000\times 50,000$ with a few million nonzero entries. + */ + bool create_sparse_plot; + /** * Default constructor. Set all elements of this structure to their * default values. */ Options(const bool show_absolute_values = false, const unsigned int block_size = 1, - const bool discontinuous = false); + const bool discontinuous = false, + const bool create_sparse_plot = false); }; /** @@ -145,7 +178,7 @@ public: void build_patches(const Matrix &matrix, const std::string &name, - const Options options = Options(false, 1, false)); + const Options options = Options(false, 1, false, false)); private: /** @@ -332,24 +365,47 @@ MatrixOut::build_patches(const Matrix &matrix, const std::string &name, const Options options) { - size_type gridpoints_x = (matrix.n() / options.block_size + - (matrix.n() % options.block_size != 0 ? 1 : 0)), - gridpoints_y = (matrix.m() / options.block_size + - (matrix.m() % options.block_size != 0 ? 1 : 0)); + size_type n_patches_x = (matrix.n() / options.block_size + + (matrix.n() % options.block_size != 0 ? 1 : 0)), + n_patches_y = (matrix.m() / options.block_size + + (matrix.m() % options.block_size != 0 ? 1 : 0)); // If continuous, the number of // plotted patches is matrix size-1 if (!options.discontinuous) { - --gridpoints_x; - --gridpoints_y; + --n_patches_x; + --n_patches_y; } + const size_type n_patches = + (options.create_sparse_plot ? + [&]() { + size_type count = 0; + for (size_type i = 0; i < n_patches_y; ++i) + { + for (size_type j = 0; j < n_patches_x; ++j) + // Use the same logic as below to determine whether we + // need to output a patch, and count if we do: + if ((((options.discontinuous == true) && + (get_gridpoint_value(matrix, i, j, options) != 0)) || + ((options.discontinuous == false) && + ((get_gridpoint_value(matrix, i, j, options) != 0) || + (get_gridpoint_value(matrix, i + 1, j, options) != 0) || + (get_gridpoint_value(matrix, i, j + 1, options) != 0) || + (get_gridpoint_value(matrix, i + 1, j + 1, options) != + 0))))) + ++count; + } + return count; + }() : + n_patches_x * n_patches_y); + // first clear old data and re-set the object to a correctly sized state: patches.clear(); try { - patches.resize(gridpoints_x * gridpoints_y); + patches.resize(n_patches); } catch (const std::bad_alloc &) { @@ -357,17 +413,32 @@ MatrixOut::build_patches(const Matrix &matrix, ExcMessage("You are trying to create a graphical " "representation of a matrix that would " "requiring outputting " + - std::to_string(gridpoints_x) + "x" + - std::to_string(gridpoints_y) + + (options.create_sparse_plot ? + std::to_string(n_patches) : + std::to_string(n_patches_x) + "x" + + std::to_string(n_patches_y)) + " patches. There is not enough memory to output " + "this many patches.")); } // now build the patches size_type index = 0; - for (size_type i = 0; i < gridpoints_y; ++i) - for (size_type j = 0; j < gridpoints_x; ++j, ++index) + for (size_type i = 0; i < n_patches_y; ++i) + for (size_type j = 0; j < n_patches_x; ++j) { + // If we are creating a sparse plot, check whether this patch + // would have any nonzero values. If not, we can skip the + // patch: + if (options.create_sparse_plot && + (((options.discontinuous == true) && + (get_gridpoint_value(matrix, i, j, options) == 0)) || + ((options.discontinuous == false) && + (get_gridpoint_value(matrix, i, j, options) == 0) && + (get_gridpoint_value(matrix, i + 1, j, options) == 0) && + (get_gridpoint_value(matrix, i, j + 1, options) == 0) && + (get_gridpoint_value(matrix, i + 1, j + 1, options) == 0)))) + continue; + patches[index].n_subdivisions = 1; patches[index].reference_cell = ReferenceCells::Quadrilateral; @@ -418,7 +489,9 @@ MatrixOut::build_patches(const Matrix &matrix, patches[index].data(0, 3) = get_gridpoint_value(matrix, i + 1, j + 1, options); } - }; + + ++index; + } // finally set the name this->name = name; diff --git a/source/lac/matrix_out.cc b/source/lac/matrix_out.cc index 618526b1ad..c4d55b24c1 100644 --- a/source/lac/matrix_out.cc +++ b/source/lac/matrix_out.cc @@ -23,10 +23,12 @@ DEAL_II_NAMESPACE_OPEN MatrixOut::Options::Options(const bool show_absolute_values, const unsigned int block_size, - const bool discontinuous) + const bool discontinuous, + const bool create_sparse_plot) : show_absolute_values(show_absolute_values) , block_size(block_size) , discontinuous(discontinuous) + , create_sparse_plot(create_sparse_plot) {}