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) {
855 } else {
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),
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
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
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
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
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
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
const Epetra_CrsMatrix * origMatrix
LightweightCrsMatrix * importMatrix
const Epetra_BlockMap * importColMap
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
void MyGlobalElementsPtr(int *&MyGlobalElementList) const
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.
Epetra_CrsMatrix * CreateTransposeLocal(OriginalTypeRef orig)
Local-only transpose operator. Don't use this unless you're sure you know what you're doing.
long long IndexBase64() const
long long MaxAllGID64() const
int LID(int GID) const
long long GID64(int LID) const
int MinLID() const
bool SameAs(const Epetra_BlockMap &Map) const
bool GlobalIndicesInt() const
const Epetra_Comm & Comm() const
long long MinAllGID64() const
int NumMyElements() const
int MaxLID() const
bool MyGID(int GID_in) const
bool GlobalIndicesLongLong() const
int MyGlobalElementsPtr(int *&MyGlobalElementList) const
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const=0
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)
int ExtractGlobalRowCopy(int_type Row, int Length, int &NumEntries, double *values, int_type *Indices) const
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
const Epetra_BlockMap & Map() const
virtual int Do(char *export_objs, int obj_size, int &len_import_objs, char *&import_objs)=0
virtual int CreateFromSends(const int &NumExportIDs, const int *ExportPIDs, bool Deterministic, int &NumRemoteIDs)=0
long long GlobalLength64() const
void Sort(bool SortAscending, int NumKeys, T *Keys, int NumDoubleCompanions, double **DoubleCompanions, int NumIntCompanions, int **IntCompanions, int NumLongLongCompanions, long long **LongLongCompanions)
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)