← Back to branch summary

~ubuntu-branches/ubuntu/wily/dolfin/wily-proposed

« back to all changes in this revision

Viewing changes to dolfin/la/PETScMatrix.cpp

  • Committer: Package Import Robot
  • Author(s): Johannes Ring
  • Date: 2014-09-22 14:35:34 UTC
  • mfrom: (1.1.17) (19.1.23 sid)
  • Revision ID: package-import@ubuntu.com-20140922143534-0yi89jyuqbgdxwm9
Tags: 1.4.0+dfsg-4
http://bugs.debian.org/760658
* debian/control: Disable libcgal-dev on i386, mipsel and sparc.
* debian/rules: Remove bad directives in pkg-config file dolfin.pc
  (closes: #760658).
* Remove debian/libdolfin-dev.lintian-overrides.

Show diffs side-by-side

added added

removed removed

Lines of Context:
dolfin/la/PETScMatrix.cpp
51
51
                              ("frobenius", NORM_FROBENIUS);
52
52
 
53
53
//-----------------------------------------------------------------------------
54
 
PETScMatrix::PETScMatrix(bool use_gpu) : _use_gpu(use_gpu)
 
54
PETScMatrix::PETScMatrix(bool use_gpu) : PETScBaseMatrix(NULL),
 
55
                                         _use_gpu(use_gpu)
55
56
{
56
57
#ifndef HAS_PETSC_CUSP
57
58
  if (use_gpu)
65
66
  // Do nothing else
66
67
}
67
68
//-----------------------------------------------------------------------------
68
 
PETScMatrix::PETScMatrix(boost::shared_ptr<Mat> A, bool use_gpu) :
69
 
  PETScBaseMatrix(A), _use_gpu(use_gpu)
70
 
{
71
 
#ifndef HAS_PETSC_CUSP
72
 
  if (use_gpu)
 
69
PETScMatrix::PETScMatrix(Mat A) : PETScBaseMatrix(A), _use_gpu(false)
 
70
{
 
71
  // Do nothing (reference count to A is incremented in base class)
 
72
}
 
73
//-----------------------------------------------------------------------------
 
74
PETScMatrix::PETScMatrix(const PETScMatrix& A) : PETScBaseMatrix(NULL),
 
75
                                                 _use_gpu(false)
 
76
{
 
77
  if (A.mat())
73
78
  {
74
 
    dolfin_error("PETScMatrix.cpp",
75
 
                 "create GPU matrix",
76
 
                 "PETSc not compiled with Cusp support");
 
79
    PetscErrorCode ierr = MatDuplicate(A.mat(), MAT_COPY_VALUES, &_matA);
 
80
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatDuplicate");
77
81
  }
78
 
#endif
79
 
}
80
 
//-----------------------------------------------------------------------------
81
 
PETScMatrix::PETScMatrix(const PETScMatrix& A): _use_gpu(false)
82
 
{
83
 
  *this = A;
84
82
}
85
83
//-----------------------------------------------------------------------------
86
84
PETScMatrix::~PETScMatrix()
87
85
{
88
 
  // Do nothing
 
86
  // Do nothing (PETSc matrix is destroyed in base class)
89
87
}
90
88
//-----------------------------------------------------------------------------
91
 
boost::shared_ptr<GenericMatrix> PETScMatrix::copy() const
 
89
std::shared_ptr<GenericMatrix> PETScMatrix::copy() const
92
90
{
93
 
  boost::shared_ptr<GenericMatrix> B;
94
 
  if (!_A)
95
 
    B.reset(new PETScMatrix());
96
 
  else
97
 
  {
98
 
    // Create copy of PETSc matrix
99
 
    boost::shared_ptr<Mat> _Acopy(new Mat, PETScMatrixDeleter());
100
 
    PetscErrorCode ierr = MatDuplicate(*_A, MAT_COPY_VALUES, _Acopy.get());
101
 
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatDuplicate");
102
 
 
103
 
    // Create PETScMatrix
104
 
    B.reset(new PETScMatrix(_Acopy));
105
 
  }
106
 
 
107
 
  return B;
 
91
  return std::shared_ptr<GenericMatrix>(new PETScMatrix(*this));
108
92
}
109
93
//-----------------------------------------------------------------------------
110
94
void PETScMatrix::init(const TensorLayout& tensor_layout)
121
105
    = tensor_layout.local_range(1);
122
106
  const std::size_t m = row_range.second - row_range.first;
123
107
  const std::size_t n = col_range.second - col_range.first;
124
 
  dolfin_assert(M > 0 && N > 0 && m > 0 && n > 0);
125
108
 
126
109
  // Get sparsity pattern
127
110
  dolfin_assert(tensor_layout.sparsity_pattern());
128
111
  const GenericSparsityPattern& sparsity_pattern
129
112
    = *tensor_layout.sparsity_pattern();
130
113
 
131
 
  // Create matrix (any old matrix is destroyed automatically)
132
 
  if (_A && !_A.unique())
 
114
  if (_matA)
133
115
  {
134
 
    dolfin_error("PETScMatrix.cpp",
135
 
                 "initialize PETSc matrix",
136
 
                 "More than one object points to the underlying PETSc object");
 
116
    #ifdef DOLFIN_DEPRECATION_ERROR
 
117
    error("PETScMatrix may not be initialized more than once. Remove build definiton -DDOLFIN_DEPRECATION_ERROR to change this to a warning.");
 
118
    #else
 
119
    warning("PETScMatrix may not be initialized more than once. In version > 1.4, this will become an error.");
 
120
    #endif
 
121
    MatDestroy(&_matA);
137
122
  }
138
 
  _A.reset(new Mat, PETScMatrixDeleter());
139
123
 
140
124
  // Initialize matrix
141
125
  if (row_range.first == 0 && row_range.second == M)
145
129
    std::vector<std::size_t> num_nonzeros(M);
146
130
    sparsity_pattern.num_nonzeros_diagonal(num_nonzeros);
147
131
 
148
 
    // FIXME: Debugging
149
 
    //cout << "Number of nonzeros for PETSc matrix:" << endl;
150
 
    //for (unsigned int i = 0; i < M; i++)
151
 
    //  cout << "i = " << i << ": " << num_nonzeros[i] << endl;
152
 
 
153
132
    // Create matrix
154
 
    ierr = MatCreate(PETSC_COMM_SELF, _A.get());
 
133
    ierr = MatCreate(PETSC_COMM_SELF, &_matA);
155
134
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatCreate");
156
135
 
157
136
    // Set size
158
 
    ierr = MatSetSizes(*_A, M, N, M, N);
 
137
    ierr = MatSetSizes(_matA, M, N, M, N);
159
138
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetSizes");
160
139
 
161
140
    // Set matrix type according to chosen architecture
162
141
    if (!_use_gpu)
163
142
    {
164
 
      ierr = MatSetType(*_A, MATSEQAIJ);
 
143
      ierr = MatSetType(_matA, MATSEQAIJ);
165
144
      if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetType");
166
145
    }
167
146
    #ifdef HAS_PETSC_CUSP
168
147
    else
169
148
    {
170
 
      ierr = MatSetType(*_A, MATSEQAIJCUSP);
 
149
      ierr = MatSetType(_matA, MATSEQAIJCUSP);
171
150
      if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetType");
172
151
    }
173
152
    #endif
175
154
    // Set block size
176
155
    if (tensor_layout.block_size > 1)
177
156
    {
178
 
     ierr =  MatSetBlockSize(*_A, tensor_layout.block_size);
 
157
     ierr =  MatSetBlockSize(_matA, tensor_layout.block_size);
179
158
     if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetBlockSize");
180
159
    }
181
160
 
186
165
    // Copy number of non-zeros to PetscInt type
187
166
    const std::vector<PetscInt> _num_nonzeros(num_nonzeros.begin(),
188
167
                                              num_nonzeros.end());
189
 
    ierr = MatSeqAIJSetPreallocation(*_A, 0, _num_nonzeros.data());
 
168
    ierr = MatSeqAIJSetPreallocation(_matA, 0, _num_nonzeros.data());
190
169
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatSeqAIJSetPreallocation");
191
170
 
192
171
    // Set column indices
202
181
      //  cout << "  Col: " << _column_indices[i][j] << endl;
203
182
      column_indices.insert(column_indices.end(), _column_indices[i].begin(), _column_indices[i].end());
204
183
    }
205
 
    MatSeqAIJSetColumnIndices(*_A, &column_indices[0]);
 
184
    MatSeqAIJSetColumnIndices(_matA, &column_indices[0]);
206
185
    */
207
186
  }
