← Back to branch summary

~corrado-maurini/dolfin/tao

« back to all changes in this revision

Viewing changes to dolfin/la/EpetraMatrix.cpp

  • Committer: corrado maurini
  • Date: 2012-12-18 12:16:08 UTC
  • mfrom: (6685.78.207 trunk)
  • Revision ID: corrado.maurini@upmc.fr-20121218121608-nk82ly9jgsld9u84
updating with trunk, fix uint in TAO solver and hacking the check for tao FindTAO.cmake

Show diffs side-by-side

added added

removed removed

Lines of Context:
dolfin/la/EpetraMatrix.cpp
101
101
  }
102
102
 
103
103
  // Get local range
104
 
  const std::pair<uint, uint> range = tensor_layout.local_range(0);
105
 
  const uint num_local_rows = range.second - range.first;
106
 
  const uint n0 = range.first;
 
104
  const std::pair<std::size_t, std::size_t> range = tensor_layout.local_range(0);
 
105
  const std::size_t num_local_rows = range.second - range.first;
 
106
  const std::size_t n0 = range.first;
107
107
 
108
108
  dolfin_assert(tensor_layout.sparsity_pattern());
109
 
  const SparsityPattern& _pattern = dynamic_cast<const SparsityPattern&>(*tensor_layout.sparsity_pattern());
110
 
  const std::vector<std::vector<dolfin::uint> > d_pattern = _pattern.diagonal_pattern(SparsityPattern::unsorted);
111
 
  const std::vector<std::vector<dolfin::uint> > o_pattern = _pattern.off_diagonal_pattern(SparsityPattern::unsorted);
 
109
  const SparsityPattern& _pattern
 
110
    = dynamic_cast<const SparsityPattern&>(*tensor_layout.sparsity_pattern());
 
111
  const std::vector<std::vector<std::size_t> > d_pattern
 
112
    = _pattern.diagonal_pattern(SparsityPattern::unsorted);
 
113
  const std::vector<std::vector<std::size_t> > o_pattern
 
114
      = _pattern.off_diagonal_pattern(SparsityPattern::unsorted);
112
115
 
113
116
  // Get number of non-zeroes per row
114
 
  std::vector<uint> num_nonzeros;
 
117
  std::vector<std::size_t> num_nonzeros;
115
118
  _pattern.num_local_nonzeros(num_nonzeros);
116
119
 
117
120
  // Create row map
118
121
  EpetraFactory& f = EpetraFactory::instance();
119
122
  Epetra_MpiComm comm = f.get_mpi_comm();
120
 
  Epetra_Map row_map((int) tensor_layout.size(0), (int) num_local_rows, 0, comm);
 
123
  Epetra_Map row_map((DolfinIndex) tensor_layout.size(0), (DolfinIndex) num_local_rows, 0, comm);
121
124
 
122
125
  // For rectangular matrices with more columns than rows, the columns which are
123
126
  // larger than those in row_map are marked as nonlocal (and assembly fails).
124
127
  // The domain_map fixes that problem, at least in the serial case.
125
128
  // FIXME: Needs attention in the parallel case. Maybe range_map is also req'd.
126
 
  const std::pair<uint, uint> colrange = tensor_layout.local_range(1);
 
129
  const std::pair<std::size_t, std::size_t> colrange = tensor_layout.local_range(1);
127
130
  const int num_local_cols = colrange.second - colrange.first;
128
 
  Epetra_Map domain_map((int) tensor_layout.size(1), (int) num_local_cols, 0, comm);
 
131
  Epetra_Map domain_map((DolfinIndex) tensor_layout.size(1), num_local_cols, 0, comm);
129
132
 
130
133
  // Create Epetra_FECrsGraph
131
 
  Epetra_CrsGraph matrix_map(Copy, row_map, reinterpret_cast<int*>(&num_nonzeros[0]));
 
134
  const std::vector<int> _num_nonzeros(num_nonzeros.begin(), num_nonzeros.end());
 
