EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
42 #include <EpetraExt_ConfigDefs.h>
43 #include <EpetraExt_MatrixMatrix.h>
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 
68 namespace EpetraExt {
69 
70 //=========================================================================
71 // ETI
72 template<int> int import_only(const Epetra_CrsMatrix& M,
73  const Epetra_Map& targetMap,
74  CrsMatrixStruct& Mview,
75  const Epetra_Import * prototypeImporter);
76 
77 template<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 //
89 template<typename int_type>
90 double 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
118 template<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 
321 int mult_A_Btrans(CrsMatrixStruct& Aview,
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
349 template<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 
486 int mult_Atrans_B(CrsMatrixStruct& Aview,
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
512 template<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 
650 int mult_Atrans_Btrans(CrsMatrixStruct& Aview,
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 // ==============================================================
677 template<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),
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 //=========================================================================
904 template<typename int_type>
905 int form_map_union(const Epetra_Map* map1,
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 //=========================================================================
969 template<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 //=========================================================================
1081 template<typename int_type>
1082 int MatrixMatrix::TMultiply(const Epetra_CrsMatrix& A,
1083  bool transposeA,
1084  const Epetra_CrsMatrix& B,
1085  bool transposeB,
1086  Epetra_CrsMatrix& C,
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,
1308  Epetra_CrsMatrix& C,
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 //=========================================================================
1328 template<typename int_type>
1329 int MatrixMatrix::TAdd(const Epetra_CrsMatrix& A,
1330  bool transposeA,
1331  double scalarA,
1332  Epetra_CrsMatrix& B,
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,
1426  Epetra_CrsMatrix& B,
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 
1444 template<typename int_type>
1445 int MatrixMatrix::TAdd(const Epetra_CrsMatrix& A,
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 //=========================================================================
1558 template<typename int_type>
1559 int MatrixMatrix::TJacobi(double omega,
1560  const Epetra_Vector & Dinv,
1561  const Epetra_CrsMatrix& A,
1562  const Epetra_CrsMatrix& B,
1563  Epetra_CrsMatrix& C,
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 
1695 int MatrixMatrix::Jacobi(double omega,
1696  const Epetra_Vector & Dinv,
1697  const Epetra_CrsMatrix& A,
1698  const Epetra_CrsMatrix& B,
1699  Epetra_CrsMatrix& C,
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 
LightweightCrsMatrix * importMatrix
const Epetra_Map & RangeMap() const
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)=0
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int MaxLID() const
bool GlobalIndicesLongLong() const
void MyGlobalElementsPtr(int *&MyGlobalElementList) const
bool SameAs(const Epetra_BlockMap &Map) const
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)
virtual int CreateFromSends(const int &NumExportIDs, const int *ExportPIDs, bool Deterministic, int &NumRemoteIDs)=0
const Epetra_Map & RowMatrixRowMap() const
virtual bool Filled()=0
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.
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)=0
bool GlobalIndicesInt() const
const Epetra_Map & ColMap() const
bool IndicesAreLocal() const
const Epetra_CrsMatrix * origMatrix
int MinLID() const
Transform to form the explicit transpose of a Epetra_RowMatrix.
virtual int MyPID() const =0
int FillComplete(bool OptimizeDataStorage=true)
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
int PutScalar(double ScalarConstant)
Epetra_Map * find_rows_containing_cols(const Epetra_CrsMatrix &M, const Epetra_Map &column_map)
int mult_Atrans_Btrans(CrsMatrixStruct &Aview, CrsMatrixStruct &Bview, CrsWrapper &C, bool keep_all_hard_zeros)
const Epetra_Map & RowMap() const
int NumMyElements() const
int form_map_union(const Epetra_Map *map1, const Epetra_Map *map2, const Epetra_Map *&mapunion)
int import_and_extract_views(const Epetra_CrsMatrix &M, const Epetra_Map &targetMap, CrsMatrixStruct &Mview, const Epetra_Import *prototypeImporter=0)
virtual const Epetra_Map & RowMap() const =0
const Epetra_Import * Importer() const
int ExtractGlobalRowView(int GlobalRow, int &NumEntries, double *&Values, int *&Indices) const
int MaxNumEntries() const
int ExtractGlobalRowCopy(int GlobalRow, int Length, int &NumEntries, double *Values, int *Indices) const
int mult_Atrans_B(CrsMatrixStruct &Aview, CrsMatrixStruct &Bview, CrsWrapper &C)
int NumMyRows() const
int LID(int GID) const
bool MyGID(int GID_in) const
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.
const Epetra_Comm & Comm() const
virtual int Do(char *export_objs, int obj_size, int &len_import_objs, char *&import_objs)=0
bool IndicesAreGlobal() const
long long IndexBase64() const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
double sparsedot(double *u, int_type *u_ind, int u_len, double *v, int_type *v_ind, int v_len)
Method for internal use...
virtual int NumProc() const =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...
bool Filled() const
const Epetra_Map & DomainMap() const
static void EPETRA_LIB_DLL_EXPORT Sort(bool SortAscending, int NumKeys, T *Keys, int NumDoubleCompanions, double **DoubleCompanions, int NumIntCompanions, int **IntCompanions, int NumLongLongCompanions, long long **LongLongCompanions)
Epetra_Map * Tfind_rows_containing_cols(const Epetra_CrsMatrix &M, const Epetra_Map &column_map)
int Scale(double ScalarConstant)
int mult_A_Btrans(CrsMatrixStruct &Aview, CrsMatrixStruct &Bview, CrsWrapper &C, bool keep_all_hard_zeros)
const Epetra_BlockMap * importColMap
const Epetra_Comm & Comm() const
int import_only(const Epetra_CrsMatrix &M, const Epetra_Map &targetMap, CrsMatrixStruct &Mview, const Epetra_Import *prototypeImporter)
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const