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 : #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 103899 : PetscMatrixBase<T>::PetscMatrixBase(const Parallel::Communicator & comm_in) :
40 : SparseMatrix<T>(comm_in),
41 98183 : _mat(nullptr),
42 103899 : _destroy_mat_on_exit(true)
43 103899 : {}
44 :
45 :
46 :
47 : // Constructor taking an existing Mat but not the responsibility
48 : // for destroying it
49 : template <typename T>
50 2108 : PetscMatrixBase<T>::PetscMatrixBase(Mat mat_in,
51 : const Parallel::Communicator & comm_in,
52 : const bool destroy_on_exit) :
53 : SparseMatrix<T>(comm_in),
54 2108 : _destroy_mat_on_exit(destroy_on_exit)
55 : {
56 2108 : this->_mat = mat_in;
57 2108 : this->_is_initialized = true;
58 :
59 2108 : this->set_context();
60 2108 : }
61 :
62 :
63 :
64 : // Destructor
65 : template <typename T>
66 106007 : PetscMatrixBase<T>::~PetscMatrixBase()
67 : {
68 106007 : this->clear();
69 106007 : }
70 :
71 :
72 :
73 : template <typename T>
74 145055 : void PetscMatrixBase<T>::clear () noexcept
75 : {
76 145055 : if ((this->initialized()) && (this->_destroy_mat_on_exit))
77 : {
78 1420 : 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 45142 : PetscErrorCode ierr = MatDestroy (&_mat);
83 : if (ierr)
84 : libmesh_warning("Warning: MatDestroy returned a non-zero error code which we ignored.");
85 :
86 45142 : this->_is_initialized = false;
87 : }
88 145055 : }
89 :
90 : template <typename T>
91 840 : void PetscMatrixBase<T>::set_destroy_mat_on_exit(bool destroy)
92 : {
93 840 : this->_destroy_mat_on_exit = destroy;
94 840 : }
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 46610 : void PetscMatrixBase<T>::set_context()
106 : {
107 1438 : libmesh_assert(this->_mat);
108 : PetscContainer container;
109 46610 : LibmeshPetscCall(PetscContainerCreate(this->comm().get(), &container));
110 46610 : LibmeshPetscCall(PetscContainerSetPointer(container, this));
111 46610 : LibmeshPetscCall(PetscObjectCompose((PetscObject)(Mat)this->_mat, "PetscMatrixCtx", (PetscObject)container));
112 46610 : LibmeshPetscCall(PetscContainerDestroy(&container));
113 46610 : }
114 :
115 : template <typename T>
116 152128 : PetscMatrixBase<T> * PetscMatrixBase<T>::get_context(Mat mat, const TIMPI::Communicator & comm)
117 : {
118 : void * ctx;
119 : PetscContainer container;
120 152128 : LibmeshPetscCall2(comm, PetscObjectQuery((PetscObject)mat, "PetscMatrixCtx", (PetscObject *)&container));
121 152128 : if (!container)
122 408 : return nullptr;
123 :
124 137974 : LibmeshPetscCall2(comm, PetscContainerGetPointer(container, &ctx));
125 3712 : libmesh_assert(ctx);
126 137974 : return static_cast<PetscMatrixBase<T> *>(ctx);
127 : }
128 :
129 : template <typename T>
130 34002 : numeric_index_type PetscMatrixBase<T>::m () const
131 : {
132 32710 : libmesh_assert (this->initialized());
133 :
134 34002 : PetscInt petsc_m=0, petsc_n=0;
135 :
136 34002 : LibmeshPetscCall(MatGetSize (this->_mat, &petsc_m, &petsc_n));
137 :
138 34002 : return static_cast<numeric_index_type>(petsc_m);
139 : }
140 :
141 : template <typename T>
142 420 : numeric_index_type PetscMatrixBase<T>::local_m () const
143 : {
144 12 : libmesh_assert (this->initialized());
145 :
146 420 : PetscInt m = 0;
147 :
148 420 : LibmeshPetscCall(MatGetLocalSize (this->_mat, &m, NULL));
149 :
150 420 : return static_cast<numeric_index_type>(m);
151 : }
152 :
153 : template <typename T>
154 34352 : numeric_index_type PetscMatrixBase<T>::n () const
155 : {
156 32720 : libmesh_assert (this->initialized());
157 :
158 34352 : PetscInt petsc_m=0, petsc_n=0;
159 :
160 34352 : LibmeshPetscCall(MatGetSize (this->_mat, &petsc_m, &petsc_n));
161 :
162 34352 : return static_cast<numeric_index_type>(petsc_n);
163 : }
164 :
165 : template <typename T>
166 140 : numeric_index_type PetscMatrixBase<T>::local_n () const
167 : {
168 4 : libmesh_assert (this->initialized());
169 :
170 140 : PetscInt n = 0;
171 :
172 140 : LibmeshPetscCall(MatGetLocalSize (this->_mat, NULL, &n));
173 :
174 140 : return static_cast<numeric_index_type>(n);
175 : }
176 :
177 : template <typename T>
178 14989 : numeric_index_type PetscMatrixBase<T>::row_start () const
179 : {
180 428 : libmesh_assert (this->initialized());
181 :
182 14989 : PetscInt start=0, stop=0;
183 :
184 14989 : LibmeshPetscCall(MatGetOwnershipRange(this->_mat, &start, &stop));
185 :
186 14989 : return static_cast<numeric_index_type>(start);
187 : }
188 :
189 : template <typename T>
190 3751 : numeric_index_type PetscMatrixBase<T>::row_stop () const
191 : {
192 162 : libmesh_assert (this->initialized());
193 :
194 3751 : PetscInt start=0, stop=0;
195 :
196 3751 : LibmeshPetscCall(MatGetOwnershipRange(this->_mat, &start, &stop));
197 :
198 3751 : return static_cast<numeric_index_type>(stop);
199 : }
200 :
201 : template <typename T>
202 560 : numeric_index_type PetscMatrixBase<T>::col_start () const
203 : {
204 16 : libmesh_assert (this->initialized());
205 :
206 560 : PetscInt start=0, stop=0;
207 :
208 560 : LibmeshPetscCall(MatGetOwnershipRangeColumn(this->_mat, &start, &stop));
209 :
210 560 : return static_cast<numeric_index_type>(start);
211 : }
212 :
213 : template <typename T>
214 70 : numeric_index_type PetscMatrixBase<T>::col_stop () const
215 : {
216 2 : libmesh_assert (this->initialized());
217 :
218 70 : PetscInt start=0, stop=0;
219 :
220 70 : LibmeshPetscCall(MatGetOwnershipRangeColumn(this->_mat, &start, &stop));
221 :
222 70 : return static_cast<numeric_index_type>(stop);
223 : }
224 :
225 : template <typename T>
226 1062126 : void PetscMatrixBase<T>::close ()
227 : {
228 28112 : 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 1062126 : MatAssemblyBeginEnd(this->comm(), this->_mat, MAT_FINAL_ASSEMBLY);
238 1062126 : }
239 :
240 : template <typename T>
241 2999200 : bool PetscMatrixBase<T>::closed() const
242 : {
243 102244 : libmesh_assert (this->initialized());
244 :
245 : PetscBool assembled;
246 :
247 2999200 : LibmeshPetscCall(MatAssembled(this->_mat, &assembled));
248 :
249 2999200 : 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
|