208
187
  else
215
194
 
216
195
    // FIXME: Try using MatStashSetInitialSize to optimise performance
217
196
 
218
 
    //info("Initializing parallel PETSc matrix (MPIAIJ) of size %d x %d.", M, N);
219
 
    //info("Local range is [%d, %d] x [%d, %d].",
220
 
    //     row_range.first, row_range.second, col_range.first, col_range.second);
221
 
 
222
197
    // Get number of nonzeros for each row from sparsity pattern
223
198
    std::vector<std::size_t> num_nonzeros_diagonal;
224
199
    std::vector<std::size_t> num_nonzeros_off_diagonal;
226
201
    sparsity_pattern.num_nonzeros_off_diagonal(num_nonzeros_off_diagonal);
227
202
 
228
203
    // Create matrix
229
 
    ierr = MatCreate(PETSC_COMM_WORLD, _A.get());
 
204
    ierr = MatCreate(PETSC_COMM_WORLD, &_matA);
230
205
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatCreate");
231
206
 
232
207
    // Set size
233
 
    ierr = MatSetSizes(*_A, m, n, M, N);
 
208
    ierr = MatSetSizes(_matA, m, n, M, N);
234
209
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetSizes");
235
210
 
236
211
    // Set matrix type
237
 
    ierr = MatSetType(*_A, MATMPIAIJ);
 
