EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
EpetraExt_MatrixMatrix.cpp
Go to the documentation of this file.
1//@HEADER
2// ***********************************************************************
3//
4// EpetraExt: Epetra Extended - Linear Algebra Services Package
5// Copyright (2011) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41
44
45#include <EpetraExt_MMHelpers.h>
46
48
49#include <Epetra_Export.h>
50#include <Epetra_Import.h>
51#include <Epetra_Util.h>
52#include <Epetra_Map.h>
53#include <Epetra_Comm.h>
54#include <Epetra_CrsMatrix.h>
55#include <Epetra_Vector.h>
56#include <Epetra_Directory.h>
57#include <Epetra_HashTable.h>
58#include <Epetra_Distributor.h>
59
60#include <Teuchos_TimeMonitor.hpp>
61
62#ifdef HAVE_VECTOR
63#include <vector>
64#endif
65
66
67
68namespace EpetraExt {
69
70//=========================================================================
71// ETI
72template<int> int import_only(const Epetra_CrsMatrix& M,
73 const Epetra_Map& targetMap,
74 CrsMatrixStruct& Mview,
75 const Epetra_Import * prototypeImporter);
76
77template<long long> int import_only(const Epetra_CrsMatrix& M,
78 const Epetra_Map& targetMap,
79 CrsMatrixStruct& Mview,
80 const Epetra_Import * prototypeImporter);
81
82
83//=========================================================================
84//
85//Method for internal use... sparsedot forms a dot-product between two
86//sparsely-populated 'vectors'.
87//Important assumption: assumes the indices in u_ind and v_ind are sorted.
88//
89template<typename int_type>
90double sparsedot(double* u, int_type* u_ind, int u_len,
91 double* v, int_type* v_ind, int v_len)
92{
93 double result = 0.0;
94
95 int v_idx = 0;
96 int u_idx = 0;
97
98 while(v_idx < v_len && u_idx < u_len) {
99 int_type ui = u_ind[u_idx];
100 int_type vi = v_ind[v_idx];
101
102 if (ui < vi) {
103 ++u_idx;
104 }
105 else if (ui > vi) {
106 ++v_idx;
107 }
108 else {
109 result += u[u_idx++]*v[v_idx++];
110 }
111 }
112
113 return(result);
114}
115
116//=========================================================================
117//kernel method for computing the local portion of C = A*B^T
118template<typename int_type>
120 CrsMatrixStruct& Bview,
121 CrsWrapper& C,
122 bool keep_all_hard_zeros)
123{
124 int i, j, k;
125 int returnValue = 0;
126
127 int maxlen = 0;
128 for(i=0; i<Aview.numRows; ++i) {
129 if (Aview.numEntriesPerRow[i] > maxlen) maxlen = Aview.numEntriesPerRow[i];
130 }
131 for(i=0; i<Bview.numRows; ++i) {
132 if (Bview.numEntriesPerRow[i] > maxlen) maxlen = Bview.numEntriesPerRow[i];
133 }
134
135 //std::cout << "Aview: " << std::endl;
136 //dumpCrsMatrixStruct(Aview);
137
138 //std::cout << "Bview: " << std::endl;
139 //dumpCrsMatrixStruct(Bview);
140
141 int numBcols = Bview.colMap->NumMyElements();
142 int numBrows = Bview.numRows;
143
144 int iworklen = maxlen*2 + numBcols;
145 int_type* iwork = new int_type[iworklen];
146
147 int_type * bcols = iwork+maxlen*2;
148 int_type* bgids = 0;
149 Bview.colMap->MyGlobalElementsPtr(bgids);
150 double* bvals = new double[maxlen*2];
151 double* avals = bvals+maxlen;
152
153 int_type max_all_b = (int_type) Bview.colMap->MaxAllGID64();
154 int_type min_all_b = (int_type) Bview.colMap->MinAllGID64();
155
156 //bcols will hold the GIDs from B's column-map for fast access
157 //during the computations below
158 for(i=0; i<numBcols; ++i) {
159 int blid = Bview.colMap->LID(bgids[i]);
160 bcols[blid] = bgids[i];
161 }
162
163 //next create arrays indicating the first and last column-index in
164 //each row of B, so that we can know when to skip certain rows below.
165 //This will provide a large performance gain for banded matrices, and
166 //a somewhat smaller gain for *most* other matrices.
167 int_type* b_firstcol = new int_type[2*numBrows];
168 int_type* b_lastcol = b_firstcol+numBrows;
169 int_type temp;
170 for(i=0; i<numBrows; ++i) {
171 b_firstcol[i] = max_all_b;
172 b_lastcol[i] = min_all_b;
173
174 int Blen_i = Bview.numEntriesPerRow[i];
175 if (Blen_i < 1) continue;
176 int* Bindices_i = Bview.indices[i];
177
178 if (Bview.remote[i]) {
179 for(k=0; k<Blen_i; ++k) {
180 temp = (int_type) Bview.importColMap->GID64(Bindices_i[k]);
181 if (temp < b_firstcol[i]) b_firstcol[i] = temp;
182 if (temp > b_lastcol[i]) b_lastcol[i] = temp;
183 }
184 }
185 else {
186 for(k=0; k<Blen_i; ++k) {
187 temp = bcols[Bindices_i[k]];
188 if (temp < b_firstcol[i]) b_firstcol[i] = temp;
189 if (temp > b_lastcol[i]) b_lastcol[i] = temp;
190 }
191 }
192 }
193
194 Epetra_Util util;
195
196 int_type* Aind = iwork;
197 int_type* Bind = iwork+maxlen;
198
199 bool C_filled = C.Filled();
200
201 //To form C = A*B^T, we're going to execute this expression:
202 //
203 // C(i,j) = sum_k( A(i,k)*B(j,k) )
204 //
205 //This is the easiest case of all to code (easier than A*B, A^T*B, A^T*B^T).
206 //But it requires the use of a 'sparsedot' function (we're simply forming
207 //dot-products with row A_i and row B_j for all i and j).
208
209 //loop over the rows of A.
210 for(i=0; i<Aview.numRows; ++i) {
211 if (Aview.remote[i]) {
212 continue;
213 }
214
215 int* Aindices_i = Aview.indices[i];
216 double* Aval_i = Aview.values[i];
217 int A_len_i = Aview.numEntriesPerRow[i];
218 if (A_len_i < 1) {
219 continue;
220 }
221
222 for(k=0; k<A_len_i; ++k) {
223 Aind[k] = (int_type) Aview.colMap->GID64(Aindices_i[k]);
224 avals[k] = Aval_i[k];
225 }
226
227 util.Sort<int_type>(true, A_len_i, Aind, 1, &avals, 0, NULL, 0, 0);
228
229 int_type mina = Aind[0];
230 int_type maxa = Aind[A_len_i-1];
231
232 if (mina > max_all_b || maxa < min_all_b) {
233 continue;
234 }
235
236 int_type global_row = (int_type) Aview.rowMap->GID64(i);
237
238 //loop over the rows of B and form results C_ij = dot(A(i,:),B(j,:))
239 for(j=0; j<Bview.numRows; ++j) {
240 if (b_firstcol[j] > maxa || b_lastcol[j] < mina) {
241 continue;
242 }
243
244 int* Bindices_j = Bview.indices[j];
245 int B_len_j = Bview.numEntriesPerRow[j];
246 if (B_len_j < 1) {
247 continue;
248 }
249
250 int_type tmp;
251 int Blen = 0;
252
253 if (Bview.remote[j]) {
254 for(k=0; k<B_len_j; ++k) {
255 tmp = (int_type) Bview.importColMap->GID64(Bindices_j[k]);
256 if (tmp < mina || tmp > maxa) {
257 continue;
258 }
259
260 bvals[Blen] = Bview.values[j][k];
261 Bind[Blen++] = tmp;
262 }
263 }
264 else {
265 for(k=0; k<B_len_j; ++k) {
266 tmp = bcols[Bindices_j[k]];
267 if (tmp < mina || tmp > maxa) {
268 continue;
269 }
270
271 bvals[Blen] = Bview.values[j][k];
272 Bind[Blen++] = tmp;
273 }
274 }
275
276 if (Blen < 1) {
277 continue;
278 }
279
280 util.Sort<int_type>(true, Blen, Bind, 1, &bvals, 0, NULL, 0, 0);
281
282 double C_ij = sparsedot(avals, Aind, A_len_i,
283 bvals, Bind, Blen);
284
285 if (!keep_all_hard_zeros && C_ij == 0.0)
286 continue;
287
288 int_type global_col = (int_type) Bview.rowMap->GID64(j);
289
290 int err = C_filled ?
291 C.SumIntoGlobalValues(global_row, 1, &C_ij, &global_col)
292 :
293 C.InsertGlobalValues(global_row, 1, &C_ij, &global_col);
294
295 if (err < 0) {
296 return(err);
297 }
298 if (err > 0) {
299 if (C_filled) {
300 //C.Filled()==true, and C doesn't have all the necessary nonzero
301 //locations, or global_row or global_col is out of range (less
302 //than 0 or non local).
303 std::cerr << "EpetraExt::MatrixMatrix::Multiply Warning: failed "
304 << "to insert value in result matrix at position "<<global_row
305 <<"," <<global_col<<", possibly because result matrix has a "
306 << "column-map that doesn't include column "<<global_col
307 <<" on this proc." <<std::endl;
308 return(err);
309 }
310 }
311 }
312 }
313
314 delete [] iwork;
315 delete [] bvals;
316 delete [] b_firstcol;
317
318 return(returnValue);
319}
320
322 CrsMatrixStruct& Bview,
323 CrsWrapper& C,
324 bool keep_all_hard_zeros)
325{
326#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
327 if(Aview.rowMap->GlobalIndicesInt() &&
328 Aview.colMap->GlobalIndicesInt() &&
329 Aview.rowMap->GlobalIndicesInt() &&
330 Aview.colMap->GlobalIndicesInt()) {
331 return mult_A_Btrans<int>(Aview, Bview, C, keep_all_hard_zeros);
332 }
333 else
334#endif
335#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
336 if(Aview.rowMap->GlobalIndicesLongLong() &&
337 Aview.colMap->GlobalIndicesLongLong() &&
338 Aview.rowMap->GlobalIndicesLongLong() &&
339 Aview.colMap->GlobalIndicesLongLong()) {
340 return mult_A_Btrans<long long>(Aview, Bview, C, keep_all_hard_zeros);
341 }
342 else
343#endif
344 throw std::runtime_error("EpetraExt::mult_A_Btrans: GlobalIndices type unknown");
345}
346
347//=========================================================================
348//kernel method for computing the local portion of C = A^T*B
349template<typename int_type>
351 CrsMatrixStruct& Bview,
352 CrsWrapper& C)
353{
354
355 int C_firstCol = Bview.colMap->MinLID();
356 int C_lastCol = Bview.colMap->MaxLID();
357
358 int C_firstCol_import = 0;
359 int C_lastCol_import = -1;
360
361 if (Bview.importColMap != NULL) {
362 C_firstCol_import = Bview.importColMap->MinLID();
363 C_lastCol_import = Bview.importColMap->MaxLID();
364 }
365
366 int C_numCols = C_lastCol - C_firstCol + 1;
367 int C_numCols_import = C_lastCol_import - C_firstCol_import + 1;
368
369 if (C_numCols_import > C_numCols) C_numCols = C_numCols_import;
370
371 double* C_row_i = new double[C_numCols];
372 int_type* C_colInds = new int_type[C_numCols];
373
374 int i, j, k;
375
376 for(j=0; j<C_numCols; ++j) {
377 C_row_i[j] = 0.0;
378 C_colInds[j] = 0;
379 }
380
381 //To form C = A^T*B, compute a series of outer-product updates.
382 //
383 // for (ith column of A^T) {
384 // C_i = outer product of A^T(:,i) and B(i,:)
385 // Where C_i is the ith matrix update,
386 // A^T(:,i) is the ith column of A^T, and
387 // B(i,:) is the ith row of B.
388 //
389
390 //dumpCrsMatrixStruct(Aview);
391 //dumpCrsMatrixStruct(Bview);
392 int localProc = Bview.colMap->Comm().MyPID();
393
394 int_type* Arows = 0;
395 Aview.rowMap->MyGlobalElementsPtr(Arows);
396
397 bool C_filled = C.Filled();
398
399 //loop over the rows of A (which are the columns of A^T).
400 for(i=0; i<Aview.numRows; ++i) {
401
402 int* Aindices_i = Aview.indices[i];
403 double* Aval_i = Aview.values[i];
404
405 //we'll need to get the row of B corresponding to Arows[i],
406 //where Arows[i] is the GID of A's ith row.
407 int Bi = Bview.rowMap->LID(Arows[i]);
408 if (Bi<0) {
409 std::cout << "mult_Atrans_B ERROR, proc "<<localProc<<" needs row "
410 <<Arows[i]<<" of matrix B, but doesn't have it."<<std::endl;
411 return(-1);
412 }
413
414 int* Bcol_inds = Bview.indices[Bi];
415 double* Bvals_i = Bview.values[Bi];
416
417 //for each column-index Aj in the i-th row of A, we'll update
418 //global-row GID(Aj) of the result matrix C. In that row of C,
419 //we'll update column-indices given by the column-indices in the
420 //ith row of B that we're now holding (Bcol_inds).
421
422 //First create a list of GIDs for the column-indices
423 //that we'll be updating.
424
425 int Blen = Bview.numEntriesPerRow[Bi];
426 if (Bview.remote[Bi]) {
427 for(j=0; j<Blen; ++j) {
428 C_colInds[j] = (int_type) Bview.importColMap->GID64(Bcol_inds[j]);
429 }
430 }
431 else {
432 for(j=0; j<Blen; ++j) {
433 C_colInds[j] = (int_type) Bview.colMap->GID64(Bcol_inds[j]);
434 }
435 }
436
437 //loop across the i-th row of A (column of A^T)
438 for(j=0; j<Aview.numEntriesPerRow[i]; ++j) {
439
440 int Aj = Aindices_i[j];
441 double Aval = Aval_i[j];
442
443 int_type global_row;
444 if (Aview.remote[i]) {
445 global_row = (int_type) Aview.importColMap->GID64(Aj);
446 }
447 else {
448 global_row = (int_type) Aview.colMap->GID64(Aj);
449 }
450
451 if (!C.RowMap().MyGID(global_row)) {
452 continue;
453 }
454
455 for(k=0; k<Blen; ++k) {
456 C_row_i[k] = Aval*Bvals_i[k];
457 }
458
459 //
460 //Now add this row-update to C.
461 //
462
463 int err = C_filled ?
464 C.SumIntoGlobalValues(global_row, Blen, C_row_i, C_colInds)
465 :
466 C.InsertGlobalValues(global_row, Blen, C_row_i, C_colInds);
467
468 if (err < 0) {
469 return(err);
470 }
471 if (err > 0) {
472 if (C_filled) {
473 //C is Filled, and doesn't have all the necessary nonzero locations.
474 return(err);
475 }
476 }
477 }
478 }
479
480 delete [] C_row_i;
481 delete [] C_colInds;
482
483 return(0);
484}
485
487 CrsMatrixStruct& Bview,
488 CrsWrapper& C)
489{
490#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
491 if(Aview.rowMap->GlobalIndicesInt() &&
492 Aview.colMap->GlobalIndicesInt() &&
493 Aview.rowMap->GlobalIndicesInt() &&
494 Aview.colMap->GlobalIndicesInt()) {
495 return mult_Atrans_B<int>(Aview, Bview, C);
496 }
497 else
498#endif
499#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
500 if(Aview.rowMap->GlobalIndicesLongLong() &&
501 Aview.colMap->GlobalIndicesLongLong() &&
502 Aview.rowMap->GlobalIndicesLongLong() &&
503 Aview.colMap->GlobalIndicesLongLong()) {
504 return mult_Atrans_B<long long>(Aview, Bview, C);
505 }
506 else
507#endif
508 throw std::runtime_error("EpetraExt::mult_Atrans_B: GlobalIndices type unknown");
509}
510
511//kernel method for computing the local portion of C = A^T*B^T
512template<typename int_type>
514 CrsMatrixStruct& Bview,
515 CrsWrapper& C,
516 bool keep_all_hard_zeros)
517{
518 int C_firstCol = Aview.rowMap->MinLID();
519 int C_lastCol = Aview.rowMap->MaxLID();
520
521 int C_firstCol_import = 0;
522 int C_lastCol_import = -1;
523
524 if (Aview.importColMap != NULL) {
525 C_firstCol_import = Aview.importColMap->MinLID();
526 C_lastCol_import = Aview.importColMap->MaxLID();
527 }
528
529 int C_numCols = C_lastCol - C_firstCol + 1;
530 int C_numCols_import = C_lastCol_import - C_firstCol_import + 1;
531
532 if (C_numCols_import > C_numCols) C_numCols = C_numCols_import;
533
534 double* dwork = new double[C_numCols];
535 int_type* iwork = new int_type[C_numCols];
536
537 double* C_col_j = dwork;
538 int_type* C_inds = iwork;
539
540 int i, j, k;
541
542 for(j=0; j<C_numCols; ++j) {
543 C_col_j[j] = 0.0;
544 C_inds[j] = -1;
545 }
546
547 int_type* A_col_inds = 0;
548 Aview.colMap->MyGlobalElementsPtr(A_col_inds);
549 int_type* A_col_inds_import = 0;
550 if(Aview.importColMap)
551 Aview.importColMap->MyGlobalElementsPtr(A_col_inds_import);
552
553 const Epetra_Map* Crowmap = &(C.RowMap());
554
555 //To form C = A^T*B^T, we're going to execute this expression:
556 //
557 // C(i,j) = sum_k( A(k,i)*B(j,k) )
558 //
559 //Our goal, of course, is to navigate the data in A and B once, without
560 //performing searches for column-indices, etc. In other words, we avoid
561 //column-wise operations like the plague...
562
563 int_type* Brows = 0;
564 Bview.rowMap->MyGlobalElementsPtr(Brows);
565
566 //loop over the rows of B
567 for(j=0; j<Bview.numRows; ++j) {
568 int* Bindices_j = Bview.indices[j];
569 double* Bvals_j = Bview.values[j];
570
571 int_type global_col = Brows[j];
572
573 //loop across columns in the j-th row of B and for each corresponding
574 //row in A, loop across columns and accumulate product
575 //A(k,i)*B(j,k) into our partial sum quantities in C_col_j. In other
576 //words, as we stride across B(j,:), we use selected rows in A to
577 //calculate updates for column j of the result matrix C.
578
579 for(k=0; k<Bview.numEntriesPerRow[j]; ++k) {
580
581 int bk = Bindices_j[k];
582 double Bval = Bvals_j[k];
583
584 int_type global_k;
585 if (Bview.remote[j]) {
586 global_k = (int_type) Bview.importColMap->GID64(bk);
587 }
588 else {
589 global_k = (int_type) Bview.colMap->GID64(bk);
590 }
591
592 //get the corresponding row in A
593 int ak = Aview.rowMap->LID(global_k);
594 if (ak<0) {
595 continue;
596 }
597
598 int* Aindices_k = Aview.indices[ak];
599 double* Avals_k = Aview.values[ak];
600
601 int C_len = 0;
602
603 if (Aview.remote[ak]) {
604 for(i=0; i<Aview.numEntriesPerRow[ak]; ++i) {
605 C_col_j[C_len] = Avals_k[i]*Bval;
606 C_inds[C_len++] = A_col_inds_import[Aindices_k[i]];
607 }
608 }
609 else {
610 for(i=0; i<Aview.numEntriesPerRow[ak]; ++i) {
611 C_col_j[C_len] = Avals_k[i]*Bval;
612 C_inds[C_len++] = A_col_inds[Aindices_k[i]];
613 }
614 }
615
616 //Now loop across the C_col_j values and put non-zeros into C.
617
618 for(i=0; i < C_len ; ++i) {
619 if (!keep_all_hard_zeros && C_col_j[i] == 0.0) continue;
620
621 int_type global_row = C_inds[i];
622 if (!Crowmap->MyGID(global_row)) {
623 continue;
624 }
625
626 int err = C.SumIntoGlobalValues(global_row, 1, &(C_col_j[i]), &global_col);
627
628 if (err < 0) {
629 return(err);
630 }
631 else {
632 if (err > 0) {
633 err = C.InsertGlobalValues(global_row, 1, &(C_col_j[i]), &global_col);
634 if (err < 0) {
635 return(err);
636 }
637 }
638 }
639 }
640
641 }
642 }
643
644 delete [] dwork;
645 delete [] iwork;
646
647 return(0);
648}
649
651 CrsMatrixStruct& Bview,
652 CrsWrapper& C,
653 bool keep_all_hard_zeros)
654{
655#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
656 if(Aview.rowMap->GlobalIndicesInt() &&
657 Aview.colMap->GlobalIndicesInt() &&
658 Aview.rowMap->GlobalIndicesInt() &&
659 Aview.colMap->GlobalIndicesInt()) {
660 return mult_Atrans_Btrans<int>(Aview, Bview, C, keep_all_hard_zeros);
661 }
662 else
663#endif
664#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
665 if(Aview.rowMap->GlobalIndicesLongLong() &&
666 Aview.colMap->GlobalIndicesLongLong() &&
667 Aview.rowMap->GlobalIndicesLongLong() &&
668 Aview.colMap->GlobalIndicesLongLong()) {
669 return mult_Atrans_Btrans<long long>(Aview, Bview, C, keep_all_hard_zeros);
670 }
671 else
672#endif
673 throw std::runtime_error("EpetraExt::mult_Atrans_Btrans: GlobalIndices type unknown");
674}
675
676// ==============================================================
677template<typename int_type>
679 const Epetra_Map& targetMap,
680 CrsMatrixStruct& Mview,
681 const Epetra_Import * prototypeImporter=0)
682{
683 //The goal of this method is to populate the 'Mview' struct with views of the
684 //rows of M, including all rows that correspond to elements in 'targetMap'.
685 //
686 //If targetMap includes local elements that correspond to remotely-owned rows
687 //of M, then those remotely-owned rows will be imported into
688 //'Mview.importMatrix', and views of them will be included in 'Mview'.
689
690 // The prototype importer, if used, has to know who owns all of the PIDs in M's rowmap.
691 // The typical use of this is to provide A's column map when I&XV is called for B, for
692 // a C = A * B multiply.
693
694#ifdef ENABLE_MMM_TIMINGS
695 using Teuchos::TimeMonitor;
696 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: MMM I&X Alloc")));
697#endif
698 Mview.deleteContents();
699
700 Mview.origMatrix = &M;
701 const Epetra_Map& Mrowmap = M.RowMap();
702 int numProcs = Mrowmap.Comm().NumProc();
703 Mview.numRows = targetMap.NumMyElements();
704 int_type* Mrows = 0;
705 targetMap.MyGlobalElementsPtr(Mrows);
706
707 if (Mview.numRows > 0) {
708 Mview.numEntriesPerRow = new int[Mview.numRows];
709 Mview.indices = new int*[Mview.numRows];
710 Mview.values = new double*[Mview.numRows];
711 Mview.remote = new bool[Mview.numRows];
712 }
713
714 Mview.origRowMap = &(M.RowMap());
715 Mview.rowMap = &targetMap;
716 Mview.colMap = &(M.ColMap());
717 Mview.domainMap = &(M.DomainMap());
718 Mview.importColMap = NULL;
719 Mview.numRemote = 0;
720
721
722#ifdef ENABLE_MMM_TIMINGS
723 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: MMM I&X Extract")));
724#endif
725
726 int *rowptr=0, *colind=0;
727 double *vals=0;
728
729 EPETRA_CHK_ERR( M.ExtractCrsDataPointers(rowptr,colind,vals) );
730
731 if(Mrowmap.SameAs(targetMap)) {
732 // Short Circuit: The Row and Target Maps are the Same
733 for(int i=0; i<Mview.numRows; ++i) {
734 Mview.numEntriesPerRow[i] = rowptr[i+1]-rowptr[i];
735 Mview.indices[i] = colind + rowptr[i];
736 Mview.values[i] = vals + rowptr[i];
737 Mview.remote[i] = false;
738 }
739 return 0;
740 }
741 else if(prototypeImporter && prototypeImporter->SourceMap().SameAs(M.RowMap()) && prototypeImporter->TargetMap().SameAs(targetMap)){
742 // The prototypeImporter can tell me what is local and what is remote
743 int * PermuteToLIDs = prototypeImporter->PermuteToLIDs();
744 int * PermuteFromLIDs = prototypeImporter->PermuteFromLIDs();
745 int * RemoteLIDs = prototypeImporter->RemoteLIDs();
746 for(int i=0; i<prototypeImporter->NumSameIDs();i++){
747 Mview.numEntriesPerRow[i] = rowptr[i+1]-rowptr[i];
748 Mview.indices[i] = colind + rowptr[i];
749 Mview.values[i] = vals + rowptr[i];
750 Mview.remote[i] = false;
751 }
752 for(int i=0; i<prototypeImporter->NumPermuteIDs();i++){
753 int to = PermuteToLIDs[i];
754 int from = PermuteFromLIDs[i];
755 Mview.numEntriesPerRow[to] = rowptr[from+1]-rowptr[from];
756 Mview.indices[to] = colind + rowptr[from];
757 Mview.values[to] = vals + rowptr[from];
758 Mview.remote[to] = false;
759 }
760 for(int i=0; i<prototypeImporter->NumRemoteIDs();i++){
761 Mview.remote[RemoteLIDs[i]] = true;
762 ++Mview.numRemote;
763 }
764 }
765 else {
766 // Only LID can tell me who is local and who is remote
767 for(int i=0; i<Mview.numRows; ++i) {
768 int mlid = Mrowmap.LID(Mrows[i]);
769 if (mlid < 0) {
770 Mview.remote[i] = true;
771 ++Mview.numRemote;
772 }
773 else {
774 Mview.numEntriesPerRow[i] = rowptr[mlid+1]-rowptr[mlid];
775 Mview.indices[i] = colind + rowptr[mlid];
776 Mview.values[i] = vals + rowptr[mlid];
777 Mview.remote[i] = false;
778 }
779 }
780 }
781
782#ifdef ENABLE_MMM_TIMINGS
783 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: MMM I&X Collective-0")));
784#endif
785
786 if (numProcs < 2) {
787 if (Mview.numRemote > 0) {
788 std::cerr << "EpetraExt::MatrixMatrix::Multiply ERROR, numProcs < 2 but "
789 << "attempting to import remote matrix rows."<<std::endl;
790 return(-1);
791 }
792
793 //If only one processor we don't need to import any remote rows, so return.
794 return(0);
795 }
796
797 //
798 //Now we will import the needed remote rows of M, if the global maximum
799 //value of numRemote is greater than 0.
800 //
801
802 int globalMaxNumRemote = 0;
803 Mrowmap.Comm().MaxAll(&Mview.numRemote, &globalMaxNumRemote, 1);
804
805 if (globalMaxNumRemote > 0) {
806#ifdef ENABLE_MMM_TIMINGS
807 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: MMM I&X Import-1")));
808#endif
809 //Create a map that describes the remote rows of M that we need.
810
811 int_type* MremoteRows = Mview.numRemote>0 ? new int_type[Mview.numRemote] : NULL;
812 int offset = 0;
813 for(int i=0; i<Mview.numRows; ++i) {
814 if (Mview.remote[i]) {
815 MremoteRows[offset++] = Mrows[i];
816 }
817 }
818
819 LightweightMap MremoteRowMap((int_type) -1, Mview.numRemote, MremoteRows, (int_type) Mrowmap.IndexBase64());
820
821#ifdef ENABLE_MMM_TIMINGS
822 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: MMM I&X Import-2")));
823#endif
824 //Create an importer with target-map MremoteRowMap and
825 //source-map Mrowmap.
826 Epetra_Import * importer=0;
827 RemoteOnlyImport * Rimporter=0;
828 if(prototypeImporter && prototypeImporter->SourceMap().SameAs(M.RowMap()) && prototypeImporter->TargetMap().SameAs(targetMap)) {
829 Rimporter = new RemoteOnlyImport(*prototypeImporter,MremoteRowMap);
830 }
831 else if(!prototypeImporter) {
832 Epetra_Map MremoteRowMap2((int_type) -1, Mview.numRemote, MremoteRows, (int_type) Mrowmap.IndexBase64(), Mrowmap.Comm());
833 importer=new Epetra_Import(MremoteRowMap2, Mrowmap);
834 }
835 else
836 throw std::runtime_error("prototypeImporter->SourceMap() does not match M.RowMap()!");
837
838
839#ifdef ENABLE_MMM_TIMINGS
840 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: MMM I&X Import-3")));
841#endif
842
843 if(Mview.importMatrix) delete Mview.importMatrix;
844 if(Rimporter) Mview.importMatrix = new LightweightCrsMatrix(M,*Rimporter);
845 else Mview.importMatrix = new LightweightCrsMatrix(M,*importer);
846
847#ifdef ENABLE_MMM_TIMINGS
848 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: MMM I&X Import-4")));
849#endif
850
851 //Finally, use the freshly imported data to fill in the gaps in our views
852 int N;
853 if (Mview.importMatrix->use_lw) {
854 N = Mview.importMatrix->RowMapLW_->NumMyElements();
855 } else {
856 N = Mview.importMatrix->RowMapEP_->NumMyElements();
857 }
858
859 if(N > 0) {
860 rowptr = &Mview.importMatrix->rowptr_[0];
861 colind = &Mview.importMatrix->colind_[0];
862 vals = &Mview.importMatrix->vals_[0];
863 }
864
865
866 for(int i=0; i<Mview.numRows; ++i) {
867 if (Mview.remote[i]) {
868 int importLID = MremoteRowMap.LID(Mrows[i]);
869 Mview.numEntriesPerRow[i] = rowptr[importLID+1]-rowptr[importLID];
870 Mview.indices[i] = colind + rowptr[importLID];
871 Mview.values[i] = vals + rowptr[importLID];
872 }
873 }
874
875
876 int_type * MyColGIDs = 0;
877 if(Mview.importMatrix->ColMap_.NumMyElements()>0) {
878 Mview.importMatrix->ColMap_.MyGlobalElementsPtr(MyColGIDs);
879 }
880 Mview.importColMap =
881 new Epetra_Map (static_cast<int_type> (-1),
882 Mview.importMatrix->ColMap_.NumMyElements (),
883 MyColGIDs,
884 static_cast<int_type> (Mview.importMatrix->ColMap_.IndexBase64 ()),
885 M.Comm ());
886 delete [] MremoteRows;
887#ifdef ENABLE_MMM_TIMINGS
888 MM=Teuchos::null;
889#endif
890
891 // Cleanup
892 delete Rimporter;
893 delete importer;
894 }
895 return(0);
896}
897
898
899
900
901
902
903//=========================================================================
904template<typename int_type>
906 const Epetra_Map* map2,
907 const Epetra_Map*& mapunion)
908{
909 //form the union of two maps
910
911 if (map1 == NULL) {
912 mapunion = new Epetra_Map(*map2);
913 return(0);
914 }
915
916 if (map2 == NULL) {
917 mapunion = new Epetra_Map(*map1);
918 return(0);
919 }
920
921 int map1_len = map1->NumMyElements();
922 int_type* map1_elements = 0;
923 map1->MyGlobalElementsPtr(map1_elements);
924 int map2_len = map2->NumMyElements();
925 int_type* map2_elements = 0;
926 map2->MyGlobalElementsPtr(map2_elements);
927
928 int_type* union_elements = new int_type[map1_len+map2_len];
929
930 int map1_offset = 0, map2_offset = 0, union_offset = 0;
931
932 while(map1_offset < map1_len && map2_offset < map2_len) {
933 int_type map1_elem = map1_elements[map1_offset];
934 int_type map2_elem = map2_elements[map2_offset];
935
936 if (map1_elem < map2_elem) {
937 union_elements[union_offset++] = map1_elem;
938 ++map1_offset;
939 }
940 else if (map1_elem > map2_elem) {
941 union_elements[union_offset++] = map2_elem;
942 ++map2_offset;
943 }
944 else {
945 union_elements[union_offset++] = map1_elem;
946 ++map1_offset;
947 ++map2_offset;
948 }
949 }
950
951 int i;
952 for(i=map1_offset; i<map1_len; ++i) {
953 union_elements[union_offset++] = map1_elements[i];
954 }
955
956 for(i=map2_offset; i<map2_len; ++i) {
957 union_elements[union_offset++] = map2_elements[i];
958 }
959
960 mapunion = new Epetra_Map((int_type) -1, union_offset, union_elements,
961 (int_type) map1->IndexBase64(), map1->Comm());
962
963 delete [] union_elements;
964
965 return(0);
966}
967
968//=========================================================================
969template<typename int_type>
971 const Epetra_Map& column_map)
972{
973 //The goal of this function is to find all rows in the matrix M that contain
974 //column-indices which are in 'column_map'. A map containing those rows is
975 //returned. (The returned map will then be used as the source row-map for
976 //importing those rows into an overlapping distribution.)
977
978 std::pair<int_type,int_type> my_col_range = get_col_range<int_type>(column_map);
979
980 const Epetra_Comm& Comm = M.RowMap().Comm();
981 int num_procs = Comm.NumProc();
982 int my_proc = Comm.MyPID();
983 std::vector<int> send_procs;
984 send_procs.reserve(num_procs-1);
985 std::vector<int_type> col_ranges;
986 col_ranges.reserve(2*(num_procs-1));
987 for(int p=0; p<num_procs; ++p) {
988 if (p == my_proc) continue;
989 send_procs.push_back(p);
990 col_ranges.push_back(my_col_range.first);
991 col_ranges.push_back(my_col_range.second);
992 }
993
994 Epetra_Distributor* distor = Comm.CreateDistributor();
995
996 int num_recv_procs = 0;
997 int num_send_procs = send_procs.size();
998 bool deterministic = true;
999 int* send_procs_ptr = send_procs.size()>0 ? &send_procs[0] : NULL;
1000 distor->CreateFromSends(num_send_procs, send_procs_ptr, deterministic, num_recv_procs);
1001
1002 int len_import_chars = 0;
1003 char* import_chars = NULL;
1004
1005 char* export_chars = col_ranges.size()>0 ? reinterpret_cast<char*>(&col_ranges[0]) : NULL;
1006 int err = distor->Do(export_chars, 2*sizeof(int_type), len_import_chars, import_chars);
1007 if (err != 0) {
1008 std::cout << "ERROR returned from Distributor::Do."<<std::endl;
1009 }
1010
1011 int_type* IntImports = reinterpret_cast<int_type*>(import_chars);
1012 int num_import_pairs = len_import_chars/(2*sizeof(int_type));
1013 int offset = 0;
1014 std::vector<int> send_procs2;
1015 send_procs2.reserve(num_procs);
1016 std::vector<int_type> proc_col_ranges;
1017 std::pair<int_type,int_type> M_col_range = get_col_range<int_type>(M);
1018 for(int i=0; i<num_import_pairs; ++i) {
1019 int_type first_col = IntImports[offset++];
1020 int_type last_col = IntImports[offset++];
1021 if (first_col <= M_col_range.second && last_col >= M_col_range.first) {
1022 send_procs2.push_back(send_procs[i]);
1023 proc_col_ranges.push_back(first_col);
1024 proc_col_ranges.push_back(last_col);
1025 }
1026 }
1027
1028 std::vector<int_type> send_rows;
1029 std::vector<int> rows_per_send_proc;
1030 pack_outgoing_rows(M, proc_col_ranges, send_rows, rows_per_send_proc);
1031
1032 Epetra_Distributor* distor2 = Comm.CreateDistributor();
1033
1034 int* send_procs2_ptr = send_procs2.size()>0 ? &send_procs2[0] : NULL;
1035 num_recv_procs = 0;
1036 err = distor2->CreateFromSends(send_procs2.size(), send_procs2_ptr, deterministic, num_recv_procs);
1037
1038 export_chars = send_rows.size()>0 ? reinterpret_cast<char*>(&send_rows[0]) : NULL;
1039 int* rows_per_send_proc_ptr = rows_per_send_proc.size()>0 ? &rows_per_send_proc[0] : NULL;
1040 len_import_chars = 0;
1041 err = distor2->Do(export_chars, (int)sizeof(int_type), rows_per_send_proc_ptr, len_import_chars, import_chars);
1042
1043 int_type* recvd_row_ints = reinterpret_cast<int_type*>(import_chars);
1044 int num_recvd_rows = len_import_chars/sizeof(int_type);
1045
1046 Epetra_Map recvd_rows((int_type) -1, num_recvd_rows, recvd_row_ints, (int_type) column_map.IndexBase64(), Comm);
1047
1048 delete distor;
1049 delete distor2;
1050 delete [] import_chars;
1051
1052 Epetra_Map* result_map = NULL;
1053
1054 err = form_map_union<int_type>(&(M.RowMap()), &recvd_rows, (const Epetra_Map*&)result_map);
1055 if (err != 0) {
1056 return(NULL);
1057 }
1058
1059 return(result_map);
1060}
1061
1063 const Epetra_Map& column_map)
1064{
1065#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1066 if(M.RowMatrixRowMap().GlobalIndicesInt() && column_map.GlobalIndicesInt()) {
1067 return Tfind_rows_containing_cols<int>(M, column_map);
1068 }
1069 else
1070#endif
1071#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1072 if(M.RowMatrixRowMap().GlobalIndicesLongLong() && column_map.GlobalIndicesLongLong()) {
1073 return Tfind_rows_containing_cols<long long>(M, column_map);
1074 }
1075 else
1076#endif
1077 throw std::runtime_error("EpetraExt::find_rows_containing_cols: GlobalIndices type unknown");
1078}
1079
1080//=========================================================================
1081template<typename int_type>
1083 bool transposeA,
1084 const Epetra_CrsMatrix& B,
1085 bool transposeB,
1087 bool call_FillComplete_on_result,
1088 bool keep_all_hard_zeros)
1089{
1090 // DEBUG
1091 // bool NewFlag=!C.IndicesAreLocal() && !C.IndicesAreGlobal();
1092#ifdef ENABLE_MMM_TIMINGS
1093 using Teuchos::TimeMonitor;
1094 Teuchos::RCP<TimeMonitor> MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: MMM All Setup")));
1095#endif
1096 // end DEBUG
1097
1098 //
1099 //This method forms the matrix-matrix product C = op(A) * op(B), where
1100 //op(A) == A if transposeA is false,
1101 //op(A) == A^T if transposeA is true,
1102 //and similarly for op(B).
1103 //
1104
1105 //A and B should already be Filled.
1106 //(Should we go ahead and call FillComplete() on them if necessary?
1107 // or error out? For now, we choose to error out.)
1108 if (!A.Filled() || !B.Filled()) {
1109 EPETRA_CHK_ERR(-1);
1110 }
1111
1112 // Is the C matrix new?
1113 bool NewFlag=!C.IndicesAreLocal() && !C.IndicesAreGlobal();
1114
1115 //We're going to refer to the different combinations of op(A) and op(B)
1116 //as scenario 1 through 4.
1117
1118 int scenario = 1;//A*B
1119 if (transposeB && !transposeA) scenario = 2;//A*B^T
1120 if (transposeA && !transposeB) scenario = 3;//A^T*B
1121 if (transposeA && transposeB) scenario = 4;//A^T*B^T
1122 if(NewFlag && call_FillComplete_on_result && transposeA && !transposeB) scenario = 5; // A^T*B, newmatrix
1123
1124 //now check size compatibility
1125 long long Aouter = transposeA ? A.NumGlobalCols64() : A.NumGlobalRows64();
1126 long long Bouter = transposeB ? B.NumGlobalRows64() : B.NumGlobalCols64();
1127 long long Ainner = transposeA ? A.NumGlobalRows64() : A.NumGlobalCols64();
1128 long long Binner = transposeB ? B.NumGlobalCols64() : B.NumGlobalRows64();
1129 if (Ainner != Binner) {
1130 std::cerr << "MatrixMatrix::Multiply: ERROR, inner dimensions of op(A) and op(B) "
1131 << "must match for matrix-matrix product. op(A) is "
1132 <<Aouter<<"x"<<Ainner << ", op(B) is "<<Binner<<"x"<<Bouter<<std::endl;
1133 return(-1);
1134 }
1135
1136 //The result matrix C must at least have a row-map that reflects the
1137 //correct row-size. Don't check the number of columns because rectangular
1138 //matrices which were constructed with only one map can still end up
1139 //having the correct capacity and dimensions when filled.
1140 if (Aouter > C.NumGlobalRows64()) {
1141 std::cerr << "MatrixMatrix::Multiply: ERROR, dimensions of result C must "
1142 << "match dimensions of op(A) * op(B). C has "<<C.NumGlobalRows64()
1143 << " rows, should have at least "<<Aouter << std::endl;
1144 return(-1);
1145 }
1146
1147 //It doesn't matter whether C is already Filled or not. If it is already
1148 //Filled, it must have space allocated for the positions that will be
1149 //referenced in forming C = op(A)*op(B). If it doesn't have enough space,
1150 //we'll error out later when trying to store result values.
1151
1152 //We're going to need to import remotely-owned sections of A and/or B
1153 //if more than 1 processor is performing this run, depending on the scenario.
1154 int numProcs = A.Comm().NumProc();
1155
1156 //If we are to use the transpose of A and/or B, we'll need to be able to
1157 //access, on the local processor, all rows that contain column-indices in
1158 //the domain-map.
1159 const Epetra_Map* domainMap_A = &(A.DomainMap());
1160 const Epetra_Map* domainMap_B = &(B.DomainMap());
1161
1162 const Epetra_Map* rowmap_A = &(A.RowMap());
1163 const Epetra_Map* rowmap_B = &(B.RowMap());
1164
1165 //Declare some 'work-space' maps which may be created depending on
1166 //the scenario, and which will be deleted before exiting this function.
1167 const Epetra_Map* workmap1 = NULL;
1168 const Epetra_Map* workmap2 = NULL;
1169 const Epetra_Map* mapunion1 = NULL;
1170
1171 //Declare a couple of structs that will be used to hold views of the data
1172 //of A and B, to be used for fast access during the matrix-multiplication.
1173 CrsMatrixStruct Aview;
1174 CrsMatrixStruct Atransview;
1175 CrsMatrixStruct Bview;
1176 Teuchos::RCP<Epetra_CrsMatrix> Atrans;
1177
1178 const Epetra_Map* targetMap_A = rowmap_A;
1179 const Epetra_Map* targetMap_B = rowmap_B;
1180
1181#ifdef ENABLE_MMM_TIMINGS
1182 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: MMM All I&X")));
1183#endif
1184 if (numProcs > 1) {
1185 //If op(A) = A^T, find all rows of A that contain column-indices in the
1186 //local portion of the domain-map. (We'll import any remote rows
1187 //that fit this criteria onto the local processor.)
1188 if (scenario == 3 || scenario == 4) {
1189 workmap1 = Tfind_rows_containing_cols<int_type>(A, *domainMap_A);
1190 targetMap_A = workmap1;
1191 }
1192 }
1193 if (scenario == 5) {
1194 targetMap_A = &(A.ColMap());
1195 }
1196
1197 //Now import any needed remote rows and populate the Aview struct.
1198 if(scenario==1 && call_FillComplete_on_result) {
1199 EPETRA_CHK_ERR(import_only<int_type>(A,*targetMap_A,Aview));
1200 }
1201 else if (scenario == 5){
1202 // Perform a local transpose of A and store that in Atransview
1203 EpetraExt::RowMatrix_Transpose at(const_cast<Epetra_Map *>(targetMap_A),false);
1204 Epetra_CrsMatrix * Anonconst = const_cast<Epetra_CrsMatrix *>(&A);
1205 Atrans = Teuchos::rcp(at.CreateTransposeLocal(*Anonconst));
1206 import_only<int_type>(*Atrans,*targetMap_A,Atransview);
1207 }
1208 else {
1209 EPETRA_CHK_ERR( import_and_extract_views<int_type>(A, *targetMap_A, Aview));
1210 }
1211
1212
1213 // NOTE: Next up is to switch to import_only for B as well, and then modify the THREE SerialCores
1214 // to add a Acol2Brow and Acol2Bimportrow array for in-algorithm lookups.
1215
1216
1217 // Make sure B's views are consistent with A even in serial.
1218 const Epetra_Map* colmap_op_A = NULL;
1219 if(scenario==1 || numProcs > 1){
1220 if (transposeA && scenario == 3) {
1221 colmap_op_A = targetMap_A;
1222 }
1223 else {
1224 colmap_op_A = &(A.ColMap());
1225 }
1226 targetMap_B = colmap_op_A;
1227 }
1228 if(scenario==5) targetMap_B = &(B.RowMap());
1229
1230
1231 if (numProcs > 1) {
1232 //If op(B) = B^T, find all rows of B that contain column-indices in the
1233 //local-portion of the domain-map, or in the column-map of op(A).
1234 //We'll import any remote rows that fit this criteria onto the
1235 //local processor.
1236 if (transposeB) {
1237 EPETRA_CHK_ERR( form_map_union<int_type>(colmap_op_A, domainMap_B, mapunion1) );
1238 workmap2 = Tfind_rows_containing_cols<int_type>(B, *mapunion1);
1239 targetMap_B = workmap2;
1240 }
1241 }
1242
1243 //Now import any needed remote rows and populate the Bview struct.
1244 if((scenario==1 && call_FillComplete_on_result) || scenario==5) {
1245 EPETRA_CHK_ERR(import_only<int_type>(B,*targetMap_B,Bview,A.Importer()));
1246 }
1247 else {
1248 EPETRA_CHK_ERR( import_and_extract_views<int_type>(B, *targetMap_B, Bview) );
1249 }
1250
1251#ifdef ENABLE_MMM_TIMINGS
1252 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: MMM All Multiply")));
1253#endif
1254
1255 // Zero if filled
1256 if(C.Filled()) C.PutScalar(0.0);
1257
1258 //Now call the appropriate method to perform the actual multiplication.
1259 CrsWrapper_Epetra_CrsMatrix ecrsmat(C);
1260
1261 switch(scenario) {
1262 case 1: EPETRA_CHK_ERR( mult_A_B(A,Aview,B,Bview,C,call_FillComplete_on_result, keep_all_hard_zeros) );
1263 break;
1264 case 2: EPETRA_CHK_ERR( mult_A_Btrans(Aview, Bview, ecrsmat, keep_all_hard_zeros) );
1265 break;
1266 case 3: EPETRA_CHK_ERR( mult_Atrans_B(Aview, Bview, ecrsmat) );
1267 break;
1268 case 4: EPETRA_CHK_ERR( mult_Atrans_Btrans(Aview, Bview, ecrsmat, keep_all_hard_zeros) );
1269 break;
1270 case 5: EPETRA_CHK_ERR( mult_AT_B_newmatrix(Atransview, Bview, C, keep_all_hard_zeros) );
1271 break;
1272 }
1273
1274
1275 if (scenario != 1 && call_FillComplete_on_result && scenario != 5) {
1276 //We'll call FillComplete on the C matrix before we exit, and give
1277 //it a domain-map and a range-map.
1278 //The domain-map will be the domain-map of B, unless
1279 //op(B)==transpose(B), in which case the range-map of B will be used.
1280 //The range-map will be the range-map of A, unless
1281 //op(A)==transpose(A), in which case the domain-map of A will be used.
1282
1283 const Epetra_Map* domainmap =
1284 transposeB ? &(B.RangeMap()) : &(B.DomainMap());
1285
1286 const Epetra_Map* rangemap =
1287 transposeA ? &(A.DomainMap()) : &(A.RangeMap());
1288
1289 if (!C.Filled()) {
1290 EPETRA_CHK_ERR( C.FillComplete(*domainmap, *rangemap) );
1291 }
1292 }
1293
1294 //Finally, delete the objects that were potentially created
1295 //during the course of importing remote sections of A and B.
1296
1297 delete mapunion1; mapunion1 = NULL;
1298 delete workmap1; workmap1 = NULL;
1299 delete workmap2; workmap2 = NULL;
1300
1301 return(0);
1302}
1303
1305 bool transposeA,
1306 const Epetra_CrsMatrix& B,
1307 bool transposeB,
1309 bool call_FillComplete_on_result,
1310 bool keep_all_hard_zeros)
1311{
1312#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1313 if(A.RowMap().GlobalIndicesInt() && B.RowMap().GlobalIndicesInt()) {
1314 return TMultiply<int>(A, transposeA, B, transposeB, C, call_FillComplete_on_result, keep_all_hard_zeros);
1315 }
1316 else
1317#endif
1318#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1319 if(A.RowMap().GlobalIndicesLongLong() && B.RowMap().GlobalIndicesLongLong()) {
1320 return TMultiply<long long>(A, transposeA, B, transposeB, C, call_FillComplete_on_result, keep_all_hard_zeros);
1321 }
1322 else
1323#endif
1324 throw std::runtime_error("EpetraExt::MatrixMatrix::Add: GlobalIndices type unknown");
1325}
1326
1327//=========================================================================
1328template<typename int_type>
1330 bool transposeA,
1331 double scalarA,
1333 double scalarB )
1334{
1335 //
1336 //This method forms the matrix-matrix sum B = scalarA * op(A) + scalarB * B, where
1337
1338 //A should already be Filled. It doesn't matter whether B is
1339 //already Filled, but if it is, then its graph must already contain
1340 //all nonzero locations that will be referenced in forming the
1341 //sum.
1342
1343 if (!A.Filled() ) {
1344 std::cerr << "EpetraExt::MatrixMatrix::Add ERROR, input matrix A.Filled() is false, it is required to be true. (Result matrix B is not required to be Filled())."<<std::endl;
1345 EPETRA_CHK_ERR(-1);
1346 }
1347
1348 //explicit tranpose A formed as necessary
1349 Epetra_CrsMatrix * Aprime = 0;
1350 EpetraExt::RowMatrix_Transpose * Atrans = 0;
1351 if( transposeA )
1352 {
1353 Atrans = new EpetraExt::RowMatrix_Transpose();
1354 Aprime = &(dynamic_cast<Epetra_CrsMatrix&>(((*Atrans)(const_cast<Epetra_CrsMatrix&>(A)))));
1355 }
1356 else
1357 Aprime = const_cast<Epetra_CrsMatrix*>(&A);
1358
1359 int MaxNumEntries = EPETRA_MAX( A.MaxNumEntries(), B.MaxNumEntries() );
1360 int A_NumEntries, B_NumEntries;
1361 int_type * A_Indices = new int_type[MaxNumEntries];
1362 double * A_Values = new double[MaxNumEntries];
1363 int* B_Indices_local;
1364 int_type* B_Indices_global;
1365 double* B_Values;
1366
1367 int NumMyRows = B.NumMyRows();
1368 int_type Row;
1369 int err;
1370 int ierr = 0;
1371
1372 if( scalarA )
1373 {
1374 //Loop over B's rows and sum into
1375 for( int i = 0; i < NumMyRows; ++i )
1376 {
1377 Row = (int_type) B.GRID64(i);
1378 EPETRA_CHK_ERR( Aprime->ExtractGlobalRowCopy( Row, MaxNumEntries, A_NumEntries, A_Values, A_Indices ) );
1379
1380 if (scalarB != 1.0) {
1381 if (!B.Filled()) {
1382 EPETRA_CHK_ERR( B.ExtractGlobalRowView( Row, B_NumEntries,
1383 B_Values, B_Indices_global));
1384 }
1385 else {
1386 EPETRA_CHK_ERR( B.ExtractMyRowView( i, B_NumEntries,
1387 B_Values, B_Indices_local));
1388 }
1389
1390 for(int jj=0; jj<B_NumEntries; ++jj) {
1391 B_Values[jj] = scalarB*B_Values[jj];
1392 }
1393 }
1394
1395 if( scalarA != 1.0 ) {
1396 for( int j = 0; j < A_NumEntries; ++j ) A_Values[j] *= scalarA;
1397 }
1398
1399 if( B.Filled() ) {//Sum In Values
1400 err = B.SumIntoGlobalValues( Row, A_NumEntries, A_Values, A_Indices );
1401 assert( err >= 0 );
1402 if (err < 0) ierr = err;
1403 }
1404 else {
1405 err = B.InsertGlobalValues( Row, A_NumEntries, A_Values, A_Indices );
1406 assert( err == 0 || err == 1 || err == 3 );
1407 if (err < 0) ierr = err;
1408 }
1409 }
1410 }
1411 else {
1412 EPETRA_CHK_ERR( B.Scale(scalarB) );
1413 }
1414
1415 delete [] A_Indices;
1416 delete [] A_Values;
1417
1418 if( Atrans ) delete Atrans;
1419
1420 return(ierr);
1421}
1422
1424 bool transposeA,
1425 double scalarA,
1427 double scalarB )
1428{
1429#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1430 if(A.RowMap().GlobalIndicesInt() && B.RowMap().GlobalIndicesInt()) {
1431 return TAdd<int>(A, transposeA, scalarA, B, scalarB);
1432 }
1433 else
1434#endif
1435#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1436 if(A.RowMap().GlobalIndicesLongLong() && B.RowMap().GlobalIndicesLongLong()) {
1437 return TAdd<long long>(A, transposeA, scalarA, B, scalarB);
1438 }
1439 else
1440#endif
1441 throw std::runtime_error("EpetraExt::MatrixMatrix::Add: GlobalIndices type unknown");
1442}
1443
1444template<typename int_type>
1446 bool transposeA,
1447 double scalarA,
1448 const Epetra_CrsMatrix & B,
1449 bool transposeB,
1450 double scalarB,
1451 Epetra_CrsMatrix * & C)
1452{
1453 //
1454 //This method forms the matrix-matrix sum C = scalarA * op(A) + scalarB * op(B), where
1455
1456 //A and B should already be Filled. C should be an empty pointer.
1457
1458 if (!A.Filled() || !B.Filled() ) {
1459 std::cerr << "EpetraExt::MatrixMatrix::Add ERROR, input matrix A.Filled() or B.Filled() is false,"
1460 << "they are required to be true. (Result matrix C should be an empty pointer)" << std::endl;
1461 EPETRA_CHK_ERR(-1);
1462 }
1463
1464 Epetra_CrsMatrix * Aprime = 0, * Bprime=0;
1465 EpetraExt::RowMatrix_Transpose * Atrans = 0,* Btrans = 0;
1466
1467 //explicit tranpose A formed as necessary
1468 if( transposeA ) {
1469 Atrans = new EpetraExt::RowMatrix_Transpose();
1470 Aprime = &(dynamic_cast<Epetra_CrsMatrix&>(((*Atrans)(const_cast<Epetra_CrsMatrix&>(A)))));
1471 }
1472 else
1473 Aprime = const_cast<Epetra_CrsMatrix*>(&A);
1474
1475 //explicit tranpose B formed as necessary
1476 if( transposeB ) {
1477 Btrans = new EpetraExt::RowMatrix_Transpose();
1478 Bprime = &(dynamic_cast<Epetra_CrsMatrix&>(((*Btrans)(const_cast<Epetra_CrsMatrix&>(B)))));
1479 }
1480 else
1481 Bprime = const_cast<Epetra_CrsMatrix*>(&B);
1482
1483 // allocate or zero the new matrix
1484 if(C!=0)
1485 C->PutScalar(0.0);
1486 else
1487 C = new Epetra_CrsMatrix(Copy,Aprime->RowMap(),0);
1488
1489 // build arrays for easy resuse
1490 int ierr = 0;
1491 Epetra_CrsMatrix * Mat[] = { Aprime,Bprime};
1492 double scalar[] = { scalarA, scalarB};
1493
1494 // do a loop over each matrix to add: A reordering might be more efficient
1495 for(int k=0;k<2;k++) {
1496 int MaxNumEntries = Mat[k]->MaxNumEntries();
1497 int NumEntries;
1498 int_type * Indices = new int_type[MaxNumEntries];
1499 double * Values = new double[MaxNumEntries];
1500
1501 int NumMyRows = Mat[k]->NumMyRows();
1502 int err;
1503 int_type Row;
1504
1505 //Loop over rows and sum into C
1506 for( int i = 0; i < NumMyRows; ++i ) {
1507 Row = (int_type) Mat[k]->GRID64(i);
1508 EPETRA_CHK_ERR( Mat[k]->ExtractGlobalRowCopy( Row, MaxNumEntries, NumEntries, Values, Indices));
1509
1510 if( scalar[k] != 1.0 )
1511 for( int j = 0; j < NumEntries; ++j ) Values[j] *= scalar[k];
1512
1513 if(C->Filled()) { // Sum in values
1514 err = C->SumIntoGlobalValues( Row, NumEntries, Values, Indices );
1515 if (err < 0) ierr = err;
1516 } else { // just add it to the unfilled CRS Matrix
1517 err = C->InsertGlobalValues( Row, NumEntries, Values, Indices );
1518 if (err < 0) ierr = err;
1519 }
1520 }
1521
1522 delete [] Indices;
1523 delete [] Values;
1524 }
1525
1526 if( Atrans ) delete Atrans;
1527 if( Btrans ) delete Btrans;
1528
1529 return(ierr);
1530}
1531
1533 bool transposeA,
1534 double scalarA,
1535 const Epetra_CrsMatrix & B,
1536 bool transposeB,
1537 double scalarB,
1538 Epetra_CrsMatrix * & C)
1539{
1540#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1541 if(A.RowMap().GlobalIndicesInt() && B.RowMap().GlobalIndicesInt()) {
1542 return TAdd<int>(A, transposeA, scalarA, B, transposeB, scalarB, C);
1543 }
1544 else
1545#endif
1546#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1547 if(A.RowMap().GlobalIndicesLongLong() && B.RowMap().GlobalIndicesLongLong()) {
1548 return TAdd<long long>(A, transposeA, scalarA, B, transposeB, scalarB, C);
1549 }
1550 else
1551#endif
1552 throw std::runtime_error("EpetraExt::MatrixMatrix::Add: GlobalIndices type unknown");
1553}
1554
1555
1556
1557//=========================================================================
1558template<typename int_type>
1559int MatrixMatrix::TJacobi(double omega,
1560 const Epetra_Vector & Dinv,
1561 const Epetra_CrsMatrix& A,
1562 const Epetra_CrsMatrix& B,
1564 bool call_FillComplete_on_result)
1565{
1566#ifdef ENABLE_MMM_TIMINGS
1567 using Teuchos::TimeMonitor;
1568 Teuchos::RCP<TimeMonitor> MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: Jacobi All Setup")));
1569#endif
1570
1571 //A and B should already be Filled.
1572 if (!A.Filled() || !B.Filled()) {
1573 EPETRA_CHK_ERR(-1);
1574 }
1575
1576 //now check size compatibility
1577 long long Aouter = A.NumGlobalRows64();
1578 long long Bouter = B.NumGlobalCols64();
1579 long long Ainner = A.NumGlobalCols64();
1580 long long Binner = B.NumGlobalRows64();
1581 long long Dlen = Dinv.GlobalLength64();
1582 if (Ainner != Binner) {
1583 std::cerr << "MatrixMatrix::Jacobi: ERROR, inner dimensions of A and B "
1584 << "must match for matrix-matrix product. A is "
1585 <<Aouter<<"x"<<Ainner << ", B is "<<Binner<<"x"<<Bouter<<std::endl;
1586 return(-1);
1587 }
1588
1589 //The result matrix C must at least have a row-map that reflects the
1590 //correct row-size. Don't check the number of columns because rectangular
1591 //matrices which were constructed with only one map can still end up
1592 //having the correct capacity and dimensions when filled.
1593 if (Aouter > C.NumGlobalRows64()) {
1594 std::cerr << "MatrixMatrix::Jacobi: ERROR, dimensions of result C must "
1595 << "match dimensions of A * B. C has "<<C.NumGlobalRows64()
1596 << " rows, should have at least "<<Aouter << std::endl;
1597 return(-1);
1598 }
1599
1600 // Check against the D matrix
1601 if(Dlen != Aouter) {
1602 std::cerr << "MatrixMatrix::Jacboi: ERROR, dimensions of result D must "
1603 << "match dimensions of A's rows. D has "<< Dlen
1604 << " rows, should have " << Aouter << std::endl;
1605 return(-1);
1606 }
1607
1608 if(!A.RowMap().SameAs(B.RowMap()) || !A.RowMap().SameAs(Dinv.Map())) {
1609 std::cerr << "MatrixMatrix::Jacboi: ERROR, RowMap of A must match RowMap of B "
1610 << "and Map of D."<<std::endl;
1611 return(-1);
1612 }
1613
1614 //It doesn't matter whether C is already Filled or not. If it is already
1615 //Filled, it must have space allocated for the positions that will be
1616 //referenced in forming C. If it doesn't have enough space,
1617 //we'll error out later when trying to store result values.
1618
1619 //We're going to need to import remotely-owned sections of A and/or B
1620 //if more than 1 processor is performing this run, depending on the scenario.
1621 int numProcs = A.Comm().NumProc();
1622
1623 // Maps
1624 const Epetra_Map* rowmap_A = &(A.RowMap());
1625 const Epetra_Map* rowmap_B = &(B.RowMap());
1626
1627
1628
1629 //Declare some 'work-space' maps which may be created depending on
1630 //the scenario, and which will be deleted before exiting this function.
1631 const Epetra_Map* workmap1 = NULL;
1632 const Epetra_Map* workmap2 = NULL;
1633 const Epetra_Map* mapunion1 = NULL;
1634
1635 //Declare a couple of structs that will be used to hold views of the data
1636 //of A and B, to be used for fast access during the matrix-multiplication.
1637 CrsMatrixStruct Aview;
1638 CrsMatrixStruct Bview;
1639
1640 const Epetra_Map* targetMap_A = rowmap_A;
1641 const Epetra_Map* targetMap_B = rowmap_B;
1642
1643#ifdef ENABLE_MMM_TIMINGS
1644 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: Jacobi All I&X")));
1645#endif
1646
1647 //Now import any needed remote rows and populate the Aview struct.
1648 if(call_FillComplete_on_result) {
1649 EPETRA_CHK_ERR(import_only<int_type>(A,*targetMap_A,Aview));
1650 }
1651 else {
1652 EPETRA_CHK_ERR( import_and_extract_views<int_type>(A, *targetMap_A, Aview));
1653 }
1654
1655 // NOTE: Next up is to switch to import_only for B as well, and then modify the THREE SerialCores
1656 // to add a Acol2Brow and Acol2Bimportrow array for in-algorithm lookups.
1657
1658 // Make sure B's views are consistent with A even in serial.
1659 const Epetra_Map* colmap_op_A = NULL;
1660 if(numProcs > 1){
1661 colmap_op_A = &(A.ColMap());
1662 targetMap_B = colmap_op_A;
1663 }
1664
1665 //Now import any needed remote rows and populate the Bview struct.
1666 if(call_FillComplete_on_result) {
1667 EPETRA_CHK_ERR(import_only<int_type>(B,*targetMap_B,Bview,A.Importer()));
1668 }
1669 else {
1670 EPETRA_CHK_ERR( import_and_extract_views<int_type>(B, *targetMap_B, Bview) );
1671 }
1672
1673#ifdef ENABLE_MMM_TIMINGS
1674 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer("EpetraExt: Jacobi All Multiply")));
1675#endif
1676
1677 // Zero if filled
1678 if(C.Filled()) C.PutScalar(0.0);
1679
1680 //Now call the appropriate method to perform the actual multiplication.
1681 CrsWrapper_Epetra_CrsMatrix ecrsmat(C);
1682 EPETRA_CHK_ERR( jacobi_A_B(omega,Dinv,A,Aview,B,Bview,C,call_FillComplete_on_result) );
1683
1684 //Finally, delete the objects that were potentially created
1685 //during the course of importing remote sections of A and B.
1686 delete mapunion1; mapunion1 = NULL;
1687 delete workmap1; workmap1 = NULL;
1688 delete workmap2; workmap2 = NULL;
1689
1690 return(0);
1691}
1692
1693
1694
1695int MatrixMatrix::Jacobi(double omega,
1696 const Epetra_Vector & Dinv,
1697 const Epetra_CrsMatrix& A,
1698 const Epetra_CrsMatrix& B,
1700 bool call_FillComplete_on_result)
1701{
1702#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
1703 if(A.RowMap().GlobalIndicesInt() && B.RowMap().GlobalIndicesInt()) {
1704 return TJacobi<int>(omega, Dinv, A, B, C, call_FillComplete_on_result);
1705 }
1706 else
1707#endif
1708#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
1709 if(A.RowMap().GlobalIndicesLongLong() && B.RowMap().GlobalIndicesLongLong()) {
1710 return TJacobi<long long>(omega, Dinv, A, B, C, call_FillComplete_on_result);
1711 }
1712 else
1713#endif
1714 throw std::runtime_error("EpetraExt::MatrixMatrix::Jacobi: GlobalIndices type unknown");
1715}
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726} // namespace EpetraExt
1727
#define EPETRA_MAX(x, y)
#define EPETRA_CHK_ERR(a)
Copy
virtual bool Filled()=0
virtual const Epetra_Map & RowMap() const =0
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)=0
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)=0
static int Jacobi(double omega, const Epetra_Vector &Dinv, const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B, Epetra_CrsMatrix &C, bool call_FillComplete_on_result=true)
Given Epetra_CrsMatrix objects A, B and C, and Epetra_Vector Dinv, form the product C = (I-omega * Di...
static int TAdd(const Epetra_CrsMatrix &A, bool transposeA, double scalarA, Epetra_CrsMatrix &B, double scalarB)
static int jacobi_A_B(double omega, const Epetra_Vector &Dinv, const Epetra_CrsMatrix &A, CrsMatrixStruct &Aview, const Epetra_CrsMatrix &B, CrsMatrixStruct &Bview, Epetra_CrsMatrix &C, bool call_FillComplete_on_result)
static int mult_AT_B_newmatrix(const CrsMatrixStruct &Atransview, const CrsMatrixStruct &Bview, Epetra_CrsMatrix &C, bool keep_all_hard_zeros)
static int Add(const Epetra_CrsMatrix &A, bool transposeA, double scalarA, Epetra_CrsMatrix &B, double scalarB)
Given Epetra_CrsMatrix objects A and B, form the sum B = a*A + b*B.
static int Multiply(const Epetra_CrsMatrix &A, bool transposeA, const Epetra_CrsMatrix &B, bool transposeB, Epetra_CrsMatrix &C, bool call_FillComplete_on_result=true, bool keep_all_hard_zeros=false)
Given Epetra_CrsMatrix objects A, B and C, form the product C = A*B.
static int TMultiply(const Epetra_CrsMatrix &A, bool transposeA, const Epetra_CrsMatrix &B, bool transposeB, Epetra_CrsMatrix &C, bool call_FillComplete_on_result, bool keep_all_hard_zeros)
static int mult_A_B(const Epetra_CrsMatrix &A, CrsMatrixStruct &Aview, const Epetra_CrsMatrix &B, CrsMatrixStruct &Bview, Epetra_CrsMatrix &C, bool call_FillComplete_on_result, bool keep_all_hard_zeros)
static int TJacobi(double omega, const Epetra_Vector &Dinv, const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B, Epetra_CrsMatrix &C, bool call_FillComplete_on_result)
Transform to form the explicit transpose of a Epetra_RowMatrix.
int MyGlobalElementsPtr(int *&MyGlobalElementList) const
virtual int NumProc() const=0
virtual int MyPID() const=0
virtual Epetra_Distributor * CreateDistributor() const=0
bool Filled() const
const Epetra_Map & RowMap() const
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
int FillComplete(bool OptimizeDataStorage=true)
int PutScalar(double ScalarConstant)
long long NumGlobalRows64() const
int Scale(double ScalarConstant)
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
const Epetra_Comm & Comm() const
const Epetra_Map & RangeMap() const
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
const Epetra_Map & ColMap() const
bool IndicesAreLocal() const
int MaxNumEntries() const
const Epetra_Map & DomainMap() const
int NumMyRows() const
const Epetra_Import * Importer() const
int ExtractGlobalRowView(int_type Row, int &NumEntries, double *&values, int_type *&Indices) const
long long GRID64(int LRID_in) const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
bool IndicesAreGlobal() const
const Epetra_Map & RowMatrixRowMap() const
long long NumGlobalCols64() const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
int mult_A_Btrans(CrsMatrixStruct &Aview, CrsMatrixStruct &Bview, CrsWrapper &C, bool keep_all_hard_zeros)
int import_and_extract_views(const Epetra_CrsMatrix &M, const Epetra_Map &targetMap, CrsMatrixStruct &Mview, const Epetra_Import *prototypeImporter=0)
double sparsedot(double *u, int_type *u_ind, int u_len, double *v, int_type *v_ind, int v_len)
Method for internal use... sparsedot forms a dot-product between two sparsely-populated 'vectors'.
int mult_Atrans_Btrans(CrsMatrixStruct &Aview, CrsMatrixStruct &Bview, CrsWrapper &C, bool keep_all_hard_zeros)
void pack_outgoing_rows(const Epetra_CrsMatrix &mtx, const std::vector< int > &proc_col_ranges, std::vector< int > &send_rows, std::vector< int > &rows_per_send_proc)
int import_only(const Epetra_CrsMatrix &M, const Epetra_Map &targetMap, CrsMatrixStruct &Mview, const Epetra_Import *prototypeImporter)
Epetra_Map * find_rows_containing_cols(const Epetra_CrsMatrix &M, const Epetra_Map &column_map)
const int N
Epetra_Map * Tfind_rows_containing_cols(const Epetra_CrsMatrix &M, const Epetra_Map &column_map)
int mult_Atrans_B(CrsMatrixStruct &Aview, CrsMatrixStruct &Bview, CrsWrapper &C)
int form_map_union(const Epetra_Map *map1, const Epetra_Map *map2, const Epetra_Map *&mapunion)