]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make it possible to create sparse outputs from MatrixOut.
authorWolfgang Bangerth <bangerth@colostate.edu>
Fri, 18 Apr 2025 20:00:20 +0000 (14:00 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Fri, 18 Apr 2025 20:00:20 +0000 (14:00 -0600)
include/deal.II/lac/matrix_out.h
source/lac/matrix_out.cc

index ee9b1d5f05fc090d1b6c9287612bf6321aff3eac..f82ac2f1b0b8a906471046cf42113864c61430b6 100644 (file)
@@ -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;
index 618526b1ad5e0377de9436cd23ba638801c49858..c4d55b24c18bca38f78adb333af259c1e517198f 100644 (file)
@@ -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)
 {}
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.