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