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 : // C++ includes
43 : #include <memory>
44 : #include <fstream>
45 :
46 : #ifdef LIBMESH_HAVE_CXX11_REGEX
47 : #include <regex>
48 : #endif
49 :
50 :
51 : namespace libMesh
52 : {
53 :
54 :
55 : //------------------------------------------------------------------
56 : // SparseMatrix Methods
57 :
58 :
59 : // Constructor
60 : template <typename T>
61 107648 : SparseMatrix<T>::SparseMatrix (const Parallel::Communicator & comm_in) :
62 : ParallelObject(comm_in),
63 101764 : _dof_map(nullptr),
64 101764 : _sp(nullptr),
65 101764 : _is_initialized(false),
66 107648 : _use_hash_table(false)
67 107648 : {}
68 :
69 :
70 :
71 : template <typename T>
72 47966 : void SparseMatrix<T>::attach_dof_map (const DofMap & dof_map)
73 : {
74 47966 : _dof_map = &dof_map;
75 47966 : if (!_sp)
76 22635 : _sp = dof_map.get_sparsity_pattern();
77 47966 : }
78 :
79 :
80 :
81 : template <typename T>
82 41068 : void SparseMatrix<T>::attach_sparsity_pattern (const SparsityPattern::Build & sp)
83 : {
84 41068 : _sp = &sp;
85 41068 : }
86 :
87 :
88 :
89 : // default implementation is to fall back to non-blocked method
90 : template <typename T>
91 0 : void SparseMatrix<T>::add_block_matrix (const DenseMatrix<T> & dm,
92 : const std::vector<numeric_index_type> & brows,
93 : const std::vector<numeric_index_type> & bcols)
94 : {
95 0 : libmesh_assert_equal_to (dm.m() / brows.size(), dm.n() / bcols.size());
96 :
97 : const numeric_index_type blocksize = cast_int<numeric_index_type>
98 0 : (dm.m() / brows.size());
99 :
100 0 : libmesh_assert_equal_to (dm.m()%blocksize, 0);
101 0 : libmesh_assert_equal_to (dm.n()%blocksize, 0);
102 :
103 0 : std::vector<numeric_index_type> rows, cols;
104 :
105 0 : rows.reserve(blocksize*brows.size());
106 0 : cols.reserve(blocksize*bcols.size());
107 :
108 0 : for (auto & row : brows)
109 : {
110 0 : numeric_index_type i = row * blocksize;
111 :
112 0 : for (unsigned int v=0; v<blocksize; v++)
113 0 : rows.push_back(i++);
114 : }
115 :
116 0 : for (auto & col : bcols)
117 : {
118 0 : numeric_index_type j = col * blocksize;
119 :
120 0 : for (unsigned int v=0; v<blocksize; v++)
121 0 : cols.push_back(j++);
122 : }
123 :
124 0 : this->add_matrix (dm, rows, cols);
125 0 : }
126 :
127 :
128 :
129 : // Full specialization of print method for Complex datatypes
130 : template <>
131 0 : void SparseMatrix<Complex>::print(std::ostream & os, const bool sparse) const
132 : {
133 : // std::complex<>::operator<<() is defined, but use this form
134 :
135 0 : if (sparse)
136 : {
137 0 : libmesh_not_implemented();
138 : }
139 :
140 0 : os << "Real part:" << std::endl;
141 0 : for (auto i : make_range(this->m()))
142 : {
143 0 : for (auto j : make_range(this->n()))
144 0 : os << std::setw(8) << (*this)(i,j).real() << " ";
145 0 : os << std::endl;
146 : }
147 :
148 0 : os << std::endl << "Imaginary part:" << std::endl;
149 0 : for (auto i : make_range(this->m()))
150 : {
151 0 : for (auto j : make_range(this->n()))
152 0 : os << std::setw(8) << (*this)(i,j).imag() << " ";
153 0 : os << std::endl;
154 : }
155 0 : }
156 :
157 :
158 :
159 :
160 :
161 :
162 : // Full specialization for Real datatypes
163 : template <typename T>
164 : std::unique_ptr<SparseMatrix<T>>
165 24092 : SparseMatrix<T>::build(const Parallel::Communicator & comm,
166 : const SolverPackage solver_package,
167 : const MatrixBuildType matrix_build_type /* = AUTOMATIC */)
168 : {
169 : // Avoid unused parameter warnings when no solver packages are enabled.
170 684 : libmesh_ignore(comm);
171 :
172 24092 : if (matrix_build_type == MatrixBuildType::DIAGONAL)
173 0 : return std::make_unique<DiagonalMatrix<T>>(comm);
174 :
175 : // Build the appropriate vector
176 24092 : switch (solver_package)
177 : {
178 :
179 : #ifdef LIBMESH_HAVE_LASPACK
180 0 : case LASPACK_SOLVERS:
181 0 : return std::make_unique<LaspackMatrix<T>>(comm);
182 : #endif
183 :
184 :
185 : #ifdef LIBMESH_HAVE_PETSC
186 23845 : case PETSC_SOLVERS:
187 23845 : return std::make_unique<PetscMatrix<T>>(comm);
188 : #endif
189 :
190 :
191 : #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
192 : case TRILINOS_SOLVERS:
193 : return std::make_unique<EpetraMatrix<T>>(comm);
194 : #endif
195 :
196 :
197 : #ifdef LIBMESH_HAVE_EIGEN
198 247 : case EIGEN_SOLVERS:
199 247 : return std::make_unique<EigenSparseMatrix<T>>(comm);
200 : #endif
201 :
202 0 : default:
203 0 : libmesh_error_msg("ERROR: Unrecognized solver package: " << solver_package);
204 : }
205 : }
206 :
207 :
208 : template <typename T>
209 2958121 : void SparseMatrix<T>::vector_mult (NumericVector<T> & dest,
210 : const NumericVector<T> & arg) const
211 : {
212 2958121 : dest.zero();
213 167712 : this->vector_mult_add(dest,arg);
214 2958121 : }
215 :
216 :
217 :
218 : template <typename T>
219 168700 : void SparseMatrix<T>::vector_mult_add (NumericVector<T> & dest,
220 : const NumericVector<T> & arg) const
221 : {
222 : /* This functionality is actually implemented in the \p
223 : NumericVector class. */
224 2959109 : dest.add_vector(arg,*this);
225 168700 : }
226 :
227 :
228 :
229 : template <typename T>
230 0 : void SparseMatrix<T>::zero_rows (std::vector<numeric_index_type> &, T)
231 : {
232 : /* This functionality isn't implemented or stubbed in every subclass yet */
233 0 : libmesh_not_implemented();
234 : }
235 :
236 :
237 : template <typename T>
238 14 : std::size_t SparseMatrix<T>::n_nonzeros() const
239 : {
240 14 : if (!_sp)
241 2 : return 0;
242 0 : return _sp->n_nonzeros();
243 : }
244 :
245 :
246 : template <typename T>
247 0 : void SparseMatrix<T>::print(std::ostream & os, const bool sparse) const
248 : {
249 0 : parallel_object_only();
250 :
251 0 : libmesh_assert (this->initialized());
252 :
253 0 : const numeric_index_type first_dof = this->row_start(),
254 0 : end_dof = this->row_stop();
255 :
256 : // We'll print the matrix from processor 0 to make sure
257 : // it's serialized properly
258 0 : if (this->processor_id() == 0)
259 : {
260 0 : libmesh_assert_equal_to (first_dof, 0);
261 0 : for (numeric_index_type i : make_range(end_dof))
262 : {
263 0 : if (sparse)
264 : {
265 0 : for (auto j : make_range(this->n()))
266 : {
267 0 : T c = (*this)(i,j);
268 0 : if (c != static_cast<T>(0.0))
269 : {
270 0 : os << i << " " << j << " " << c << std::endl;
271 : }
272 : }
273 : }
274 : else
275 : {
276 0 : for (auto j : make_range(this->n()))
277 0 : os << (*this)(i,j) << " ";
278 0 : os << std::endl;
279 : }
280 : }
281 :
282 0 : std::vector<numeric_index_type> ibuf, jbuf;
283 0 : std::vector<T> cbuf;
284 0 : numeric_index_type currenti = end_dof;
285 0 : for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
286 : {
287 0 : this->comm().receive(p, ibuf);
288 0 : this->comm().receive(p, jbuf);
289 0 : this->comm().receive(p, cbuf);
290 0 : libmesh_assert_equal_to (ibuf.size(), jbuf.size());
291 0 : libmesh_assert_equal_to (ibuf.size(), cbuf.size());
292 :
293 0 : if (ibuf.empty())
294 0 : continue;
295 0 : libmesh_assert_greater_equal (ibuf.front(), currenti);
296 0 : libmesh_assert_greater_equal (ibuf.back(), ibuf.front());
297 :
298 0 : std::size_t currentb = 0;
299 0 : for (;currenti <= ibuf.back(); ++currenti)
300 : {
301 0 : if (sparse)
302 : {
303 0 : for (numeric_index_type j=0; j<this->n(); j++)
304 : {
305 0 : if (currentb < ibuf.size() &&
306 0 : ibuf[currentb] == currenti &&
307 0 : jbuf[currentb] == j)
308 : {
309 0 : os << currenti << " " << j << " " << cbuf[currentb] << std::endl;
310 0 : currentb++;
311 : }
312 : }
313 : }
314 : else
315 : {
316 0 : for (auto j : make_range(this->n()))
317 : {
318 0 : if (currentb < ibuf.size() &&
319 0 : ibuf[currentb] == currenti &&
320 0 : jbuf[currentb] == j)
321 : {
322 0 : os << cbuf[currentb] << " ";
323 0 : currentb++;
324 : }
325 : else
326 0 : os << static_cast<T>(0.0) << " ";
327 : }
328 0 : os << std::endl;
329 : }
330 : }
331 : }
332 0 : if (!sparse)
333 : {
334 0 : for (; currenti != this->m(); ++currenti)
335 : {
336 0 : for (numeric_index_type j=0; j<this->n(); j++)
337 0 : os << static_cast<T>(0.0) << " ";
338 0 : os << std::endl;
339 : }
340 : }
341 : }
342 : else
343 : {
344 0 : std::vector<numeric_index_type> ibuf, jbuf;
345 0 : std::vector<T> cbuf;
346 :
347 : // We'll assume each processor has access to entire
348 : // matrix rows, so (*this)(i,j) is valid if i is a local index.
349 0 : for (numeric_index_type i : make_range(first_dof, end_dof))
350 : {
351 0 : for (auto j : make_range(this->n()))
352 : {
353 0 : T c = (*this)(i,j);
354 0 : if (c != static_cast<T>(0.0))
355 : {
356 0 : ibuf.push_back(i);
357 0 : jbuf.push_back(j);
358 0 : cbuf.push_back(c);
359 : }
360 : }
361 : }
362 0 : this->comm().send(0,ibuf);
363 0 : this->comm().send(0,jbuf);
364 0 : this->comm().send(0,cbuf);
365 : }
366 0 : }
367 :
368 :
369 : template <typename T>
370 14 : void SparseMatrix<T>::print_matlab(const std::string & name) const
371 : {
372 2 : parallel_object_only();
373 :
374 2 : libmesh_assert (this->initialized());
375 :
376 14 : const numeric_index_type first_dof = this->row_start(),
377 14 : end_dof = this->row_stop();
378 :
379 : // We'll print the matrix from processor 0 to make sure
380 : // it's serialized properly
381 16 : if (this->processor_id() == 0)
382 : {
383 12 : std::unique_ptr<std::ofstream> file;
384 :
385 14 : if (name != "")
386 24 : file = std::make_unique<std::ofstream>(name.c_str());
387 :
388 14 : std::ostream & os = (name == "") ? libMesh::out : *file;
389 :
390 14 : std::size_t sparsity_nonzeros = this->n_nonzeros();
391 :
392 2 : std::size_t real_nonzeros = 0;
393 :
394 2 : libmesh_assert_equal_to(first_dof, 0);
395 70 : for (numeric_index_type i : make_range(end_dof))
396 : {
397 320 : for (auto j : make_range(this->n()))
398 : {
399 224 : T c = (*this)(i,j);
400 208 : if (c != static_cast<T>(0.0))
401 224 : ++real_nonzeros;
402 : }
403 : }
404 :
405 :
406 14 : for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
407 : {
408 0 : std::size_t nonzeros_on_p = 0;
409 0 : this->comm().receive(p, nonzeros_on_p);
410 0 : real_nonzeros += nonzeros_on_p;
411 : }
412 :
413 : if (sparsity_nonzeros &&
414 : sparsity_nonzeros != real_nonzeros)
415 : libmesh_warning(sparsity_nonzeros <<
416 : " nonzeros allocated, but " <<
417 : real_nonzeros << " used.");
418 :
419 : // We probably want to be more consistent than that, if our
420 : // sparsity is overallocated.
421 :
422 : // Print a header similar to PETSc's mat_view ascii_matlab
423 14 : os << "%Mat Object: () " << this->n_processors() << " MPI processes\n"
424 6 : << "% type: " << (this->n_processors() > 1 ? "mpi" : "seq") << "aij\n"
425 40 : << "% Size = " << this->m() << ' ' << this->n() << '\n'
426 14 : << "% Nonzeros = " << real_nonzeros << '\n'
427 14 : << "zzz = zeros(" << real_nonzeros << ",3);\n"
428 14 : << "zzz = [\n";
429 :
430 70 : for (numeric_index_type i : make_range(end_dof))
431 : {
432 : // FIXME - we need a base class way to iterate over a
433 : // SparseMatrix row.
434 320 : for (auto j : make_range(this->n()))
435 : {
436 224 : T c = (*this)(i,j);
437 208 : if (c != static_cast<T>(0.0))
438 : {
439 : // Convert from 0-based to 1-based indexing
440 432 : os << (i+1) << ' ' << (j+1) << " " << c << '\n';
441 : }
442 : }
443 : }
444 :
445 4 : std::vector<numeric_index_type> ibuf, jbuf;
446 4 : std::vector<T> cbuf;
447 14 : for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
448 : {
449 0 : this->comm().receive(p, ibuf);
450 0 : this->comm().receive(p, jbuf);
451 0 : this->comm().receive(p, cbuf);
452 0 : libmesh_assert_equal_to (ibuf.size(), jbuf.size());
453 0 : libmesh_assert_equal_to (ibuf.size(), cbuf.size());
454 :
455 0 : for (auto n : index_range(ibuf))
456 0 : os << ibuf[n] << ' ' << jbuf[n] << " " << cbuf[n] << '\n';
457 : }
458 :
459 12 : os << "];\n" << "Mat_sparse = spconvert(zzz);" << std::endl;
460 10 : }
461 : else
462 : {
463 0 : std::vector<numeric_index_type> ibuf, jbuf;
464 0 : std::vector<T> cbuf;
465 0 : std::size_t my_nonzeros = 0;
466 :
467 : // We'll assume each processor has access to entire
468 : // matrix rows, so (*this)(i,j) is valid if i is a local index.
469 0 : for (numeric_index_type i : make_range(first_dof, end_dof))
470 : {
471 0 : for (auto j : make_range(this->n()))
472 : {
473 0 : T c = (*this)(i,j);
474 0 : if (c != static_cast<T>(0.0))
475 : {
476 0 : ibuf.push_back(i);
477 0 : jbuf.push_back(j);
478 0 : cbuf.push_back(c);
479 0 : ++my_nonzeros;
480 : }
481 : }
482 : }
483 0 : this->comm().send(0,my_nonzeros);
484 0 : this->comm().send(0,ibuf);
485 0 : this->comm().send(0,jbuf);
486 0 : this->comm().send(0,cbuf);
487 : }
488 14 : }
489 :
490 :
491 :
492 : template <typename T>
493 0 : void SparseMatrix<T>::print_petsc_binary(const std::string &)
494 : {
495 0 : libmesh_not_implemented_msg
496 : ("libMesh cannot write PETSc binary-format files from non-PETSc matrices");
497 : }
498 :
499 :
500 :
501 : template <typename T>
502 0 : void SparseMatrix<T>::print_petsc_hdf5(const std::string &)
503 : {
504 0 : libmesh_not_implemented_msg
505 : ("libMesh cannot write PETSc HDF5-format files from non-PETSc matrices");
506 : }
507 :
508 :
509 :
510 : template <typename T>
511 1060 : void SparseMatrix<T>::read(const std::string & filename)
512 : {
513 : {
514 1124 : std::ifstream in (filename.c_str());
515 1060 : libmesh_error_msg_if
516 : (!in.good(), "ERROR: cannot read file:\n\t" <<
517 : filename);
518 996 : }
519 :
520 1060 : std::string_view basename = Utility::basename_of(filename);
521 :
522 1060 : const bool gzipped_file = Utility::ends_with(filename, ".gz");
523 :
524 1060 : if (gzipped_file)
525 12 : basename.remove_suffix(3);
526 :
527 2120 : if (Utility::ends_with(basename, ".matlab") ||
528 1160 : Utility::ends_with(basename, ".m"))
529 990 : this->read_matlab(filename);
530 70 : else if (Utility::ends_with(basename, ".petsc64"))
531 : {
532 : #ifndef LIBMESH_HAVE_PETSC
533 0 : libmesh_error_msg("Cannot load PETSc matrix file " <<
534 : filename << " without PETSc-enabled libMesh.");
535 : #endif
536 : #if LIBMESH_DOF_ID_BYTES != 8
537 0 : libmesh_error_msg("Cannot load 64-bit PETSc matrix file " <<
538 : filename << " with non-64-bit libMesh.");
539 : #endif
540 66 : this->read_petsc_binary(filename);
541 : }
542 4 : else if (Utility::ends_with(basename, ".petsc32"))
543 : {
544 : #ifndef LIBMESH_HAVE_PETSC
545 0 : libmesh_error_msg("Cannot load PETSc matrix file " <<
546 : filename << " without PETSc-enabled libMesh.");
547 : #endif
548 : #if LIBMESH_DOF_ID_BYTES != 4
549 0 : libmesh_error_msg("Cannot load 32-bit PETSc matrix file " <<
550 : filename << " with non-32-bit libMesh.");
551 : #endif
552 4 : this->read_petsc_binary(filename);
553 : }
554 : else
555 0 : libmesh_error_msg(" ERROR: Unrecognized matrix file extension on: "
556 : << basename
557 : << "\n I understand the following:\n\n"
558 : << " *.matlab -- Matlab sparse matrix format\n"
559 : << " *.m -- Matlab sparse matrix format\n"
560 : << " *.petsc32 -- PETSc binary format, 32-bit\n"
561 : << " *.petsc64 -- PETSc binary format, 64-bit\n"
562 : );
563 1060 : }
564 :
565 :
566 : template <typename T>
567 1411 : void SparseMatrix<T>::read_matlab(const std::string & filename)
568 : {
569 84 : LOG_SCOPE("read_matlab()", "SparseMatrix");
570 :
571 : #ifndef LIBMESH_HAVE_CXX11_REGEX
572 : libmesh_not_implemented(); // What is your compiler?!? Email us!
573 : libmesh_ignore(filename);
574 : #else
575 42 : parallel_object_only();
576 :
577 1411 : const bool gzipped_file = Utility::ends_with(filename, ".gz");
578 :
579 : // The sizes we get from the file
580 1411 : std::size_t m = 0,
581 1411 : n = 0;
582 :
583 : // If we don't already have this size, we'll need to reinit, and
584 : // determine which rows+columns each processor is in charge of.
585 84 : std::vector<numeric_index_type> new_row_starts, new_row_stops,
586 84 : new_col_starts, new_col_stops;
587 :
588 : numeric_index_type new_row_start, new_row_stop,
589 : new_col_start, new_col_stop;
590 :
591 : // We'll read through the file three times: once to get a reliable
592 : // value for the matrix size (so we can divvy it up among
593 : // processors), then again to get the sparsity to send to each
594 : // processor, then a final time to get the entries to send to each
595 : // processor.
596 : //
597 : // We'll use an istream here; it might be an ifstream if we're
598 : // opening a raw ASCII file or a gzstream if we're opening a
599 : // compressed one.
600 1369 : std::unique_ptr<std::istream> file;
601 :
602 : // We'll need a temporary structure to cache matrix entries, because
603 : // we need to read through the whole file before we know the size
604 : // and sparsity structure with which we can init().
605 : //
606 : // Reading through the file three times via `seekg` doesn't work
607 : // with our gzstream wrapper, and seems to take three times as long
608 : // even with a plain ifstream. What happened to disk caching!?
609 84 : std::vector<std::tuple<numeric_index_type, numeric_index_type, T>> entries;
610 :
611 : // First read through the file, saving size and entry data
612 : {
613 : // We'll read the matrix on processor 0 rather than try to juggle
614 : // parallel I/O.
615 1453 : if (this->processor_id() == 0)
616 : {
617 : // We'll be using regular expressions to make ourselves slightly
618 : // more robust to formatting.
619 702 : const std::regex start_regex // assignment like "zzz = ["
620 : ("\\s*\\w+\\s*=\\s*\\[");
621 702 : const std::regex end_regex // end of assignment
622 : ("^[^%]*\\]");
623 :
624 644 : if (gzipped_file)
625 : {
626 : #ifdef LIBMESH_HAVE_GZSTREAM
627 255 : auto inf = std::make_unique<igzstream>();
628 9 : libmesh_assert(inf);
629 9 : inf->open(filename.c_str(), std::ios::in);
630 237 : file = std::move(inf);
631 : #else
632 : libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
633 : #endif
634 228 : }
635 : else
636 : {
637 418 : auto inf = std::make_unique<std::ifstream>();
638 20 : libmesh_assert(inf);
639 :
640 418 : std::string new_name = Utility::unzip_file(filename);
641 :
642 398 : inf->open(new_name.c_str(), std::ios::in);
643 378 : file = std::move(inf);
644 358 : }
645 :
646 : // If we have a matrix with all-zero trailing rows, the only
647 : // way to get the size is if it ended up in a comment
648 702 : const std::regex size_regex // comment like "% size = 8 8"
649 : ("%\\s*[Ss][Ii][Zz][Ee]\\s*=\\s*(\\d+)\\s+(\\d+)");
650 673 : const std::string whitespace = " \t";
651 :
652 29 : bool have_started = false;
653 29 : bool have_ended = false;
654 644 : std::size_t largest_i_seen = 0, largest_j_seen = 0;
655 :
656 : // Data for the row we're working on
657 : // Use 1-based indexing for current_row, as in the file
658 644 : std::size_t current_row = 1;
659 :
660 333873 : for (std::string line; std::getline(*file, line);)
661 : {
662 7616 : std::smatch sm;
663 :
664 : // First, try to match an entry. This is the most common
665 : // case so we won't rely on slow std::regex for it.
666 : // stringstream is at least an improvement over that.
667 :
668 : // Look for row/col/val like "1 1 -2.0e-4"
669 :
670 182125 : std::istringstream l(line);
671 :
672 : std::size_t i, j;
673 1822 : T value;
674 :
675 9438 : l >> i >> j >> value;
676 :
677 174538 : if (!l.fail())
678 : {
679 170030 : libmesh_error_msg_if
680 : (!have_started, "Confused by premature entries in matrix file " << filename);
681 :
682 323461 : entries.emplace_back(cast_int<numeric_index_type>(i),
683 177443 : cast_int<numeric_index_type>(j),
684 : value);
685 :
686 170030 : libmesh_error_msg_if
687 : (!i || !j, "Expected 1-based indexing in matrix file "
688 : << filename);
689 :
690 170030 : current_row = std::max(current_row, i);
691 :
692 170030 : libmesh_error_msg_if
693 : (i < current_row,
694 : "Can't handle out-of-order entries in matrix file "
695 : << filename);
696 :
697 170030 : largest_i_seen = std::max(i, largest_i_seen);
698 170030 : largest_j_seen = std::max(j, largest_j_seen);
699 : }
700 :
701 4508 : else if (std::regex_search(line, sm, size_regex))
702 : {
703 673 : const std::string msize = sm[1];
704 644 : const std::string nsize = sm[2];
705 644 : m = std::stoull(msize);
706 644 : n = std::stoull(nsize);
707 : }
708 :
709 3864 : else if (std::regex_search(line, start_regex))
710 29 : have_started = true;
711 :
712 3220 : else if (std::regex_search(line, end_regex))
713 : {
714 29 : have_ended = true;
715 58 : break;
716 : }
717 : }
718 :
719 644 : libmesh_error_msg_if
720 : (!have_started, "Confused by missing assignment beginning in matrix file " << filename);
721 :
722 644 : libmesh_error_msg_if
723 : (!have_ended, "Confused by missing assignment ending in matrix file " << filename);
724 :
725 644 : libmesh_error_msg_if
726 : (m > largest_i_seen, "Confused by missing final row(s) in matrix file " << filename);
727 :
728 644 : libmesh_error_msg_if
729 : (m > 0 && m < largest_i_seen, "Confused by extra final row(s) in matrix file " << filename);
730 :
731 644 : if (!m)
732 0 : m = largest_i_seen;
733 :
734 644 : libmesh_error_msg_if
735 : (n > largest_j_seen, "Confused by missing final column(s) in matrix file " << filename);
736 :
737 644 : libmesh_error_msg_if
738 : (n > 0 && n < largest_j_seen, "Confused by extra final column(s) in matrix file " << filename);
739 :
740 644 : if (!n)
741 0 : n = largest_j_seen;
742 :
743 644 : this->comm().broadcast(m);
744 644 : this->comm().broadcast(n);
745 586 : }
746 : else
747 : {
748 754 : this->comm().broadcast(m);
749 767 : this->comm().broadcast(n);
750 : }
751 :
752 1411 : if (this->initialized() &&
753 1413 : m == this->m() &&
754 70 : n == this->n())
755 : {
756 70 : new_row_start = this->row_start(),
757 70 : new_row_stop = this->row_stop();
758 :
759 72 : new_col_start = this->col_start(),
760 70 : new_col_stop = this->col_stop();
761 : }
762 : else
763 : {
764 : // Determine which rows/columns each processor will be in charge of
765 1381 : new_row_start = this->processor_id() * m / this->n_processors(),
766 1341 : new_row_stop = (this->processor_id()+1) * m / this->n_processors();
767 :
768 1381 : new_col_start = this->processor_id() * n / this->n_processors(),
769 1341 : new_col_stop = (this->processor_id()+1) * n / this->n_processors();
770 : }
771 :
772 1411 : this->comm().gather(0, new_row_start, new_row_starts);
773 1411 : this->comm().gather(0, new_row_stop, new_row_stops);
774 1411 : this->comm().gather(0, new_col_start, new_col_starts);
775 1411 : this->comm().gather(0, new_col_stop, new_col_stops);
776 :
777 : } // Done reading entry data and broadcasting matrix size
778 :
779 : // Calculate the matrix sparsity and initialize it second
780 : {
781 : // Deduce the sparsity pattern, or at least the maximum number of
782 : // on- and off- diagonal non-zeros per row.
783 1411 : numeric_index_type on_diagonal_nonzeros =0,
784 1411 : off_diagonal_nonzeros =0;
785 :
786 1453 : if (this->processor_id() == 0)
787 : {
788 : // Data for the row we're working on
789 : // Use 1-based indexing for current_row, as in the file
790 29 : numeric_index_type current_row = 1;
791 29 : processor_id_type current_proc = 0;
792 644 : numeric_index_type current_on_diagonal_nonzeros = 0;
793 644 : numeric_index_type current_off_diagonal_nonzeros = 0;
794 :
795 170674 : for (auto [i, j, value] : entries)
796 : {
797 170030 : if (i > current_row)
798 : {
799 876 : current_row = i;
800 : // +1 for 1-based indexing in file
801 22252 : while (current_row >= (new_row_stops[current_proc]+1))
802 767 : ++current_proc;
803 20596 : current_on_diagonal_nonzeros = 0;
804 20596 : current_off_diagonal_nonzeros = 0;
805 : }
806 :
807 : // +1 for 1-based indexing in file
808 183980 : if (j >= (new_col_starts[current_proc]+1) &&
809 165008 : j < (new_col_stops[current_proc]+1))
810 : {
811 143474 : ++current_on_diagonal_nonzeros;
812 143474 : on_diagonal_nonzeros =
813 : std::max(on_diagonal_nonzeros,
814 5555 : current_on_diagonal_nonzeros);
815 : }
816 : else
817 : {
818 26556 : ++current_off_diagonal_nonzeros;
819 26556 : off_diagonal_nonzeros =
820 : std::max(off_diagonal_nonzeros,
821 1858 : current_off_diagonal_nonzeros);
822 : }
823 : }
824 : }
825 :
826 1411 : this->comm().broadcast(on_diagonal_nonzeros);
827 1411 : this->comm().broadcast(off_diagonal_nonzeros);
828 :
829 1411 : this->init(m, n,
830 : new_row_stop-new_row_start,
831 : new_col_stop-new_col_start,
832 : on_diagonal_nonzeros,
833 : off_diagonal_nonzeros);
834 : }
835 :
836 : // Set the matrix values last.
837 : // Convert from 1-based to 0-based indexing
838 1453 : if (this->processor_id() == 0)
839 170674 : for (auto [i, j, value] : entries)
840 170030 : this->set(i-1, j-1, value);
841 :
842 1411 : this->close();
843 : #endif
844 2738 : }
845 :
846 :
847 :
848 : template <typename T>
849 0 : void SparseMatrix<T>::read_petsc_binary(const std::string &)
850 : {
851 0 : libmesh_not_implemented_msg
852 : ("libMesh cannot read PETSc binary-format files into non-PETSc matrices");
853 : }
854 :
855 :
856 :
857 : template <typename T>
858 0 : void SparseMatrix<T>::read_petsc_hdf5(const std::string &)
859 : {
860 0 : libmesh_not_implemented_msg
861 : ("libMesh cannot read PETSc HDF5-format files into non-PETSc matrices");
862 : }
863 :
864 :
865 :
866 : template <typename T>
867 75 : void SparseMatrix<T>::scale(const T scale)
868 : {
869 4 : libmesh_assert(this->closed());
870 :
871 375 : for (const auto i : make_range(this->row_start(), this->row_stop()))
872 1500 : for (const auto j : make_range(this->col_start(), this->col_stop()))
873 1200 : this->set(i, j, (*this)(i, j) * scale);
874 75 : }
875 :
876 :
877 :
878 : template <typename T>
879 290 : Real SparseMatrix<T>::l1_norm_diff(const SparseMatrix<T> & other_mat) const
880 : {
881 290 : auto diff_mat = this->clone();
882 290 : diff_mat->add(-1.0, other_mat);
883 580 : return diff_mat->l1_norm();
884 266 : }
885 :
886 : //------------------------------------------------------------------
887 : // Explicit instantiations
888 : template class LIBMESH_EXPORT SparseMatrix<Number>;
889 :
890 : } // namespace libMesh
|