212
    ierr = MatSetType(_matA, MATMPIAIJ);
238
213
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetType");
239
214
 
240
215
    // Set block size
241
216
    if (tensor_layout.block_size > 1)
242
217
    {
243
 
      ierr = MatSetBlockSize(*_A, tensor_layout.block_size);
 
218
      ierr = MatSetBlockSize(_matA, tensor_layout.block_size);
244
219
      if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetBlockSize");
245
220
    }
246
221
    // Allocate space (using data from sparsity pattern)
250
225
    const std::vector<PetscInt>
251
226
      _num_nonzeros_off_diagonal(num_nonzeros_off_diagonal.begin(),
252
227
                                 num_nonzeros_off_diagonal.end());
253
 
    ierr = MatMPIAIJSetPreallocation(*_A, 0, _num_nonzeros_diagonal.data(),
 
228
    ierr = MatMPIAIJSetPreallocation(_matA, 0, _num_nonzeros_diagonal.data(),
254
229
                                     0, _num_nonzeros_off_diagonal.data());
255
230
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatMPIAIJSetPreallocation");
256
231
  }
258
233
  // Set some options
259
234
 
260
235
  // Do not allow more entries than have been pre-allocated
261
 
  ierr = MatSetOption(*_A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
262
 
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetOption");
263
 
 
264
 
  ierr = MatSetOption(*_A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
265
 
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetOption");
266
 
 
267
 
  ierr = MatSetFromOptions(*_A);
 
236
  ierr = MatSetOption(_matA, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
 
237
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetOption");
 
238
 
 
239
  ierr = MatSetOption(_matA, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
 
240
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetOption");
 
241
 
 
242
  ierr = MatSetFromOptions(_matA);
268
243
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetFromOptions");
269
244
 
270
245
  #if PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR > 2
271
 
  ierr = MatSetUp(*_A.get());
 
246
  ierr = MatSetUp(_matA);
272
247
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetUp");
273
248
  #endif
274
249
}
275
250
//-----------------------------------------------------------------------------
 
251
bool PETScMatrix::empty() const
 
252
{
 
253
  return _matA ? false : true;
 
254
}
 
255
//-----------------------------------------------------------------------------
276
256
void PETScMatrix::get(double* block,
277
257
                      std::size_t m, const dolfin::la_index* rows,
278
258
                      std::size_t n, const dolfin::la_index* cols) const
279
259
{
280
260
  // Get matrix entries (must be on this process)
281
 
  dolfin_assert(_A);
282
 
  PetscErrorCode ierr = MatGetValues(*_A, m, rows, n, cols, block);
 
261
  dolfin_assert(_matA);
 
262
  PetscErrorCode ierr = MatGetValues(_matA, m, rows, n, cols, block);
283
263
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatGetValues");
284
264
}
285
265
//-----------------------------------------------------------------------------
287
267
                      std::size_t m, const dolfin::la_index* rows,
288
268
                      std::size_t n, const dolfin::la_index* cols)
289
269
{
290
 
  dolfin_assert(_A);
291
 
  PetscErrorCode ierr = MatSetValues(*_A, m, rows, n, cols, block,
 
270
  dolfin_assert(_matA);
 
271
  PetscErrorCode ierr = MatSetValues(_matA, m, rows, n, cols, block,
292
272
                                    INSERT_VALUES);
293
273
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetValues");
294
274
}
297
277
                      std::size_t m, const dolfin::la_index* rows,
298
278
                      std::size_t n, const dolfin::la_index* cols)
299
279
{
300
 
  dolfin_assert(_A);
301
 
  PetscErrorCode ierr = MatSetValues(*_A, m, rows, n, cols, block, ADD_VALUES);
 
280
  dolfin_assert(_matA);
 
281
  PetscErrorCode ierr = MatSetValues(_matA, m, rows, n, cols, block, ADD_VALUES);
302
282
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatSetValues");
303
283
}
304
284
//-----------------------------------------------------------------------------
308
288
  PetscErrorCode ierr;
309
289
 
310
290
  const PETScMatrix* AA = &as_type<const PETScMatrix>(A);
311
 
  dolfin_assert(_A);
 
291
  dolfin_assert(_matA);
312
292
  dolfin_assert(AA->mat());
313
293
  if (same_nonzero_pattern)
314
294
  {
315
 
    ierr = MatAXPY(*_A, a, *AA->mat(), SAME_NONZERO_PATTERN);
 
295
    ierr = MatAXPY(_matA, a, AA->mat(), SAME_NONZERO_PATTERN);
316
296
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatAXPY");
317
297
  }
318
298
  else
319
299
  {
320
 
    ierr = MatAXPY(*_A, a, *AA->mat(), DIFFERENT_NONZERO_PATTERN);
 
300
    ierr = MatAXPY(_matA, a, AA->mat(), DIFFERENT_NONZERO_PATTERN);
321
301
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatAXPY");
322
302
  }
323
303
}
325
305
void PETScMatrix::getrow(std::size_t row, std::vector<std::size_t>& columns,
326
306
                         std::vector<double>& values) const
327
307
{
328
 
  dolfin_assert(_A);
 
308
  dolfin_assert(_matA);
329
309
 
330
310
  PetscErrorCode ierr;
331
311
  const PetscInt *cols = 0;
332
312
  const double *vals = 0;
333
313
  PetscInt ncols = 0;
334
 
  ierr = MatGetRow(*_A, row, &ncols, &cols, &vals);
 
314
  ierr = MatGetRow(_matA, row, &ncols, &cols, &vals);
335
315
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatGetRow");
336
316
 
337
317
  // Assign values to std::vectors
338
318
  columns.assign(cols, cols + ncols);
339
319
  values.assign(vals, vals + ncols);
340
320
 
341
 
  ierr = MatRestoreRow(*_A, row, &ncols, &cols, &vals);
 
321
  ierr = MatRestoreRow(_matA, row, &ncols, &cols, &vals);
342
322
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatRestorRow");
343
323
}
344
324
//-----------------------------------------------------------------------------
346
326
                         const std::vector<std::size_t>& columns,
347
327
                         const std::vector<double>& values)
