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 109649 : SparseMatrix<T>::SparseMatrix (const Parallel::Communicator & comm_in) :
67 : ParallelObject(comm_in),
68 103641 : _dof_map(nullptr),
69 103641 : _sp(nullptr),
70 103641 : _is_initialized(false),
71 109649 : _use_hash_table(false)
72 109649 : {}
73 :
74 :
75 :
76 : template <typename T>
77 49528 : void SparseMatrix<T>::attach_dof_map (const DofMap & dof_map)
78 : {
79 49528 : _dof_map = &dof_map;
80 49528 : if (!_sp)
81 24197 : _sp = dof_map.get_sparsity_pattern();
82 49528 : }
83 :
84 :
85 :
86 : template <typename T>
87 42630 : void SparseMatrix<T>::attach_sparsity_pattern (const SparsityPattern::Build & sp)
88 : {
89 42630 : _sp = &sp;
90 42630 : }
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 25654 : 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 728 : libmesh_ignore(comm);
176 :
177 25654 : if (matrix_build_type == MatrixBuildType::DIAGONAL)
178 0 : return std::make_unique<DiagonalMatrix<T>>(comm);
179 :
180 : // Build the appropriate vector
181 25654 : 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 25385 : case PETSC_SOLVERS:
192 25385 : 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 269 : case EIGEN_SOLVERS:
204 269 : 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 : #elif LIBMESH_DOF_ID_BYTES != 8
541 0 : libmesh_error_msg("Cannot load 64-bit PETSc matrix file " <<
542 : filename << " with non-64-bit libMesh.");
543 : #endif
544 66 : this->read_petsc_binary(filename);
545 : }
546 8 : else if (Utility::ends_with(basename, ".petsc32"))
547 : {
548 : #ifndef LIBMESH_HAVE_PETSC
549 0 : libmesh_error_msg("Cannot load PETSc matrix file " <<
550 : filename << " without PETSc-enabled libMesh.");
551 : #elif LIBMESH_DOF_ID_BYTES != 4
552 0 : libmesh_error_msg("Cannot load 32-bit PETSc matrix file " <<
553 : filename << " with non-32-bit libMesh.");
554 : #endif
555 4 : this->read_petsc_binary(filename);
556 : }
557 4 : else if (Utility::ends_with(basename, ".h5"))
558 : {
559 8 : this->read_coreform_hdf5(filename);
560 : }
561 : else
562 0 : libmesh_error_msg(" ERROR: Unrecognized matrix file extension on: "
563 : << basename
564 : << "\n I understand the following:\n\n"
565 : << " *.h5 -- CoreForm HDF5 sparse matrix format\n"
566 : << " *.matlab -- Matlab sparse matrix format\n"
567 : << " *.m -- Matlab sparse matrix format\n"
568 : << " *.petsc32 -- PETSc binary format, 32-bit\n"
569 : << " *.petsc64 -- PETSc binary format, 64-bit\n"
570 : );
571 1068 : }
572 :
573 :
574 :
575 : template <typename T>
576 4 : void SparseMatrix<T>::read_coreform_hdf5(const std::string & filename,
577 : const std::string & groupname)
578 : {
579 : #ifndef LIBMESH_HAVE_HDF5
580 0 : libmesh_ignore(filename, groupname);
581 0 : libmesh_error_msg("ERROR: need HDF5 support to handle .h5 files!!!");
582 : #else
583 : LOG_SCOPE("read_coreform_hdf5()", "SparseMatrix");
584 :
585 : std::size_t num_rows, num_cols;
586 :
587 : hid_t group;
588 :
589 18 : auto check_open = [&filename](hid_t id, const std::string & objname)
590 : {
591 18 : if (id == H5I_INVALID_HID)
592 0 : libmesh_error_msg("Couldn't open " + objname + " in " + filename);
593 : };
594 :
595 52 : auto check_hdf5 = [&filename](auto hdf5val, const std::string & objname)
596 : {
597 48 : if (hdf5val < 0)
598 0 : libmesh_error_msg("HDF5 error from " + objname + " in " + filename);
599 : };
600 :
601 :
602 4 : if (this->processor_id() == 0)
603 : {
604 3 : const hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
605 :
606 3 : if (file == H5I_INVALID_HID)
607 0 : libmesh_file_error(filename);
608 :
609 3 : group = H5Gopen(file, groupname.c_str(), H5P_DEFAULT);
610 3 : check_open(group, groupname);
611 :
612 27 : auto read_size_attribute = [&filename, &group, &check_open, &check_hdf5]
613 : (const std::string & attribute_name)
614 : {
615 6 : unsigned long long returnval = 0;
616 :
617 6 : const hid_t attr = H5Aopen(group, attribute_name.c_str(), H5P_DEFAULT);
618 6 : check_open(attr, attribute_name);
619 :
620 6 : const hid_t attr_type = H5Aget_type(attr);
621 6 : check_hdf5(attr_type, attribute_name + " type");
622 :
623 : // HDF5 will convert between the file's integer type and ours, but
624 : // we do expect an integer type.
625 6 : if (H5Tget_class(attr_type) != H5T_INTEGER)
626 0 : libmesh_error_msg("Non-integer type for " + attribute_name + " in " + filename);
627 :
628 6 : H5Tclose(attr_type);
629 :
630 : // HDF5 is supposed to handle both upscaling and endianness
631 : // conversions here
632 6 : herr_t errval = H5Aread(attr, H5T_NATIVE_ULLONG, &returnval);
633 6 : check_hdf5(errval, attribute_name + " read");
634 :
635 6 : H5Aclose(attr);
636 :
637 6 : return returnval;
638 : };
639 :
640 3 : num_cols = read_size_attribute("num_cols");
641 3 : num_rows = read_size_attribute("num_rows");
642 :
643 3 : this->comm().broadcast(num_cols);
644 3 : this->comm().broadcast(num_rows);
645 : }
646 : else
647 : {
648 1 : this->comm().broadcast(num_cols);
649 1 : this->comm().broadcast(num_rows);
650 : }
651 :
652 : numeric_index_type new_row_start, new_row_stop,
653 : new_col_start, new_col_stop;
654 :
655 : // If we need to reinit, we need to determine which rows+columns
656 : // each processor is in charge of.
657 : std::vector<numeric_index_type> new_row_starts, new_row_stops,
658 : new_col_starts, new_col_stops;
659 :
660 4 : if (this->initialized() &&
661 4 : num_cols == this->m() &&
662 0 : num_rows == this->n())
663 : {
664 0 : new_row_start = this->row_start(),
665 0 : new_row_stop = this->row_stop();
666 :
667 0 : new_col_start = this->col_start(),
668 0 : new_col_stop = this->col_stop();
669 : }
670 : else
671 : {
672 : // Determine which rows/columns each processor will be in charge of
673 4 : new_row_start = this->processor_id() * num_rows / this->n_processors(),
674 4 : new_row_stop = (this->processor_id()+1) * num_rows / this->n_processors();
675 :
676 4 : new_col_start = this->processor_id() * num_cols / this->n_processors(),
677 4 : new_col_stop = (this->processor_id()+1) * num_cols / this->n_processors();
678 : }
679 :
680 4 : this->comm().gather(0, new_row_start, new_row_starts);
681 4 : this->comm().gather(0, new_row_stop, new_row_stops);
682 4 : this->comm().gather(0, new_col_start, new_col_starts);
683 4 : this->comm().gather(0, new_col_stop, new_col_stops);
684 :
685 4 : numeric_index_type on_diagonal_nonzeros = 0,
686 4 : off_diagonal_nonzeros = 0;
687 :
688 : std::vector<std::size_t> cols, row_offsets;
689 : std::vector<double> vals;
690 :
691 4 : if (this->processor_id() == 0)
692 : {
693 12 : auto read_vector = [&filename, &group, &check_open, &check_hdf5]
694 : (const std::string & dataname, auto hdf5_class,
695 : auto hdf5_type, auto & datavec)
696 : {
697 9 : const hid_t data = H5Dopen1(group, dataname.c_str());
698 9 : check_open(data, dataname.c_str());
699 :
700 9 : const hid_t data_type = H5Dget_type(data);
701 9 : check_hdf5(data_type, dataname + " type");
702 :
703 : // HDF5 will convert between the file's integer type and ours, but
704 : // we do expect an integer type.
705 9 : if (H5Tget_class(data_type) != hdf5_class)
706 0 : libmesh_error_msg("Unexpected type for " + dataname + " in " + filename);
707 :
708 9 : H5Tclose(data_type);
709 :
710 9 : const hid_t dataspace = H5Dget_space(data);
711 9 : check_hdf5(dataspace, dataname + " space");
712 :
713 9 : int ndims = H5Sget_simple_extent_ndims(dataspace);
714 9 : if (ndims != 1)
715 0 : libmesh_error_msg("Non-vector space for " + dataname + " in " + filename);
716 :
717 : hsize_t len, maxlen;
718 9 : herr_t errval = H5Sget_simple_extent_dims(dataspace, &len, &maxlen);
719 9 : check_hdf5(errval, dataname + " dims");
720 :
721 9 : datavec.resize(len);
722 :
723 9 : errval = H5Dread(data, hdf5_type, H5S_ALL, H5S_ALL, H5P_DEFAULT, datavec.data());
724 18 : check_hdf5(errval, dataname + " read");
725 : };
726 :
727 3 : read_vector("cols", H5T_INTEGER, H5T_NATIVE_ULLONG, cols);
728 3 : read_vector("row_offsets", H5T_INTEGER, H5T_NATIVE_ULLONG, row_offsets);
729 6 : read_vector("vals", H5T_FLOAT, H5T_NATIVE_DOUBLE, vals);
730 :
731 3 : if (cols.size() != vals.size())
732 0 : libmesh_error_msg("Inconsistent cols/vals sizes in " + filename);
733 :
734 3 : if (row_offsets.size() != num_rows + 1)
735 0 : libmesh_error_msg("Inconsistent row_offsets size in " + filename);
736 :
737 : // Data for the row we're working on
738 : numeric_index_type current_row = 0;
739 : processor_id_type current_proc = 0;
740 3 : numeric_index_type current_on_diagonal_nonzeros = 0;
741 3 : numeric_index_type current_off_diagonal_nonzeros = 0;
742 3 : if (row_offsets[0] != 0)
743 0 : libmesh_error_msg("Unexpected row_offsets[0] in " + filename);
744 :
745 381 : for (auto i : index_range(cols))
746 : {
747 456 : while (row_offsets[current_row+1] <= i)
748 : {
749 : ++current_row;
750 78 : if (row_offsets[current_row] < row_offsets[current_row-1])
751 0 : libmesh_error_msg("Non-monotonic row_offsets in " + filename);
752 78 : current_on_diagonal_nonzeros = 0;
753 78 : current_off_diagonal_nonzeros = 0;
754 : }
755 :
756 379 : while (current_row >= new_row_stops[current_proc])
757 1 : ++current_proc;
758 :
759 : // 0-based indexing in file
760 378 : if (cols[i] >= new_col_starts[current_proc] &&
761 352 : cols[i] < new_col_stops[current_proc])
762 : {
763 317 : ++current_on_diagonal_nonzeros;
764 317 : on_diagonal_nonzeros =
765 : std::max(on_diagonal_nonzeros,
766 : current_on_diagonal_nonzeros);
767 : }
768 : else
769 : {
770 61 : ++current_off_diagonal_nonzeros;
771 61 : off_diagonal_nonzeros =
772 : std::max(off_diagonal_nonzeros,
773 : current_off_diagonal_nonzeros);
774 : }
775 : }
776 : }
777 :
778 4 : this->comm().broadcast(on_diagonal_nonzeros);
779 4 : this->comm().broadcast(off_diagonal_nonzeros);
780 :
781 4 : this->init(num_rows, num_cols,
782 : new_row_stop-new_row_start,
783 : new_col_stop-new_col_start,
784 : on_diagonal_nonzeros,
785 : off_diagonal_nonzeros);
786 :
787 : // Set the matrix values last.
788 4 : if (this->processor_id() == 0)
789 : {
790 : numeric_index_type current_row = 0;
791 381 : for (auto i : index_range(cols))
792 : {
793 456 : while (row_offsets[current_row+1] <= i)
794 : {
795 : ++current_row;
796 : libmesh_assert_greater_equal (row_offsets[current_row],
797 : row_offsets[current_row-1]);
798 : }
799 378 : this->set(current_row, cols[i], vals[i]);
800 : }
801 : }
802 :
803 4 : this->close();
804 : #endif // LIBMESH_HAVE_HDF5
805 4 : }
806 :
807 :
808 :
809 : template <typename T>
810 1415 : void SparseMatrix<T>::read_matlab(const std::string & filename)
811 : {
812 84 : LOG_SCOPE("read_matlab()", "SparseMatrix");
813 :
814 : #ifndef LIBMESH_HAVE_CXX11_REGEX
815 : libmesh_not_implemented(); // What is your compiler?!? Email us!
816 : libmesh_ignore(filename);
817 : #else
818 42 : parallel_object_only();
819 :
820 1415 : const bool gzipped_file = Utility::ends_with(filename, ".gz");
821 :
822 : // The sizes we get from the file
823 1415 : std::size_t m = 0,
824 1415 : n = 0;
825 :
826 : // If we don't already have this size, we'll need to reinit, and
827 : // determine which rows+columns each processor is in charge of.
828 84 : std::vector<numeric_index_type> new_row_starts, new_row_stops,
829 84 : new_col_starts, new_col_stops;
830 :
831 : numeric_index_type new_row_start, new_row_stop,
832 : new_col_start, new_col_stop;
833 :
834 : // We'll read through the file three times: once to get a reliable
835 : // value for the matrix size (so we can divvy it up among
836 : // processors), then again to get the sparsity to send to each
837 : // processor, then a final time to get the entries to send to each
838 : // processor.
839 : //
840 : // We'll use an istream here; it might be an ifstream if we're
841 : // opening a raw ASCII file or a gzstream if we're opening a
842 : // compressed one.
843 1373 : std::unique_ptr<std::istream> file;
844 :
845 : // We'll need a temporary structure to cache matrix entries, because
846 : // we need to read through the whole file before we know the size
847 : // and sparsity structure with which we can init().
848 : //
849 : // Reading through the file three times via `seekg` doesn't work
850 : // with our gzstream wrapper, and seems to take three times as long
851 : // even with a plain ifstream. What happened to disk caching!?
852 84 : std::vector<std::tuple<numeric_index_type, numeric_index_type, T>> entries;
853 :
854 : // First read through the file, saving size and entry data
855 : {
856 : // We'll read the matrix on processor 0 rather than try to juggle
857 : // parallel I/O.
858 1457 : if (this->processor_id() == 0)
859 : {
860 : // We'll be using regular expressions to make ourselves slightly
861 : // more robust to formatting.
862 705 : const std::regex start_regex // assignment like "zzz = ["
863 : ("\\s*\\w+\\s*=\\s*\\[");
864 705 : const std::regex end_regex // end of assignment
865 : ("^[^%]*\\]");
866 :
867 647 : if (gzipped_file)
868 : {
869 : #ifdef LIBMESH_HAVE_GZSTREAM
870 255 : auto inf = std::make_unique<igzstream>();
871 9 : libmesh_assert(inf);
872 9 : inf->open(filename.c_str(), std::ios::in);
873 237 : file = std::move(inf);
874 : #else
875 : libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
876 : #endif
877 228 : }
878 : else
879 : {
880 421 : auto inf = std::make_unique<std::ifstream>();
881 20 : libmesh_assert(inf);
882 :
883 421 : std::string new_name = Utility::unzip_file(filename);
884 :
885 401 : inf->open(new_name.c_str(), std::ios::in);
886 381 : file = std::move(inf);
887 361 : }
888 :
889 : // If we have a matrix with all-zero trailing rows, the only
890 : // way to get the size is if it ended up in a comment
891 705 : const std::regex size_regex // comment like "% size = 8 8"
892 : ("%\\s*[Ss][Ii][Zz][Ee]\\s*=\\s*(\\d+)\\s+(\\d+)");
893 676 : const std::string whitespace = " \t";
894 :
895 29 : bool have_started = false;
896 29 : bool have_ended = false;
897 647 : std::size_t largest_i_seen = 0, largest_j_seen = 0;
898 :
899 : // Data for the row we're working on
900 : // Use 1-based indexing for current_row, as in the file
901 647 : std::size_t current_row = 1;
902 :
903 334671 : for (std::string line; std::getline(*file, line);)
904 : {
905 7616 : std::smatch sm;
906 :
907 : // First, try to match an entry. This is the most common
908 : // case so we won't rely on slow std::regex for it.
909 : // stringstream is at least an improvement over that.
910 :
911 : // Look for row/col/val like "1 1 -2.0e-4"
912 :
913 182524 : std::istringstream l(line);
914 :
915 : std::size_t i, j;
916 1822 : T value;
917 :
918 9438 : l >> i >> j >> value;
919 :
920 174937 : if (!l.fail())
921 : {
922 170408 : libmesh_error_msg_if
923 : (!have_started, "Confused by premature entries in matrix file " << filename);
924 :
925 324217 : entries.emplace_back(cast_int<numeric_index_type>(i),
926 177821 : cast_int<numeric_index_type>(j),
927 : value);
928 :
929 170408 : libmesh_error_msg_if
930 : (!i || !j, "Expected 1-based indexing in matrix file "
931 : << filename);
932 :
933 170408 : current_row = std::max(current_row, i);
934 :
935 170408 : libmesh_error_msg_if
936 : (i < current_row,
937 : "Can't handle out-of-order entries in matrix file "
938 : << filename);
939 :
940 170408 : largest_i_seen = std::max(i, largest_i_seen);
941 170408 : largest_j_seen = std::max(j, largest_j_seen);
942 : }
943 :
944 4529 : else if (std::regex_search(line, sm, size_regex))
945 : {
946 676 : const std::string msize = sm[1];
947 647 : const std::string nsize = sm[2];
948 647 : m = std::stoull(msize);
949 647 : n = std::stoull(nsize);
950 : }
951 :
952 3882 : else if (std::regex_search(line, start_regex))
953 29 : have_started = true;
954 :
955 3235 : else if (std::regex_search(line, end_regex))
956 : {
957 29 : have_ended = true;
958 58 : break;
959 : }
960 : }
961 :
962 647 : libmesh_error_msg_if
963 : (!have_started, "Confused by missing assignment beginning in matrix file " << filename);
964 :
965 647 : libmesh_error_msg_if
966 : (!have_ended, "Confused by missing assignment ending in matrix file " << filename);
967 :
968 647 : libmesh_error_msg_if
969 : (m > largest_i_seen, "Confused by missing final row(s) in matrix file " << filename);
970 :
971 647 : libmesh_error_msg_if
972 : (m > 0 && m < largest_i_seen, "Confused by extra final row(s) in matrix file " << filename);
973 :
974 647 : if (!m)
975 0 : m = largest_i_seen;
976 :
977 647 : libmesh_error_msg_if
978 : (n > largest_j_seen, "Confused by missing final column(s) in matrix file " << filename);
979 :
980 647 : libmesh_error_msg_if
981 : (n > 0 && n < largest_j_seen, "Confused by extra final column(s) in matrix file " << filename);
982 :
983 647 : if (!n)
984 0 : n = largest_j_seen;
985 :
986 647 : this->comm().broadcast(m);
987 647 : this->comm().broadcast(n);
988 589 : }
989 : else
990 : {
991 755 : this->comm().broadcast(m);
992 768 : this->comm().broadcast(n);
993 : }
994 :
995 1415 : if (this->initialized() &&
996 1417 : m == this->m() &&
997 70 : n == this->n())
998 : {
999 70 : new_row_start = this->row_start(),
1000 70 : new_row_stop = this->row_stop();
1001 :
1002 72 : new_col_start = this->col_start(),
1003 70 : new_col_stop = this->col_stop();
1004 : }
1005 : else
1006 : {
1007 : // Determine which rows/columns each processor will be in charge of
1008 1385 : new_row_start = this->processor_id() * m / this->n_processors(),
1009 1345 : new_row_stop = (this->processor_id()+1) * m / this->n_processors();
1010 :
1011 1385 : new_col_start = this->processor_id() * n / this->n_processors(),
1012 1345 : new_col_stop = (this->processor_id()+1) * n / this->n_processors();
1013 : }
1014 :
1015 1415 : this->comm().gather(0, new_row_start, new_row_starts);
1016 1415 : this->comm().gather(0, new_row_stop, new_row_stops);
1017 1415 : this->comm().gather(0, new_col_start, new_col_starts);
1018 1415 : this->comm().gather(0, new_col_stop, new_col_stops);
1019 :
1020 : } // Done reading entry data and broadcasting matrix size
1021 :
1022 : // Calculate the matrix sparsity and initialize it second
1023 : {
1024 : // Deduce the sparsity pattern, or at least the maximum number of
1025 : // on- and off- diagonal non-zeros per row.
1026 1415 : numeric_index_type on_diagonal_nonzeros =0,
1027 1415 : off_diagonal_nonzeros =0;
1028 :
1029 1457 : if (this->processor_id() == 0)
1030 : {
1031 : // Data for the row we're working on
1032 : // Use 1-based indexing for current_row, as in the file
1033 29 : numeric_index_type current_row = 1;
1034 29 : processor_id_type current_proc = 0;
1035 647 : numeric_index_type current_on_diagonal_nonzeros = 0;
1036 647 : numeric_index_type current_off_diagonal_nonzeros = 0;
1037 :
1038 171055 : for (auto [i, j, value] : entries)
1039 : {
1040 170408 : if (i > current_row)
1041 : {
1042 876 : current_row = i;
1043 : // +1 for 1-based indexing in file
1044 22331 : while (current_row >= (new_row_stops[current_proc]+1))
1045 768 : ++current_proc;
1046 20674 : current_on_diagonal_nonzeros = 0;
1047 20674 : current_off_diagonal_nonzeros = 0;
1048 : }
1049 :
1050 : // +1 for 1-based indexing in file
1051 184358 : if (j >= (new_col_starts[current_proc]+1) &&
1052 165360 : j < (new_col_stops[current_proc]+1))
1053 : {
1054 143791 : ++current_on_diagonal_nonzeros;
1055 143791 : on_diagonal_nonzeros =
1056 : std::max(on_diagonal_nonzeros,
1057 5555 : current_on_diagonal_nonzeros);
1058 : }
1059 : else
1060 : {
1061 26617 : ++current_off_diagonal_nonzeros;
1062 26617 : off_diagonal_nonzeros =
1063 : std::max(off_diagonal_nonzeros,
1064 1858 : current_off_diagonal_nonzeros);
1065 : }
1066 : }
1067 : }
1068 :
1069 1415 : this->comm().broadcast(on_diagonal_nonzeros);
1070 1415 : this->comm().broadcast(off_diagonal_nonzeros);
1071 :
1072 1415 : this->init(m, n,
1073 : new_row_stop-new_row_start,
1074 : new_col_stop-new_col_start,
1075 : on_diagonal_nonzeros,
1076 : off_diagonal_nonzeros);
1077 : }
1078 :
1079 : // Set the matrix values last.
1080 : // Convert from 1-based to 0-based indexing
1081 1457 : if (this->processor_id() == 0)
1082 171055 : for (auto [i, j, value] : entries)
1083 170408 : this->set(i-1, j-1, value);
1084 :
1085 1415 : this->close();
1086 : #endif
1087 2746 : }
1088 :
1089 :
1090 :
1091 : template <typename T>
1092 0 : void SparseMatrix<T>::read_petsc_binary(const std::string &)
1093 : {
1094 0 : libmesh_not_implemented_msg
1095 : ("libMesh cannot read PETSc binary-format files into non-PETSc matrices");
1096 : }
1097 :
1098 :
1099 :
1100 : template <typename T>
1101 0 : void SparseMatrix<T>::read_petsc_hdf5(const std::string &)
1102 : {
1103 0 : libmesh_not_implemented_msg
1104 : ("libMesh cannot read PETSc HDF5-format files into non-PETSc matrices");
1105 : }
1106 :
1107 :
1108 :
1109 : template <typename T>
1110 75 : void SparseMatrix<T>::scale(const T scale)
1111 : {
1112 4 : libmesh_assert(this->closed());
1113 :
1114 375 : for (const auto i : make_range(this->row_start(), this->row_stop()))
1115 1500 : for (const auto j : make_range(this->col_start(), this->col_stop()))
1116 1200 : this->set(i, j, (*this)(i, j) * scale);
1117 75 : }
1118 :
1119 :
1120 :
1121 : template <typename T>
1122 290 : Real SparseMatrix<T>::l1_norm_diff(const SparseMatrix<T> & other_mat) const
1123 : {
1124 290 : auto diff_mat = this->clone();
1125 290 : diff_mat->add(-1.0, other_mat);
1126 580 : return diff_mat->l1_norm();
1127 266 : }
1128 :
1129 : //------------------------------------------------------------------
1130 : // Explicit instantiations
1131 : template class LIBMESH_EXPORT SparseMatrix<Number>;
1132 :
1133 : } // namespace libMesh
|
