Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 : // Anonymous namespace for helper functions
56 : namespace
57 : {
58 : #ifdef LIBMESH_HAVE_HDF5
59 38 : void check_open(const std::string & filename,
60 : hid_t id,
61 : const std::string & objname)
62 : {
63 38 : if (id == H5I_INVALID_HID)
64 0 : libmesh_error_msg("Couldn't open " + objname + " in " + filename);
65 :
66 38 : }
67 :
68 : template <typename T>
69 126 : void check_hdf5(const std::string & filename, T hdf5val, const std::string & objname)
70 : {
71 126 : if (hdf5val < 0)
72 0 : libmesh_error_msg("HDF5 error from " + objname + " in " + filename);
73 126 : }
74 : #endif
75 : }
76 :
77 :
78 :
79 : namespace libMesh
80 : {
81 :
82 :
83 : //------------------------------------------------------------------
84 : // SparseMatrix Methods
85 :
86 :
87 : // Constructor
88 : template <typename T>
89 114547 : SparseMatrix<T>::SparseMatrix (const Parallel::Communicator & comm_in) :
90 : ParallelObject(comm_in),
91 108351 : _dof_map(nullptr),
92 108351 : _sp(nullptr),
93 108351 : _is_initialized(false),
94 114547 : _use_hash_table(false)
95 114547 : {}
96 :
97 :
98 :
99 : template <typename T>
100 53339 : void SparseMatrix<T>::attach_dof_map (const DofMap & dof_map)
101 : {
102 53339 : _dof_map = &dof_map;
103 53339 : if (!_sp)
104 24508 : _sp = dof_map.get_sparsity_pattern();
105 53339 : }
106 :
107 :
108 :
109 : template <typename T>
110 42783 : void SparseMatrix<T>::attach_sparsity_pattern (const SparsityPattern::Build & sp)
111 : {
112 42783 : _sp = &sp;
113 42783 : }
114 :
115 :
116 :
117 : // default implementation is to fall back to non-blocked method
118 : template <typename T>
119 0 : void SparseMatrix<T>::add_block_matrix (const DenseMatrix<T> & dm,
120 : const std::vector<numeric_index_type> & brows,
121 : const std::vector<numeric_index_type> & bcols)
122 : {
123 0 : libmesh_assert_equal_to (dm.m() / brows.size(), dm.n() / bcols.size());
124 :
125 : const numeric_index_type blocksize = cast_int<numeric_index_type>
126 0 : (dm.m() / brows.size());
127 :
128 0 : libmesh_assert_equal_to (dm.m()%blocksize, 0);
129 0 : libmesh_assert_equal_to (dm.n()%blocksize, 0);
130 :
131 0 : std::vector<numeric_index_type> rows, cols;
132 :
133 0 : rows.reserve(blocksize*brows.size());
134 0 : cols.reserve(blocksize*bcols.size());
135 :
136 0 : for (auto & row : brows)
137 : {
138 0 : numeric_index_type i = row * blocksize;
139 :
140 0 : for (unsigned int v=0; v<blocksize; v++)
141 0 : rows.push_back(i++);
142 : }
143 :
144 0 : for (auto & col : bcols)
145 : {
146 0 : numeric_index_type j = col * blocksize;
147 :
148 0 : for (unsigned int v=0; v<blocksize; v++)
149 0 : cols.push_back(j++);
150 : }
151 :
152 0 : this->add_matrix (dm, rows, cols);
153 0 : }
154 :
155 :
156 :
157 : // Full specialization of print method for Complex datatypes
158 : template <>
159 0 : void SparseMatrix<Complex>::print(std::ostream & os, const bool sparse) const
160 : {
161 : // std::complex<>::operator<<() is defined, but use this form
162 :
163 0 : if (sparse)
164 : {
165 0 : libmesh_not_implemented();
166 : }
167 :
168 0 : os << "Real part:" << std::endl;
169 0 : for (auto i : make_range(this->m()))
170 : {
171 0 : for (auto j : make_range(this->n()))
172 0 : os << std::setw(8) << (*this)(i,j).real() << " ";
173 0 : os << std::endl;
174 : }
175 :
176 0 : os << std::endl << "Imaginary part:" << std::endl;
177 0 : for (auto i : make_range(this->m()))
178 : {
179 0 : for (auto j : make_range(this->n()))
180 0 : os << std::setw(8) << (*this)(i,j).imag() << " ";
181 0 : os << std::endl;
182 : }
183 0 : }
184 :
185 :
186 :
187 :
188 :
189 :
190 : // Full specialization for Real datatypes
191 : template <typename T>
192 : std::unique_ptr<SparseMatrix<T>>
193 29447 : SparseMatrix<T>::build(const Parallel::Communicator & comm,
194 : const SolverPackage solver_package,
195 : const MatrixBuildType matrix_build_type /* = AUTOMATIC */)
196 : {
197 : // Avoid unused parameter warnings when no solver packages are enabled.
198 814 : libmesh_ignore(comm);
199 :
200 29447 : if (matrix_build_type == MatrixBuildType::DIAGONAL)
201 0 : return std::make_unique<DiagonalMatrix<T>>(comm);
202 :
203 : // Build the appropriate vector
204 29447 : switch (solver_package)
205 : {
206 :
207 : #ifdef LIBMESH_HAVE_LASPACK
208 0 : case LASPACK_SOLVERS:
209 0 : return std::make_unique<LaspackMatrix<T>>(comm);
210 : #endif
211 :
212 :
213 : #ifdef LIBMESH_HAVE_PETSC
214 29187 : case PETSC_SOLVERS:
215 29187 : return std::make_unique<PetscMatrix<T>>(comm);
216 : #endif
217 :
218 :
219 : #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
220 : case TRILINOS_SOLVERS:
221 : return std::make_unique<EpetraMatrix<T>>(comm);
222 : #endif
223 :
224 :
225 : #ifdef LIBMESH_HAVE_EIGEN
226 260 : case EIGEN_SOLVERS:
227 260 : return std::make_unique<EigenSparseMatrix<T>>(comm);
228 : #endif
229 :
230 0 : default:
231 0 : libmesh_error_msg("ERROR: Unrecognized solver package: " << solver_package);
232 : }
233 : }
234 :
235 :
236 : template <typename T>
237 2958121 : void SparseMatrix<T>::vector_mult (NumericVector<T> & dest,
238 : const NumericVector<T> & arg) const
239 : {
240 2958121 : dest.zero();
241 167712 : this->vector_mult_add(dest,arg);
242 2958121 : }
243 :
244 :
245 :
246 : template <typename T>
247 168700 : void SparseMatrix<T>::vector_mult_add (NumericVector<T> & dest,
248 : const NumericVector<T> & arg) const
249 : {
250 : /* This functionality is actually implemented in the \p
251 : NumericVector class. */
252 2959109 : dest.add_vector(arg,*this);
253 168700 : }
254 :
255 :
256 :
257 : template <typename T>
258 0 : void SparseMatrix<T>::zero_rows (std::vector<numeric_index_type> &, T)
259 : {
260 : /* This functionality isn't implemented or stubbed in every subclass yet */
261 0 : libmesh_not_implemented();
262 : }
263 :
264 :
265 : template <typename T>
266 28 : std::size_t SparseMatrix<T>::n_nonzeros() const
267 : {
268 28 : if (!_sp)
269 4 : return 0;
270 0 : return _sp->n_nonzeros();
271 : }
272 :
273 :
274 : template <typename T>
275 0 : void SparseMatrix<T>::print(std::ostream & os, const bool sparse) const
276 : {
277 0 : parallel_object_only();
278 :
279 0 : libmesh_assert (this->initialized());
280 :
281 0 : const numeric_index_type first_dof = this->row_start(),
282 0 : end_dof = this->row_stop();
283 :
284 : // We'll print the matrix from processor 0 to make sure
285 : // it's serialized properly
286 0 : if (this->processor_id() == 0)
287 : {
288 0 : libmesh_assert_equal_to (first_dof, 0);
289 0 : for (numeric_index_type i : make_range(end_dof))
290 : {
291 0 : if (sparse)
292 : {
293 0 : for (auto j : make_range(this->n()))
294 : {
295 0 : T c = (*this)(i,j);
296 0 : if (c != static_cast<T>(0.0))
297 : {
298 0 : os << i << " " << j << " " << c << std::endl;
299 : }
300 : }
301 : }
302 : else
303 : {
304 0 : for (auto j : make_range(this->n()))
305 0 : os << (*this)(i,j) << " ";
306 0 : os << std::endl;
307 : }
308 : }
309 :
310 0 : std::vector<numeric_index_type> ibuf, jbuf;
311 0 : std::vector<T> cbuf;
312 0 : numeric_index_type currenti = end_dof;
313 0 : for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
314 : {
315 0 : this->comm().receive(p, ibuf);
316 0 : this->comm().receive(p, jbuf);
317 0 : this->comm().receive(p, cbuf);
318 0 : libmesh_assert_equal_to (ibuf.size(), jbuf.size());
319 0 : libmesh_assert_equal_to (ibuf.size(), cbuf.size());
320 :
321 0 : if (ibuf.empty())
322 0 : continue;
323 0 : libmesh_assert_greater_equal (ibuf.front(), currenti);
324 0 : libmesh_assert_greater_equal (ibuf.back(), ibuf.front());
325 :
326 0 : std::size_t currentb = 0;
327 0 : for (;currenti <= ibuf.back(); ++currenti)
328 : {
329 0 : if (sparse)
330 : {
331 0 : for (numeric_index_type j=0; j<this->n(); j++)
332 : {
333 0 : if (currentb < ibuf.size() &&
334 0 : ibuf[currentb] == currenti &&
335 0 : jbuf[currentb] == j)
336 : {
337 0 : os << currenti << " " << j << " " << cbuf[currentb] << std::endl;
338 0 : currentb++;
339 : }
340 : }
341 : }
342 : else
343 : {
344 0 : for (auto j : make_range(this->n()))
345 : {
346 0 : if (currentb < ibuf.size() &&
347 0 : ibuf[currentb] == currenti &&
348 0 : jbuf[currentb] == j)
349 : {
350 0 : os << cbuf[currentb] << " ";
351 0 : currentb++;
352 : }
353 : else
354 0 : os << static_cast<T>(0.0) << " ";
355 : }
356 0 : os << std::endl;
357 : }
358 : }
359 : }
360 0 : if (!sparse)
361 : {
362 0 : for (; currenti != this->m(); ++currenti)
363 : {
364 0 : for (numeric_index_type j=0; j<this->n(); j++)
365 0 : os << static_cast<T>(0.0) << " ";
366 0 : os << std::endl;
367 : }
368 : }
369 : }
370 : else
371 : {
372 0 : std::vector<numeric_index_type> ibuf, jbuf;
373 0 : std::vector<T> cbuf;
374 :
375 : // We'll assume each processor has access to entire
376 : // matrix rows, so (*this)(i,j) is valid if i is a local index.
377 0 : for (numeric_index_type i : make_range(first_dof, end_dof))
378 : {
379 0 : for (auto j : make_range(this->n()))
380 : {
381 0 : T c = (*this)(i,j);
382 0 : if (c != static_cast<T>(0.0))
383 : {
384 0 : ibuf.push_back(i);
385 0 : jbuf.push_back(j);
386 0 : cbuf.push_back(c);
387 : }
388 : }
389 : }
390 0 : this->comm().send(0,ibuf);
391 0 : this->comm().send(0,jbuf);
392 0 : this->comm().send(0,cbuf);
393 : }
394 0 : }
395 :
396 :
397 : template <typename T>
398 28 : void SparseMatrix<T>::print_matlab(const std::string & name) const
399 : {
400 4 : parallel_object_only();
401 :
402 4 : libmesh_assert (this->initialized());
403 :
404 28 : const numeric_index_type first_dof = this->row_start(),
405 28 : end_dof = this->row_stop();
406 :
407 : // We'll print the matrix from processor 0 to make sure
408 : // it's serialized properly
409 32 : if (this->processor_id() == 0)
410 : {
411 24 : std::unique_ptr<std::ofstream> file;
412 :
413 28 : if (name != "")
414 48 : file = std::make_unique<std::ofstream>(name.c_str());
415 :
416 28 : std::ostream & os = (name == "") ? libMesh::out : *file;
417 :
418 28 : std::size_t sparsity_nonzeros = this->n_nonzeros();
419 :
420 4 : std::size_t real_nonzeros = 0;
421 :
422 4 : libmesh_assert_equal_to(first_dof, 0);
423 140 : for (numeric_index_type i : make_range(end_dof))
424 : {
425 640 : for (auto j : make_range(this->n()))
426 : {
427 448 : T c = (*this)(i,j);
428 416 : if (c != static_cast<T>(0.0))
429 448 : ++real_nonzeros;
430 : }
431 : }
432 :
433 :
434 28 : for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
435 : {
436 0 : std::size_t nonzeros_on_p = 0;
437 0 : this->comm().receive(p, nonzeros_on_p);
438 0 : real_nonzeros += nonzeros_on_p;
439 : }
440 :
441 : if (sparsity_nonzeros &&
442 : sparsity_nonzeros != real_nonzeros)
443 : libmesh_warning(sparsity_nonzeros <<
444 : " nonzeros allocated, but " <<
445 : real_nonzeros << " used.");
446 :
447 : // We probably want to be more consistent than that, if our
448 : // sparsity is overallocated.
449 :
450 : // Print a header similar to PETSc's mat_view ascii_matlab
451 28 : os << "%Mat Object: () " << this->n_processors() << " MPI processes\n"
452 12 : << "% type: " << (this->n_processors() > 1 ? "mpi" : "seq") << "aij\n"
453 80 : << "% Size = " << this->m() << ' ' << this->n() << '\n'
454 28 : << "% Nonzeros = " << real_nonzeros << '\n'
455 28 : << "zzz = zeros(" << real_nonzeros << ",3);\n"
456 28 : << "zzz = [\n";
457 :
458 140 : for (numeric_index_type i : make_range(end_dof))
459 : {
460 : // FIXME - we need a base class way to iterate over a
461 : // SparseMatrix row.
462 640 : for (auto j : make_range(this->n()))
463 : {
464 448 : T c = (*this)(i,j);
465 416 : if (c != static_cast<T>(0.0))
466 : {
467 : // Convert from 0-based to 1-based indexing
468 864 : os << (i+1) << ' ' << (j+1) << " " << c << '\n';
469 : }
470 : }
471 : }
472 :
473 8 : std::vector<numeric_index_type> ibuf, jbuf;
474 8 : std::vector<T> cbuf;
475 28 : for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
476 : {
477 0 : this->comm().receive(p, ibuf);
478 0 : this->comm().receive(p, jbuf);
479 0 : this->comm().receive(p, cbuf);
480 0 : libmesh_assert_equal_to (ibuf.size(), jbuf.size());
481 0 : libmesh_assert_equal_to (ibuf.size(), cbuf.size());
482 :
483 0 : for (auto n : index_range(ibuf))
484 0 : os << ibuf[n] << ' ' << jbuf[n] << " " << cbuf[n] << '\n';
485 : }
486 :
487 24 : os << "];\n" << "Mat_sparse = spconvert(zzz);" << std::endl;
488 20 : }
489 : else
490 : {
491 0 : std::vector<numeric_index_type> ibuf, jbuf;
492 0 : std::vector<T> cbuf;
493 0 : std::size_t my_nonzeros = 0;
494 :
495 : // We'll assume each processor has access to entire
496 : // matrix rows, so (*this)(i,j) is valid if i is a local index.
497 0 : for (numeric_index_type i : make_range(first_dof, end_dof))
498 : {
499 0 : for (auto j : make_range(this->n()))
500 : {
501 0 : T c = (*this)(i,j);
502 0 : if (c != static_cast<T>(0.0))
503 : {
504 0 : ibuf.push_back(i);
505 0 : jbuf.push_back(j);
506 0 : cbuf.push_back(c);
507 0 : ++my_nonzeros;
508 : }
509 : }
510 : }
511 0 : this->comm().send(0,my_nonzeros);
512 0 : this->comm().send(0,ibuf);
513 0 : this->comm().send(0,jbuf);
514 0 : this->comm().send(0,cbuf);
515 : }
516 28 : }
517 :
518 :
519 :
520 : template <typename T>
521 3 : void SparseMatrix<T>::print_coreform_hdf5(const std::string & filename,
522 : const std::string & groupname) const
523 : {
524 : #if defined(LIBMESH_USE_COMPLEX_NUMBERS) || !defined(LIBMESH_HAVE_HDF5)
525 : // TODO: HDF5 version 2.0.0 and later adds native support for
526 : // complex numbers with the H5T_NATIVE_DOUBLE_COMPLEX type [0], so
527 : // we could consider supporting T==std::complex<Real> in the future.
528 : // [0]: https://forum.hdfgroup.org/t/coming-in-the-next-hdf5-release-native-support-for-complex-number-datatypes/13543
529 0 : libmesh_ignore(filename, groupname);
530 0 : libmesh_error_msg("ERROR: need HDF5 support to handle .h5 files!!!");
531 : #else
532 : LOG_SCOPE("print_coreform_hdf5()", "SparseMatrix");
533 :
534 : // In this implementation, we copy the SparseMatrix entries into a
535 : // std::vector<double>, so this won't work for any Number type for
536 : // which sizeof(Number) > sizeof(double).
537 : if constexpr (sizeof(T) > sizeof(double))
538 : libmesh_not_implemented();
539 :
540 3 : const numeric_index_type first_dof = this->row_start(),
541 3 : end_dof = this->row_stop();
542 :
543 : std::vector<std::size_t> cols, row_offsets;
544 : std::vector<double> vals;
545 :
546 3 : if (this->processor_id() == 0)
547 2 : row_offsets.push_back(0);
548 :
549 17 : for (numeric_index_type i : make_range(first_dof, end_dof))
550 : {
551 174 : for (auto j : make_range(this->n()))
552 : {
553 146 : T c = (*this)(i,j);
554 146 : if (c != static_cast<T>(0.0))
555 : {
556 81 : cols.push_back(j);
557 81 : vals.push_back(c);
558 : }
559 : }
560 : // This is a *local* row offset; proc 0 may need to adjust later
561 14 : row_offsets.push_back(cols.size());
562 : }
563 :
564 3 : if (this->processor_id() == 0)
565 : {
566 2 : const hid_t file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC,
567 : H5P_DEFAULT, H5P_DEFAULT);
568 2 : if (file == H5I_INVALID_HID)
569 0 : libmesh_file_error(filename);
570 :
571 2 : const hid_t group =
572 2 : H5Gcreate2(file, groupname.c_str(), H5P_DEFAULT, H5P_DEFAULT,
573 : H5P_DEFAULT);
574 2 : check_open(filename, group, groupname);
575 :
576 8 : auto write_size_attribute = [&filename, &group]
577 : (const std::string & attribute_name, unsigned long long writeval)
578 : {
579 4 : const hid_t fspace = H5Screate(H5S_SCALAR);
580 4 : check_hdf5(filename, fspace, attribute_name + " fspace");
581 :
582 4 : const hid_t attr = H5Acreate2(group, attribute_name.c_str(),
583 4 : H5T_STD_I64LE, fspace,
584 : H5P_DEFAULT, H5P_DEFAULT);
585 4 : check_hdf5(filename, attr, attribute_name);
586 :
587 : // HDF5 is supposed to handle both upscaling and endianness
588 : // conversions here
589 4 : const herr_t errval = H5Awrite(attr, H5T_NATIVE_ULLONG, &writeval);
590 4 : check_hdf5(filename, errval, attribute_name + " write");
591 :
592 4 : H5Aclose(attr);
593 4 : H5Sclose(fspace);
594 : };
595 :
596 2 : write_size_attribute("num_cols", this->n());
597 2 : write_size_attribute("num_rows", this->m());
598 :
599 14 : auto write_vector = [&filename, &group]
600 : (const std::string & dataname, auto hdf5_file_type,
601 : auto hdf5_native_type, auto & datavec)
602 : {
603 6 : const hsize_t len[1] = {cast_int<hsize_t>(datavec.size())};
604 :
605 6 : const hid_t space = H5Screate_simple(1, len, nullptr);
606 12 : check_hdf5(filename, space, dataname + " space");
607 :
608 : const hid_t data =
609 6 : H5Dcreate2(group, dataname.c_str(), hdf5_file_type, space,
610 : H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
611 12 : check_hdf5(filename, data, dataname + " data");
612 :
613 6 : const hid_t errval =
614 6 : H5Dwrite(data, hdf5_native_type, H5S_ALL, H5S_ALL,
615 : H5P_DEFAULT, datavec.data());
616 6 : check_hdf5(filename, errval, dataname + " write");
617 :
618 6 : H5Dclose(data);
619 6 : H5Sclose(space);
620 : };
621 :
622 : std::vector<std::size_t> vals_sizes, first_dofs;
623 2 : this->comm().gather(0, vals.size(), vals_sizes);
624 2 : this->comm().gather(0, cast_int<std::size_t>(first_dof), first_dofs);
625 2 : first_dofs.push_back(this->m());
626 :
627 2 : this->comm().allgather(cols);
628 2 : this->comm().allgather(vals);
629 2 : this->comm().allgather(row_offsets);
630 :
631 : libmesh_assert_equal_to(vals.size(),
632 : cols.size());
633 : libmesh_assert_equal_to(vals.size(),
634 : std::accumulate(vals_sizes.begin(),
635 : vals_sizes.end(),
636 : std::size_t(0)));
637 :
638 : std::size_t extra_offset = 0;
639 3 : for (auto p : make_range(processor_id_type(1), this->n_processors()))
640 : {
641 1 : extra_offset += vals_sizes[p-1];
642 6 : for (auto i : make_range(first_dofs[p]+1, first_dofs[p+1]+1))
643 5 : row_offsets[i] += extra_offset;
644 : }
645 :
646 2 : write_vector("cols", H5T_STD_U64LE, H5T_NATIVE_ULLONG, cols);
647 2 : write_vector("row_offsets", H5T_STD_U64LE, H5T_NATIVE_ULLONG, row_offsets);
648 2 : write_vector("vals", H5T_IEEE_F64LE, H5T_NATIVE_DOUBLE, vals);
649 :
650 2 : H5Gclose(group);
651 2 : H5Fclose(file);
652 : }
653 : else
654 : {
655 : std::vector<std::size_t> dummy;
656 1 : this->comm().gather(0, vals.size(), dummy);
657 1 : this->comm().gather(0, cast_int<std::size_t>(first_dof), dummy);
658 1 : this->comm().allgather(cols);
659 1 : this->comm().allgather(vals);
660 1 : this->comm().allgather(row_offsets);
661 : }
662 : #endif // LIBMESH_HAVE_HDF5
663 3 : }
664 :
665 :
666 :
667 : template <typename T>
668 0 : void SparseMatrix<T>::print_petsc_binary(const std::string &) const
669 : {
670 0 : libmesh_not_implemented_msg
671 : ("libMesh cannot write PETSc binary-format files from non-PETSc matrices");
672 : }
673 :
674 :
675 :
676 : template <typename T>
677 0 : void SparseMatrix<T>::print_petsc_hdf5(const std::string &) const
678 : {
679 0 : libmesh_not_implemented_msg
680 : ("libMesh cannot write PETSc HDF5-format files from non-PETSc matrices");
681 : }
682 :
683 :
684 :
685 : template <typename T>
686 1216 : void SparseMatrix<T>::read(const std::string & filename)
687 : {
688 : {
689 1292 : std::ifstream in (filename.c_str());
690 1216 : libmesh_error_msg_if
691 : (!in.good(), "ERROR: cannot read file:\n\t" <<
692 : filename);
693 1140 : }
694 :
695 1216 : std::string_view basename = Utility::basename_of(filename);
696 :
697 1216 : const bool gzipped_file = Utility::ends_with(filename, ".gz");
698 :
699 1216 : if (gzipped_file)
700 18 : basename.remove_suffix(3);
701 :
702 2432 : if (Utility::ends_with(basename, ".matlab") ||
703 1330 : Utility::ends_with(basename, ".m"))
704 1138 : this->read_matlab(filename);
705 78 : else if (Utility::ends_with(basename, ".petsc64"))
706 : {
707 : #ifndef LIBMESH_HAVE_PETSC
708 0 : libmesh_error_msg("Cannot load PETSc matrix file " <<
709 : filename << " without PETSc-enabled libMesh.");
710 : #elif LIBMESH_DOF_ID_BYTES != 8
711 0 : libmesh_error_msg("Cannot load 64-bit PETSc matrix file " <<
712 : filename << " with non-64-bit libMesh.");
713 : #endif
714 66 : if (gzipped_file)
715 0 : libmesh_not_implemented_msg("Gzipped PETSc matrices are not currently supported");
716 66 : this->read_petsc_binary(filename);
717 : }
718 12 : else if (Utility::ends_with(basename, ".petsc32"))
719 : {
720 : #ifndef LIBMESH_HAVE_PETSC
721 0 : libmesh_error_msg("Cannot load PETSc matrix file " <<
722 : filename << " without PETSc-enabled libMesh.");
723 : #elif LIBMESH_DOF_ID_BYTES != 4
724 0 : libmesh_error_msg("Cannot load 32-bit PETSc matrix file " <<
725 : filename << " with non-32-bit libMesh.");
726 : #endif
727 4 : if (gzipped_file)
728 0 : libmesh_not_implemented_msg("Gzipped PETSc matrices are not currently supported");
729 4 : this->read_petsc_binary(filename);
730 : }
731 8 : else if (Utility::ends_with(basename, ".h5"))
732 : {
733 8 : if (gzipped_file)
734 0 : libmesh_not_implemented_msg("Gzipped HDF5 matrices are not currently supported");
735 16 : this->read_coreform_hdf5(filename);
736 : }
737 : else
738 0 : libmesh_error_msg(" ERROR: Unrecognized matrix file extension on: "
739 : << basename
740 : << "\n I understand the following:\n\n"
741 : << " *.h5 -- CoreForm HDF5 sparse matrix format\n"
742 : << " *.matlab -- Matlab sparse matrix format\n"
743 : << " *.matlab.gz -- Matlab sparse matrix format, gzipped\n"
744 : << " *.m -- Matlab sparse matrix format\n"
745 : << " *.m.gz -- Matlab sparse matrix format, gzipped\n"
746 : << " *.petsc32 -- PETSc binary format, 32-bit\n"
747 : << " *.petsc64 -- PETSc binary format, 64-bit\n"
748 : );
749 1216 : }
750 :
751 :
752 :
753 : template <typename T>
754 171 : void SparseMatrix<T>::print(const std::string & filename) const
755 : {
756 : {
757 187 : std::ofstream outstr (filename.c_str());
758 171 : libmesh_error_msg_if
759 : (!outstr.good(), "ERROR: cannot write to file:\n\t" <<
760 : filename);
761 155 : }
762 :
763 171 : std::string_view basename = Utility::basename_of(filename);
764 :
765 171 : const bool gzipped_file = Utility::ends_with(filename, ".gz");
766 :
767 171 : if (gzipped_file)
768 4 : basename.remove_suffix(3);
769 :
770 342 : if (Utility::ends_with(basename, ".matlab") ||
771 182 : Utility::ends_with(basename, ".m"))
772 168 : this->print_matlab(filename);
773 3 : else if (Utility::ends_with(basename, ".petsc64"))
774 : {
775 : #ifndef LIBMESH_HAVE_PETSC
776 0 : libmesh_error_msg("Cannot load PETSc matrix file " <<
777 : filename << " without PETSc-enabled libMesh.");
778 : #elif LIBMESH_DOF_ID_BYTES != 8
779 0 : libmesh_error_msg("Cannot load 64-bit PETSc matrix file " <<
780 : filename << " with non-64-bit libMesh.");
781 : #endif
782 0 : if (gzipped_file)
783 0 : libmesh_not_implemented_msg("Gzipped PETSc matrices are not currently supported");
784 0 : this->print_petsc_binary(filename);
785 : }
786 3 : else if (Utility::ends_with(basename, ".petsc32"))
787 : {
788 : #ifndef LIBMESH_HAVE_PETSC
789 0 : libmesh_error_msg("Cannot load PETSc matrix file " <<
790 : filename << " without PETSc-enabled libMesh.");
791 : #elif LIBMESH_DOF_ID_BYTES != 4
792 0 : libmesh_error_msg("Cannot load 32-bit PETSc matrix file " <<
793 : filename << " with non-32-bit libMesh.");
794 : #endif
795 0 : if (gzipped_file)
796 0 : libmesh_not_implemented_msg("Gzipped PETSc matrices are not currently supported");
797 0 : this->print_petsc_binary(filename);
798 : }
799 3 : else if (Utility::ends_with(basename, ".h5"))
800 : {
801 3 : if (gzipped_file)
802 0 : libmesh_not_implemented_msg("Gzipped HDF5 matrices are not currently supported");
803 6 : this->print_coreform_hdf5(filename);
804 : }
805 : else
806 0 : libmesh_error_msg(" ERROR: Unrecognized matrix file extension on: "
807 : << basename
808 : << "\n I understand the following:\n\n"
809 : << " *.h5 -- CoreForm HDF5 sparse matrix format\n"
810 : << " *.matlab -- Matlab sparse matrix format\n"
811 : << " *.matlab.gz -- Matlab sparse matrix format, gzipped\n"
812 : << " *.m -- Matlab sparse matrix format\n"
813 : << " *.m.gz -- Matlab sparse matrix format, gzipped\n"
814 : << " *.petsc32 -- PETSc binary format, 32-bit\n"
815 : << " *.petsc64 -- PETSc binary format, 64-bit\n"
816 : );
817 171 : }
818 :
819 :
820 :
821 : template <typename T>
822 8 : void SparseMatrix<T>::read_coreform_hdf5(const std::string & filename,
823 : const std::string & groupname)
824 : {
825 : #ifndef LIBMESH_HAVE_HDF5
826 0 : libmesh_ignore(filename, groupname);
827 0 : libmesh_error_msg("ERROR: need HDF5 support to handle .h5 files!!!");
828 : #else
829 : LOG_SCOPE("read_coreform_hdf5()", "SparseMatrix");
830 :
831 8 : std::size_t num_rows = 0, num_cols = 0;
832 :
833 : // These are only used on pid 0, but avoid "uninitialized" warnings
834 8 : hid_t group = H5I_INVALID_HID;
835 : hid_t file = H5I_INVALID_HID;
836 :
837 8 : if (this->processor_id() == 0)
838 : {
839 6 : file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
840 :
841 6 : if (file == H5I_INVALID_HID)
842 0 : libmesh_file_error(filename);
843 :
844 6 : group = H5Gopen(file, groupname.c_str(), H5P_DEFAULT);
845 6 : check_open(filename, group, groupname);
846 :
847 54 : auto read_size_attribute = [&filename, &group]
848 : (const std::string & attribute_name)
849 : {
850 12 : unsigned long long returnval = 0;
851 :
852 12 : const hid_t attr = H5Aopen(group, attribute_name.c_str(), H5P_DEFAULT);
853 12 : check_open(filename, attr, attribute_name);
854 :
855 12 : const hid_t attr_type = H5Aget_type(attr);
856 12 : check_hdf5(filename, attr_type, attribute_name + " type");
857 :
858 : // HDF5 will convert between the file's integer type and ours, but
859 : // we do expect an integer type.
860 12 : if (H5Tget_class(attr_type) != H5T_INTEGER)
861 0 : libmesh_error_msg("Non-integer type for " + attribute_name + " in " + filename);
862 :
863 12 : H5Tclose(attr_type);
864 :
865 : // HDF5 is supposed to handle both upscaling and endianness
866 : // conversions here
867 12 : const herr_t errval = H5Aread(attr, H5T_NATIVE_ULLONG, &returnval);
868 12 : check_hdf5(filename, errval, attribute_name + " read");
869 :
870 12 : H5Aclose(attr);
871 :
872 12 : return returnval;
873 : };
874 :
875 6 : num_cols = read_size_attribute("num_cols");
876 6 : num_rows = read_size_attribute("num_rows");
877 :
878 6 : this->comm().broadcast(num_cols);
879 6 : this->comm().broadcast(num_rows);
880 : }
881 : else
882 : {
883 2 : this->comm().broadcast(num_cols);
884 2 : this->comm().broadcast(num_rows);
885 : }
886 :
887 : numeric_index_type new_row_start, new_row_stop,
888 : new_col_start, new_col_stop;
889 :
890 : // If we need to reinit, we need to determine which rows+columns
891 : // each processor is in charge of.
892 : std::vector<numeric_index_type> new_row_starts, new_row_stops,
893 : new_col_starts, new_col_stops;
894 :
895 8 : if (this->initialized() &&
896 8 : num_cols == this->n() &&
897 0 : num_rows == this->m())
898 : {
899 0 : new_row_start = this->row_start(),
900 0 : new_row_stop = this->row_stop();
901 :
902 0 : new_col_start = this->col_start(),
903 0 : new_col_stop = this->col_stop();
904 : }
905 : else
906 : {
907 : // Determine which rows/columns each processor will be in charge of
908 8 : new_row_start = this->processor_id() * num_rows / this->n_processors(),
909 8 : new_row_stop = (this->processor_id()+1) * num_rows / this->n_processors();
910 :
911 8 : new_col_start = this->processor_id() * num_cols / this->n_processors(),
912 8 : new_col_stop = (this->processor_id()+1) * num_cols / this->n_processors();
913 : }
914 :
915 8 : this->comm().gather(0, new_row_start, new_row_starts);
916 8 : this->comm().gather(0, new_row_stop, new_row_stops);
917 8 : this->comm().gather(0, new_col_start, new_col_starts);
918 8 : this->comm().gather(0, new_col_stop, new_col_stops);
919 :
920 8 : numeric_index_type on_diagonal_nonzeros = 0,
921 8 : off_diagonal_nonzeros = 0;
922 :
923 : std::vector<std::size_t> cols, row_offsets;
924 : std::vector<double> vals;
925 :
926 8 : if (this->processor_id() == 0)
927 : {
928 42 : auto read_vector = [&filename, &group]
929 : (const std::string & dataname, auto hdf5_class,
930 : auto hdf5_type, auto & datavec)
931 : {
932 18 : const hid_t data = H5Dopen1(group, dataname.c_str());
933 18 : check_open(filename, data, dataname.c_str());
934 :
935 18 : const hid_t data_type = H5Dget_type(data);
936 18 : check_hdf5(filename, data_type, dataname + " type");
937 :
938 : // HDF5 will convert between the file's integer type and ours, but
939 : // we do expect an integer type.
940 18 : if (H5Tget_class(data_type) != hdf5_class)
941 0 : libmesh_error_msg("Unexpected type for " + dataname + " in " + filename);
942 :
943 18 : H5Tclose(data_type);
944 :
945 18 : const hid_t dataspace = H5Dget_space(data);
946 18 : check_hdf5(filename, dataspace, dataname + " space");
947 :
948 18 : int ndims = H5Sget_simple_extent_ndims(dataspace);
949 18 : if (ndims != 1)
950 0 : libmesh_error_msg("Non-vector space for " + dataname + " in " + filename);
951 :
952 : hsize_t len, maxlen;
953 18 : herr_t errval = H5Sget_simple_extent_dims(dataspace, &len, &maxlen);
954 18 : check_hdf5(filename, errval, dataname + " dims");
955 :
956 18 : datavec.resize(len);
957 :
958 18 : errval = H5Dread(data, hdf5_type, H5S_ALL, H5S_ALL, H5P_DEFAULT, datavec.data());
959 18 : check_hdf5(filename, errval, dataname + " read");
960 :
961 18 : H5Dclose(data);
962 : };
963 :
964 6 : read_vector("cols", H5T_INTEGER, H5T_NATIVE_ULLONG, cols);
965 6 : read_vector("row_offsets", H5T_INTEGER, H5T_NATIVE_ULLONG, row_offsets);
966 12 : read_vector("vals", H5T_FLOAT, H5T_NATIVE_DOUBLE, vals);
967 :
968 6 : if (cols.size() != vals.size())
969 0 : libmesh_error_msg("Inconsistent cols/vals sizes in " + filename);
970 :
971 6 : if (row_offsets.size() != num_rows + 1)
972 0 : libmesh_error_msg("Inconsistent row_offsets size in " + filename);
973 :
974 : // Data for the row we're working on
975 : numeric_index_type current_row = 0;
976 : processor_id_type current_proc = 0;
977 6 : numeric_index_type current_on_diagonal_nonzeros = 0;
978 6 : numeric_index_type current_off_diagonal_nonzeros = 0;
979 6 : if (row_offsets[0] != 0)
980 0 : libmesh_error_msg("Unexpected row_offsets[0] in " + filename);
981 :
982 481 : for (auto i : index_range(cols))
983 : {
984 568 : while (row_offsets[current_row+1] <= i)
985 : {
986 : ++current_row;
987 93 : if (row_offsets[current_row] < row_offsets[current_row-1])
988 0 : libmesh_error_msg("Non-monotonic row_offsets in " + filename);
989 93 : current_on_diagonal_nonzeros = 0;
990 93 : current_off_diagonal_nonzeros = 0;
991 : }
992 :
993 477 : while (current_row >= new_row_stops[current_proc])
994 2 : ++current_proc;
995 :
996 : // 0-based indexing in file
997 475 : if (cols[i] >= new_col_starts[current_proc] &&
998 449 : cols[i] < new_col_stops[current_proc])
999 : {
1000 414 : ++current_on_diagonal_nonzeros;
1001 414 : on_diagonal_nonzeros =
1002 : std::max(on_diagonal_nonzeros,
1003 : current_on_diagonal_nonzeros);
1004 : }
1005 : else
1006 : {
1007 61 : ++current_off_diagonal_nonzeros;
1008 61 : off_diagonal_nonzeros =
1009 : std::max(off_diagonal_nonzeros,
1010 : current_off_diagonal_nonzeros);
1011 : }
1012 : }
1013 : }
1014 :
1015 8 : this->comm().broadcast(on_diagonal_nonzeros);
1016 8 : this->comm().broadcast(off_diagonal_nonzeros);
1017 :
1018 8 : this->init(num_rows, num_cols,
1019 : new_row_stop-new_row_start,
1020 : new_col_stop-new_col_start,
1021 : on_diagonal_nonzeros,
1022 : off_diagonal_nonzeros);
1023 :
1024 : // Set the matrix values last.
1025 8 : if (this->processor_id() == 0)
1026 : {
1027 : numeric_index_type current_row = 0;
1028 481 : for (auto i : index_range(cols))
1029 : {
1030 568 : while (row_offsets[current_row+1] <= i)
1031 : {
1032 : ++current_row;
1033 : libmesh_assert_greater_equal (row_offsets[current_row],
1034 : row_offsets[current_row-1]);
1035 : }
1036 475 : this->set(current_row, cols[i], vals[i]);
1037 : }
1038 :
1039 6 : H5Gclose(group);
1040 6 : H5Fclose(file);
1041 : }
1042 :
1043 8 : this->close();
1044 : #endif // LIBMESH_HAVE_HDF5
1045 8 : }
1046 :
1047 :
1048 :
1049 : template <typename T>
1050 1559 : void SparseMatrix<T>::read_matlab(const std::string & filename)
1051 : {
1052 96 : LOG_SCOPE("read_matlab()", "SparseMatrix");
1053 :
1054 : #ifndef LIBMESH_HAVE_CXX11_REGEX
1055 : libmesh_not_implemented(); // What is your compiler?!? Email us!
1056 : libmesh_ignore(filename);
1057 : #else
1058 48 : parallel_object_only();
1059 :
1060 1559 : const bool gzipped_file = Utility::ends_with(filename, ".gz");
1061 :
1062 : // The sizes we get from the file
1063 1559 : std::size_t m = 0,
1064 1559 : n = 0;
1065 :
1066 : // If we don't already have this size, we'll need to reinit, and
1067 : // determine which rows+columns each processor is in charge of.
1068 96 : std::vector<numeric_index_type> new_row_starts, new_row_stops,
1069 96 : new_col_starts, new_col_stops;
1070 :
1071 : numeric_index_type new_row_start, new_row_stop,
1072 : new_col_start, new_col_stop;
1073 :
1074 : // We'll read through the file three times: once to get a reliable
1075 : // value for the matrix size (so we can divvy it up among
1076 : // processors), then again to get the sparsity to send to each
1077 : // processor, then a final time to get the entries to send to each
1078 : // processor.
1079 : //
1080 : // We'll use an istream here; it might be an ifstream if we're
1081 : // opening a raw ASCII file or a gzstream if we're opening a
1082 : // compressed one.
1083 1511 : std::unique_ptr<std::istream> file;
1084 :
1085 : // We'll need a temporary structure to cache matrix entries, because
1086 : // we need to read through the whole file before we know the size
1087 : // and sparsity structure with which we can init().
1088 : //
1089 : // Reading through the file three times via `seekg` doesn't work
1090 : // with our gzstream wrapper, and seems to take three times as long
1091 : // even with a plain ifstream. What happened to disk caching!?
1092 96 : std::vector<std::tuple<numeric_index_type, numeric_index_type, T>> entries;
1093 :
1094 : // First read through the file, saving size and entry data
1095 : {
1096 : // We'll read the matrix on processor 0 rather than try to juggle
1097 : // parallel I/O.
1098 1607 : if (this->processor_id() == 0)
1099 : {
1100 : // We'll be using regular expressions to make ourselves slightly
1101 : // more robust to formatting.
1102 800 : const std::regex start_regex // assignment like "zzz = ["
1103 : ("\\s*\\w+\\s*=\\s*\\[");
1104 800 : const std::regex end_regex // end of assignment
1105 : ("^[^%]*\\]");
1106 :
1107 732 : if (gzipped_file)
1108 : {
1109 : #ifdef LIBMESH_HAVE_GZSTREAM
1110 345 : auto inf = std::make_unique<igzstream>();
1111 14 : libmesh_assert(inf);
1112 14 : inf->open(filename.c_str(), std::ios::in);
1113 317 : file = std::move(inf);
1114 : #else
1115 : libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
1116 : #endif
1117 303 : }
1118 : else
1119 : {
1120 421 : auto inf = std::make_unique<std::ifstream>();
1121 20 : libmesh_assert(inf);
1122 :
1123 421 : std::string new_name = Utility::unzip_file(filename);
1124 :
1125 401 : inf->open(new_name.c_str(), std::ios::in);
1126 381 : file = std::move(inf);
1127 361 : }
1128 :
1129 : // If we have a matrix with all-zero trailing rows, the only
1130 : // way to get the size is if it ended up in a comment
1131 800 : const std::regex size_regex // comment like "% size = 8 8"
1132 : ("%\\s*[Ss][Ii][Zz][Ee]\\s*=\\s*(\\d+)\\s+(\\d+)");
1133 766 : const std::string whitespace = " \t";
1134 :
1135 34 : bool have_started = false;
1136 34 : bool have_ended = false;
1137 732 : std::size_t largest_i_seen = 0, largest_j_seen = 0;
1138 :
1139 : // Data for the row we're working on
1140 : // Use 1-based indexing for current_row, as in the file
1141 732 : std::size_t current_row = 1;
1142 :
1143 343420 : for (std::string line; std::getline(*file, line);)
1144 : {
1145 7780 : std::smatch sm;
1146 :
1147 : // First, try to match an entry. This is the most common
1148 : // case so we won't rely on slow std::regex for it.
1149 : // stringstream is at least an improvement over that.
1150 :
1151 : // Look for row/col/val like "1 1 -2.0e-4"
1152 :
1153 187219 : std::istringstream l(line);
1154 :
1155 : std::size_t i, j;
1156 1822 : T value;
1157 :
1158 9602 : l >> i >> j >> value;
1159 :
1160 179473 : if (!l.fail())
1161 : {
1162 174349 : libmesh_error_msg_if
1163 : (!have_started, "Confused by premature entries in matrix file " << filename);
1164 :
1165 331841 : entries.emplace_back(cast_int<numeric_index_type>(i),
1166 181891 : cast_int<numeric_index_type>(j),
1167 : value);
1168 :
1169 174349 : libmesh_error_msg_if
1170 : (!i || !j, "Expected 1-based indexing in matrix file "
1171 : << filename);
1172 :
1173 174349 : current_row = std::max(current_row, i);
1174 :
1175 174349 : libmesh_error_msg_if
1176 : (i < current_row,
1177 : "Can't handle out-of-order entries in matrix file "
1178 : << filename);
1179 :
1180 174349 : largest_i_seen = std::max(i, largest_i_seen);
1181 174349 : largest_j_seen = std::max(j, largest_j_seen);
1182 : }
1183 :
1184 5124 : else if (std::regex_search(line, sm, size_regex))
1185 : {
1186 766 : const std::string msize = sm[1];
1187 732 : const std::string nsize = sm[2];
1188 732 : m = std::stoull(msize);
1189 732 : n = std::stoull(nsize);
1190 : }
1191 :
1192 4392 : else if (std::regex_search(line, start_regex))
1193 34 : have_started = true;
1194 :
1195 3660 : else if (std::regex_search(line, end_regex))
1196 : {
1197 34 : have_ended = true;
1198 68 : break;
1199 : }
1200 : }
1201 :
1202 732 : libmesh_error_msg_if
1203 : (!have_started, "Confused by missing assignment beginning in matrix file " << filename);
1204 :
1205 732 : libmesh_error_msg_if
1206 : (!have_ended, "Confused by missing assignment ending in matrix file " << filename);
1207 :
1208 732 : libmesh_error_msg_if
1209 : (m > largest_i_seen, "Confused by missing final row(s) in matrix file " << filename);
1210 :
1211 732 : libmesh_error_msg_if
1212 : (m > 0 && m < largest_i_seen, "Confused by extra final row(s) in matrix file " << filename);
1213 :
1214 732 : if (!m)
1215 0 : m = largest_i_seen;
1216 :
1217 732 : libmesh_error_msg_if
1218 : (n > largest_j_seen, "Confused by missing final column(s) in matrix file " << filename);
1219 :
1220 732 : libmesh_error_msg_if
1221 : (n > 0 && n < largest_j_seen, "Confused by extra final column(s) in matrix file " << filename);
1222 :
1223 732 : if (!n)
1224 0 : n = largest_j_seen;
1225 :
1226 732 : this->comm().broadcast(m);
1227 732 : this->comm().broadcast(n);
1228 664 : }
1229 : else
1230 : {
1231 813 : this->comm().broadcast(m);
1232 827 : this->comm().broadcast(n);
1233 : }
1234 :
1235 1559 : if (this->initialized() &&
1236 1561 : m == this->m() &&
1237 70 : n == this->n())
1238 : {
1239 70 : new_row_start = this->row_start(),
1240 70 : new_row_stop = this->row_stop();
1241 :
1242 72 : new_col_start = this->col_start(),
1243 70 : new_col_stop = this->col_stop();
1244 : }
1245 : else
1246 : {
1247 : // Determine which rows/columns each processor will be in charge of
1248 1535 : new_row_start = this->processor_id() * m / this->n_processors(),
1249 1489 : new_row_stop = (this->processor_id()+1) * m / this->n_processors();
1250 :
1251 1535 : new_col_start = this->processor_id() * n / this->n_processors(),
1252 1489 : new_col_stop = (this->processor_id()+1) * n / this->n_processors();
1253 : }
1254 :
1255 1559 : this->comm().gather(0, new_row_start, new_row_starts);
1256 1559 : this->comm().gather(0, new_row_stop, new_row_stops);
1257 1559 : this->comm().gather(0, new_col_start, new_col_starts);
1258 1559 : this->comm().gather(0, new_col_stop, new_col_stops);
1259 :
1260 : } // Done reading entry data and broadcasting matrix size
1261 :
1262 : // Calculate the matrix sparsity and initialize it second
1263 : {
1264 : // Deduce the sparsity pattern, or at least the maximum number of
1265 : // on- and off- diagonal non-zeros per row.
1266 1559 : numeric_index_type on_diagonal_nonzeros =0,
1267 1559 : off_diagonal_nonzeros =0;
1268 :
1269 1607 : if (this->processor_id() == 0)
1270 : {
1271 : // Data for the row we're working on
1272 : // Use 1-based indexing for current_row, as in the file
1273 34 : numeric_index_type current_row = 1;
1274 34 : processor_id_type current_proc = 0;
1275 732 : numeric_index_type current_on_diagonal_nonzeros = 0;
1276 732 : numeric_index_type current_off_diagonal_nonzeros = 0;
1277 :
1278 175081 : for (auto [i, j, value] : entries)
1279 : {
1280 174349 : if (i > current_row)
1281 : {
1282 897 : current_row = i;
1283 : // +1 for 1-based indexing in file
1284 23009 : while (current_row >= (new_row_stops[current_proc]+1))
1285 827 : ++current_proc;
1286 21271 : current_on_diagonal_nonzeros = 0;
1287 21271 : current_off_diagonal_nonzeros = 0;
1288 : }
1289 :
1290 : // +1 for 1-based indexing in file
1291 188557 : if (j >= (new_col_starts[current_proc]+1) &&
1292 169430 : j < (new_col_stops[current_proc]+1))
1293 : {
1294 147732 : ++current_on_diagonal_nonzeros;
1295 147732 : on_diagonal_nonzeros =
1296 : std::max(on_diagonal_nonzeros,
1297 5684 : current_on_diagonal_nonzeros);
1298 : }
1299 : else
1300 : {
1301 26617 : ++current_off_diagonal_nonzeros;
1302 26617 : off_diagonal_nonzeros =
1303 : std::max(off_diagonal_nonzeros,
1304 1858 : current_off_diagonal_nonzeros);
1305 : }
1306 : }
1307 : }
1308 :
1309 1559 : this->comm().broadcast(on_diagonal_nonzeros);
1310 1559 : this->comm().broadcast(off_diagonal_nonzeros);
1311 :
1312 1559 : this->init(m, n,
1313 : new_row_stop-new_row_start,
1314 : new_col_stop-new_col_start,
1315 : on_diagonal_nonzeros,
1316 : off_diagonal_nonzeros);
1317 : }
1318 :
1319 : // Set the matrix values last.
1320 : // Convert from 1-based to 0-based indexing
1321 1607 : if (this->processor_id() == 0)
1322 175081 : for (auto [i, j, value] : entries)
1323 174349 : this->set(i-1, j-1, value);
1324 :
1325 1559 : this->close();
1326 : #endif
1327 3022 : }
1328 :
1329 :
1330 :
1331 : template <typename T>
1332 0 : void SparseMatrix<T>::read_petsc_binary(const std::string &)
1333 : {
1334 0 : libmesh_not_implemented_msg
1335 : ("libMesh cannot read PETSc binary-format files into non-PETSc matrices");
1336 : }
1337 :
1338 :
1339 :
1340 : template <typename T>
1341 0 : void SparseMatrix<T>::read_petsc_hdf5(const std::string &)
1342 : {
1343 0 : libmesh_not_implemented_msg
1344 : ("libMesh cannot read PETSc HDF5-format files into non-PETSc matrices");
1345 : }
1346 :
1347 :
1348 :
1349 : template <typename T>
1350 75 : void SparseMatrix<T>::scale(const T scale)
1351 : {
1352 4 : libmesh_assert(this->closed());
1353 :
1354 375 : for (const auto i : make_range(this->row_start(), this->row_stop()))
1355 1500 : for (const auto j : make_range(this->col_start(), this->col_stop()))
1356 1200 : this->set(i, j, (*this)(i, j) * scale);
1357 75 : }
1358 :
1359 :
1360 :
1361 : template <typename T>
1362 290 : Real SparseMatrix<T>::l1_norm_diff(const SparseMatrix<T> & other_mat) const
1363 : {
1364 290 : auto diff_mat = this->clone();
1365 290 : diff_mat->add(-1.0, other_mat);
1366 580 : return diff_mat->l1_norm();
1367 266 : }
1368 :
1369 : //------------------------------------------------------------------
1370 : // Explicit instantiations
1371 : template class LIBMESH_EXPORT SparseMatrix<Number>;
1372 :
1373 : } // namespace libMesh
|
