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-2022 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
28namespace 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
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 //============================================================================
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 =
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 * 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
int * column_index()
Access to C-style column index array.
Definition: matrices.h:1072
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
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
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
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.
bool built() const
double * values_pt()
access function to the underlying values
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....
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
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....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...