348
328
{
349
 
  dolfin_assert(_A);
 
329
  dolfin_assert(_matA);
350
330
 
351
331
  // Check size of arrays
352
332
  if (columns.size() != values.size())
369
349
//-----------------------------------------------------------------------------
370
350
void PETScMatrix::zero(std::size_t m, const dolfin::la_index* rows)
371
351
{
372
 
  dolfin_assert(_A);
 
352
  dolfin_assert(_matA);
373
353
 
374
354
  PetscErrorCode ierr;
375
355
  IS is = 0;
377
357
  const PetscInt _m = m;
378
358
  ierr = ISCreateGeneral(PETSC_COMM_SELF, _m, rows, PETSC_COPY_VALUES, &is);
379
359
  if (ierr != 0) petsc_error(ierr, __FILE__, "ISCreateGeneral");
380
 
  ierr = MatZeroRowsIS(*_A, is, null, NULL, NULL);
 
360
  ierr = MatZeroRowsIS(_matA, is, null, NULL, NULL);
381
361
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatZeroRowsIS");
382
362
  ierr = ISDestroy(&is);
383
363
  if (ierr != 0) petsc_error(ierr, __FILE__, "ISDestroy");
385
365
//-----------------------------------------------------------------------------
386
366
void PETScMatrix::ident(std::size_t m, const dolfin::la_index* rows)
387
367
{
388
 
  dolfin_assert(_A);
 
368
  dolfin_assert(_matA);
389
369
 
390
370
  PetscErrorCode ierr;
391
371
  IS is = 0;
393
373
  const PetscInt _m = m;
394
374
  ierr = ISCreateGeneral(PETSC_COMM_SELF, _m, rows, PETSC_COPY_VALUES, &is);
395
375
  if (ierr != 0) petsc_error(ierr, __FILE__, "ISCreateGeneral");
396
 
  ierr = MatZeroRowsIS(*_A, is, one, NULL, NULL);
 
376
  ierr = MatZeroRowsIS(_matA, is, one, NULL, NULL);
397
377
  if (ierr == PETSC_ERR_ARG_WRONGSTATE)
398
378
  {
399
379
    dolfin_error("PETScMatrix.cpp",
409
389
//-----------------------------------------------------------------------------
410
390
void PETScMatrix::mult(const GenericVector& x, GenericVector& y) const
411
391
{
412
 
  dolfin_assert(_A);
 
392
  dolfin_assert(_matA);
413
393
 
414
394
  const PETScVector& xx = as_type<const PETScVector>(x);
415
395
  PETScVector& yy = as_type<PETScVector>(y);
423
403
 
424
404
  // Resize RHS if empty
425
405
  if (yy.size() == 0)
426
 
    resize(yy, 0);
 
406
    init_vector(yy, 0);
427
407
 
428
408
  if (size(0) != yy.size())
429
409
  {
432
412
                 "Vector for matrix-vector result has wrong size");
433
413
  }
434
414
 
435
 
  PetscErrorCode ierr = MatMult(*_A, *xx.vec(), *yy.vec());
 
415
  PetscErrorCode ierr = MatMult(_matA, xx.vec(), yy.vec());
436
416
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatMult");
437
417
}
438
418
//-----------------------------------------------------------------------------
439
419
void PETScMatrix::transpmult(const GenericVector& x, GenericVector& y) const
440
420
{
441
 
  dolfin_assert(_A);
 
421
  dolfin_assert(_matA);
442
422
 
443
423
  const PETScVector& xx = as_type<const PETScVector>(x);
444
424
  PETScVector& yy = as_type<PETScVector>(y);
452
432
 
453
433
  // Resize RHS if empty
454
434
  if (yy.size() == 0)
455
 
    resize(yy, 1);
 
435
    init_vector(yy, 1);
456
436
 
457
437
  if (size(1) != yy.size())
458
438
  {
461
441
                 "Vector for transpose matrix-vector result has wrong size");
462
442
  }
463
443
 
464
 
  PetscErrorCode ierr = MatMultTranspose(*_A, *xx.vec(), *yy.vec());
 
444
  PetscErrorCode ierr = MatMultTranspose(_matA, xx.vec(), yy.vec());
465
445
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatMultTranspose");
466
446
}
467
447
//-----------------------------------------------------------------------------
 
448
void PETScMatrix::set_diagonal(const GenericVector& x)
 
449
{
 
450
  dolfin_assert(_matA);
 
451
 
 
452
  const PETScVector& xx = x.down_cast<PETScVector>();
 
453
  if (size(1) != size(0) || size(0) != xx.size())
 
454
  {
 
455
    dolfin_error("PETScMatrix.cpp",
 
456
                 "set diagonal of a PETSc matrix",
 
457
                 "Matrix and vector dimensions don't match for matrix-vector set");
 
458
  }
 
459
 
 
460
  PetscErrorCode ierr = MatDiagonalSet(_matA, xx.vec(), INSERT_VALUES);
 
461
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatDiagonalSet");
 
462
  apply("insert");
 
463
}
 
464
//-----------------------------------------------------------------------------
468
465
double PETScMatrix::norm(std::string norm_type) const
469
466
{
470
 
  dolfin_assert(_A);
 
467
  dolfin_assert(_matA);
471
468
 
472
469
  // Check that norm is known
473
470
  if (norm_types.count(norm_type) == 0)
478
475
  }
479
476
 
480
477
  double value = 0.0;
481
 
  PetscErrorCode ierr = MatNorm(*_A, norm_types.find(norm_type)->second,
 
478
  PetscErrorCode ierr = MatNorm(_matA, norm_types.find(norm_type)->second,
482
479
                                &value);
483
480
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatNorm");
484
481
  return value;
488
485
{
489
486
  Timer timer("Apply (PETScMatrix)");
490
487
 
491
 
  dolfin_assert(_A);
 
488
  dolfin_assert(_matA);
492
489
  PetscErrorCode ierr;
493
490
  if (mode == "add")
494
491
  {
495
 
    ierr = MatAssemblyBegin(*_A, MAT_FINAL_ASSEMBLY);
 
492
    ierr = MatAssemblyBegin(_matA, MAT_FINAL_ASSEMBLY);
496
493
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatAssemblyBegin");
497
 
    ierr = MatAssemblyEnd(*_A, MAT_FINAL_ASSEMBLY);
 
494
    ierr = MatAssemblyEnd(_matA, MAT_FINAL_ASSEMBLY);
498
495
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatAssemblyEnd");
499
496
  }
500
497
  else if (mode == "insert")
501
498
  {
502
 
    ierr = MatAssemblyBegin(*_A, MAT_FINAL_ASSEMBLY);
 
499
    ierr = MatAssemblyBegin(_matA, MAT_FINAL_ASSEMBLY);
503
500
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatAssemblyBegin");
504
 
    ierr = MatAssemblyEnd(*_A, MAT_FINAL_ASSEMBLY);
 
501
    ierr = MatAssemblyEnd(_matA, MAT_FINAL_ASSEMBLY);
505
502
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatAssemblyEnd");
506
503
  }
507
504
  else if (mode == "flush")
508
505
  {
509
 
    ierr = MatAssemblyBegin(*_A, MAT_FLUSH_ASSEMBLY);
 
506
    ierr = MatAssemblyBegin(_matA, MAT_FLUSH_ASSEMBLY);
510
507
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatAssemblyBegin");
511
 
    ierr = MatAssemblyEnd(*_A, MAT_FLUSH_ASSEMBLY);
 
508
    ierr = MatAssemblyEnd(_matA, MAT_FLUSH_ASSEMBLY);
512
509
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatAssemblyEnd");
513
510
  }
514
511
  else
519
516
  }
520
517
}
521
518
//-----------------------------------------------------------------------------
 
519
MPI_Comm PETScMatrix::mpi_comm() const
 
520
{
 
521
  dolfin_assert(_matA);
 
522
  MPI_Comm mpi_comm = MPI_COMM_NULL;
 
523
  PetscObjectGetComm((PetscObject)_matA, &mpi_comm);
 
524
  return mpi_comm;
 
525
}
 
526
//-----------------------------------------------------------------------------
522
527
void PETScMatrix::zero()
523
528
{
524
 
  dolfin_assert(_A);
525
 
  PetscErrorCode ierr = MatZeroEntries(*_A);
 
529
  dolfin_assert(_matA);
 
530
  PetscErrorCode ierr = MatZeroEntries(_matA);
526
531
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatZeroEntries");
527
532
}
528
533
//-----------------------------------------------------------------------------
529
534
const PETScMatrix& PETScMatrix::operator*= (double a)
530
535
{
531
 
  dolfin_assert(_A);
532
 
  PetscErrorCode ierr = MatScale(*_A, a);
 
536
  dolfin_assert(_matA);
 
537
  PetscErrorCode ierr = MatScale(_matA, a);
533
538
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatScale");
534
539
  return *this;
535
540
}
536
541
//-----------------------------------------------------------------------------
537
542
const PETScMatrix& PETScMatrix::operator/= (double a)
538
543
{
539
 
  dolfin_assert(_A);
540
 
  MatScale(*_A, 1.0/a);
 
544
  dolfin_assert(_matA);
 
545
  MatScale(_matA, 1.0/a);
541
546
  return *this;
542
547
}
543
548
//-----------------------------------------------------------------------------
549
554
//-----------------------------------------------------------------------------
550
555
bool PETScMatrix::is_symmetric(double tol) const
551
556
{
552
 
  dolfin_assert(_A);
 
557
  dolfin_assert(_matA);
553
558
  PetscBool symmetric = PETSC_FALSE;
554
 
  PetscErrorCode ierr = MatIsSymmetric(*_A, tol, &symmetric);
 
559
  PetscErrorCode ierr = MatIsSymmetric(_matA, tol, &symmetric);
555
560
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatIsSymmetric");
556
561
  return symmetric == PETSC_TRUE ? true : false;
557
562
}
565
570
    return PETScCuspFactory::instance();
566
571
  #endif
567
572
 
568
 
  // Return something to keep the compiler happy. Code will never be reached.
 
573
  // Return something to keep the compiler happy. Code will never be
 
574
  // reached.
569
575
  return PETScFactory::instance();
570
576
}
571
577
//-----------------------------------------------------------------------------
572
578
const PETScMatrix& PETScMatrix::operator= (const PETScMatrix& A)
573
579
{
574
580
  if (!A.mat())
575
 
    _A.reset();
 
581
  {
 
582
    if (_matA)
 
583
    {
 
584
      #ifdef DOLFIN_DEPRECATION_ERROR
 
585
      error("PETScVector may not be initialized more than once. Remove build definiton -DDOLFIN_DEPRECATION_ERROR to change this to a warning. Error is in PETScMatrix::operator=.");
 
586
      #else
 
587
      warning("PETScVector may not be initialized more than once. In version > 1.4, this will become an error. Warning is in PETScMatrix::operator=.");
 
588
      #endif
 
589
      MatDestroy(&_matA);
 
590
    }
 
591
    _matA = NULL;
 
592
  }
576
593
  else if (this != &A) // Check for self-assignment
577
594
  {
578
 
    if (_A && !_A.unique())
 
595
    if (_matA)
579
596
    {
580
 
      dolfin_error("PETScMatrix.cpp",
581
 
                 "assign to PETSc matrix",
582
 
                 "More than one object points to the underlying PETSc object");
 
597
      // Get reference count to _matA
 
598
      PetscInt ref_count = 0;
 
599
      PetscObjectGetReference((PetscObject)_matA, &ref_count);
 
600
      if (ref_count > 1)
 
601
      {
 
602
        dolfin_error("PETScMatrix.cpp",
 
603
                     "assign to PETSc matrix",
 
604
                     "More than one object points to the underlying PETSc object");
 
605
      }
 
606
      #ifdef DOLFIN_DEPRECATION_ERROR
 
607
      error("PETScMatrix may not be initialized more than once. Remove build definiton -DDOLFIN_DEPRECATION_ERROR to change this to a warning. Error is in PETScMatrix::operator=.");
 
608
      #else
 
609
      warning("PETScMatrix may not be initialized more than once. In version > 1.4, this will become an error. Warning is in PETScMatrix::operator=.");
 
610
      #endif
 
611
      MatDestroy(&_matA);
583
612
    }
584
 
    _A.reset(new Mat, PETScMatrixDeleter());
585
613
 
586
614
    // Duplicate with the same pattern as A.A
587
 
    PetscErrorCode ierr = MatDuplicate(*A.mat(), MAT_COPY_VALUES, _A.get());
 
615
    PetscErrorCode ierr = MatDuplicate(A.mat(), MAT_COPY_VALUES, &_matA);
588
616
    if (ierr != 0) petsc_error(ierr, __FILE__, "MatDuplicate");
589
617
  }
590
618
  return *this;
599
627
                               FILE_MODE_WRITE, &view_out);
600
628
  if (ierr != 0) petsc_error(ierr, __FILE__, "PetscViewerBinaryOpen");
601
629
 
602
 
  ierr = MatView(*(_A.get()), view_out);
 
630
  ierr = MatView(_matA, view_out);
603
631
  if (ierr != 0) petsc_error(ierr, __FILE__, "MatView");
604
632
 
605
633
  ierr = PetscViewerDestroy(&view_out);
608
636
//-----------------------------------------------------------------------------
609
637
std::string PETScMatrix::str(bool verbose) const
610
638
{
611
 
  if (!_A)
 
639
  if (!_matA)
612
640
    return "<Uninitialized PETScMatrix>";
613
641
 
614
642
  std::stringstream s;
617
645
    warning("Verbose output for PETScMatrix not implemented, calling PETSc MatView directly.");
618
646
 
619
647
    // FIXME: Maybe this could be an option?
620
 
    dolfin_assert(_A);
 
648
    dolfin_assert(_matA);
621
649
    PetscErrorCode ierr;
622
 
    if (MPI::num_processes() > 1)
 
650
    if (MPI::size(MPI_COMM_WORLD) > 1)
623
651
    {
624
 
      ierr = MatView(*_A, PETSC_VIEWER_STDOUT_WORLD);
 
652
      ierr = MatView(_matA, PETSC_VIEWER_STDOUT_WORLD);
625
653
      if (ierr != 0) petsc_error(ierr, __FILE__, "MatView");
626
654
    }
627
655
    else
628
656
    {
629
 
      ierr = MatView(*_A, PETSC_VIEWER_STDOUT_SELF);
 
657
      ierr = MatView(_matA, PETSC_VIEWER_STDOUT_SELF);
630
658
      if (ierr != 0) petsc_error(ierr, __FILE__, "MatView");
631
659
    }
632
660
  }