Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : #include "libmesh/libmesh_config.h"
21 :
22 : #ifdef LIBMESH_HAVE_PETSC
23 :
24 : // Local includes
25 : #include "libmesh/petsc_matrix_base.h"
26 : #include "libmesh/petsc_matrix_shell_matrix.h"
27 : #include "libmesh/petsc_matrix.h"
28 :
29 : namespace libMesh
30 : {
31 :
32 :
33 : //-----------------------------------------------------------------------
34 : // PetscMatrixBase members
35 :
36 :
37 : // Constructor
38 : template <typename T>
39 10262 : PetscMatrixBase<T>::PetscMatrixBase(const Parallel::Communicator & comm_in) :
40 : SparseMatrix<T>(comm_in),
41 4258 : _mat(nullptr),
42 10262 : _destroy_mat_on_exit(true)
43 10262 : {}
44 :
45 :
46 :
47 : // Constructor taking an existing Mat but not the responsibility
48 : // for destroying it
49 : template <typename T>
50 92 : PetscMatrixBase<T>::PetscMatrixBase(Mat mat_in,
51 : const Parallel::Communicator & comm_in,
52 : const bool destroy_on_exit) :
53 : SparseMatrix<T>(comm_in),
54 92 : _destroy_mat_on_exit(destroy_on_exit)
55 : {
56 92 : this->_mat = mat_in;
57 92 : this->_is_initialized = true;
58 :
59 92 : this->set_context();
60 92 : }
61 :
62 :
63 :
64 : // Destructor
65 : template <typename T>
66 10354 : PetscMatrixBase<T>::~PetscMatrixBase()
67 : {
68 10354 : this->clear();
69 10354 : }
70 :
71 :
72 :
73 : template <typename T>
74 14408 : void PetscMatrixBase<T>::clear () noexcept
75 : {
76 14408 : if ((this->initialized()) && (this->_destroy_mat_on_exit))
77 : {
78 1662 : exceptionless_semiparallel_only();
79 :
80 : // If we encounter an error here, print a warning but otherwise
81 : // keep going since we may be recovering from an exception.
82 4942 : PetscErrorCode ierr = MatDestroy (&_mat);
83 : if (ierr)
84 : libmesh_warning("Warning: MatDestroy returned a non-zero error code which we ignored.");
85 :
86 4942 : this->_is_initialized = false;
87 : }
88 14408 : }
89 :
90 : template <typename T>
91 72 : void PetscMatrixBase<T>::set_destroy_mat_on_exit(bool destroy)
92 : {
93 72 : this->_destroy_mat_on_exit = destroy;
94 72 : }
95 :
96 :
97 : template <typename T>
98 0 : void PetscMatrixBase<T>::swap(PetscMatrixBase<T> & m_in)
99 : {
100 0 : std::swap(_mat, m_in._mat);
101 0 : std::swap(_destroy_mat_on_exit, m_in._destroy_mat_on_exit);
102 0 : }
103 :
104 : template <typename T>
105 4990 : void PetscMatrixBase<T>::set_context()
106 : {
107 1678 : libmesh_assert(this->_mat);
108 : PetscContainer container;
109 4990 : LibmeshPetscCall(PetscContainerCreate(this->comm().get(), &container));
110 4990 : LibmeshPetscCall(PetscContainerSetPointer(container, this));
111 4990 : LibmeshPetscCall(PetscObjectCompose((PetscObject)(Mat)this->_mat, "PetscMatrixCtx", (PetscObject)container));
112 4990 : LibmeshPetscCall(PetscContainerDestroy(&container));
113 4990 : }
114 :
115 : template <typename T>
116 14704 : PetscMatrixBase<T> * PetscMatrixBase<T>::get_context(Mat mat, const TIMPI::Communicator & comm)
117 : {
118 : void * ctx;
119 : PetscContainer container;
120 14704 : LibmeshPetscCall2(comm, PetscObjectQuery((PetscObject)mat, "PetscMatrixCtx", (PetscObject *)&container));
121 14704 : if (!container)
122 408 : return nullptr;
123 :
124 13478 : LibmeshPetscCall2(comm, PetscContainerGetPointer(container, &ctx));
125 3712 : libmesh_assert(ctx);
126 13478 : return static_cast<PetscMatrixBase<T> *>(ctx);
127 : }
128 :
129 : template <typename T>
130 32788 : numeric_index_type PetscMatrixBase<T>::m () const
131 : {
132 32710 : libmesh_assert (this->initialized());
133 :
134 32788 : PetscInt petsc_m=0, petsc_n=0;
135 :
136 32788 : LibmeshPetscCall(MatGetSize (this->_mat, &petsc_m, &petsc_n));
137 :
138 32788 : return static_cast<numeric_index_type>(petsc_m);
139 : }
140 :
141 : template <typename T>
142 36 : numeric_index_type PetscMatrixBase<T>::local_m () const
143 : {
144 12 : libmesh_assert (this->initialized());
145 :
146 36 : PetscInt m = 0;
147 :
148 36 : LibmeshPetscCall(MatGetLocalSize (this->_mat, &m, NULL));
149 :
150 36 : return static_cast<numeric_index_type>(m);
151 : }
152 :
153 : template <typename T>
154 32827 : numeric_index_type PetscMatrixBase<T>::n () const
155 : {
156 32720 : libmesh_assert (this->initialized());
157 :
158 32827 : PetscInt petsc_m=0, petsc_n=0;
159 :
160 32827 : LibmeshPetscCall(MatGetSize (this->_mat, &petsc_m, &petsc_n));
161 :
162 32827 : return static_cast<numeric_index_type>(petsc_n);
163 : }
164 :
165 : template <typename T>
166 12 : numeric_index_type PetscMatrixBase<T>::local_n () const
167 : {
168 4 : libmesh_assert (this->initialized());
169 :
170 12 : PetscInt n = 0;
171 :
172 12 : LibmeshPetscCall(MatGetLocalSize (this->_mat, NULL, &n));
173 :
174 12 : return static_cast<numeric_index_type>(n);
175 : }
176 :
177 : template <typename T>
178 1836 : numeric_index_type PetscMatrixBase<T>::row_start () const
179 : {
180 566 : libmesh_assert (this->initialized());
181 :
182 1836 : PetscInt start=0, stop=0;
183 :
184 1836 : LibmeshPetscCall(MatGetOwnershipRange(this->_mat, &start, &stop));
185 :
186 1836 : return static_cast<numeric_index_type>(start);
187 : }
188 :
189 : template <typename T>
190 504 : numeric_index_type PetscMatrixBase<T>::row_stop () const
191 : {
192 166 : libmesh_assert (this->initialized());
193 :
194 504 : PetscInt start=0, stop=0;
195 :
196 504 : LibmeshPetscCall(MatGetOwnershipRange(this->_mat, &start, &stop));
197 :
198 504 : return static_cast<numeric_index_type>(stop);
199 : }
200 :
201 : template <typename T>
202 78 : numeric_index_type PetscMatrixBase<T>::col_start () const
203 : {
204 24 : libmesh_assert (this->initialized());
205 :
206 78 : PetscInt start=0, stop=0;
207 :
208 78 : LibmeshPetscCall(MatGetOwnershipRangeColumn(this->_mat, &start, &stop));
209 :
210 78 : return static_cast<numeric_index_type>(start);
211 : }
212 :
213 : template <typename T>
214 6 : numeric_index_type PetscMatrixBase<T>::col_stop () const
215 : {
216 2 : libmesh_assert (this->initialized());
217 :
218 6 : PetscInt start=0, stop=0;
219 :
220 6 : LibmeshPetscCall(MatGetOwnershipRangeColumn(this->_mat, &start, &stop));
221 :
222 6 : return static_cast<numeric_index_type>(stop);
223 : }
224 :
225 : template <typename T>
226 92520 : void PetscMatrixBase<T>::close ()
227 : {
228 28672 : semiparallel_only();
229 :
230 : // BSK - 1/19/2004
231 : // strictly this check should be OK, but it seems to
232 : // fail on matrix-free matrices. Do they falsely
233 : // state they are assembled? Check with the developers...
234 : // if (this->closed())
235 : // return;
236 :
237 92520 : MatAssemblyBeginEnd(this->comm(), this->_mat, MAT_FINAL_ASSEMBLY);
238 92520 : }
239 :
240 : template <typename T>
241 272858 : bool PetscMatrixBase<T>::closed() const
242 : {
243 102378 : libmesh_assert (this->initialized());
244 :
245 : PetscBool assembled;
246 :
247 272858 : LibmeshPetscCall(MatAssembled(this->_mat, &assembled));
248 :
249 272858 : return (assembled == PETSC_TRUE);
250 : }
251 :
252 : //------------------------------------------------------------------
253 : // Explicit instantiations
254 : template class LIBMESH_EXPORT PetscMatrixBase<Number>;
255 :
256 : } // namespace libMesh
257 :
258 :
259 : #endif // #ifdef LIBMESH_HAVE_PETSC
|