135
  Epetra_CrsGraph matrix_map(Copy, row_map, _num_nonzeros.data());
132
136
 
133
137
  // Add diagonal block indices
134
 
  for (uint local_row = 0; local_row < d_pattern.size(); local_row++)
 
138
  for (std::size_t local_row = 0; local_row < d_pattern.size(); local_row++)
135
139
  {
136
 
    const uint global_row = local_row + n0;
137
 
    std::vector<uint>& entries = const_cast<std::vector<uint>&>(d_pattern[local_row]);
138
 
    matrix_map.InsertGlobalIndices(global_row, entries.size(),
139
 
                                   reinterpret_cast<int*>(&entries[0]));
 
140
    const DolfinIndex global_row = local_row + n0;
 
141
    std::vector<DolfinIndex> entries(d_pattern[local_row].begin(), d_pattern[local_row].end());
 
142
    matrix_map.InsertGlobalIndices(global_row, (DolfinIndex) entries.size(), entries.data());
140
143
  }
141
144
 
142
145
  // Add off-diagonal block indices (parallel only)
143
 
  for (uint local_row = 0; local_row < o_pattern.size(); local_row++)
 
146
  for (std::size_t local_row = 0; local_row < o_pattern.size(); local_row++)
144
147
  {
145
 
    const uint global_row = local_row + n0;
146
 
    std::vector<uint>& entries = const_cast<std::vector<uint>&>(o_pattern[local_row]);
147
 
    matrix_map.InsertGlobalIndices(global_row, entries.size(),
148
 
                                   reinterpret_cast<int*>(&entries[0]));
 
148
    const std::size_t global_row = local_row + n0;
 
149
    std::vector<DolfinIndex> entries(o_pattern[local_row].begin(),o_pattern[local_row].end());
 
150
    matrix_map.InsertGlobalIndices(global_row, entries.size(), entries.data());
149
151
  }
150
152
 
151
153
  try
174
176
  return B;
175
177
}
176
178
//-----------------------------------------------------------------------------
177
 
dolfin::uint EpetraMatrix::size(uint dim) const
 
