Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local Includes
21 : #include "libmesh/sparse_matrix.h"
22 :
23 : // libMesh includes
24 : #include "libmesh/dof_map.h"
25 : #include "libmesh/dense_matrix.h"
26 : #include "libmesh/diagonal_matrix.h"
27 : #include "libmesh/laspack_matrix.h"
28 : #include "libmesh/libmesh_logging.h"
29 : #include "libmesh/eigen_sparse_matrix.h"
30 : #include "libmesh/parallel.h"
31 : #include "libmesh/petsc_matrix.h"
32 : #include "libmesh/trilinos_epetra_matrix.h"
33 : #include "libmesh/numeric_vector.h"
34 : #include "libmesh/enum_solver_package.h"
35 :
36 :
37 : // gzstream for reading compressed files as a stream
38 : #ifdef LIBMESH_HAVE_GZSTREAM
39 : # include "gzstream.h"
40 : #endif
41 :
42 : #ifdef LIBMESH_HAVE_HDF5
43 : // Sticking to the C API to be safest on portability
44 : # include "hdf5.h"
45 : #endif
46 :
47 : // C++ includes
48 : #include <memory>
49 : #include <fstream>
50 :
51 : #ifdef LIBMESH_HAVE_CXX11_REGEX
52 : #include <regex>
53 : #endif
54 :
55 :
56 : namespace libMesh
57 : {
58 :
59 :
60 : //------------------------------------------------------------------
61 : // SparseMatrix Methods
62 :
63 :
64 : // Constructor
65 : template <typename T>
66 108294 : SparseMatrix<T>::SparseMatrix (const Parallel::Communicator & comm_in) :
67 : ParallelObject(comm_in),
68 102370 : _dof_map(nullptr),
69 102370 : _sp(nullptr),
70 102370 : _is_initialized(false),
71 108294 : _use_hash_table(false)
72 108294 : {}
73 :
74 :
75 :
76 : template <typename T>
77 48463 : void SparseMatrix<T>::attach_dof_map (const DofMap & dof_map)
78 : {
79 48463 : _dof_map = &dof_map;
80 48463 : if (!_sp)
81 23132 : _sp = dof_map.get_sparsity_pattern();
82 48463 : }
83 :
84 :
85 :
86 : template <typename T>
87 41565 : void SparseMatrix<T>::attach_sparsity_pattern (const SparsityPattern::Build & sp)
88 : {
89 41565 : _sp = &sp;
90 41565 : }
91 :
92 :
93 :
94 : // default implementation is to fall back to non-blocked method
95 : template <typename T>
96 0 : void SparseMatrix<T>::add_block_matrix (const DenseMatrix<T> & dm,
97 : const std::vector<numeric_index_type> & brows,
98 : const std::vector<numeric_index_type> & bcols)
99 : {
100 0 : libmesh_assert_equal_to (dm.m() / brows.size(), dm.n() / bcols.size());
101 :
102 : const numeric_index_type blocksize = cast_int<numeric_index_type>
103 0 : (dm.m() / brows.size());
104 :
105 0 : libmesh_assert_equal_to (dm.m()%blocksize, 0);
106 0 : libmesh_assert_equal_to (dm.n()%blocksize, 0);
107 :
108 0 : std::vector<numeric_index_type> rows, cols;
109 :
110 0 : rows.reserve(blocksize*brows.size());
111 0 : cols.reserve(blocksize*bcols.size());
112 :
113 0 : for (auto & row : brows)
114 : {
115 0 : numeric_index_type i = row * blocksize;
116 :
117 0 : for (unsigned int v=0; v<blocksize; v++)
118 0 : rows.push_back(i++);
119 : }
120 :
121 0 : for (auto & col : bcols)
122 : {
123 0 : numeric_index_type j = col * blocksize;
124 :
125 0 : for (unsigned int v=0; v<blocksize; v++)
126 0 : cols.push_back(j++);
127 : }
128 :
129 0 : this->add_matrix (dm, rows, cols);
130 0 : }
131 :
132 :
133 :
134 : // Full specialization of print method for Complex datatypes
135 : template <>
136 0 : void SparseMatrix<Complex>::print(std::ostream & os, const bool sparse) const
137 : {
138 : // std::complex<>::operator<<() is defined, but use this form
139 :
140 0 : if (sparse)
141 : {
142 0 : libmesh_not_implemented();
143 : }
144 :
145 0 : os << "Real part:" << std::endl;
146 0 : for (auto i : make_range(this->m()))
147 : {
148 0 : for (auto j : make_range(this->n()))
149 0 : os << std::setw(8) << (*this)(i,j).real() << " ";
150 0 : os << std::endl;
151 : }
152 :
153 0 : os << std::endl << "Imaginary part:" << std::endl;
154 0 : for (auto i : make_range(this->m()))
155 : {
156 0 : for (auto j : make_range(this->n()))
157 0 : os << std::setw(8) << (*this)(i,j).imag() << " ";
158 0 : os << std::endl;
159 : }
160 0 : }
161 :
162 :
163 :
164 :
165 :
166 :
167 : // Full specialization for Real datatypes
168 : template <typename T>
169 : std::unique_ptr<SparseMatrix<T>>
170 24589 : SparseMatrix<T>::build(const Parallel::Communicator & comm,
171 : const SolverPackage solver_package,
172 : const MatrixBuildType matrix_build_type /* = AUTOMATIC */)
173 : {
174 : // Avoid unused parameter warnings when no solver packages are enabled.
175 698 : libmesh_ignore(comm);
176 :
177 24589 : if (matrix_build_type == MatrixBuildType::DIAGONAL)
178 0 : return std::make_unique<DiagonalMatrix<T>>(comm);
179 :
180 : // Build the appropriate vector
181 24589 : switch (solver_package)
182 : {
183 :
184 : #ifdef LIBMESH_HAVE_LASPACK
185 0 : case LASPACK_SOLVERS:
186 0 : return std::make_unique<LaspackMatrix<T>>(comm);
187 : #endif
188 :
189 :
190 : #ifdef LIBMESH_HAVE_PETSC
191 24335 : case PETSC_SOLVERS:
192 24335 : return std::make_unique<PetscMatrix<T>>(comm);
193 : #endif
194 :
195 :
196 : #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
197 : case TRILINOS_SOLVERS:
198 : return std::make_unique<EpetraMatrix<T>>(comm);
199 : #endif
200 :
201 :
202 : #ifdef LIBMESH_HAVE_EIGEN
203 254 : case EIGEN_SOLVERS:
204 254 : return std::make_unique<EigenSparseMatrix<T>>(comm);
205 : #endif
206 :
207 0 : default:
208 0 : libmesh_error_msg("ERROR: Unrecognized solver package: " << solver_package);
209 : }
210 : }
211 :
212 :
213 : template <typename T>
214 2958121 : void SparseMatrix<T>::vector_mult (NumericVector<T> & dest,
215 : const NumericVector<T> & arg) const
216 : {
217 2958121 : dest.zero();
218 167712 : this->vector_mult_add(dest,arg);
219 2958121 : }
220 :
221 :
222 :
223 : template <typename T>
224 168700 : void SparseMatrix<T>::vector_mult_add (NumericVector<T> & dest,
225 : const NumericVector<T> & arg) const
226 : {
227 : /* This functionality is actually implemented in the \p
228 : NumericVector class. */
229 2959109 : dest.add_vector(arg,*this);
230 168700 : }
231 :
232 :
233 :
234 : template <typename T>
235 0 : void SparseMatrix<T>::zero_rows (std::vector<numeric_index_type> &, T)
236 : {
237 : /* This functionality isn't implemented or stubbed in every subclass yet */
238 0 : libmesh_not_implemented();
239 : }
240 :
241 :
242 : template <typename T>
243 14 : std::size_t SparseMatrix<T>::n_nonzeros() const
244 : {
245 14 : if (!_sp)
246 2 : return 0;
247 0 : return _sp->n_nonzeros();
248 : }
249 :
250 :
251 : template <typename T>
252 0 : void SparseMatrix<T>::print(std::ostream & os, const bool sparse) const
253 : {
254 0 : parallel_object_only();
255 :
256 0 : libmesh_assert (this->initialized());
257 :
258 0 : const numeric_index_type first_dof = this->row_start(),
259 0 : end_dof = this->row_stop();
260 :
261 : // We'll print the matrix from processor 0 to make sure
262 : // it's serialized properly
263 0 : if (this->processor_id() == 0)
264 : {
265 0 : libmesh_assert_equal_to (first_dof, 0);
266 0 : for (numeric_index_type i : make_range(end_dof))
267 : {
268 0 : if (sparse)
269 : {
270 0 : for (auto j : make_range(this->n()))
271 : {
272 0 : T c = (*this)(i,j);
273 0 : if (c != static_cast<T>(0.0))
274 : {
275 0 : os << i << " " << j << " " << c << std::endl;
276 : }
277 : }
278 : }
279 : else
280 : {
281 0 : for (auto j : make_range(this->n()))
282 0 : os << (*this)(i,j) << " ";
283 0 : os << std::endl;
284 : }
285 : }
286 :
287 0 : std::vector<numeric_index_type> ibuf, jbuf;
288 0 : std::vector<T> cbuf;
289 0 : numeric_index_type currenti = end_dof;
290 0 : for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
291 : {
292 0 : this->comm().receive(p, ibuf);
293 0 : this->comm().receive(p, jbuf);
294 0 : this->comm().receive(p, cbuf);
295 0 : libmesh_assert_equal_to (ibuf.size(), jbuf.size());
296 0 : libmesh_assert_equal_to (ibuf.size(), cbuf.size());
297 :
298 0 : if (ibuf.empty())
299 0 : continue;
300 0 : libmesh_assert_greater_equal (ibuf.front(), currenti);
301 0 : libmesh_assert_greater_equal (ibuf.back(), ibuf.front());
302 :
303 0 : std::size_t currentb = 0;
304 0 : for (;currenti <= ibuf.back(); ++currenti)
305 : {
306 0 : if (sparse)
307 : {
308 0 : for (numeric_index_type j=0; j<this->n(); j++)
309 : {
310 0 : if (currentb < ibuf.size() &&
311 0 : ibuf[currentb] == currenti &&
312 0 : jbuf[currentb] == j)
313 : {
314 0 : os << currenti << " " << j << " " << cbuf[currentb] << std::endl;
315 0 : currentb++;
316 : }
317 : }
318 : }
319 : else
320 : {
321 0 : for (auto j : make_range(this->n()))
322 : {
323 0 : if (currentb < ibuf.size() &&
324 0 : ibuf[currentb] == currenti &&
325 0 : jbuf[currentb] == j)
326 : {
327 0 : os << cbuf[currentb] << " ";
328 0 : currentb++;
329 : }
330 : else
331 0 : os << static_cast<T>(0.0) << " ";
332 : }
333 0 : os << std::endl;
334 : }
335 : }
336 : }
337 0 : if (!sparse)
338 : {
339 0 : for (; currenti != this->m(); ++currenti)
340 : {
341 0 : for (numeric_index_type j=0; j<this->n(); j++)
342 0 : os << static_cast<T>(0.0) << " ";
343 0 : os << std::endl;
344 : }
345 : }
346 : }
347 : else
348 : {
349 0 : std::vector<numeric_index_type> ibuf, jbuf;
350 0 : std::vector<T> cbuf;
351 :
352 : // We'll assume each processor has access to entire
353 : // matrix rows, so (*this)(i,j) is valid if i is a local index.
354 0 : for (numeric_index_type i : make_range(first_dof, end_dof))
355 : {
356 0 : for (auto j : make_range(this->n()))
357 : {
358 0 : T c = (*this)(i,j);
359 0 : if (c != static_cast<T>(0.0))
360 : {
361 0 : ibuf.push_back(i);
362 0 : jbuf.push_back(j);
363 0 : cbuf.push_back(c);
364 : }
365 : }
366 : }
367 0 : this->comm().send(0,ibuf);
368 0 : this->comm().send(0,jbuf);
369 0 : this->comm().send(0,cbuf);
370 : }
371 0 : }
372 :
373 :
374 : template <typename T>
375 14 : void SparseMatrix<T>::print_matlab(const std::string & name) const
376 : {
377 2 : parallel_object_only();
378 :
379 2 : libmesh_assert (this->initialized());
380 :
381 14 : const numeric_index_type first_dof = this->row_start(),
382 14 : end_dof = this->row_stop();
383 :
384 : // We'll print the matrix from processor 0 to make sure
385 : // it's serialized properly
386 16 : if (this->processor_id() == 0)
387 : {
388 12 : std::unique_ptr<std::ofstream> file;
389 :
390 14 : if (name != "")
391 24 : file = std::make_unique<std::ofstream>(name.c_str());
392 :
393 14 : std::ostream & os = (name == "") ? libMesh::out : *file;
394 :
395 14 : std::size_t sparsity_nonzeros = this->n_nonzeros();
396 :
397 2 : std::size_t real_nonzeros = 0;
398 :
399 2 : libmesh_assert_equal_to(first_dof, 0);
400 70 : for (numeric_index_type i : make_range(end_dof))
401 : {
402 320 : for (auto j : make_range(this->n()))
403 : {
404 224 : T c = (*this)(i,j);
405 208 : if (c != static_cast<T>(0.0))
406 224 : ++real_nonzeros;
407 : }
408 : }
409 :
410 :
411 14 : for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
412 : {
413 0 : std::size_t nonzeros_on_p = 0;
414 0 : this->comm().receive(p, nonzeros_on_p);
415 0 : real_nonzeros += nonzeros_on_p;
416 : }
417 :
418 : if (sparsity_nonzeros &&
419 : sparsity_nonzeros != real_nonzeros)
420 : libmesh_warning(sparsity_nonzeros <<
421 : " nonzeros allocated, but " <<
422 : real_nonzeros << " used.");
423 :
424 : // We probably want to be more consistent than that, if our
425 : // sparsity is overallocated.
426 :
427 : // Print a header similar to PETSc's mat_view ascii_matlab
428 14 : os << "%Mat Object: () " << this->n_processors() << " MPI processes\n"
429 6 : << "% type: " << (this->n_processors() > 1 ? "mpi" : "seq") << "aij\n"
430 40 : << "% Size = " << this->m() << ' ' << this->n() << '\n'
431 14 : << "% Nonzeros = " << real_nonzeros << '\n'
432 14 : << "zzz = zeros(" << real_nonzeros << ",3);\n"
433 14 : << "zzz = [\n";
434 :
435 70 : for (numeric_index_type i : make_range(end_dof))
436 : {
437 : // FIXME - we need a base class way to iterate over a
438 : // SparseMatrix row.
439 320 : for (auto j : make_range(this->n()))
440 : {
441 224 : T c = (*this)(i,j);
442 208 : if (c != static_cast<T>(0.0))
443 : {
444 : // Convert from 0-based to 1-based indexing
445 432 : os << (i+1) << ' ' << (j+1) << " " << c << '\n';
446 : }
447 : }
448 : }
449 :
450 4 : std::vector<numeric_index_type> ibuf, jbuf;
451 4 : std::vector<T> cbuf;
452 14 : for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
453 : {
454 0 : this->comm().receive(p, ibuf);
455 0 : this->comm().receive(p, jbuf);
456 0 : this->comm().receive(p, cbuf);
457 0 : libmesh_assert_equal_to (ibuf.size(), jbuf.size());
458 0 : libmesh_assert_equal_to (ibuf.size(), cbuf.size());
459 :
460 0 : for (auto n : index_range(ibuf))
461 0 : os << ibuf[n] << ' ' << jbuf[n] << " " << cbuf[n] << '\n';
462 : }
463 :
464 12 : os << "];\n" << "Mat_sparse = spconvert(zzz);" << std::endl;
465 10 : }
466 : else
467 : {
468 0 : std::vector<numeric_index_type> ibuf, jbuf;
469 0 : std::vector<T> cbuf;
470 0 : std::size_t my_nonzeros = 0;
471 :
472 : // We'll assume each processor has access to entire
473 : // matrix rows, so (*this)(i,j) is valid if i is a local index.
474 0 : for (numeric_index_type i : make_range(first_dof, end_dof))
475 : {
476 0 : for (auto j : make_range(this->n()))
477 : {
478 0 : T c = (*this)(i,j);
479 0 : if (c != static_cast<T>(0.0))
480 : {
481 0 : ibuf.push_back(i);
482 0 : jbuf.push_back(j);
483 0 : cbuf.push_back(c);
484 0 : ++my_nonzeros;
485 : }
486 : }
487 : }
488 0 : this->comm().send(0,my_nonzeros);
489 0 : this->comm().send(0,ibuf);
490 0 : this->comm().send(0,jbuf);
491 0 : this->comm().send(0,cbuf);
492 : }
493 14 : }
494 :
495 :
496 :
497 : template <typename T>
498 0 : void SparseMatrix<T>::print_petsc_binary(const std::string &)
499 : {
500 0 : libmesh_not_implemented_msg
501 : ("libMesh cannot write PETSc binary-format files from non-PETSc matrices");
502 : }
503 :
504 :
505 :
506 : template <typename T>
507 0 : void SparseMatrix<T>::print_petsc_hdf5(const std::string &)
508 : {
509 0 : libmesh_not_implemented_msg
510 : ("libMesh cannot write PETSc HDF5-format files from non-PETSc matrices");
511 : }
512 :
513 :
514 :
515 : template <typename T>
516 1068 : void SparseMatrix<T>::read(const std::string & filename)
517 : {
518 : {
519 1132 : std::ifstream in (filename.c_str());
520 1068 : libmesh_error_msg_if
521 : (!in.good(), "ERROR: cannot read file:\n\t" <<
522 : filename);
523 1004 : }
524 :
525 1068 : std::string_view basename = Utility::basename_of(filename);
526 :
527 1068 : const bool gzipped_file = Utility::ends_with(filename, ".gz");
528 :
529 1068 : if (gzipped_file)
530 12 : basename.remove_suffix(3);
531 :
532 2136 : if (Utility::ends_with(basename, ".matlab") ||
533 1172 : Utility::ends_with(basename, ".m"))
534 994 : this->read_matlab(filename);
535 74 : else if (Utility::ends_with(basename, ".petsc64"))
536 : {
537 : #ifndef LIBMESH_HAVE_PETSC
538 0 : libmesh_error_msg("Cannot load PETSc matrix file " <<
539 : filename << " without PETSc-enabled libMesh.");
540 : #endif
541 : #if LIBMESH_DOF_ID_BYTES != 8
542 0 : libmesh_error_msg("Cannot load 64-bit PETSc matrix file " <<
543 : filename << " with non-64-bit libMesh.");
544 : #endif
545 66 : this->read_petsc_binary(filename);
546 : }
547 8 : else if (Utility::ends_with(basename, ".petsc32"))
548 : {
549 : #ifndef LIBMESH_HAVE_PETSC
550 0 : libmesh_error_msg("Cannot load PETSc matrix file " <<
551 : filename << " without PETSc-enabled libMesh.");
552 : #endif
553 : #if LIBMESH_DOF_ID_BYTES != 4
554 0 : libmesh_error_msg("Cannot load 32-bit PETSc matrix file " <<
555 : filename << " with non-32-bit libMesh.");
556 : #endif
557 4 : this->read_petsc_binary(filename);
558 : }
559 4 : else if (Utility::ends_with(basename, ".h5"))
560 : {
561 8 : this->read_coreform_hdf5(filename);
562 : }
563 : else
564 0 : libmesh_error_msg(" ERROR: Unrecognized matrix file extension on: "
565 : << basename
566 : << "\n I understand the following:\n\n"
567 : << " *.h5 -- CoreForm HDF5 sparse matrix format\n"
568 : << " *.matlab -- Matlab sparse matrix format\n"
569 : << " *.m -- Matlab sparse matrix format\n"
570 : << " *.petsc32 -- PETSc binary format, 32-bit\n"
571 : << " *.petsc64 -- PETSc binary format, 64-bit\n"
572 : );
573 1068 : }
574 :
575 :
576 :
577 : template <typename T>
578 4 : void SparseMatrix<T>::read_coreform_hdf5(const std::string & filename,
579 : const std::string & groupname)
580 : {
581 : #ifndef LIBMESH_HAVE_HDF5
582 0 : libmesh_ignore(filename, groupname);
583 0 : libmesh_error_msg("ERROR: need HDF5 support to handle .h5 files!!!");
584 : #else
585 : LOG_SCOPE("read_coreform_hdf5()", "SparseMatrix");
586 :
587 : std::size_t num_rows, num_cols;
588 :
589 : hid_t group;
590 :
591 18 : auto check_open = [&filename](hid_t id, const std::string & objname)
592 : {
593 18 : if (id == H5I_INVALID_HID)
594 0 : libmesh_error_msg("Couldn't open " + objname + " in " + filename);
595 : };
596 :
597 52 : auto check_hdf5 = [&filename](auto hdf5val, const std::string & objname)
598 : {
599 48 : if (hdf5val < 0)
600 0 : libmesh_error_msg("HDF5 error from " + objname + " in " + filename);
601 : };
602 :
603 :
604 4 : if (this->processor_id() == 0)
605 : {
606 3 : const hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
607 :
608 3 : if (file == H5I_INVALID_HID)
609 0 : libmesh_file_error(filename);
610 :
611 3 : group = H5Gopen(file, groupname.c_str(), H5P_DEFAULT);
612 3 : check_open(group, groupname);
613 :
614 27 : auto read_size_attribute = [&filename, &group, &check_open, &check_hdf5]
615 : (const std::string & attribute_name)
616 : {
617 6 : unsigned long long returnval = 0;
618 :
619 6 : const hid_t attr = H5Aopen(group, attribute_name.c_str(), H5P_DEFAULT);
620 6 : check_open(attr, attribute_name);
621 :
622 6 : const hid_t attr_type = H5Aget_type(attr);
623 6 : check_hdf5(attr_type, attribute_name + " type");
624 :
625 : // HDF5 will convert between the file's integer type and ours, but
626 : // we do expect an integer type.
627 6 : if (H5Tget_class(attr_type) != H5T_INTEGER)
628 0 : libmesh_error_msg("Non-integer type for " + attribute_name + " in " + filename);
629 :
630 6 : H5Tclose(attr_type);
631 :
632 : // HDF5 is supposed to handle both upscaling and endianness
633 : // conversions here
634 6 : herr_t errval = H5Aread(attr, H5T_NATIVE_ULLONG, &returnval);
635 6 : check_hdf5(errval, attribute_name + " read");
636 :
637 6 : H5Aclose(attr);
638 :
639 6 : return returnval;
640 : };
641 :
642 3 : num_cols = read_size_attribute("num_cols");
643 3 : num_rows = read_size_attribute("num_rows");
644 :
645 3 : this->comm().broadcast(num_cols);
646 3 : this->comm().broadcast(num_rows);
647 : }
648 : else
649 : {
650 1 : this->comm().broadcast(num_cols);
651 1 : this->comm().broadcast(num_rows);
652 : }
653 :
654 : numeric_index_type new_row_start, new_row_stop,
655 : new_col_start, new_col_stop;
656 :
657 : // If we need to reinit, we need to determine which rows+columns
658 : // each processor is in charge of.
659 : std::vector<numeric_index_type> new_row_starts, new_row_stops,
660 : new_col_starts, new_col_stops;
661 :
662 4 : if (this->initialized() &&
663 4 : num_cols == this->m() &&
664 0 : num_rows == this->n())
665 : {
666 0 : new_row_start = this->row_start(),
667 0 : new_row_stop = this->row_stop();
668 :
669 0 : new_col_start = this->col_start(),
670 0 : new_col_stop = this->col_stop();
671 : }
672 : else
673 : {
674 : // Determine which rows/columns each processor will be in charge of
675 4 : new_row_start = this->processor_id() * num_rows / this->n_processors(),
676 4 : new_row_stop = (this->processor_id()+1) * num_rows / this->n_processors();
677 :
678 4 : new_col_start = this->processor_id() * num_cols / this->n_processors(),
679 4 : new_col_stop = (this->processor_id()+1) * num_cols / this->n_processors();
680 : }
681 :
682 4 : this->comm().gather(0, new_row_start, new_row_starts);
683 4 : this->comm().gather(0, new_row_stop, new_row_stops);
684 4 : this->comm().gather(0, new_col_start, new_col_starts);
685 4 : this->comm().gather(0, new_col_stop, new_col_stops);
686 :
687 4 : numeric_index_type on_diagonal_nonzeros = 0,
688 4 : off_diagonal_nonzeros = 0;
689 :
690 : std::vector<std::size_t> cols, row_offsets;
691 : std::vector<double> vals;
692 :
693 4 : if (this->processor_id() == 0)
694 : {
695 12 : auto read_vector = [&filename, &group, &check_open, &check_hdf5]
696 : (const std::string & dataname, auto hdf5_class,
697 : auto hdf5_type, auto & datavec)
698 : {
699 9 : const hid_t data = H5Dopen1(group, dataname.c_str());
700 9 : check_open(data, dataname.c_str());
701 :
702 9 : const hid_t data_type = H5Dget_type(data);
703 9 : check_hdf5(data_type, dataname + " type");
704 :
705 : // HDF5 will convert between the file's integer type and ours, but
706 : // we do expect an integer type.
707 9 : if (H5Tget_class(data_type) != hdf5_class)
708 0 : libmesh_error_msg("Unexpected type for " + dataname + " in " + filename);
709 :
710 9 : H5Tclose(data_type);
711 :
712 9 : const hid_t dataspace = H5Dget_space(data);
713 9 : check_hdf5(dataspace, dataname + " space");
714 :
715 9 : int ndims = H5Sget_simple_extent_ndims(dataspace);
716 9 : if (ndims != 1)
717 0 : libmesh_error_msg("Non-vector space for " + dataname + " in " + filename);
718 :
719 : hsize_t len, maxlen;
720 9 : herr_t errval = H5Sget_simple_extent_dims(dataspace, &len, &maxlen);
721 9 : check_hdf5(errval, dataname + " dims");
722 :
723 9 : datavec.resize(len);
724 :
725 9 : errval = H5Dread(data, hdf5_type, H5S_ALL, H5S_ALL, H5P_DEFAULT, datavec.data());
726 18 : check_hdf5(errval, dataname + " read");
727 : };
728 :
729 3 : read_vector("cols", H5T_INTEGER, H5T_NATIVE_ULLONG, cols);
730 3 : read_vector("row_offsets", H5T_INTEGER, H5T_NATIVE_ULLONG, row_offsets);
731 6 : read_vector("vals", H5T_FLOAT, H5T_NATIVE_DOUBLE, vals);
732 :
733 3 : if (cols.size() != vals.size())
734 0 : libmesh_error_msg("Inconsistent cols/vals sizes in " + filename);
735 :
736 3 : if (row_offsets.size() != num_rows + 1)
737 0 : libmesh_error_msg("Inconsistent row_offsets size in " + filename);
738 :
739 : // Data for the row we're working on
740 : numeric_index_type current_row = 0;
741 : processor_id_type current_proc = 0;
742 3 : numeric_index_type current_on_diagonal_nonzeros = 0;
743 3 : numeric_index_type current_off_diagonal_nonzeros = 0;
744 3 : if (row_offsets[0] != 0)
745 0 : libmesh_error_msg("Unexpected row_offsets[0] in " + filename);
746 :
747 381 : for (auto i : index_range(cols))
748 : {
749 456 : while (row_offsets[current_row+1] <= i)
750 : {
751 : ++current_row;
752 78 : if (row_offsets[current_row] < row_offsets[current_row-1])
753 0 : libmesh_error_msg("Non-monotonic row_offsets in " + filename);
754 78 : current_on_diagonal_nonzeros = 0;
755 78 : current_off_diagonal_nonzeros = 0;
756 : }
757 :
758 379 : while (current_row >= (new_row_stops[current_proc]+1))
759 1 : ++current_proc;
760 :
761 : // 0-based indexing in file
762 378 : if (cols[i] >= new_col_starts[current_proc] &&
763 352 : cols[i] < new_col_stops[current_proc])
764 : {
765 316 : ++current_on_diagonal_nonzeros;
766 316 : on_diagonal_nonzeros =
767 : std::max(on_diagonal_nonzeros,
768 : current_on_diagonal_nonzeros);
769 : }
770 : else
771 : {
772 62 : ++current_off_diagonal_nonzeros;
773 62 : off_diagonal_nonzeros =
774 : std::max(off_diagonal_nonzeros,
775 : current_off_diagonal_nonzeros);
776 : }
777 : }
778 : }
779 :
780 4 : this->comm().broadcast(on_diagonal_nonzeros);
781 4 : this->comm().broadcast(off_diagonal_nonzeros);
782 :
783 4 : this->init(num_rows, num_cols,
784 : new_row_stop-new_row_start,
785 : new_col_stop-new_col_start,
786 : on_diagonal_nonzeros,
787 : off_diagonal_nonzeros);
788 :
789 : // Set the matrix values last.
790 4 : if (this->processor_id() == 0)
791 : {
792 : numeric_index_type current_row = 0;
793 381 : for (auto i : index_range(cols))
794 : {
795 456 : while (row_offsets[current_row+1] <= i)
796 : {
797 : ++current_row;
798 : libmesh_assert_greater_equal (row_offsets[current_row],
799 : row_offsets[current_row-1]);
800 : }
801 378 : this->set(current_row, cols[i], vals[i]);
802 : }
803 : }
804 :
805 4 : this->close();
806 : #endif // LIBMESH_HAVE_HDF5
807 4 : }
808 :
809 :
810 :
811 : template <typename T>
812 1415 : void SparseMatrix<T>::read_matlab(const std::string & filename)
813 : {
814 84 : LOG_SCOPE("read_matlab()", "SparseMatrix");
815 :
816 : #ifndef LIBMESH_HAVE_CXX11_REGEX
817 : libmesh_not_implemented(); // What is your compiler?!? Email us!
818 : libmesh_ignore(filename);
819 : #else
820 42 : parallel_object_only();
821 :
822 1415 : const bool gzipped_file = Utility::ends_with(filename, ".gz");
823 :
824 : // The sizes we get from the file
825 1415 : std::size_t m = 0,
826 1415 : n = 0;
827 :
828 : // If we don't already have this size, we'll need to reinit, and
829 : // determine which rows+columns each processor is in charge of.
830 84 : std::vector<numeric_index_type> new_row_starts, new_row_stops,
831 84 : new_col_starts, new_col_stops;
832 :
833 : numeric_index_type new_row_start, new_row_stop,
834 : new_col_start, new_col_stop;
835 :
836 : // We'll read through the file three times: once to get a reliable
837 : // value for the matrix size (so we can divvy it up among
838 : // processors), then again to get the sparsity to send to each
839 : // processor, then a final time to get the entries to send to each
840 : // processor.
841 : //
842 : // We'll use an istream here; it might be an ifstream if we're
843 : // opening a raw ASCII file or a gzstream if we're opening a
844 : // compressed one.
845 1373 : std::unique_ptr<std::istream> file;
846 :
847 : // We'll need a temporary structure to cache matrix entries, because
848 : // we need to read through the whole file before we know the size
849 : // and sparsity structure with which we can init().
850 : //
851 : // Reading through the file three times via `seekg` doesn't work
852 : // with our gzstream wrapper, and seems to take three times as long
853 : // even with a plain ifstream. What happened to disk caching!?
854 84 : std::vector<std::tuple<numeric_index_type, numeric_index_type, T>> entries;
855 :
856 : // First read through the file, saving size and entry data
857 : {
858 : // We'll read the matrix on processor 0 rather than try to juggle
859 : // parallel I/O.
860 1457 : if (this->processor_id() == 0)
861 : {
862 : // We'll be using regular expressions to make ourselves slightly
863 : // more robust to formatting.
864 705 : const std::regex start_regex // assignment like "zzz = ["
865 : ("\\s*\\w+\\s*=\\s*\\[");
866 705 : const std::regex end_regex // end of assignment
867 : ("^[^%]*\\]");
868 :
869 647 : if (gzipped_file)
870 : {
871 : #ifdef LIBMESH_HAVE_GZSTREAM
872 255 : auto inf = std::make_unique<igzstream>();
873 9 : libmesh_assert(inf);
874 9 : inf->open(filename.c_str(), std::ios::in);
875 237 : file = std::move(inf);
876 : #else
877 : libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
878 : #endif
879 228 : }
880 : else
881 : {
882 421 : auto inf = std::make_unique<std::ifstream>();
883 20 : libmesh_assert(inf);
884 :
885 421 : std::string new_name = Utility::unzip_file(filename);
886 :
887 401 : inf->open(new_name.c_str(), std::ios::in);
888 381 : file = std::move(inf);
889 361 : }
890 :
891 : // If we have a matrix with all-zero trailing rows, the only
892 : // way to get the size is if it ended up in a comment
893 705 : const std::regex size_regex // comment like "% size = 8 8"
894 : ("%\\s*[Ss][Ii][Zz][Ee]\\s*=\\s*(\\d+)\\s+(\\d+)");
895 676 : const std::string whitespace = " \t";
896 :
897 29 : bool have_started = false;
898 29 : bool have_ended = false;
899 647 : std::size_t largest_i_seen = 0, largest_j_seen = 0;
900 :
901 : // Data for the row we're working on
902 : // Use 1-based indexing for current_row, as in the file
903 647 : std::size_t current_row = 1;
904 :
905 334671 : for (std::string line; std::getline(*file, line);)
906 : {
907 7616 : std::smatch sm;
908 :
909 : // First, try to match an entry. This is the most common
910 : // case so we won't rely on slow std::regex for it.
911 : // stringstream is at least an improvement over that.
912 :
913 : // Look for row/col/val like "1 1 -2.0e-4"
914 :
915 182524 : std::istringstream l(line);
916 :
917 : std::size_t i, j;
918 1822 : T value;
919 :
920 9438 : l >> i >> j >> value;
921 :
922 174937 : if (!l.fail())
923 : {
924 170408 : libmesh_error_msg_if
925 : (!have_started, "Confused by premature entries in matrix file " << filename);
926 :
927 324217 : entries.emplace_back(cast_int<numeric_index_type>(i),
928 177821 : cast_int<numeric_index_type>(j),
929 : value);
930 :
931 170408 : libmesh_error_msg_if
932 : (!i || !j, "Expected 1-based indexing in matrix file "
933 : << filename);
934 :
935 170408 : current_row = std::max(current_row, i);
936 :
937 170408 : libmesh_error_msg_if
938 : (i < current_row,
939 : "Can't handle out-of-order entries in matrix file "
940 : << filename);
941 :
942 170408 : largest_i_seen = std::max(i, largest_i_seen);
943 170408 : largest_j_seen = std::max(j, largest_j_seen);
944 : }
945 :
946 4529 : else if (std::regex_search(line, sm, size_regex))
947 : {
948 676 : const std::string msize = sm[1];
949 647 : const std::string nsize = sm[2];
950 647 : m = std::stoull(msize);
951 647 : n = std::stoull(nsize);
952 : }
953 :
954 3882 : else if (std::regex_search(line, start_regex))
955 29 : have_started = true;
956 :
957 3235 : else if (std::regex_search(line, end_regex))
958 : {
959 29 : have_ended = true;
960 58 : break;
961 : }
962 : }
963 :
964 647 : libmesh_error_msg_if
965 : (!have_started, "Confused by missing assignment beginning in matrix file " << filename);
966 :
967 647 : libmesh_error_msg_if
968 : (!have_ended, "Confused by missing assignment ending in matrix file " << filename);
969 :
970 647 : libmesh_error_msg_if
971 : (m > largest_i_seen, "Confused by missing final row(s) in matrix file " << filename);
972 :
973 647 : libmesh_error_msg_if
974 : (m > 0 && m < largest_i_seen, "Confused by extra final row(s) in matrix file " << filename);
975 :
976 647 : if (!m)
977 0 : m = largest_i_seen;
978 :
979 647 : libmesh_error_msg_if
980 : (n > largest_j_seen, "Confused by missing final column(s) in matrix file " << filename);
981 :
982 647 : libmesh_error_msg_if
983 : (n > 0 && n < largest_j_seen, "Confused by extra final column(s) in matrix file " << filename);
984 :
985 647 : if (!n)
986 0 : n = largest_j_seen;
987 :
988 647 : this->comm().broadcast(m);
989 647 : this->comm().broadcast(n);
990 589 : }
991 : else
992 : {
993 755 : this->comm().broadcast(m);
994 768 : this->comm().broadcast(n);
995 : }
996 :
997 1415 : if (this->initialized() &&
998 1417 : m == this->m() &&
999 70 : n == this->n())
1000 : {
1001 70 : new_row_start = this->row_start(),
1002 70 : new_row_stop = this->row_stop();
1003 :
1004 72 : new_col_start = this->col_start(),
1005 70 : new_col_stop = this->col_stop();
1006 : }
1007 : else
1008 : {
1009 : // Determine which rows/columns each processor will be in charge of
1010 1385 : new_row_start = this->processor_id() * m / this->n_processors(),
1011 1345 : new_row_stop = (this->processor_id()+1) * m / this->n_processors();
1012 :
1013 1385 : new_col_start = this->processor_id() * n / this->n_processors(),
1014 1345 : new_col_stop = (this->processor_id()+1) * n / this->n_processors();
1015 : }
1016 :
1017 1415 : this->comm().gather(0, new_row_start, new_row_starts);
1018 1415 : this->comm().gather(0, new_row_stop, new_row_stops);
1019 1415 : this->comm().gather(0, new_col_start, new_col_starts);
1020 1415 : this->comm().gather(0, new_col_stop, new_col_stops);
1021 :
1022 : } // Done reading entry data and broadcasting matrix size
1023 :
1024 : // Calculate the matrix sparsity and initialize it second
1025 : {
1026 : // Deduce the sparsity pattern, or at least the maximum number of
1027 : // on- and off- diagonal non-zeros per row.
1028 1415 : numeric_index_type on_diagonal_nonzeros =0,
1029 1415 : off_diagonal_nonzeros =0;
1030 :
1031 1457 : if (this->processor_id() == 0)
1032 : {
1033 : // Data for the row we're working on
1034 : // Use 1-based indexing for current_row, as in the file
1035 29 : numeric_index_type current_row = 1;
1036 29 : processor_id_type current_proc = 0;
1037 647 : numeric_index_type current_on_diagonal_nonzeros = 0;
1038 647 : numeric_index_type current_off_diagonal_nonzeros = 0;
1039 :
1040 171055 : for (auto [i, j, value] : entries)
1041 : {
1042 170408 : if (i > current_row)
1043 : {
1044 876 : current_row = i;
1045 : // +1 for 1-based indexing in file
1046 22331 : while (current_row >= (new_row_stops[current_proc]+1))
1047 768 : ++current_proc;
1048 20674 : current_on_diagonal_nonzeros = 0;
1049 20674 : current_off_diagonal_nonzeros = 0;
1050 : }
1051 :
1052 : // +1 for 1-based indexing in file
1053 184358 : if (j >= (new_col_starts[current_proc]+1) &&
1054 165360 : j < (new_col_stops[current_proc]+1))
1055 : {
1056 143791 : ++current_on_diagonal_nonzeros;
1057 143791 : on_diagonal_nonzeros =
1058 : std::max(on_diagonal_nonzeros,
1059 5555 : current_on_diagonal_nonzeros);
1060 : }
1061 : else
1062 : {
1063 26617 : ++current_off_diagonal_nonzeros;
1064 26617 : off_diagonal_nonzeros =
1065 : std::max(off_diagonal_nonzeros,
1066 1858 : current_off_diagonal_nonzeros);
1067 : }
1068 : }
1069 : }
1070 :
1071 1415 : this->comm().broadcast(on_diagonal_nonzeros);
1072 1415 : this->comm().broadcast(off_diagonal_nonzeros);
1073 :
1074 1415 : this->init(m, n,
1075 : new_row_stop-new_row_start,
1076 : new_col_stop-new_col_start,
1077 : on_diagonal_nonzeros,
1078 : off_diagonal_nonzeros);
1079 : }
1080 :
1081 : // Set the matrix values last.
1082 : // Convert from 1-based to 0-based indexing
1083 1457 : if (this->processor_id() == 0)
1084 171055 : for (auto [i, j, value] : entries)
1085 170408 : this->set(i-1, j-1, value);
1086 :
1087 1415 : this->close();
1088 : #endif
1089 2746 : }
1090 :
1091 :
1092 :
1093 : template <typename T>
1094 0 : void SparseMatrix<T>::read_petsc_binary(const std::string &)
1095 : {
1096 0 : libmesh_not_implemented_msg
1097 : ("libMesh cannot read PETSc binary-format files into non-PETSc matrices");
1098 : }
1099 :
1100 :
1101 :
1102 : template <typename T>
1103 0 : void SparseMatrix<T>::read_petsc_hdf5(const std::string &)
1104 : {
1105 0 : libmesh_not_implemented_msg
1106 : ("libMesh cannot read PETSc HDF5-format files into non-PETSc matrices");
1107 : }
1108 :
1109 :
1110 :
1111 : template <typename T>
1112 75 : void SparseMatrix<T>::scale(const T scale)
1113 : {
1114 4 : libmesh_assert(this->closed());
1115 :
1116 375 : for (const auto i : make_range(this->row_start(), this->row_stop()))
1117 1500 : for (const auto j : make_range(this->col_start(), this->col_stop()))
1118 1200 : this->set(i, j, (*this)(i, j) * scale);
1119 75 : }
1120 :
1121 :
1122 :
1123 : template <typename T>
1124 290 : Real SparseMatrix<T>::l1_norm_diff(const SparseMatrix<T> & other_mat) const
1125 : {
1126 290 : auto diff_mat = this->clone();
1127 290 : diff_mat->add(-1.0, other_mat);
1128 580 : return diff_mat->l1_norm();
1129 266 : }
1130 :
1131 : //------------------------------------------------------------------
1132 : // Explicit instantiations
1133 : template class LIBMESH_EXPORT SparseMatrix<Number>;
1134 :
1135 : } // namespace libMesh
|
