trilinos_helpers.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #include "trilinos_helpers.h"
27 
28 namespace oomph
29 {
30  // VECTOR METHODS
31  // =============================================================
32 
33 
34  //=============================================================================
35  /// create an Epetra_Vector from an oomph-lib DoubleVector.
36  /// If oomph_vec is NOT distributed (i.e. locally replicated) and
37  /// on more than one processor, then the returned Epetra_Vector will be
38  /// uniformly distributed. If the oomph_vec is distributed then the
39  /// Epetra_Vector returned will have the same distribution as oomph_vec.
40  //=============================================================================
42  const DoubleVector& oomph_vec)
43  {
44 #ifdef PARANOID
45  // check the the oomph lib vector is setup
46  if (!oomph_vec.built())
47  {
48  std::ostringstream error_message;
49  error_message << "The oomph-lib vector (oomph_v) must be setup.";
50  throw OomphLibError(
51  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
52  }
53 #endif
54 
55  // create the corresponding Epetra_Map
56  LinearAlgebraDistribution* dist_pt = 0;
57  if (oomph_vec.distributed())
58  {
59  dist_pt = new LinearAlgebraDistribution(oomph_vec.distribution_pt());
60  }
61  else
62  {
63  dist_pt = new LinearAlgebraDistribution(
64  oomph_vec.distribution_pt()->communicator_pt(), oomph_vec.nrow(), true);
65  }
66  Epetra_Map* epetra_map_pt = create_epetra_map(dist_pt);
67 
68  // first first coefficient of the oomph vector to be inserted into the
69  // Epetra_Vector
70  unsigned offset = 0;
71  if (!oomph_vec.distributed())
72  {
73  offset = dist_pt->first_row();
74  }
75 
76  // copy the values into the oomph-lib vector
77  // const_cast OK because Epetra_Vector construction is Copying values and
78  // therefore does not modify data.
79  double* v_pt = const_cast<double*>(oomph_vec.values_pt());
80  Epetra_Vector* epetra_vec_pt =
81  new Epetra_Vector(Copy, *epetra_map_pt, v_pt + offset);
82 
83  // clean up
84  delete epetra_map_pt;
85  delete dist_pt;
86 
87  // return
88  return epetra_vec_pt;
89  }
90 
91  //=============================================================================
92  /// create an Epetra_Vector based on the argument oomph-lib
93  /// LinearAlgebraDistribution
94  /// If dist is NOT distributed and
95  /// on more than one processor, then the returned Epetra_Vector will be
96  /// uniformly distributed. If dist is distributed then the Epetra_Vector
97  /// returned will have the same distribution as dist.
98  /// The coefficient values are not set.
99  //=============================================================================
101  const LinearAlgebraDistribution* dist_pt)
102  {
103 #ifdef PARANOID
104  // check the the oomph lib vector is setup
105  if (!dist_pt->built())
106  {
107  std::ostringstream error_message;
108  error_message << "The LinearAlgebraDistribution dist_pt must be setup.";
109  throw OomphLibError(
110  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
111  }
112 #endif
113 
114  // create the corresponding Epetra_Map
115  LinearAlgebraDistribution* target_dist_pt = 0;
116  if (dist_pt->distributed())
117  {
118  target_dist_pt = new LinearAlgebraDistribution(dist_pt);
119  }
120  else
121  {
122  target_dist_pt = new LinearAlgebraDistribution(
123  dist_pt->communicator_pt(), dist_pt->nrow(), true);
124  }
125  Epetra_Map* epetra_map_pt = create_epetra_map(target_dist_pt);
126 
127  // create epetra_vector
128  Epetra_Vector* epetra_vec_pt = new Epetra_Vector(*epetra_map_pt, false);
129 
130  // clean up
131  delete epetra_map_pt;
132  delete target_dist_pt;
133 
134  // return
135  return epetra_vec_pt;
136  }
137 
138  //=============================================================================
139  /// Create an Epetra_Vector equivalent of DoubleVector
140  /// The argument DoubleVector must be built.
141  /// The Epetra_Vector will point to, and NOT COPY the underlying data in the
142  /// DoubleVector.
143  /// The oomph-lib DoubleVector and the returned Epetra_Vector will have the
144  /// the same distribution.
145  //=============================================================================
147  DoubleVector& oomph_vec)
148  {
149 #ifdef PARANOID
150  // check the the oomph lib vector is setup
151  if (!oomph_vec.built())
152  {
153  std::ostringstream error_message;
154  error_message << "The oomph-lib vector (oomph_v) must be setup.";
155  throw OomphLibError(
156  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
157  }
158 #endif
159 
160  // create the corresponding Epetra_Map
161  Epetra_Map* epetra_map_pt = create_epetra_map(oomph_vec.distribution_pt());
162 
163  // copy the values into the oomph-lib vector
164  double* v_pt = oomph_vec.values_pt();
165  Epetra_Vector* epetra_vec_pt =
166  new Epetra_Vector(View, *epetra_map_pt, v_pt);
167 
168  // clean up
169  delete epetra_map_pt;
170 
171  // return
172  return epetra_vec_pt;
173  }
174 
175  //=============================================================================
176  /// Helper function to copy the contents of a Trilinos vector to an
177  /// oomph-lib distributed vector. The distribution of the two vectors must
178  /// be identical
179  //=============================================================================
181  const Epetra_Vector* epetra_vec_pt, DoubleVector& oomph_vec)
182  {
183 #ifdef PARANOID
184  // check the the oomph lib vector is setup
185  if (!oomph_vec.built())
186  {
187  std::ostringstream error_message;
188  error_message << "The oomph-lib vector (oomph_v) must be setup.";
189  throw OomphLibError(
190  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
191  }
192 #endif
193 
194  // if the oomph-lib vector is distributed
195  if (oomph_vec.distributed())
196  {
197  // extract values from epetra_v_pt
198  double* v_values;
199  epetra_vec_pt->ExtractView(&v_values);
200 
201  // copy the values
202  unsigned nrow_local = oomph_vec.nrow_local();
203  for (unsigned i = 0; i < nrow_local; i++)
204  {
205  oomph_vec[i] = v_values[i];
206  }
207  }
208 
209  // else teh oomph-lib vector is not distributed
210  else
211  {
212  // get the values vector
213 #ifdef OOMPH_HAS_MPI
214  int nproc = epetra_vec_pt->Map().Comm().NumProc();
215  if (nproc == 1)
216  {
217  epetra_vec_pt->ExtractCopy(oomph_vec.values_pt());
218  }
219  else
220  {
221  // get the local values
222  double* local_values;
223  epetra_vec_pt->ExtractView(&local_values);
224 
225  // my rank
226  int my_rank = epetra_vec_pt->Map().Comm().MyPID();
227 
228  // number of local rows
229  Vector<int> nrow_local(nproc);
230  nrow_local[my_rank] = epetra_vec_pt->MyLength();
231 
232  // gather the First_row vector
233  int my_nrow_local_copy = nrow_local[my_rank];
234  MPI_Allgather(
235  &my_nrow_local_copy,
236  1,
237  MPI_INT,
238  &nrow_local[0],
239  1,
240  MPI_INT,
241  oomph_vec.distribution_pt()->communicator_pt()->mpi_comm());
242 
243  // number of local rows
244  Vector<int> first_row(nproc);
245  first_row[my_rank] = epetra_vec_pt->Map().MyGlobalElements()[0];
246 
247  // gather the First_row vector
248  int my_first_row = first_row[my_rank];
249  MPI_Allgather(
250  &my_first_row,
251  1,
252  MPI_INT,
253  &first_row[0],
254  1,
255  MPI_INT,
256  oomph_vec.distribution_pt()->communicator_pt()->mpi_comm());
257 
258  // gather the local solution values
259  MPI_Allgatherv(
260  local_values,
261  nrow_local[my_rank],
262  MPI_DOUBLE,
263  oomph_vec.values_pt(),
264  &nrow_local[0],
265  &first_row[0],
266  MPI_DOUBLE,
267  oomph_vec.distribution_pt()->communicator_pt()->mpi_comm());
268  }
269 #else
270  epetra_vec_pt->ExtractCopy(oomph_vec.values_pt());
271 #endif
272  }
273  }
274 
275  // MATRIX METHODS
276  // =============================================================
277 
278 
279  //=============================================================================
280  /// create an Epetra_CrsMatrix from an oomph-lib CRDoubleMatrix.
281  /// If oomph_matrix_pt is NOT distributed (i.e. locally replicated) and
282  /// on more than one processor, then the returned Epetra_Vector will be
283  /// uniformly distributed. If the oomph_matrix_pt is distributed then the
284  /// Epetra_CrsMatrix returned will have the same distribution as
285  /// oomph_matrix_pt.
286  /// The LinearAlgebraDistribution argument dist_pt should specify the
287  /// distribution of the object this matrix will operate on.
288  //=============================================================================
290  const CRDoubleMatrix* oomph_matrix_pt,
291  const LinearAlgebraDistribution* dist_pt)
292  {
293 #ifdef PARANOID
294  if (!oomph_matrix_pt->built())
295  {
296  std::ostringstream error_message;
297  error_message << "The oomph-lib matrix must be built.";
298  throw OomphLibError(
299  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
300  }
301  if (!oomph_matrix_pt->built())
302  {
303  std::ostringstream error_message;
304  error_message << "The oomph-lib matrix must be built.";
305  throw OomphLibError(
306  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
307  }
308  if (!oomph_matrix_pt->built())
309  {
310  std::ostringstream error_message;
311  error_message << "The oomph-lib matrix must be built.";
312  throw OomphLibError(
313  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
314  }
315 #endif
316 
317  // get pointers to the matrix values, column indices etc
318  // const_cast is safe because we use the Epetra_Vector "Copy" construction
319  // method
320  int* column = const_cast<int*>(oomph_matrix_pt->column_index());
321  double* value = const_cast<double*>(oomph_matrix_pt->value());
322  int* row_start = const_cast<int*>(oomph_matrix_pt->row_start());
323 
324  // create the corresponding Epetra_Map
325  LinearAlgebraDistribution* target_dist_pt = 0;
326  if (oomph_matrix_pt->distributed())
327  {
328  target_dist_pt =
329  new LinearAlgebraDistribution(oomph_matrix_pt->distribution_pt());
330  }
331  else
332  {
333  target_dist_pt = new LinearAlgebraDistribution(
334  oomph_matrix_pt->distribution_pt()->communicator_pt(),
335  oomph_matrix_pt->nrow(),
336  true);
337  }
338  Epetra_Map* epetra_map_pt = create_epetra_map(target_dist_pt);
339 
340  // first first coefficient of the oomph vector to be inserted into the
341  // Epetra_Vector
342  unsigned offset = 0;
343  if (!oomph_matrix_pt->distributed())
344  {
345  offset = target_dist_pt->first_row();
346  }
347 
348  // get my nrow_local and first_row
349  unsigned nrow_local = target_dist_pt->nrow_local();
350  unsigned first_row = target_dist_pt->first_row();
351 
352  // store the number of non zero entries per row
353  int* nnz_per_row = new int[nrow_local];
354  for (unsigned row = 0; row < nrow_local; row++)
355  {
356  nnz_per_row[row] = row_start[row + offset + 1] - row_start[offset + row];
357  }
358 
359  // create the matrix
360  Epetra_CrsMatrix* epetra_matrix_pt =
361  new Epetra_CrsMatrix(Copy, *epetra_map_pt, nnz_per_row, true);
362 
363  // insert the values
364  for (unsigned row = 0; row < nrow_local; row++)
365  {
366  // get pointer to this row in values/columns
367  int ptr = row_start[row + offset];
368 #ifdef PARANOID
369  int err = 0;
370  err =
371 #endif
372  epetra_matrix_pt->InsertGlobalValues(
373  first_row + row, nnz_per_row[row], value + ptr, column + ptr);
374 #ifdef PARANOID
375  if (err != 0)
376  {
377  std::ostringstream error_message;
378  error_message
379  << "Epetra Matrix Insert Global Values : epetra_error_flag = " << err;
380  throw OomphLibError(error_message.str(),
381  OOMPH_CURRENT_FUNCTION,
382  OOMPH_EXCEPTION_LOCATION);
383  }
384 #endif
385  }
386 
387  // complete the build of the trilinos matrix
388  LinearAlgebraDistribution* target_col_dist_pt = 0;
389  if (dist_pt->distributed())
390  {
391  target_col_dist_pt = new LinearAlgebraDistribution(dist_pt);
392  }
393  else
394  {
395  target_col_dist_pt = new LinearAlgebraDistribution(
396  dist_pt->communicator_pt(), dist_pt->nrow(), true);
397  }
398  Epetra_Map* epetra_domain_map_pt = create_epetra_map(target_col_dist_pt);
399 #ifdef PARANOID
400  int err = 0;
401  err =
402 #endif
403  epetra_matrix_pt->FillComplete(*epetra_domain_map_pt, *epetra_map_pt);
404 #ifdef PARANOID
405  if (err != 0)
406  {
407  std::ostringstream error_message;
408  error_message
409  << "Epetra Matrix Fill Complete Error : epetra_error_flag = " << err;
410  throw OomphLibError(
411  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
412  }
413 #endif
414 
415  // tidy up memory
416  delete[] nnz_per_row;
417  delete epetra_map_pt;
418  delete epetra_domain_map_pt;
419  delete target_dist_pt;
420  delete target_col_dist_pt;
421 
422  // return
423  return epetra_matrix_pt;
424  }
425 
426 
427  //=============================================================================
428  /// Class to allow sorting of column indices in conversion to epetra matrix
429  //=============================================================================
431  {
432  public:
433  /// Constructor: Pass number of first column and the number of local columns
434  DistributionPredicate(const int& first_col, const int& ncol_local)
435  : First_col(first_col), Last_col(first_col + ncol_local - 1)
436  {
437  }
438 
439  /// Comparison operator: is column col in the range
440  /// between (including) First_col and Last_col
441  bool operator()(const int& col)
442  {
443  if (col >= First_col && col <= Last_col)
444  {
445  return true;
446  }
447  else
448  {
449  return false;
450  }
451  }
452 
453  private:
454  /// First column held locally
456 
457  /// Last colum held locally
458  int Last_col;
459  };
460 
461 
462  //=============================================================================
463  /// create and Epetra_CrsMatrix from an oomph-lib CRDoubleMatrix.
464  /// Specialisation for Trilinos AztecOO.
465  /// If oomph_matrix_pt is NOT distributed (i.e. locally replicated) and
466  /// on more than one processor, then the returned Epetra_Vector will be
467  /// uniformly distributed. If the oomph_matrix_pt is distributed then the
468  /// Epetra_CrsMatrix returned will have the same distribution as
469  /// oomph_matrix_pt.
470  /// For AztecOO, the column map is ordered such that the local rows are
471  /// first.
472  //=============================================================================
473  Epetra_CrsMatrix* TrilinosEpetraHelpers::
475  CRDoubleMatrix* oomph_matrix_pt)
476  {
477 #ifdef PARANOID
478  if (!oomph_matrix_pt->built())
479  {
480  std::ostringstream error_message;
481  error_message << "The oomph-lib matrix must be built.";
482  throw OomphLibError(
483  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
484  }
485 #endif
486 
487  // get pointers to the matrix values, column indices etc
488  // const_cast is safe because we use the Epetra_Vector "Copy" construction
489  // method
490  int* column = const_cast<int*>(oomph_matrix_pt->column_index());
491  double* value = const_cast<double*>(oomph_matrix_pt->value());
492  int* row_start = const_cast<int*>(oomph_matrix_pt->row_start());
493 
494  // create the corresponding Epetra_Map
495  LinearAlgebraDistribution* target_dist_pt = 0;
496  if (oomph_matrix_pt->distributed())
497  {
498  target_dist_pt =
499  new LinearAlgebraDistribution(oomph_matrix_pt->distribution_pt());
500  }
501  else
502  {
503  target_dist_pt = new LinearAlgebraDistribution(
504  oomph_matrix_pt->distribution_pt()->communicator_pt(),
505  oomph_matrix_pt->nrow(),
506  true);
507  }
508  Epetra_Map* epetra_map_pt = create_epetra_map(target_dist_pt);
509 
510  // create the epetra column map
511 
512 #ifdef OOMPH_HAS_MPI
513  int first_col = oomph_matrix_pt->first_row();
514  int ncol_local = oomph_matrix_pt->nrow_local();
515 
516  // Build colum map
517  Epetra_Map* epetra_col_map_pt = 0;
518  {
519  // Vector of column indices; on processor goes first
520  std::vector<int> col_index_vector;
521  col_index_vector.reserve(oomph_matrix_pt->nnz() + ncol_local);
522  col_index_vector.resize(ncol_local);
523 
524  // Global column indices corresponding to on-processor rows
525  for (int c = 0; c < ncol_local; ++c)
526  {
527  col_index_vector[c] = c + first_col;
528  }
529 
530  // Remember where the on-processor rows (columns) end
531  std::vector<int>::iterator mid = col_index_vector.end();
532 
533  // Now insert ALL column indices of ALL entries
534  col_index_vector.insert(mid, column, column + oomph_matrix_pt->nnz());
535 
536  // Loop over the newly added entries and remove them if they
537  // refer to on-processor columns
538  std::vector<int>::iterator end =
539  std::remove_if(mid,
540  col_index_vector.end(),
541  DistributionPredicate(first_col, ncol_local));
542 
543  // Now sort the newly added entries
544  std::sort(mid, end);
545 
546  //...and remove duplicates
547  end = std::unique(mid, end);
548 
549  // Make the map
550  epetra_col_map_pt = new Epetra_Map(
551  -1,
552  end - col_index_vector.begin(),
553  &col_index_vector[0],
554  0,
555  Epetra_MpiComm(
556  oomph_matrix_pt->distribution_pt()->communicator_pt()->mpi_comm()));
557 
558  // Hack to clear memory
559  std::vector<int>().swap(col_index_vector);
560  }
561 
562 #else
563 
564  int ncol = oomph_matrix_pt->ncol();
565  Epetra_Map* epetra_col_map_pt =
566  new Epetra_LocalMap(ncol, 0, Epetra_SerialComm());
567 
568 #endif
569 
570  // first first coefficient of the oomph vector to be inserted into the
571  // Epetra_Vector
572  unsigned offset = 0;
573  if (!oomph_matrix_pt->distributed())
574  {
575  offset = target_dist_pt->first_row();
576  }
577 
578  // get my nrow_local and first_row
579  unsigned nrow_local = target_dist_pt->nrow_local();
580  unsigned first_row = target_dist_pt->first_row();
581 
582  // store the number of non zero entries per row
583  int* nnz_per_row = new int[nrow_local];
584  for (unsigned row = 0; row < nrow_local; ++row)
585  {
586  nnz_per_row[row] = row_start[row + offset + 1] - row_start[offset + row];
587  }
588 
589  // create the matrix
590  Epetra_CrsMatrix* epetra_matrix_pt = new Epetra_CrsMatrix(
591  Copy, *epetra_map_pt, *epetra_col_map_pt, nnz_per_row, true);
592 
593  // insert the values
594  for (unsigned row = 0; row < nrow_local; row++)
595  {
596  // get pointer to this row in values/columns
597  int ptr = row_start[row + offset];
598 #ifdef PARANOID
599  int err = 0;
600  err =
601 #endif
602  epetra_matrix_pt->InsertGlobalValues(
603  first_row + row, nnz_per_row[row], value + ptr, column + ptr);
604 #ifdef PARANOID
605  if (err != 0)
606  {
607  std::ostringstream error_message;
608  error_message
609  << "Epetra Matrix Insert Global Values : epetra_error_flag = " << err;
610  throw OomphLibError(error_message.str(),
611  OOMPH_CURRENT_FUNCTION,
612  OOMPH_EXCEPTION_LOCATION);
613  }
614 #endif
615  }
616 
617  // complete the build of the trilinos matrix
618 #ifdef PARANOID
619  int err = 0;
620  err =
621 #endif
622  epetra_matrix_pt->FillComplete();
623 
624 #ifdef PARANOID
625  if (err != 0)
626  {
627  std::ostringstream error_message;
628  error_message
629  << "Epetra Matrix Fill Complete Error : epetra_error_flag = " << err;
630  throw OomphLibError(
631  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
632  }
633 #endif
634 
635  // tidy up memory
636  delete[] nnz_per_row;
637  delete epetra_map_pt;
638  delete epetra_col_map_pt;
639  delete target_dist_pt;
640 
641  // return
642  return epetra_matrix_pt;
643  }
644 
645 
646  // MATRIX OPERATION METHODS ==================================================
647 
648 
649  //============================================================================
650  /// Function to perform a matrix-vector multiplication on a
651  /// oomph-lib matrix and vector using Trilinos functionality.
652  /// NOTE 1. the matrix and the vectors must have the same communicator.
653  /// NOTE 2. The vector will be returned with the same distribution
654  /// as the matrix, unless a distribution is predefined in the solution
655  /// vector in which case the vector will be returned with that distribution.
656  //============================================================================
657  void TrilinosEpetraHelpers::multiply(const CRDoubleMatrix* oomph_matrix_pt,
658  const DoubleVector& oomph_x,
659  DoubleVector& oomph_y)
660  {
661 #ifdef PARANOID
662  // check that this matrix is built
663  if (!oomph_matrix_pt->built())
664  {
665  std::ostringstream error_message_stream;
666  error_message_stream << "This matrix has not been built";
667  throw OomphLibError(error_message_stream.str(),
668  OOMPH_CURRENT_FUNCTION,
669  OOMPH_EXCEPTION_LOCATION);
670  }
671 
672  // check that the distribution of the matrix and the soln are the same
673  if (oomph_y.built())
674  {
675  if (!(*oomph_matrix_pt->distribution_pt() == *oomph_y.distribution_pt()))
676  {
677  std::ostringstream error_message_stream;
678  error_message_stream
679  << "The soln vector and this matrix must have the same distribution.";
680  throw OomphLibError(error_message_stream.str(),
681  OOMPH_CURRENT_FUNCTION,
682  OOMPH_EXCEPTION_LOCATION);
683  }
684  }
685 
686  // check that the distribution of the oomph-lib vector x is setup
687  if (!oomph_x.built())
688  {
689  std::ostringstream error_message_stream;
690  error_message_stream << "The x vector must be setup";
691  throw OomphLibError(error_message_stream.str(),
692  OOMPH_CURRENT_FUNCTION,
693  OOMPH_EXCEPTION_LOCATION);
694  }
695 #endif
696 
697  // setup the distribution
698  if (!oomph_y.distribution_pt()->built())
699  {
700  oomph_y.build(oomph_matrix_pt->distribution_pt(), 0.0);
701  }
702 
703  // convert matrix1 to epetra matrix
704  Epetra_CrsMatrix* epetra_matrix_pt = create_distributed_epetra_matrix(
705  oomph_matrix_pt, oomph_x.distribution_pt());
706 
707  // convert x to Trilinos vector
708  Epetra_Vector* epetra_x_pt = create_distributed_epetra_vector(oomph_x);
709 
710  // create Trilinos vector for soln ( 'viewing' the contents of the oomph-lib
711  // matrix)
712  Epetra_Vector* epetra_y_pt = create_distributed_epetra_vector(oomph_y);
713 
714  // do the multiply
715 #ifdef PARANOID
716  int epetra_error_flag = 0;
717  epetra_error_flag =
718 #endif
719  epetra_matrix_pt->Multiply(false, *epetra_x_pt, *epetra_y_pt);
720 
721  // return the solution
722  copy_to_oomphlib_vector(epetra_y_pt, oomph_y);
723 
724  // throw error if there is an epetra error
725 #ifdef PARANOID
726  if (epetra_error_flag != 0)
727  {
728  std::ostringstream error_message;
729  error_message
730  << "Epetra Matrix Vector Multiply Error : epetra_error_flag = "
731  << epetra_error_flag;
732  throw OomphLibError(
733  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
734  }
735 #endif
736 
737  // clean up
738  delete epetra_matrix_pt;
739  delete epetra_x_pt;
740  delete epetra_y_pt;
741  }
742 
743  //=============================================================================
744  /// Function to perform a matrix-matrix multiplication on oomph-lib
745  /// matrices by using Trilinos functionality.
746  /// \b NOTE 1. There are two Trilinos matrix-matrix multiplication methods
747  /// available, using either the EpetraExt::MatrixMatrix class (if use_ml ==
748  /// false) or using ML (Epetra_MatrixMult method)
749  /// \b NOTE 2. the solution matrix (matrix_soln) will be returned with the
750  /// same distribution as matrix1
751  /// \b NOTE 3. All matrices must share the same communicator.
752  //=============================================================================
754  const CRDoubleMatrix& matrix_2,
755  CRDoubleMatrix& matrix_soln,
756  const bool& use_ml)
757  {
758 #ifdef PARANOID
759  // check that matrix 1 is built
760  if (!matrix_1.built())
761  {
762  std::ostringstream error_message_stream;
763  error_message_stream << "This matrix matrix_1 has not been built";
764  throw OomphLibError(error_message_stream.str(),
765  OOMPH_CURRENT_FUNCTION,
766  OOMPH_EXCEPTION_LOCATION);
767  }
768  // check that matrix 2 is built
769  if (!matrix_2.built())
770  {
771  std::ostringstream error_message_stream;
772  error_message_stream << "This matrix matrix_2 has not been built";
773  throw OomphLibError(error_message_stream.str(),
774  OOMPH_CURRENT_FUNCTION,
775  OOMPH_EXCEPTION_LOCATION);
776  }
777  // check matrix dimensions are compatable
778  if (matrix_1.ncol() != matrix_2.nrow())
779  {
780  std::ostringstream error_message;
781  error_message
782  << "Matrix dimensions incompatible for matrix-matrix multiplication"
783  << "ncol() for first matrix: " << matrix_1.ncol()
784  << "nrow() for second matrix: " << matrix_2.nrow();
785  throw OomphLibError(
786  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
787  }
788 
789  // check that the have the same communicator
790  OomphCommunicator temp_comm(matrix_1.distribution_pt()->communicator_pt());
791  if (temp_comm != *matrix_2.distribution_pt()->communicator_pt())
792  {
793  std::ostringstream error_message;
794  error_message
795  << "Matrix dimensions incompatible for matrix-matrix multiplication"
796  << "ncol() for first matrix: " << matrix_1.ncol()
797  << "nrow() for second matrix: " << matrix_2.nrow();
798  throw OomphLibError(
799  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
800  }
801  // check that the distribution of the matrix and the soln are the same
802  if (matrix_soln.distribution_pt()->built())
803  {
804  if (!(*matrix_soln.distribution_pt() == *matrix_1.distribution_pt()))
805  {
806  std::ostringstream error_message_stream;
807  error_message_stream << "The solution matrix and matrix_1 must have "
808  "the same distribution.";
809  throw OomphLibError(error_message_stream.str(),
810  OOMPH_CURRENT_FUNCTION,
811  OOMPH_EXCEPTION_LOCATION);
812  }
813  }
814 #endif
815 
816  // setup the distribution
817  if (!matrix_soln.distribution_pt()->built())
818  {
819  matrix_soln.build(matrix_1.distribution_pt());
820  }
821 
822  // temporary fix
823  // ML MM method only appears to work for square matrices
824  // Should be investigated further.
825  bool temp_use_ml = false;
826  if ((*matrix_1.distribution_pt() == *matrix_2.distribution_pt()) &&
827  (matrix_1.ncol() == matrix_2.ncol()))
828  {
829  temp_use_ml = use_ml;
830  }
831 
832  // create matrix 1
833  Epetra_CrsMatrix* epetra_matrix_1_pt =
834  create_distributed_epetra_matrix(&matrix_1, matrix_2.distribution_pt());
835 
836  // create matrix 2
837  LinearAlgebraDistribution matrix_2_column_dist(
838  matrix_2.distribution_pt()->communicator_pt(), matrix_2.ncol(), true);
839  Epetra_CrsMatrix* epetra_matrix_2_pt =
840  create_distributed_epetra_matrix(&matrix_2, &matrix_2_column_dist);
841 
842  // create the Trilinos epetra matrix to hold solution - will have same map
843  // (and number of rows) as matrix_1
844  Epetra_CrsMatrix* solution_pt;
845 
846  // do the multiplication
847  // ---------------------
848  if (temp_use_ml)
849  {
850  // there is a problem using this function, many pages of
851  // warning messages are issued....
852  // "tmpresult->InsertGlobalValues returned 3"
853  // and
854  // "Result_epet->InsertGlobalValues returned 3"
855  // from function ML_back_to_epetraCrs(...) in
856  // file ml_epetra_utils.cpp unless the
857  // relevant lines are commented out. However this function
858  // is much faster (at least on Biowulf) than the alternative
859  // below.
860  solution_pt = Epetra_MatrixMult(epetra_matrix_1_pt, epetra_matrix_2_pt);
861  }
862  else
863  {
864  // this method requires us to pass in the solution matrix
865  solution_pt = new Epetra_CrsMatrix(Copy, epetra_matrix_1_pt->RowMap(), 0);
866 #ifdef PARANOID
867  int epetra_error_flag = 0;
868  epetra_error_flag =
869 #endif
870  EpetraExt::MatrixMatrix::Multiply(
871  *epetra_matrix_1_pt, false, *epetra_matrix_2_pt, false, *solution_pt);
872 #ifdef PARANOID
873  if (epetra_error_flag != 0)
874  {
875  std::ostringstream error_message;
876  error_message << "error flag from Multiply(): " << epetra_error_flag
877  << " from TrilinosHelpers::multiply" << std::endl;
878  throw OomphLibError(error_message.str(),
879  OOMPH_CURRENT_FUNCTION,
880  OOMPH_EXCEPTION_LOCATION);
881  }
882 #endif
883  }
884 
885  // extract values and put into solution
886  // ------------------------------------
887 
888  // find
889  int nnz_local = solution_pt->NumMyNonzeros();
890  int nrow_local = matrix_1.nrow_local();
891 
892  // do some checks
893 #ifdef PARANOID
894  // check number of global rows in soluton matches that in matrix_1
895  if ((int)matrix_1.nrow() != solution_pt->NumGlobalRows())
896  {
897  std::ostringstream error_message;
898  error_message << "Incorrect number of global rows in solution matrix. "
899  << "nrow() for first input matrix: " << matrix_1.nrow()
900  << " nrow() for solution: " << solution_pt->NumGlobalRows();
901  throw OomphLibError(
902  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
903  }
904 
905  // check number of local rows in soluton matches that in matrix_1
906  if (static_cast<int>(matrix_1.nrow_local()) != solution_pt->NumMyRows())
907  {
908  std::ostringstream error_message;
909  error_message << "Incorrect number of local rows in solution matrix. "
910  << "nrow_local() for first input matrix: "
911  << matrix_1.nrow_local() << " nrow_local() for solution: "
912  << solution_pt->NumMyRows();
913  throw OomphLibError(
914  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
915  }
916 
917  // check number of global columns in soluton matches that in matrix_2
918  if ((int)matrix_2.ncol() != solution_pt->NumGlobalCols())
919  {
920  std::ostringstream error_message;
921  error_message << "Incorrect number of global columns in solution matrix. "
922  << "ncol() for second input matrix: " << matrix_2.ncol()
923  << " ncol() for solution: " << solution_pt->NumGlobalCols();
924  throw OomphLibError(
925  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
926  }
927 
928  // check global index of the first row matches
929  if (static_cast<int>(matrix_1.first_row()) != solution_pt->GRID(0))
930  {
931  std::ostringstream error_message;
932  error_message
933  << "Incorrect global index for first row of solution matrix. "
934  << "first_row() for first input matrix : " << matrix_1.first_row()
935  << " first_row() for solution: " << solution_pt->GRID(0);
936  throw OomphLibError(
937  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
938  }
939 #endif
940 
941  // extract values from Epetra matrix row by row
942  double* value = new double[nnz_local];
943  int* column_index = new int[nnz_local];
944  int* row_start = new int[nrow_local + 1];
945  int ptr = 0;
946  int num_entries = 0;
947  int first = matrix_soln.first_row();
948  int last = first + matrix_soln.nrow_local();
949  for (int row = first; row < last; row++)
950  {
951  row_start[row - first] = ptr;
952  solution_pt->ExtractGlobalRowCopy(
953  row, nnz_local, num_entries, value + ptr, column_index + ptr);
954  ptr += num_entries;
955  }
956  row_start[nrow_local] = ptr;
957 
958  // delete Trilinos objects
959  delete epetra_matrix_1_pt;
960  delete epetra_matrix_2_pt;
961  delete solution_pt;
962 
963  // Build the Oomph-lib solution matrix using build function
964  matrix_soln.build(matrix_1.distribution_pt());
965  matrix_soln.build_without_copy(
966  matrix_2.ncol(), nnz_local, value, column_index, row_start);
967  }
968 
969 
970  // HELPER METHODS
971  // =============================================================
972 
973 
974  //=============================================================================
975  /// create an Epetra_Map corresponding to the LinearAlgebraDistribution
976  //=============================================================================
978  const LinearAlgebraDistribution* const dist_pt)
979  {
980 #ifdef OOMPH_HAS_MPI
981  if (dist_pt->distributed())
982  {
983  unsigned first_row = dist_pt->first_row();
984  unsigned nrow_local = dist_pt->nrow_local();
985  int* my_global_rows = new int[nrow_local];
986  for (unsigned i = 0; i < nrow_local; ++i)
987  {
988  my_global_rows[i] = first_row + i;
989  }
990  Epetra_Map* epetra_map_pt =
991  new Epetra_Map(dist_pt->nrow(),
992  nrow_local,
993  my_global_rows,
994  0,
995  Epetra_MpiComm(dist_pt->communicator_pt()->mpi_comm()));
996  delete[] my_global_rows;
997  return epetra_map_pt;
998  }
999  else
1000  {
1001  return new Epetra_LocalMap(
1002  int(dist_pt->nrow()),
1003  int(0),
1004  Epetra_MpiComm(dist_pt->communicator_pt()->mpi_comm()));
1005  }
1006 #else
1007  return new Epetra_LocalMap(
1008  int(dist_pt->nrow()), int(0), Epetra_SerialComm());
1009 #endif
1010  }
1011 
1012 } // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
int * column_index()
Access to C-style column index array.
Definition: matrices.h:1072
int * row_start()
Access to C-style row_start array.
Definition: matrices.h:1060
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:1008
unsigned long nnz() const
Return the number of nonzero entries (the local nnz)
Definition: matrices.h:1096
bool built() const
access function to the Built flag - indicates whether the matrix has been build - i....
Definition: matrices.h:1210
void build_without_copy(const unsigned &ncol, const unsigned &nnz, double *value, int *column_index, int *row_start)
keeps the existing distribution and just matrix that is stored without copying the matrix data
Definition: matrices.cc:1710
double * value()
Access to C-style value array.
Definition: matrices.h:1084
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:1002
void build(const LinearAlgebraDistribution *distribution_pt, const unsigned &ncol, const Vector< double > &value, const Vector< int > &column_index, const Vector< int > &row_start)
build method: vector of values, vector of column indices, vector of row starts and number of rows and...
Definition: matrices.cc:1672
bool distributed() const
distribution is serial or distributed
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
unsigned nrow() const
access function to the number of global rows.
unsigned nrow_local() const
access function for the num of local rows on this processor.
unsigned first_row() const
access function for the first row on this processor
Class to allow sorting of column indices in conversion to epetra matrix.
DistributionPredicate(const int &first_col, const int &ncol_local)
Constructor: Pass number of first column and the number of local columns.
bool operator()(const int &col)
Comparison operator: is column col in the range between (including) First_col and Last_col.
int First_col
First column held locally.
int Last_col
Last colum held locally.
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
double * values_pt()
access function to the underlying values
bool built() const
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
bool distributed() const
access function to the distributed - indicates whether the distribution is serial or distributed
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero.
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
bool built() const
if the communicator_pt is null then the distribution is not setup then false is returned,...
unsigned nrow() const
access function to the number of global rows.
unsigned nrow_local() const
access function for the num of local rows on this processor. If no MPI then Nrow is returned.
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
An OomphLibError object which should be thrown when an run-time error is encountered....
Epetra_CrsMatrix * create_distributed_epetra_matrix(const CRDoubleMatrix *oomph_matrix_pt, const LinearAlgebraDistribution *dist_pt)
create an Epetra_CrsMatrix from an oomph-lib CRDoubleMatrix. If oomph_matrix_pt is NOT distributed (i...
Epetra_CrsMatrix * create_distributed_epetra_matrix_for_aztecoo(CRDoubleMatrix *oomph_matrix_pt)
create and Epetra_CrsMatrix from an oomph-lib CRDoubleMatrix. Specialisation for Trilinos AztecOO....
void copy_to_oomphlib_vector(const Epetra_Vector *epetra_vec_pt, DoubleVector &oomph_vec)
Helper function to copy the contents of a Trilinos vector to an oomph-lib distributed vector....
void multiply(const CRDoubleMatrix *matrix, const DoubleVector &x, DoubleVector &soln)
Function to perform a matrix-vector multiplication on a oomph-lib matrix and vector using Trilinos fu...
Epetra_Map * create_epetra_map(const LinearAlgebraDistribution *const dist)
create an Epetra_Map corresponding to the LinearAlgebraDistribution
Epetra_Vector * create_distributed_epetra_vector(const DoubleVector &oomph_vec)
create an Epetra_Vector from an oomph-lib DoubleVector. If oomph_vec is NOT distributed (i....
Epetra_Vector * create_epetra_vector_view_data(DoubleVector &oomph_vec)
create an Epetra_Vector equivalent of DoubleVector The argument DoubleVector must be built....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...