179
std::size_t EpetraMatrix::size(std::size_t dim) const
178
180
{
179
181
  if (dim > 1)
180
182
  {
185
187
 
186
188
  if (A)
187
189
  {
188
 
    const int M = A->NumGlobalRows();
189
 
    const int N = A->NumGlobalCols();
 
190
    const std::size_t M = A->NumGlobalRows64();
 
191
    const std::size_t N = A->NumGlobalCols64();
190
192
    return (dim == 0 ? M : N);
191
193
  }
192
194
  else
193
195
    return 0;
194
196
}
195
197
//-----------------------------------------------------------------------------
196
 
std::pair<dolfin::uint, dolfin::uint> EpetraMatrix::local_range(uint dim) const
 
198
std::pair<std::size_t, std::size_t> EpetraMatrix::local_range(std::size_t dim) const
197
199
{
198
200
  dolfin_assert(dim < 2);
199
201
  if (dim == 1)
206
208
  const Epetra_BlockMap& row_map = A->RowMap();
207
209
  dolfin_assert(row_map.LinearMap());
208
210
 
209
 
  return std::make_pair(row_map.MinMyGID(), row_map.MaxMyGID() + 1);
 
211
  return std::make_pair(row_map.MinMyGID64(), row_map.MaxMyGID64() + 1);
210
212
}
211
213
//-----------------------------------------------------------------------------
212
 
void EpetraMatrix::resize(GenericVector& z, uint dim) const
 
214
void EpetraMatrix::resize(GenericVector& z, std::size_t dim) const
213
215
{
214
216
  dolfin_assert(A);
215
217
 
231
233
  _z.reset(*map);
232
234
}
233
235
//-----------------------------------------------------------------------------
234
 
void EpetraMatrix::get(double* block, uint m, const uint* rows,
235
 
                       uint n, const uint* cols) const
 
236
void EpetraMatrix::get(double* block, std::size_t m, const DolfinIndex* rows,
 
237
                       std::size_t n, const DolfinIndex* cols) const
236
238
{
237
239
  dolfin_assert(A);
238
240
 
241
243
  double* values;
242
244
 
243
245
  // For each row in rows
244
 
  for(uint i = 0; i < m; ++i)
 
246
  for(std::size_t i = 0; i < m; ++i)
245
247
  {
246
248
    // Extract the values and indices from row: rows[i]
247
249
    if (A->IndicesAreLocal())
268
270
 
269
271
    // Step the indices to the start of cols
270
272
    int k = 0;
271
 
    while (indices[k] < static_cast<int>(cols[0]))
 
273
    while (indices[k] < (int) cols[0])
272
274
      k++;
273
275
 
274
276
    // Fill the collumns in the block
275
 
    for (uint j = 0; j < n; j++)
 
277
    for (std::size_t j = 0; j < n; j++)
276
278
    {
277
 
      if (k < num_entities and indices[k] == static_cast<int>(cols[j]))
 
279
      if (k < num_entities and indices[k] == (int) cols[j])
278
280
      {
279
281
        block[i*n + j] = values[k];
280
282
        k++;
286
288
}
287
289
//-----------------------------------------------------------------------------
288
290
void EpetraMatrix::set(const double* block,
289
 
                       uint m, const uint* rows,
290
 
                       uint n, const uint* cols)
 
291
                       std::size_t m, const DolfinIndex* rows,
 
292
                       std::size_t n, const DolfinIndex* cols)
291
293
{
292
294
  // This function is awkward and somewhat restrictive because of the
293
295
  // poor support for setting off-process values in Epetra
294
296
 
295
297
  dolfin_assert(A);
296
298
 
297
 
  // Check that all rows are local to this process
298
 
  /*
299
 
  for (uint i = 0; i < m; ++i)
300
 
  {
301
 
    if (!A->MyGRID(rows[i]))
302
 
    {
303
 
      dolfin_error("EpetraMatrix.cpp",
304
 
                   "set block of values for Epetra matrix",
305
 
                   "Only values in row belonging to process can be set");
306
 
    }
307
 
  }
308
 
  */
309
 
 
310
 
  const int err = A->ReplaceGlobalValues(m,
311
 
      reinterpret_cast<const int*>(rows), n,
312
 
      reinterpret_cast<const int*>(cols), block, Epetra_FECrsMatrix::ROW_MAJOR);
 
299
  const int err = A->ReplaceGlobalValues(m, rows, n, cols, block,
 
300
                                         Epetra_FECrsMatrix::ROW_MAJOR);
313
301
  if (err != 0)
314
302
  {
315
303
    dolfin_error("EpetraMatrix.cpp",
319
307
}
320
308
//-----------------------------------------------------------------------------
321
309
void EpetraMatrix::add(const double* block,
322
 
                       uint m, const uint* rows,
323
 
                       uint n, const uint* cols)
 
310
                       std::size_t m, const DolfinIndex* rows,
 
311
                       std::size_t n, const DolfinIndex* cols)
324
312
{
325
313
  dolfin_assert(A);
326
 
  const std::pair<uint, uint> local_row_range = local_range(0);
327
 
  for (uint i = 0; i < m; ++i)
328
 
  {
329
 
    const uint row = rows[i];
330
 
    if (row >= local_row_range.first && row < local_row_range.second)
331
 
    {
332
 
      A->Epetra_CrsMatrix::SumIntoGlobalValues((int) row, n, block + i*n,
333
 
                                          reinterpret_cast<const int*>(cols));
334
 
    }
335
 
    else
336
 
    {
337
 
      A->SumIntoGlobalValues(1, reinterpret_cast<const int*>(rows + i),
338
 
                             n, reinterpret_cast<const int*>(cols), block + i*n,
339
 
                             Epetra_FECrsMatrix::ROW_MAJOR);
340
 
    }
341
 
  }
342
 
 
343
 
  /*
344
 
  const int err = A->SumIntoGlobalValues(m, reinterpret_cast<const int*>(rows),
345
 
                                         n, reinterpret_cast<const int*>(cols), block,
 
314
  const int err = A->SumIntoGlobalValues(m, rows, n, cols, block,
346
315
                                         Epetra_FECrsMatrix::ROW_MAJOR);
347
316
  if (err != 0)
348
317
  {
350
319
                 "add block of values to Epetra matrix",
351
320
                 "Did not manage to perform Epetra_FECrsMatrix::SumIntoGlobalValues");
352
321
  }
353
 
  */
354
322
}
355
323
//-----------------------------------------------------------------------------
356
 
void EpetraMatrix::axpy(double a, const GenericMatrix& A, bool same_nonzero_pattern)
 
324
void EpetraMatrix::axpy(double a, const GenericMatrix& A,
 
325
                        bool same_nonzero_pattern)
357
326
{
358
327
  const EpetraMatrix* AA = &as_type<const EpetraMatrix>(A);
359
328
  if (!AA->mat()->Filled())
363
332
                 "Epetra matrix is not in the correct state");
364
333
  }
365
334
 
366
 
  const int err = EpetraExt::MatrixMatrix::Add(*(AA->mat()), false, a, *(this->A), 1.0);
 
335
  const int err = EpetraExt::MatrixMatrix::Add(*(AA->mat()), false, a,
 
336
                                               *(this->A), 1.0);
367
337
  if (err != 0)
368
338
  {
369
339
    dolfin_error("EpetraMatrix.cpp",
409
379
  if (mode == "add")
410
380
    err = A->GlobalAssemble(true, Add);
411
381
  else if (mode == "insert")
412
 
    //err = A->GlobalAssemble(true, Insert);
413
382
    err = A->GlobalAssemble(true);
414
383
  else if (mode == "flush")
415
384
  {
448
417
  return s.str();
449
418
}
450
419
//-----------------------------------------------------------------------------
451
 
void EpetraMatrix::ident(uint m, const uint* rows)
 
420
void EpetraMatrix::ident(std::size_t m, const DolfinIndex* rows)
452
421
{
453
422
  dolfin_assert(A);
454
423
  dolfin_assert(A->Filled() == true);
459
428
  // but is due to the sparsity pattern computation). This function only
460
429
  // work for locally owned rows (the PETSc version works for any row).
461
430
 
462
 
  typedef boost::unordered_set<uint> MySet;
 
431
  typedef boost::unordered_set<std::size_t> MySet;
463
432
 
464
433
  // Build lists of local and nonlocal rows
465
434
  MySet local_rows;
466
 
  std::vector<uint> non_local_rows;
467
 
  for (uint i = 0; i < m; ++i)
 
435
  std::vector<std::size_t> non_local_rows;
 
436
  for (std::size_t i = 0; i < m; ++i)
468
437
  {
469
 
    if (A->MyGlobalRow(static_cast<int>(rows[i])))
 
438
    if (A->MyGlobalRow(rows[i]))
470
439
      local_rows.insert(rows[i]);
471
440
    else
472
441
      non_local_rows.push_back(rows[i]);
476
445
  if (MPI::num_processes() > 1)
477
446
  {
478
447
    // Send list of nonlocal rows to all processes
479
 
    std::vector<uint> destinations;
480
 
    std::vector<uint> send_data;
481
 
    for (uint i = 0; i < MPI::num_processes(); ++i)
 
448
    std::vector<std::size_t> destinations;
 
449
    std::vector<std::size_t> send_data;
 
450
    for (std::size_t i = 0; i < MPI::num_processes(); ++i)
482
451
    {
483
452
      if (i != MPI::process_number())
484
453
      {
488
457
      }
489
458
    }
490
459
 
491
 
    std::vector<uint> received_data;
 
460
    std::vector<std::size_t> received_data;
492
461
    MPI::distribute(send_data, destinations, received_data);
493
462
 
494
463
    // Unpack data
495
 
    for (uint i = 0; i < received_data.size(); ++i)
 
464
    for (std::size_t i = 0; i < received_data.size(); ++i)
496
465
    {
497
466
      // Insert row into set if it's local
498
 
      const uint new_index = received_data[i];
499
 
      if (A->MyGlobalRow(static_cast<int>(new_index)))
 
467
      const int new_index = received_data[i];
 
468
      if (A->MyGlobalRow(new_index))
500
469
        local_rows.insert(new_index);
501
470
    }
502
471
  }
507
476
  for (global_row = local_rows.begin(); global_row != local_rows.end(); ++global_row)
508
477
  {
509
478
    // Get local row index
510
 
    const int local_row = A->LRID(static_cast<int>(*global_row));
 
479
    const DolfinIndex _global_row = *global_row;
 
480
    const int local_row = A->LRID(_global_row);
511
481
 
512
482
    // If this process owns row, then zero row
513
483
    if (local_row >= 0)
540
510
  }
541
511
}
542
512
//-----------------------------------------------------------------------------
543
 
void EpetraMatrix::zero(uint m, const uint* rows)
 
513
void EpetraMatrix::zero(std::size_t m, const DolfinIndex* rows)
544
514
{
545
515
  // FIXME: This can be made more efficient by eliminating creation of
546
516
  //        some obejcts inside the loop
547
517
 
548
518
  dolfin_assert(A);
549
519
  const Epetra_CrsGraph& graph = A->Graph();
550
 
  for (uint i = 0; i < m; ++i)
 
520
  for (std::size_t i = 0; i < m; ++i)
551
521
  {
552
522
    const int row = rows[i];
553
523
    const int num_nz = graph.NumGlobalIndices(row);
639
609
  }
640
610
}
641
611
//-----------------------------------------------------------------------------
642
 
void EpetraMatrix::getrow(uint row, std::vector<uint>& columns,
 
612
void EpetraMatrix::getrow(std::size_t row, std::vector<std::size_t>& columns,
643
613
                          std::vector<double>& values) const
644
614
{
645
615
  dolfin_assert(A);
646
616
 
647
617
  // Get local row index
648
 
  const int local_row_index = A->LRID(static_cast<int>(row));
 
618
  const DolfinIndex _row = row;
 
619
  const int local_row_index = A->LRID(_row);
649
620
 
650
621
  // If this process has part of the row, get values
651
622
  if (local_row_index >= 0)
680
651
  }
681
652
}
682
653
//-----------------------------------------------------------------------------
683
 
void EpetraMatrix::setrow(uint row, const std::vector<uint>& columns,
 
654
void EpetraMatrix::setrow(std::size_t row,
 
655
                          const std::vector<std::size_t>& columns,
684
656
                          const std::vector<double>& values)
685
657
{
686
 
  static bool print_msg_once=true;
 
658
  static bool print_msg_once = true;
687
659
  if (print_msg_once)
688
660
  {
689
661
    info("EpetraMatrix::setrow is implemented inefficiently");
690
662
    print_msg_once = false;
691
663
  }
692
664
 
693
 
  for (uint i=0; i < columns.size(); i++)
694
 
    set(&values[i], 1, &row, 1, &columns[i]);
 
665
  for (std::size_t i = 0; i < columns.size(); i++)
 
666
  {
 
667
    DolfinIndex _row = row;
 
668
    DolfinIndex _col = columns[i];
 
669
    set(&values[i], 1, &_row, 1, &_col);
 
670
  }
695
671
}
696
672
//-----------------------------------------------------------------------------
697
673
GenericLinearAlgebraFactory& EpetraMatrix::factory() const