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 : // 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 109657 : SparseMatrix<T>::SparseMatrix (const Parallel::Communicator & comm_in) :
90 : ParallelObject(comm_in),
91 103649 : _dof_map(nullptr),
92 103649 : _sp(nullptr),
93 103649 : _is_initialized(false),
94 109657 : _use_hash_table(false)
95 109657 : {}
96 :
97 :
98 :
99 : template <typename T>
100 49528 : void SparseMatrix<T>::attach_dof_map (const DofMap & dof_map)
101 : {
102 49528 : _dof_map = &dof_map;
103 49528 : if (!_sp)
104 24197 : _sp = dof_map.get_sparsity_pattern();
105 49528 : }
106 :
107 :
108 :
109 : template <typename T>
110 42630 : void SparseMatrix<T>::attach_sparsity_pattern (const SparsityPattern::Build & sp)
111 : {
112 42630 : _sp = &sp;
113 42630 : }
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 25654 : 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 728 : libmesh_ignore(comm);
199 :
200 25654 : if (matrix_build_type == MatrixBuildType::DIAGONAL)
201 0 : return std::make_unique<DiagonalMatrix<T>>(comm);
202 :
203 : // Build the appropriate vector
204 25654 : 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 25385 : case PETSC_SOLVERS:
215 25385 : 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 269 : case EIGEN_SOLVERS:
227 269 : 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 : #ifndef LIBMESH_HAVE_HDF5
525 0 : libmesh_ignore(filename, groupname);
526 0 : libmesh_error_msg("ERROR: need HDF5 support to handle .h5 files!!!");
527 : #else
528 : LOG_SCOPE("print_coreform_hdf5()", "SparseMatrix");
529 :
530 3 : const numeric_index_type first_dof = this->row_start(),
531 3 : end_dof = this->row_stop();
532 :
533 : std::vector<std::size_t> cols, row_offsets;
534 : std::vector<double> vals;
535 :
536 3 : if (this->processor_id() == 0)
537 2 : row_offsets.push_back(0);
538 :
539 17 : for (numeric_index_type i : make_range(first_dof, end_dof))
540 : {
541 174 : for (auto j : make_range(this->n()))
542 : {
543 146 : T c = (*this)(i,j);
544 146 : if (c != static_cast<T>(0.0))
545 : {
546 81 : cols.push_back(j);
547 81 : vals.push_back(c);
548 : }
549 : }
550 : // This is a *local* row offset; proc 0 may need to adjust later
551 14 : row_offsets.push_back(cols.size());
552 : }
553 :
554 3 : if (this->processor_id() == 0)
555 : {
556 2 : const hid_t file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC,
557 : H5P_DEFAULT, H5P_DEFAULT);
558 2 : if (file == H5I_INVALID_HID)
559 0 : libmesh_file_error(filename);
560 :
561 2 : const hid_t group =
562 2 : H5Gcreate2(file, groupname.c_str(), H5P_DEFAULT, H5P_DEFAULT,
563 : H5P_DEFAULT);
564 2 : check_open(filename, group, groupname);
565 :
566 8 : auto write_size_attribute = [&filename, &group]
567 : (const std::string & attribute_name, unsigned long long writeval)
568 : {
569 4 : const hid_t fspace = H5Screate(H5S_SCALAR);
570 4 : check_hdf5(filename, fspace, attribute_name + " fspace");
571 :
572 4 : const hid_t attr = H5Acreate2(group, attribute_name.c_str(),
573 4 : H5T_STD_I64LE, fspace,
574 : H5P_DEFAULT, H5P_DEFAULT);
575 4 : check_hdf5(filename, attr, attribute_name);
576 :
577 : // HDF5 is supposed to handle both upscaling and endianness
578 : // conversions here
579 4 : const herr_t errval = H5Awrite(attr, H5T_NATIVE_ULLONG, &writeval);
580 4 : check_hdf5(filename, errval, attribute_name + " write");
581 :
582 4 : H5Aclose(attr);
583 4 : H5Sclose(fspace);
584 : };
585 :
586 2 : write_size_attribute("num_cols", this->n());
587 2 : write_size_attribute("num_rows", this->m());
588 :
589 14 : auto write_vector = [&filename, &group]
590 : (const std::string & dataname, auto hdf5_file_type,
591 : auto hdf5_native_type, auto & datavec)
592 : {
593 6 : const hsize_t len[1] = {cast_int<hsize_t>(datavec.size())};
594 :
595 6 : const hid_t space = H5Screate_simple(1, len, nullptr);
596 12 : check_hdf5(filename, space, dataname + " space");
597 :
598 : const hid_t data =
599 6 : H5Dcreate2(group, dataname.c_str(), hdf5_file_type, space,
600 : H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
601 12 : check_hdf5(filename, data, dataname + " data");
602 :
603 6 : const hid_t errval =
604 6 : H5Dwrite(data, hdf5_native_type, H5S_ALL, H5S_ALL,
605 : H5P_DEFAULT, datavec.data());
606 6 : check_hdf5(filename, errval, dataname + " write");
607 :
608 6 : H5Dclose(data);
609 6 : H5Sclose(space);
610 : };
611 :
612 : std::vector<std::size_t> vals_sizes, first_dofs;
613 2 : this->comm().gather(0, vals.size(), vals_sizes);
614 2 : this->comm().gather(0, cast_int<std::size_t>(first_dof), first_dofs);
615 2 : first_dofs.push_back(this->m());
616 :
617 2 : this->comm().allgather(cols);
618 2 : this->comm().allgather(vals);
619 2 : this->comm().allgather(row_offsets);
620 :
621 : libmesh_assert_equal_to(vals.size(),
622 : cols.size());
623 : libmesh_assert_equal_to(vals.size(),
624 : std::accumulate(vals_sizes.begin(),
625 : vals_sizes.end(),
626 : std::size_t(0)));
627 :
628 : std::size_t extra_offset = 0;
629 3 : for (auto p : make_range(processor_id_type(1), this->n_processors()))
630 : {
631 1 : extra_offset += vals_sizes[p-1];
632 6 : for (auto i : make_range(first_dofs[p]+1, first_dofs[p+1]+1))
633 5 : row_offsets[i] += extra_offset;
634 : }
635 :
636 2 : write_vector("cols", H5T_STD_U64LE, H5T_NATIVE_ULLONG, cols);
637 2 : write_vector("row_offsets", H5T_STD_U64LE, H5T_NATIVE_ULLONG, row_offsets);
638 2 : write_vector("vals", H5T_IEEE_F64LE, H5T_NATIVE_DOUBLE, vals);
639 :
640 2 : H5Gclose(group);
641 2 : H5Fclose(file);
642 : }
643 : else
644 : {
645 : std::vector<std::size_t> dummy;
646 1 : this->comm().gather(0, vals.size(), dummy);
647 1 : this->comm().gather(0, cast_int<std::size_t>(first_dof), dummy);
648 1 : this->comm().allgather(cols);
649 1 : this->comm().allgather(vals);
650 1 : this->comm().allgather(row_offsets);
651 : }
652 : #endif // LIBMESH_HAVE_HDF5
653 3 : }
654 :
655 :
656 :
657 : template <typename T>
658 0 : void SparseMatrix<T>::print_petsc_binary(const std::string &) const
659 : {
660 0 : libmesh_not_implemented_msg
661 : ("libMesh cannot write PETSc binary-format files from non-PETSc matrices");
662 : }
663 :
664 :
665 :
666 : template <typename T>
667 0 : void SparseMatrix<T>::print_petsc_hdf5(const std::string &) const
668 : {
669 0 : libmesh_not_implemented_msg
670 : ("libMesh cannot write PETSc HDF5-format files from non-PETSc matrices");
671 : }
672 :
673 :
674 :
675 : template <typename T>
676 1216 : void SparseMatrix<T>::read(const std::string & filename)
677 : {
678 : {
679 1292 : std::ifstream in (filename.c_str());
680 1216 : libmesh_error_msg_if
681 : (!in.good(), "ERROR: cannot read file:\n\t" <<
682 : filename);
683 1140 : }
684 :
685 1216 : std::string_view basename = Utility::basename_of(filename);
686 :
687 1216 : const bool gzipped_file = Utility::ends_with(filename, ".gz");
688 :
689 1216 : if (gzipped_file)
690 18 : basename.remove_suffix(3);
691 :
692 2432 : if (Utility::ends_with(basename, ".matlab") ||
693 1330 : Utility::ends_with(basename, ".m"))
694 1138 : this->read_matlab(filename);
695 78 : else if (Utility::ends_with(basename, ".petsc64"))
696 : {
697 : #ifndef LIBMESH_HAVE_PETSC
698 0 : libmesh_error_msg("Cannot load PETSc matrix file " <<
699 : filename << " without PETSc-enabled libMesh.");
700 : #elif LIBMESH_DOF_ID_BYTES != 8
701 0 : libmesh_error_msg("Cannot load 64-bit PETSc matrix file " <<
702 : filename << " with non-64-bit libMesh.");
703 : #endif
704 66 : if (gzipped_file)
705 0 : libmesh_not_implemented_msg("Gzipped PETSc matrices are not currently supported");
706 66 : this->read_petsc_binary(filename);
707 : }
708 12 : else if (Utility::ends_with(basename, ".petsc32"))
709 : {
710 : #ifndef LIBMESH_HAVE_PETSC
711 0 : libmesh_error_msg("Cannot load PETSc matrix file " <<
712 : filename << " without PETSc-enabled libMesh.");
713 : #elif LIBMESH_DOF_ID_BYTES != 4
714 0 : libmesh_error_msg("Cannot load 32-bit PETSc matrix file " <<
715 : filename << " with non-32-bit libMesh.");
716 : #endif
717 4 : if (gzipped_file)
718 0 : libmesh_not_implemented_msg("Gzipped PETSc matrices are not currently supported");
719 4 : this->read_petsc_binary(filename);
720 : }
721 8 : else if (Utility::ends_with(basename, ".h5"))
722 : {
723 8 : if (gzipped_file)
724 0 : libmesh_not_implemented_msg("Gzipped HDF5 matrices are not currently supported");
725 16 : this->read_coreform_hdf5(filename);
726 : }
727 : else
728 0 : libmesh_error_msg(" ERROR: Unrecognized matrix file extension on: "
729 : << basename
730 : << "\n I understand the following:\n\n"
731 : << " *.h5 -- CoreForm HDF5 sparse matrix format\n"
732 : << " *.matlab -- Matlab sparse matrix format\n"
733 : << " *.matlab.gz -- Matlab sparse matrix format, gzipped\n"
734 : << " *.m -- Matlab sparse matrix format\n"
735 : << " *.m.gz -- Matlab sparse matrix format, gzipped\n"
736 : << " *.petsc32 -- PETSc binary format, 32-bit\n"
737 : << " *.petsc64 -- PETSc binary format, 64-bit\n"
738 : );
739 1216 : }
740 :
741 :
742 :
743 : template <typename T>
744 171 : void SparseMatrix<T>::print(const std::string & filename) const
745 : {
746 : {
747 187 : std::ofstream outstr (filename.c_str());
748 171 : libmesh_error_msg_if
749 : (!outstr.good(), "ERROR: cannot write to file:\n\t" <<
750 : filename);
751 155 : }
752 :
753 171 : std::string_view basename = Utility::basename_of(filename);
754 :
755 171 : const bool gzipped_file = Utility::ends_with(filename, ".gz");
756 :
757 171 : if (gzipped_file)
758 4 : basename.remove_suffix(3);
759 :
760 342 : if (Utility::ends_with(basename, ".matlab") ||
761 182 : Utility::ends_with(basename, ".m"))
762 168 : this->print_matlab(filename);
763 3 : else if (Utility::ends_with(basename, ".petsc64"))
764 : {
765 : #ifndef LIBMESH_HAVE_PETSC
766 0 : libmesh_error_msg("Cannot load PETSc matrix file " <<
767 : filename << " without PETSc-enabled libMesh.");
768 : #elif LIBMESH_DOF_ID_BYTES != 8
769 0 : libmesh_error_msg("Cannot load 64-bit PETSc matrix file " <<
770 : filename << " with non-64-bit libMesh.");
771 : #endif
772 0 : if (gzipped_file)
773 0 : libmesh_not_implemented_msg("Gzipped PETSc matrices are not currently supported");
774 0 : this->print_petsc_binary(filename);
775 : }
776 3 : else if (Utility::ends_with(basename, ".petsc32"))
777 : {
778 : #ifndef LIBMESH_HAVE_PETSC
779 0 : libmesh_error_msg("Cannot load PETSc matrix file " <<
780 : filename << " without PETSc-enabled libMesh.");
781 : #elif LIBMESH_DOF_ID_BYTES != 4
782 0 : libmesh_error_msg("Cannot load 32-bit PETSc matrix file " <<
783 : filename << " with non-32-bit libMesh.");
784 : #endif
785 0 : if (gzipped_file)
786 0 : libmesh_not_implemented_msg("Gzipped PETSc matrices are not currently supported");
787 0 : this->print_petsc_binary(filename);
788 : }
789 3 : else if (Utility::ends_with(basename, ".h5"))
790 : {
791 3 : if (gzipped_file)
792 0 : libmesh_not_implemented_msg("Gzipped HDF5 matrices are not currently supported");
793 6 : this->print_coreform_hdf5(filename);
794 : }
795 : else
796 0 : libmesh_error_msg(" ERROR: Unrecognized matrix file extension on: "
797 : << basename
798 : << "\n I understand the following:\n\n"
799 : << " *.h5 -- CoreForm HDF5 sparse matrix format\n"
800 : << " *.matlab -- Matlab sparse matrix format\n"
801 : << " *.matlab.gz -- Matlab sparse matrix format, gzipped\n"
802 : << " *.m -- Matlab sparse matrix format\n"
803 : << " *.m.gz -- Matlab sparse matrix format, gzipped\n"
804 : << " *.petsc32 -- PETSc binary format, 32-bit\n"
805 : << " *.petsc64 -- PETSc binary format, 64-bit\n"
806 : );
807 171 : }
808 :
809 :
810 :
811 : template <typename T>
812 8 : void SparseMatrix<T>::read_coreform_hdf5(const std::string & filename,
813 : const std::string & groupname)
814 : {
815 : #ifndef LIBMESH_HAVE_HDF5
816 0 : libmesh_ignore(filename, groupname);
817 0 : libmesh_error_msg("ERROR: need HDF5 support to handle .h5 files!!!");
818 : #else
819 : LOG_SCOPE("read_coreform_hdf5()", "SparseMatrix");
820 :
821 : std::size_t num_rows, num_cols;
822 :
823 : // These are only used on pid 0, but avoid "uninitialized" warnings
824 8 : hid_t group = H5I_INVALID_HID;
825 : hid_t file = H5I_INVALID_HID;
826 :
827 8 : if (this->processor_id() == 0)
828 : {
829 6 : file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
830 :
831 6 : if (file == H5I_INVALID_HID)
832 0 : libmesh_file_error(filename);
833 :
834 6 : group = H5Gopen(file, groupname.c_str(), H5P_DEFAULT);
835 6 : check_open(filename, group, groupname);
836 :
837 54 : auto read_size_attribute = [&filename, &group]
838 : (const std::string & attribute_name)
839 : {
840 12 : unsigned long long returnval = 0;
841 :
842 12 : const hid_t attr = H5Aopen(group, attribute_name.c_str(), H5P_DEFAULT);
843 12 : check_open(filename, attr, attribute_name);
844 :
845 12 : const hid_t attr_type = H5Aget_type(attr);
846 12 : check_hdf5(filename, attr_type, attribute_name + " type");
847 :
848 : // HDF5 will convert between the file's integer type and ours, but
849 : // we do expect an integer type.
850 12 : if (H5Tget_class(attr_type) != H5T_INTEGER)
851 0 : libmesh_error_msg("Non-integer type for " + attribute_name + " in " + filename);
852 :
853 12 : H5Tclose(attr_type);
854 :
855 : // HDF5 is supposed to handle both upscaling and endianness
856 : // conversions here
857 12 : const herr_t errval = H5Aread(attr, H5T_NATIVE_ULLONG, &returnval);
858 12 : check_hdf5(filename, errval, attribute_name + " read");
859 :
860 12 : H5Aclose(attr);
861 :
862 12 : return returnval;
863 : };
864 :
865 6 : num_cols = read_size_attribute("num_cols");
866 6 : num_rows = read_size_attribute("num_rows");
867 :
868 6 : this->comm().broadcast(num_cols);
869 6 : this->comm().broadcast(num_rows);
870 : }
871 : else
872 : {
873 2 : this->comm().broadcast(num_cols);
874 2 : this->comm().broadcast(num_rows);
875 : }
876 :
877 : numeric_index_type new_row_start, new_row_stop,
878 : new_col_start, new_col_stop;
879 :
880 : // If we need to reinit, we need to determine which rows+columns
881 : // each processor is in charge of.
882 : std::vector<numeric_index_type> new_row_starts, new_row_stops,
883 : new_col_starts, new_col_stops;
884 :
885 8 : if (this->initialized() &&
886 8 : num_cols == this->n() &&
887 0 : num_rows == this->m())
888 : {
889 0 : new_row_start = this->row_start(),
890 0 : new_row_stop = this->row_stop();
891 :
892 0 : new_col_start = this->col_start(),
893 0 : new_col_stop = this->col_stop();
894 : }
895 : else
896 : {
897 : // Determine which rows/columns each processor will be in charge of
898 8 : new_row_start = this->processor_id() * num_rows / this->n_processors(),
899 8 : new_row_stop = (this->processor_id()+1) * num_rows / this->n_processors();
900 :
901 8 : new_col_start = this->processor_id() * num_cols / this->n_processors(),
902 8 : new_col_stop = (this->processor_id()+1) * num_cols / this->n_processors();
903 : }
904 :
905 8 : this->comm().gather(0, new_row_start, new_row_starts);
906 8 : this->comm().gather(0, new_row_stop, new_row_stops);
907 8 : this->comm().gather(0, new_col_start, new_col_starts);
908 8 : this->comm().gather(0, new_col_stop, new_col_stops);
909 :
910 8 : numeric_index_type on_diagonal_nonzeros = 0,
911 8 : off_diagonal_nonzeros = 0;
912 :
913 : std::vector<std::size_t> cols, row_offsets;
914 : std::vector<double> vals;
915 :
916 8 : if (this->processor_id() == 0)
917 : {
918 42 : auto read_vector = [&filename, &group]
919 : (const std::string & dataname, auto hdf5_class,
920 : auto hdf5_type, auto & datavec)
921 : {
922 18 : const hid_t data = H5Dopen1(group, dataname.c_str());
923 18 : check_open(filename, data, dataname.c_str());
924 :
925 18 : const hid_t data_type = H5Dget_type(data);
926 18 : check_hdf5(filename, data_type, dataname + " type");
927 :
928 : // HDF5 will convert between the file's integer type and ours, but
929 : // we do expect an integer type.
930 18 : if (H5Tget_class(data_type) != hdf5_class)
931 0 : libmesh_error_msg("Unexpected type for " + dataname + " in " + filename);
932 :
933 18 : H5Tclose(data_type);
934 :
935 18 : const hid_t dataspace = H5Dget_space(data);
936 18 : check_hdf5(filename, dataspace, dataname + " space");
937 :
938 18 : int ndims = H5Sget_simple_extent_ndims(dataspace);
939 18 : if (ndims != 1)
940 0 : libmesh_error_msg("Non-vector space for " + dataname + " in " + filename);
941 :
942 : hsize_t len, maxlen;
943 18 : herr_t errval = H5Sget_simple_extent_dims(dataspace, &len, &maxlen);
944 18 : check_hdf5(filename, errval, dataname + " dims");
945 :
946 18 : datavec.resize(len);
947 :
948 18 : errval = H5Dread(data, hdf5_type, H5S_ALL, H5S_ALL, H5P_DEFAULT, datavec.data());
949 18 : check_hdf5(filename, errval, dataname + " read");
950 :
951 18 : H5Dclose(data);
952 : };
953 :
954 6 : read_vector("cols", H5T_INTEGER, H5T_NATIVE_ULLONG, cols);
955 6 : read_vector("row_offsets", H5T_INTEGER, H5T_NATIVE_ULLONG, row_offsets);
956 12 : read_vector("vals", H5T_FLOAT, H5T_NATIVE_DOUBLE, vals);
957 :
958 6 : if (cols.size() != vals.size())
959 0 : libmesh_error_msg("Inconsistent cols/vals sizes in " + filename);
960 :
961 6 : if (row_offsets.size() != num_rows + 1)
962 0 : libmesh_error_msg("Inconsistent row_offsets size in " + filename);
963 :
964 : // Data for the row we're working on
965 : numeric_index_type current_row = 0;
966 : processor_id_type current_proc = 0;
967 6 : numeric_index_type current_on_diagonal_nonzeros = 0;
968 6 : numeric_index_type current_off_diagonal_nonzeros = 0;
969 6 : if (row_offsets[0] != 0)
970 0 : libmesh_error_msg("Unexpected row_offsets[0] in " + filename);
971 :
972 481 : for (auto i : index_range(cols))
973 : {
974 568 : while (row_offsets[current_row+1] <= i)
975 : {
976 : ++current_row;
977 93 : if (row_offsets[current_row] < row_offsets[current_row-1])
978 0 : libmesh_error_msg("Non-monotonic row_offsets in " + filename);
979 93 : current_on_diagonal_nonzeros = 0;
980 93 : current_off_diagonal_nonzeros = 0;
981 : }
982 :
983 477 : while (current_row >= new_row_stops[current_proc])
984 2 : ++current_proc;
985 :
986 : // 0-based indexing in file
987 475 : if (cols[i] >= new_col_starts[current_proc] &&
988 449 : cols[i] < new_col_stops[current_proc])
989 : {
990 414 : ++current_on_diagonal_nonzeros;
991 414 : on_diagonal_nonzeros =
992 : std::max(on_diagonal_nonzeros,
993 : current_on_diagonal_nonzeros);
994 : }
995 : else
996 : {
997 61 : ++current_off_diagonal_nonzeros;
998 61 : off_diagonal_nonzeros =
999 : std::max(off_diagonal_nonzeros,
1000 : current_off_diagonal_nonzeros);
1001 : }
1002 : }
1003 : }
1004 :
1005 8 : this->comm().broadcast(on_diagonal_nonzeros);
1006 8 : this->comm().broadcast(off_diagonal_nonzeros);
1007 :
1008 8 : this->init(num_rows, num_cols,
1009 : new_row_stop-new_row_start,
1010 : new_col_stop-new_col_start,
1011 : on_diagonal_nonzeros,
1012 : off_diagonal_nonzeros);
1013 :
1014 : // Set the matrix values last.
1015 8 : if (this->processor_id() == 0)
1016 : {
1017 : numeric_index_type current_row = 0;
1018 481 : for (auto i : index_range(cols))
1019 : {
1020 568 : while (row_offsets[current_row+1] <= i)
1021 : {
1022 : ++current_row;
1023 : libmesh_assert_greater_equal (row_offsets[current_row],
1024 : row_offsets[current_row-1]);
1025 : }
1026 475 : this->set(current_row, cols[i], vals[i]);
1027 : }
1028 :
1029 6 : H5Gclose(group);
1030 6 : H5Fclose(file);
1031 : }
1032 :
1033 8 : this->close();
1034 : #endif // LIBMESH_HAVE_HDF5
1035 8 : }
1036 :
1037 :
1038 :
1039 : template <typename T>
1040 1559 : void SparseMatrix<T>::read_matlab(const std::string & filename)
1041 : {
1042 96 : LOG_SCOPE("read_matlab()", "SparseMatrix");
1043 :
1044 : #ifndef LIBMESH_HAVE_CXX11_REGEX
1045 : libmesh_not_implemented(); // What is your compiler?!? Email us!
1046 : libmesh_ignore(filename);
1047 : #else
1048 48 : parallel_object_only();
1049 :
1050 1559 : const bool gzipped_file = Utility::ends_with(filename, ".gz");
1051 :
1052 : // The sizes we get from the file
1053 1559 : std::size_t m = 0,
1054 1559 : n = 0;
1055 :
1056 : // If we don't already have this size, we'll need to reinit, and
1057 : // determine which rows+columns each processor is in charge of.
1058 96 : std::vector<numeric_index_type> new_row_starts, new_row_stops,
1059 96 : new_col_starts, new_col_stops;
1060 :
1061 : numeric_index_type new_row_start, new_row_stop,
1062 : new_col_start, new_col_stop;
1063 :
1064 : // We'll read through the file three times: once to get a reliable
1065 : // value for the matrix size (so we can divvy it up among
1066 : // processors), then again to get the sparsity to send to each
1067 : // processor, then a final time to get the entries to send to each
1068 : // processor.
1069 : //
1070 : // We'll use an istream here; it might be an ifstream if we're
1071 : // opening a raw ASCII file or a gzstream if we're opening a
1072 : // compressed one.
1073 1511 : std::unique_ptr<std::istream> file;
1074 :
1075 : // We'll need a temporary structure to cache matrix entries, because
1076 : // we need to read through the whole file before we know the size
1077 : // and sparsity structure with which we can init().
1078 : //
1079 : // Reading through the file three times via `seekg` doesn't work
1080 : // with our gzstream wrapper, and seems to take three times as long
1081 : // even with a plain ifstream. What happened to disk caching!?
1082 96 : std::vector<std::tuple<numeric_index_type, numeric_index_type, T>> entries;
1083 :
1084 : // First read through the file, saving size and entry data
1085 : {
1086 : // We'll read the matrix on processor 0 rather than try to juggle
1087 : // parallel I/O.
1088 1607 : if (this->processor_id() == 0)
1089 : {
1090 : // We'll be using regular expressions to make ourselves slightly
1091 : // more robust to formatting.
1092 800 : const std::regex start_regex // assignment like "zzz = ["
1093 : ("\\s*\\w+\\s*=\\s*\\[");
1094 800 : const std::regex end_regex // end of assignment
1095 : ("^[^%]*\\]");
1096 :
1097 732 : if (gzipped_file)
1098 : {
1099 : #ifdef LIBMESH_HAVE_GZSTREAM
1100 345 : auto inf = std::make_unique<igzstream>();
1101 14 : libmesh_assert(inf);
1102 14 : inf->open(filename.c_str(), std::ios::in);
1103 317 : file = std::move(inf);
1104 : #else
1105 : libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
1106 : #endif
1107 303 : }
1108 : else
1109 : {
1110 421 : auto inf = std::make_unique<std::ifstream>();
1111 20 : libmesh_assert(inf);
1112 :
1113 421 : std::string new_name = Utility::unzip_file(filename);
1114 :
1115 401 : inf->open(new_name.c_str(), std::ios::in);
1116 381 : file = std::move(inf);
1117 361 : }
1118 :
1119 : // If we have a matrix with all-zero trailing rows, the only
1120 : // way to get the size is if it ended up in a comment
1121 800 : const std::regex size_regex // comment like "% size = 8 8"
1122 : ("%\\s*[Ss][Ii][Zz][Ee]\\s*=\\s*(\\d+)\\s+(\\d+)");
1123 766 : const std::string whitespace = " \t";
1124 :
1125 34 : bool have_started = false;
1126 34 : bool have_ended = false;
1127 732 : std::size_t largest_i_seen = 0, largest_j_seen = 0;
1128 :
1129 : // Data for the row we're working on
1130 : // Use 1-based indexing for current_row, as in the file
1131 732 : std::size_t current_row = 1;
1132 :
1133 343420 : for (std::string line; std::getline(*file, line);)
1134 : {
1135 7780 : std::smatch sm;
1136 :
1137 : // First, try to match an entry. This is the most common
1138 : // case so we won't rely on slow std::regex for it.
1139 : // stringstream is at least an improvement over that.
1140 :
1141 : // Look for row/col/val like "1 1 -2.0e-4"
1142 :
1143 187219 : std::istringstream l(line);
1144 :
1145 : std::size_t i, j;
1146 1822 : T value;
1147 :
1148 9602 : l >> i >> j >> value;
1149 :
1150 179473 : if (!l.fail())
1151 : {
1152 174349 : libmesh_error_msg_if
1153 : (!have_started, "Confused by premature entries in matrix file " << filename);
1154 :
1155 331841 : entries.emplace_back(cast_int<numeric_index_type>(i),
1156 181891 : cast_int<numeric_index_type>(j),
1157 : value);
1158 :
1159 174349 : libmesh_error_msg_if
1160 : (!i || !j, "Expected 1-based indexing in matrix file "
1161 : << filename);
1162 :
1163 174349 : current_row = std::max(current_row, i);
1164 :
1165 174349 : libmesh_error_msg_if
1166 : (i < current_row,
1167 : "Can't handle out-of-order entries in matrix file "
1168 : << filename);
1169 :
1170 174349 : largest_i_seen = std::max(i, largest_i_seen);
1171 174349 : largest_j_seen = std::max(j, largest_j_seen);
1172 : }
1173 :
1174 5124 : else if (std::regex_search(line, sm, size_regex))
1175 : {
1176 766 : const std::string msize = sm[1];
1177 732 : const std::string nsize = sm[2];
1178 732 : m = std::stoull(msize);
1179 732 : n = std::stoull(nsize);
1180 : }
1181 :
1182 4392 : else if (std::regex_search(line, start_regex))
1183 34 : have_started = true;
1184 :
1185 3660 : else if (std::regex_search(line, end_regex))
1186 : {
1187 34 : have_ended = true;
1188 68 : break;
1189 : }
1190 : }
1191 :
1192 732 : libmesh_error_msg_if
1193 : (!have_started, "Confused by missing assignment beginning in matrix file " << filename);
1194 :
1195 732 : libmesh_error_msg_if
1196 : (!have_ended, "Confused by missing assignment ending in matrix file " << filename);
1197 :
1198 732 : libmesh_error_msg_if
1199 : (m > largest_i_seen, "Confused by missing final row(s) in matrix file " << filename);
1200 :
1201 732 : libmesh_error_msg_if
1202 : (m > 0 && m < largest_i_seen, "Confused by extra final row(s) in matrix file " << filename);
1203 :
1204 732 : if (!m)
1205 0 : m = largest_i_seen;
1206 :
1207 732 : libmesh_error_msg_if
1208 : (n > largest_j_seen, "Confused by missing final column(s) in matrix file " << filename);
1209 :
1210 732 : libmesh_error_msg_if
1211 : (n > 0 && n < largest_j_seen, "Confused by extra final column(s) in matrix file " << filename);
1212 :
1213 732 : if (!n)
1214 0 : n = largest_j_seen;
1215 :
1216 732 : this->comm().broadcast(m);
1217 732 : this->comm().broadcast(n);
1218 664 : }
1219 : else
1220 : {
1221 813 : this->comm().broadcast(m);
1222 827 : this->comm().broadcast(n);
1223 : }
1224 :
1225 1559 : if (this->initialized() &&
1226 1561 : m == this->m() &&
1227 70 : n == this->n())
1228 : {
1229 70 : new_row_start = this->row_start(),
1230 70 : new_row_stop = this->row_stop();
1231 :
1232 72 : new_col_start = this->col_start(),
1233 70 : new_col_stop = this->col_stop();
1234 : }
1235 : else
1236 : {
1237 : // Determine which rows/columns each processor will be in charge of
1238 1535 : new_row_start = this->processor_id() * m / this->n_processors(),
1239 1489 : new_row_stop = (this->processor_id()+1) * m / this->n_processors();
1240 :
1241 1535 : new_col_start = this->processor_id() * n / this->n_processors(),
1242 1489 : new_col_stop = (this->processor_id()+1) * n / this->n_processors();
1243 : }
1244 :
1245 1559 : this->comm().gather(0, new_row_start, new_row_starts);
1246 1559 : this->comm().gather(0, new_row_stop, new_row_stops);
1247 1559 : this->comm().gather(0, new_col_start, new_col_starts);
1248 1559 : this->comm().gather(0, new_col_stop, new_col_stops);
1249 :
1250 : } // Done reading entry data and broadcasting matrix size
1251 :
1252 : // Calculate the matrix sparsity and initialize it second
1253 : {
1254 : // Deduce the sparsity pattern, or at least the maximum number of
1255 : // on- and off- diagonal non-zeros per row.
1256 1559 : numeric_index_type on_diagonal_nonzeros =0,
1257 1559 : off_diagonal_nonzeros =0;
1258 :
1259 1607 : if (this->processor_id() == 0)
1260 : {
1261 : // Data for the row we're working on
1262 : // Use 1-based indexing for current_row, as in the file
1263 34 : numeric_index_type current_row = 1;
1264 34 : processor_id_type current_proc = 0;
1265 732 : numeric_index_type current_on_diagonal_nonzeros = 0;
1266 732 : numeric_index_type current_off_diagonal_nonzeros = 0;
1267 :
1268 175081 : for (auto [i, j, value] : entries)
1269 : {
1270 174349 : if (i > current_row)
1271 : {
1272 897 : current_row = i;
1273 : // +1 for 1-based indexing in file
1274 23009 : while (current_row >= (new_row_stops[current_proc]+1))
1275 827 : ++current_proc;
1276 21271 : current_on_diagonal_nonzeros = 0;
1277 21271 : current_off_diagonal_nonzeros = 0;
1278 : }
1279 :
1280 : // +1 for 1-based indexing in file
1281 188557 : if (j >= (new_col_starts[current_proc]+1) &&
1282 169430 : j < (new_col_stops[current_proc]+1))
1283 : {
1284 147732 : ++current_on_diagonal_nonzeros;
1285 147732 : on_diagonal_nonzeros =
1286 : std::max(on_diagonal_nonzeros,
1287 5684 : current_on_diagonal_nonzeros);
1288 : }
1289 : else
1290 : {
1291 26617 : ++current_off_diagonal_nonzeros;
1292 26617 : off_diagonal_nonzeros =
1293 : std::max(off_diagonal_nonzeros,
1294 1858 : current_off_diagonal_nonzeros);
1295 : }
1296 : }
1297 : }
1298 :
1299 1559 : this->comm().broadcast(on_diagonal_nonzeros);
1300 1559 : this->comm().broadcast(off_diagonal_nonzeros);
1301 :
1302 1559 : this->init(m, n,
1303 : new_row_stop-new_row_start,
1304 : new_col_stop-new_col_start,
1305 : on_diagonal_nonzeros,
1306 : off_diagonal_nonzeros);
1307 : }
1308 :
1309 : // Set the matrix values last.
1310 : // Convert from 1-based to 0-based indexing
1311 1607 : if (this->processor_id() == 0)
1312 175081 : for (auto [i, j, value] : entries)
1313 174349 : this->set(i-1, j-1, value);
1314 :
1315 1559 : this->close();
1316 : #endif
1317 3022 : }
1318 :
1319 :
1320 :
1321 : template <typename T>
1322 0 : void SparseMatrix<T>::read_petsc_binary(const std::string &)
1323 : {
1324 0 : libmesh_not_implemented_msg
1325 : ("libMesh cannot read PETSc binary-format files into non-PETSc matrices");
1326 : }
1327 :
1328 :
1329 :
1330 : template <typename T>
1331 0 : void SparseMatrix<T>::read_petsc_hdf5(const std::string &)
1332 : {
1333 0 : libmesh_not_implemented_msg
1334 : ("libMesh cannot read PETSc HDF5-format files into non-PETSc matrices");
1335 : }
1336 :
1337 :
1338 :
1339 : template <typename T>
1340 75 : void SparseMatrix<T>::scale(const T scale)
1341 : {
1342 4 : libmesh_assert(this->closed());
1343 :
1344 375 : for (const auto i : make_range(this->row_start(), this->row_stop()))
1345 1500 : for (const auto j : make_range(this->col_start(), this->col_stop()))
1346 1200 : this->set(i, j, (*this)(i, j) * scale);
1347 75 : }
1348 :
1349 :
1350 :
1351 : template <typename T>
1352 290 : Real SparseMatrix<T>::l1_norm_diff(const SparseMatrix<T> & other_mat) const
1353 : {
1354 290 : auto diff_mat = this->clone();
1355 290 : diff_mat->add(-1.0, other_mat);
1356 580 : return diff_mat->l1_norm();
1357 266 : }
1358 :
1359 : //------------------------------------------------------------------
1360 : // Explicit instantiations
1361 : template class LIBMESH_EXPORT SparseMatrix<Number>;
1362 :
1363 : } // namespace libMesh
|
