Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Epetra_CrsSingletonFilter.cpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ************************************************************************
4 //
5 // Epetra: Linear Algebra Services Package
6 // Copyright 2011 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 */
43 
44 #include "Epetra_ConfigDefs.h"
45 #include "Epetra_Map.h"
46 #include "Epetra_Util.h"
47 #include "Epetra_Export.h"
48 #include "Epetra_Import.h"
49 #include "Epetra_MultiVector.h"
50 #include "Epetra_Vector.h"
51 #include "Epetra_IntVector.h"
52 
53 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
54 #include "Epetra_LongLongVector.h"
55 #endif
56 
57 #include "Epetra_Comm.h"
58 #include "Epetra_LinearProblem.h"
59 #include "Epetra_MapColoring.h"
61 
62 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES // FIXME
63 // FIXME long long : whole file
64 
65 //==============================================================================
68 }
69 //==============================================================================
71 
72 
73  if (ReducedProblem_!=0) delete ReducedProblem_;
74  if (ReducedMatrix_!=0) delete ReducedMatrix_;
75  if (ReducedLHS_!=0) delete ReducedLHS_;
76  if (ReducedRHS_!=0) delete ReducedRHS_;
86  if (RowMapColors_ != 0) delete RowMapColors_;
87  if (ColMapColors_ != 0) delete ColMapColors_;
88 
89  if (ColSingletonRowLIDs_ != 0) delete [] ColSingletonRowLIDs_;
90  if (ColSingletonColLIDs_ != 0) delete [] ColSingletonColLIDs_;
92  if (ColSingletonPivots_ != 0) delete [] ColSingletonPivots_;
93  if (tempExportX_ != 0) delete tempExportX_;
94  if (Indices_ != 0) delete [] Indices_;
95  if (tempX_ != 0) delete tempX_;
96  if (tempB_ != 0) delete tempB_;
97 
98 }
99 //==============================================================================
101 
102 // Initialize all attributes that have trivial default values
103 
104  FullProblem_ = 0;
105  ReducedProblem_ = 0;
106  FullMatrix_ = 0;
107  ReducedMatrix_ = 0;
108  ReducedRHS_ = 0;
109  ReducedLHS_ = 0;
118 
123 
124  AbsoluteThreshold_ = 0;
125  RelativeThreshold_ = 0;
126 
127  NumMyRowSingletons_ = -1;
128  NumMyColSingletons_ = -1;
131  RatioOfDimensions_ = -1.0;
132  RatioOfNonzeros_ = -1.0;
133 
134  HaveReducedProblem_ = false;
136  AnalysisDone_ = false;
137  SymmetricElimination_ = true;
138 
139  tempExportX_ = 0;
140  tempX_ = 0;
141  tempB_ = 0;
142 
143  Indices_ = 0;
144 
145  RowMapColors_ = 0;
146  ColMapColors_ = 0;
147 
148  FullMatrixIsCrsMatrix_ = false;
149  MaxNumMyEntries_ = 0;
150  return;
151 }
152 //==============================================================================
154 
155  int i, j, jj;
156 
157  FullMatrix_ = fullMatrix;
158 
159  if (AnalysisDone_) EPETRA_CHK_ERR(-1); // Analysis already done once. Cannot do it again
160  if (fullMatrix==0) EPETRA_CHK_ERR(-2); // Input matrix pointer is zero
161  if (fullMatrix->NumGlobalRows64()==0) EPETRA_CHK_ERR(-3); // Full matrix has zero dimension.
162  if (fullMatrix->NumGlobalNonzeros64()==0) EPETRA_CHK_ERR(-4); // Full matrix has no nonzero terms.
163 
164 
165  // First check for columns with single entries and find columns with singleton rows
166  Epetra_IntVector ColProfiles(FullMatrixColMap()); ColProfiles.PutValue(0);
167  Epetra_IntVector ColHasRowWithSingleton(FullMatrixColMap()); ColHasRowWithSingleton.PutValue(0);
168 
169  // RowIDs[j] will contain the local row ID associated with the jth column,
170  // if the jth col has a single entry
171  Epetra_IntVector RowIDs(FullMatrixColMap()); RowIDs.PutValue(-1);
172 
173  // Define MapColoring objects
174  RowMapColors_ = new Epetra_MapColoring(FullMatrixRowMap()); // Initial colors are all 0
176  Epetra_MapColoring & rowMapColors = *RowMapColors_;
177  Epetra_MapColoring & colMapColors = *ColMapColors_;
178 
179 
180  int NumMyRows = fullMatrix->NumMyRows();
181  int NumMyCols = fullMatrix->NumMyCols();
182 
183  // Set up for accessing full matrix. Will do so row-by-row.
185 
186  // Scan matrix for singleton rows, build up column profiles
187  int NumIndices;
188  int * Indices;
190  for (i=0; i<NumMyRows; i++) {
191  // Get ith row
192  EPETRA_CHK_ERR(GetRow(i, NumIndices, Indices));
193  for (j=0; j<NumIndices; j++) {
194  int ColumnIndex = Indices[j];
195  ColProfiles[ColumnIndex]++; // Increment column count
196  // Record local row ID for current column
197  // will use to identify row to eliminate if column is a singleton
198  RowIDs[ColumnIndex] = i;
199  }
200  // If row has single entry, color it and associated column with color=1
201  if (NumIndices==1) {
202  j = Indices[0];
203  ColHasRowWithSingleton[j]++;
204  rowMapColors[i] = 1;
205  colMapColors[j] = 1;
207  }
208  }
209 
210  // 1) The vector ColProfiles has column nonzero counts for each processor's contribution
211  // Combine these to get total column profile information and then redistribute to processors
212  // so each can determine if it is the owner of the row associated with the singleton column
213  // 2) The vector ColHasRowWithSingleton[i] contain count of singleton rows that are associated with
214  // the ith column on this processor. Must tell other processors that they should also eliminate
215  // these columns.
216 
217  // Make a copy of ColProfiles for later use when detecting columns that disappear locally
218 
219  Epetra_IntVector NewColProfiles(ColProfiles);
220 
221  // If RowMatrixImporter is non-trivial, we need to perform a gather/scatter to accumulate results
222 
223  if (fullMatrix->RowMatrixImporter()!=0) {
224  Epetra_IntVector tmpVec(FullMatrixDomainMap()); // Use for gather/scatter of column vectors
225  EPETRA_CHK_ERR(tmpVec.Export(ColProfiles, *fullMatrix->RowMatrixImporter(), Add));
226  EPETRA_CHK_ERR(ColProfiles.Import(tmpVec, *fullMatrix->RowMatrixImporter(), Insert));
227 
228  EPETRA_CHK_ERR(tmpVec.PutValue(0));
229  EPETRA_CHK_ERR(tmpVec.Export(ColHasRowWithSingleton, *fullMatrix->RowMatrixImporter(), Add));
230  EPETRA_CHK_ERR(ColHasRowWithSingleton.Import(tmpVec, *fullMatrix->RowMatrixImporter(), Insert));
231  }
232  // ColProfiles now contains the nonzero column entry count for all columns that have
233  // an entry on this processor.
234  // ColHasRowWithSingleton now contains a count of singleton rows associated with the corresponding
235  // local column. Next we check to make sure no column is associated with more than one singleton row.
236 
237  if (ColHasRowWithSingleton.MaxValue()>1) {
238  EPETRA_CHK_ERR(-2); // At least one col is associated with two singleton rows, can't handle it.
239  }
240 
241 
242  Epetra_IntVector RowHasColWithSingleton(fullMatrix->RowMatrixRowMap()); // Use to check for errors
243  RowHasColWithSingleton.PutValue(0);
244 
246  // Count singleton columns (that were not already counted as singleton rows)
247  for (j=0; j<NumMyCols; j++) {
248  i = RowIDs[j];
249  // Check if column is a singleton
250  if (ColProfiles[j]==1) {
251  // Check to see if this column already eliminated by the row check above
252  if (rowMapColors[i]!=1) {
253  RowHasColWithSingleton[i]++; // Increment col singleton counter for ith row
254  rowMapColors[i] = 2; // Use 2 for now, to distinguish between row eliminated directly or via column singletons
255  colMapColors[j] = 1;
257  // If we delete a row, we need to keep track of associated column entries that were also deleted
258  // in case all entries in a column are eventually deleted, in which case the column should
259  // also be deleted.
260  EPETRA_CHK_ERR(GetRow(i, NumIndices, Indices));
261  for (jj=0; jj<NumIndices; jj++) NewColProfiles[Indices[jj]]--;
262 
263  }
264  }
265  // Check if some other processor eliminated this column
266  else if (ColHasRowWithSingleton[j]==1 && rowMapColors[i]!=1) {
267  colMapColors[j] = 1;
268  }
269  }
270  if (RowHasColWithSingleton.MaxValue()>1) {
271  EPETRA_CHK_ERR(-3); // At lease one row is associated with two singleton cols, can't handle it.
272  }
273 
274 
275  // Generate arrays that keep track of column singleton row, col and pivot info needed for post-solve phase
276  EPETRA_CHK_ERR(CreatePostSolveArrays(RowIDs, rowMapColors, ColProfiles, NewColProfiles,
277  ColHasRowWithSingleton));
278 
279  for (i=0; i<NumMyRows; i++) if (rowMapColors[i]==2) rowMapColors[i] = 1; // Convert all eliminated rows to same color
280 
283  AnalysisDone_ = true;
284  return(0);
285 }
286 //==============================================================================
287 
289 
290  int i, j;
291  if (HaveReducedProblem_) EPETRA_CHK_ERR(-1); // Setup already done once. Cannot do it again
292  if (Problem==0) EPETRA_CHK_ERR(-2); // Null problem pointer
293 
294  FullProblem_ = Problem;
295  FullMatrix_ = dynamic_cast<Epetra_RowMatrix *>(Problem->GetMatrix());
296  if (FullMatrix_==0) EPETRA_CHK_ERR(-3); // Need a RowMatrix
297  if (Problem->GetRHS()==0) EPETRA_CHK_ERR(-4); // Need a RHS
298  if (Problem->GetLHS()==0) EPETRA_CHK_ERR(-5); // Need a LHS
299  // Generate reduced row and column maps
300 
301  Epetra_MapColoring & rowMapColors = *RowMapColors_;
302  Epetra_MapColoring & colMapColors = *ColMapColors_;
303 
304  ReducedMatrixRowMap_ = rowMapColors.GenerateMap(0);
305  ReducedMatrixColMap_ = colMapColors.GenerateMap(0);
306 
307  // Create domain and range map colorings by exporting map coloring of column and row maps
308 
309  if (FullMatrix()->RowMatrixImporter()!=0) {
310  Epetra_MapColoring DomainMapColors(FullMatrixDomainMap());
311  EPETRA_CHK_ERR(DomainMapColors.Export(*ColMapColors_, *FullMatrix()->RowMatrixImporter(), AbsMax));
312  OrigReducedMatrixDomainMap_ = DomainMapColors.GenerateMap(0);
313  }
314  else
316 
318  if (FullCrsMatrix()->Exporter()!=0) { // Non-trivial exporter
319  Epetra_MapColoring RangeMapColors(FullMatrixRangeMap());
320  EPETRA_CHK_ERR(RangeMapColors.Export(*RowMapColors_, *FullCrsMatrix()->Exporter(),
321  AbsMax));
322  ReducedMatrixRangeMap_ = RangeMapColors.GenerateMap(0);
323  }
324  else
326  }
327  else
329 
330  // Check to see if the reduced system domain and range maps are the same.
331  // If not, we need to remap entries of the LHS multivector so that they are distributed
332  // conformally with the rows of the reduced matrix and the RHS multivector
337  else {
341  }
342 
343  // Create pointer to Full RHS, LHS
344  Epetra_MultiVector * FullRHS = FullProblem()->GetRHS();
345  Epetra_MultiVector * FullLHS = FullProblem()->GetLHS();
346  int NumVectors = FullLHS->NumVectors();
347 
348  // Create importers
351 
352  // Construct Reduced Matrix
354 
355  // Create storage for temporary X values due to explicit elimination of rows
356  tempExportX_ = new Epetra_MultiVector(FullMatrixColMap(), NumVectors);
357 
358  int NumEntries;
359  int * Indices;
360  double * Values;
361  int NumMyRows = FullMatrix()->NumMyRows();
362  int ColSingletonCounter = 0;
363  for (i=0; i<NumMyRows; i++) {
364  int curGRID = FullMatrixRowMap().GID64(i); // CJ FIXME TODO long long
365  if (ReducedMatrixRowMap()->MyGID(curGRID)) { // Check if this row should go into reduced matrix
366 
367  EPETRA_CHK_ERR(GetRowGCIDs(i, NumEntries, Values, Indices)); // Get current row (Indices are global)
368 
369  int ierr = ReducedMatrix()->InsertGlobalValues(curGRID, NumEntries,
370  Values, Indices); // Insert into reduce matrix
371  // Positive errors will occur because we are submitting col entries that are not part of
372  // reduced system. However, because we specified a column map to the ReducedMatrix constructor
373  // these extra column entries will be ignored and we will be politely reminded by a positive
374  // error code
375  if (ierr<0) EPETRA_CHK_ERR(ierr);
376  }
377  else {
378  EPETRA_CHK_ERR(GetRow(i, NumEntries, Values, Indices)); // Get current row
379  if (NumEntries==1) {
380  double pivot = Values[0];
381  if (pivot==0.0) EPETRA_CHK_ERR(-1); // Encountered zero row, unable to continue
382  int indX = Indices[0];
383  for (j=0; j<NumVectors; j++)
384  (*tempExportX_)[j][indX] = (*FullRHS)[j][i]/pivot;
385  }
386  // Otherwise, this is a singleton column and we will scan for the pivot element needed
387  // for post-solve equations
388  else {
389  int targetCol = ColSingletonColLIDs_[ColSingletonCounter];
390  for (j=0; j<NumEntries; j++) {
391  if (Indices[j]==targetCol) {
392  double pivot = Values[j];
393  if (pivot==0.0) EPETRA_CHK_ERR(-2); // Encountered zero column, unable to continue
394  ColSingletonPivotLIDs_[ColSingletonCounter] = j; // Save for later use
395  ColSingletonPivots_[ColSingletonCounter] = pivot;
396  ColSingletonCounter++;
397  break;
398  }
399  }
400  }
401  }
402  }
403 
404  // Now convert to local indexing. We have constructed things so that the domain and range of the
405  // matrix will have the same map. If the reduced matrix domain and range maps were not the same, the
406  // differences were addressed in the ConstructRedistributeExporter() method
408 
409  // Construct Reduced LHS (Puts any initial guess values into reduced system)
410 
413  FullLHS->PutScalar(0.0); // zero out Full LHS since we will inject values as we get them
414 
415  // Construct Reduced RHS
416 
417  // First compute influence of already-known values of X on RHS
418  tempX_ = new Epetra_MultiVector(FullMatrixDomainMap(), NumVectors);
419  tempB_ = new Epetra_MultiVector(FullRHS->Map(), NumVectors);
420 
421  //Inject known X values into tempX for purpose of computing tempB = FullMatrix*tempX
422  // Also inject into full X since we already know the solution
423 
424  if (FullMatrix()->RowMatrixImporter()!=0) {
425  EPETRA_CHK_ERR(tempX_->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
426  EPETRA_CHK_ERR(FullLHS->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
427  }
428  else {
429  tempX_->Update(1.0, *tempExportX_, 0.0);
430  FullLHS->Update(1.0, *tempExportX_, 0.0);
431  }
432 
433 
434  EPETRA_CHK_ERR(FullMatrix()->Multiply(false, *tempX_, *tempB_));
435 
436  EPETRA_CHK_ERR(tempB_->Update(1.0, *FullRHS, -1.0)); // tempB now has influence of already-known X values
437 
439  EPETRA_CHK_ERR(ReducedRHS_->Import(*tempB_, *Full2ReducedRHSImporter_, Insert));
440 
441  // Finally construct Reduced Linear Problem
442 
443  ReducedProblem_ = new Epetra_LinearProblem(ReducedMatrix_, ReducedLHS_, ReducedRHS_);
444 
445  double fn = (double) FullMatrix()->NumGlobalRows64();
446  double fnnz = (double) FullMatrix()->NumGlobalNonzeros64();
447  double rn = (double) ReducedMatrix()->NumGlobalRows64();
448  double rnnz = (double) ReducedMatrix()->NumGlobalNonzeros64();
449 
450  RatioOfDimensions_ = (fn-rn)/fn;
451  RatioOfNonzeros_ = (fnnz-rnnz)/fnnz;
452  HaveReducedProblem_ = true;
453 
454  return(0);
455 }
456 
457 //==============================================================================
459 
460  int i, j;
461 
462  if (Problem==0) EPETRA_CHK_ERR(-1); // Null problem pointer
463 
464  FullProblem_ = Problem;
465  FullMatrix_ = dynamic_cast<Epetra_RowMatrix *>(Problem->GetMatrix());
466  if (FullMatrix_==0) EPETRA_CHK_ERR(-2); // Need a RowMatrix
467  if (Problem->GetRHS()==0) EPETRA_CHK_ERR(-3); // Need a RHS
468  if (Problem->GetLHS()==0) EPETRA_CHK_ERR(-4); // Need a LHS
469  if (!HaveReducedProblem_) EPETRA_CHK_ERR(-5); // Must have set up reduced problem
470 
471  // Create pointer to Full RHS, LHS
472  Epetra_MultiVector * FullRHS = FullProblem()->GetRHS();
473  Epetra_MultiVector * FullLHS = FullProblem()->GetLHS();
474  int NumVectors = FullLHS->NumVectors();
475 
476  int NumEntries;
477  int * Indices;
478  double * Values;
479  int NumMyRows = FullMatrix()->NumMyRows();
480  int ColSingletonCounter = 0;
481  for (i=0; i<NumMyRows; i++) {
482  int curGRID = FullMatrixRowMap().GID64(i); // CJ FIXME TODO long long
483  if (ReducedMatrixRowMap()->MyGID(curGRID)) { // Check if this row should go into reduced matrix
484  EPETRA_CHK_ERR(GetRowGCIDs(i, NumEntries, Values, Indices)); // Get current row (indices global)
485  int ierr = ReducedMatrix()->ReplaceGlobalValues(curGRID, NumEntries,
486  Values, Indices);
487  // Positive errors will occur because we are submitting col entries that are not part of
488  // reduced system. However, because we specified a column map to the ReducedMatrix constructor
489  // these extra column entries will be ignored and we will be politely reminded by a positive
490  // error code
491  if (ierr<0) EPETRA_CHK_ERR(ierr);
492  }
493  // Otherwise if singleton row we explicitly eliminate this row and solve for corresponding X value
494  else {
495  EPETRA_CHK_ERR(GetRow(i, NumEntries, Values, Indices)); // Get current row
496  if (NumEntries==1) {
497  double pivot = Values[0];
498  if (pivot==0.0) EPETRA_CHK_ERR(-1); // Encountered zero row, unable to continue
499  int indX = Indices[0];
500  for (j=0; j<NumVectors; j++)
501  (*tempExportX_)[j][indX] = (*FullRHS)[j][i]/pivot;
502  }
503  // Otherwise, this is a singleton column and we will scan for the pivot element needed
504  // for post-solve equations
505  else {
506  j = ColSingletonPivotLIDs_[ColSingletonCounter];
507  double pivot = Values[j];
508  if (pivot==0.0) EPETRA_CHK_ERR(-2); // Encountered zero column, unable to continue
509  ColSingletonPivots_[ColSingletonCounter] = pivot;
510  ColSingletonCounter++;
511  }
512  }
513  }
514 
515  assert(ColSingletonCounter==NumMyColSingletons_); // Sanity test
516 
517  // Update Reduced LHS (Puts any initial guess values into reduced system)
518 
519  ReducedLHS_->PutScalar(0.0); // zero out Reduced LHS
521  FullLHS->PutScalar(0.0); // zero out Full LHS since we will inject values as we get them
522 
523  // Construct Reduced RHS
524 
525  // Zero out temp space
526  tempX_->PutScalar(0.0);
527  tempB_->PutScalar(0.0);
528 
529  //Inject known X values into tempX for purpose of computing tempB = FullMatrix*tempX
530  // Also inject into full X since we already know the solution
531 
532  if (FullMatrix()->RowMatrixImporter()!=0) {
533  EPETRA_CHK_ERR(tempX_->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
534  EPETRA_CHK_ERR(FullLHS->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
535  }
536  else {
537  tempX_->Update(1.0, *tempExportX_, 0.0);
538  FullLHS->Update(1.0, *tempExportX_, 0.0);
539  }
540 
541 
542  EPETRA_CHK_ERR(FullMatrix()->Multiply(false, *tempX_, *tempB_));
543 
544  EPETRA_CHK_ERR(tempB_->Update(1.0, *FullRHS, -1.0)); // tempB now has influence of already-known X values
545 
546  ReducedRHS_->PutScalar(0.0);
548  return(0);
549 }
550 
551 //==============================================================================
553  Epetra_Export * & RedistributeExporter,
554  Epetra_Map * & RedistributeMap) {
555 
556  int IndexBase = SourceMap->IndexBase(); // CJ TODO FIXME long long
557  if (IndexBase!=TargetMap->IndexBase()) EPETRA_CHK_ERR(-1);
558 
559  const Epetra_Comm & Comm = TargetMap->Comm();
560 
561  int TargetNumMyElements = TargetMap->NumMyElements();
562  int SourceNumMyElements = SourceMap->NumMyElements();
563 
564  // ContiguousTargetMap has same number of elements per PE as TargetMap, but uses contigious indexing
565  Epetra_Map ContiguousTargetMap(-1, TargetNumMyElements, IndexBase,Comm);
566 
567  // Same for ContiguousSourceMap
568  Epetra_Map ContiguousSourceMap(-1, SourceNumMyElements, IndexBase, Comm);
569 
570  assert(ContiguousSourceMap.NumGlobalElements64()==ContiguousTargetMap.NumGlobalElements64());
571 
572  // Now create a vector that contains the global indices of the Source Epetra_MultiVector
573 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
574  Epetra_IntVector *SourceIndices = 0;
575  Epetra_IntVector *TargetIndices = 0;
576 #endif
577 
578 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
579  Epetra_LongLongVector *SourceIndices_LL = 0;
580  Epetra_LongLongVector *TargetIndices_LL = 0;
581 #endif
582 
583  if(SourceMap->GlobalIndicesInt())
584 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
585  SourceIndices = new Epetra_IntVector(View, ContiguousSourceMap, SourceMap->MyGlobalElements());
586 #else
587  throw "Epetra_CrsSingletonFilter::ConstructRedistributeExporter: GlobalIndicesInt but no int API";
588 #endif
589  else if(SourceMap->GlobalIndicesLongLong())
590 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
591  SourceIndices_LL = new Epetra_LongLongVector(View, ContiguousSourceMap, SourceMap->MyGlobalElements64());
592 #else
593  throw "Epetra_CrsSingletonFilter::ConstructRedistributeExporter: GlobalIndicesLongLong but no long long API";
594 #endif
595  else
596  throw "Epetra_CrsSingletonFilter::ConstructRedistributeExporter: Unknown global index type";
597 
598  // Create an exporter to send the SourceMap global IDs to the target distribution
599  Epetra_Export Exporter(ContiguousSourceMap, ContiguousTargetMap);
600 
601  // Create a vector to catch the global IDs in the target distribution
602  if(TargetMap->GlobalIndicesInt()) {
603 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
604  TargetIndices = new Epetra_IntVector(ContiguousTargetMap);
605  TargetIndices->Export(*SourceIndices, Exporter, Insert);
606 #else
607  throw "Epetra_CrsSingletonFilter::ConstructRedistributeExporter: GlobalIndicesInt but no int API";
608 #endif
609  } else if(TargetMap->GlobalIndicesLongLong()) {
610 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
611  TargetIndices_LL = new Epetra_LongLongVector(ContiguousTargetMap);
612  TargetIndices_LL->Export(*SourceIndices_LL, Exporter, Insert);
613 #else
614  throw "Epetra_CrsSingletonFilter::ConstructRedistributeExporter: GlobalIndicesLongLong but no long long API";
615 #endif
616  }
617  else
618  throw "Epetra_CrsSingletonFilter::ConstructRedistributeExporter: Unknown global index type";
619 
620  // Create a new map that describes how the Source MultiVector should be laid out so that it has
621  // the same number of elements on each processor as the TargetMap
622 
623  if(TargetMap->GlobalIndicesInt())
624 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
625  RedistributeMap = new Epetra_Map(-1, TargetNumMyElements, TargetIndices->Values(), IndexBase, Comm);
626 #else
627  throw "Epetra_CrsSingletonFilter::ConstructRedistributeExporter: GlobalIndicesInt but no int API";
628 #endif
629  else if(TargetMap->GlobalIndicesLongLong())
630 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
631  RedistributeMap = new Epetra_Map((long long) -1, TargetNumMyElements, TargetIndices_LL->Values(), IndexBase, Comm);
632 #else
633  throw "Epetra_CrsSingletonFilter::ConstructRedistributeExporter: GlobalIndicesLongLong but no long long API";
634 #endif
635  else
636  throw "Epetra_CrsSingletonFilter::ConstructRedistributeExporter: Unknown global index type";
637 
638  // This exporter will finally redistribute the Source MultiVector to the same layout as the TargetMap
639  RedistributeExporter = new Epetra_Export(*SourceMap, *RedistributeMap);
640  return(0);
641 }
642 //==============================================================================
644 
645  int jj, k;
646 
647  Epetra_MultiVector * FullLHS = FullProblem()->GetLHS();
648  Epetra_MultiVector * FullRHS = FullProblem()->GetRHS();
649 
650  tempX_->PutScalar(0.0); tempExportX_->PutScalar(0.0);
651  // Inject values that the user computed for the reduced problem into the full solution vector
653  FullLHS->Update(1.0, *tempX_, 1.0);
654 
655  // Next we will use our full solution vector which is populated with pre-filter solution
656  // values and reduced system solution values to compute the sum of the row contributions
657  // that must be subtracted to get the post-filter solution values
658 
659  EPETRA_CHK_ERR(FullMatrix()->Multiply(false, *FullLHS, *tempB_));
660 
661 
662 
663  // Finally we loop through the local rows that were associated with column singletons and compute the
664  // solution for these equations.
665 
666  int NumVectors = tempB_->NumVectors();
667  for (k=0; k<NumMyColSingletons_; k++) {
668  int i = ColSingletonRowLIDs_[k];
669  int j = ColSingletonColLIDs_[k];
670  double pivot = ColSingletonPivots_[k];
671  for (jj=0; jj<NumVectors; jj++)
672  (*tempExportX_)[jj][j]= ((*FullRHS)[jj][i] - (*tempB_)[jj][i])/pivot;
673  }
674 
675  // Finally, insert values from post-solve step and we are done!!!!
676 
677 
678  if (FullMatrix()->RowMatrixImporter()!=0) {
679  EPETRA_CHK_ERR(tempX_->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
680  }
681  else {
682  tempX_->Update(1.0, *tempExportX_, 0.0);
683  }
684 
685  FullLHS->Update(1.0, *tempX_, 1.0);
686 
687  return(0);
688 }
689 //==============================================================================
691 
693 
694  // Cast to CrsMatrix, if possible. Can save some work.
695  FullCrsMatrix_ = dynamic_cast<Epetra_CrsMatrix *>(FullMatrix());
696  FullMatrixIsCrsMatrix_ = (FullCrsMatrix_!=0); // Pointer is non-zero if cast worked
697  Indices_ = new int[MaxNumMyEntries_];
699 
700  return(0);
701 }
702 //==============================================================================
703 int Epetra_CrsSingletonFilter::GetRow(int Row, int & NumIndices, int * & Indices) {
704 
705  if (FullMatrixIsCrsMatrix_) { // View of current row
706  EPETRA_CHK_ERR(FullCrsMatrix()->Graph().ExtractMyRowView(Row, NumIndices, Indices));
707  }
708  else { // Copy of current row (we must get the values, but we ignore them)
709  EPETRA_CHK_ERR(FullMatrix()->ExtractMyRowCopy(Row, MaxNumMyEntries_, NumIndices,
710  Values_.Values(), Indices_));
711  Indices = Indices_;
712  }
713  return(0);
714 }
715 //==============================================================================
716 int Epetra_CrsSingletonFilter::GetRow(int Row, int & NumIndices,
717  double * & Values, int * & Indices) {
718 
719  if (FullMatrixIsCrsMatrix_) { // View of current row
720  EPETRA_CHK_ERR(FullCrsMatrix_->ExtractMyRowView(Row, NumIndices, Values, Indices));
721  }
722  else { // Copy of current row (we must get the values, but we ignore them)
723  EPETRA_CHK_ERR(FullMatrix()->ExtractMyRowCopy(Row, MaxNumMyEntries_, NumIndices,
724  Values_.Values(), Indices_));
725  Values = Values_.Values();
726  Indices = Indices_;
727  }
728  return(0);
729 }
730 //==============================================================================
731 int Epetra_CrsSingletonFilter::GetRowGCIDs(int Row, int & NumIndices,
732  double * & Values, int * & GlobalIndices) {
733 
734  EPETRA_CHK_ERR(FullMatrix()->ExtractMyRowCopy(Row, MaxNumMyEntries_, NumIndices,
735  Values_.Values(), Indices_));
736  for (int j=0; j<NumIndices; j++) Indices_[j] = FullMatrixColMap().GID64(Indices_[j]); // FIXME long long
737  Values = Values_.Values();
738  GlobalIndices = Indices_;
739  return(0);
740 }
741 //==============================================================================
743  const Epetra_MapColoring & rowMapColors,
744  const Epetra_IntVector & ColProfiles,
745  const Epetra_IntVector & NewColProfiles,
746  const Epetra_IntVector & ColHasRowWithSingleton) {
747 
748  int j;
749 
750  if (NumMyColSingletons_==0) return(0); // Nothing to do
751 
752  Epetra_MapColoring & colMapColors = *ColMapColors_;
753 
754  int NumMyCols = FullMatrix()->NumMyCols();
755 
756  // We will need these arrays for the post-solve phase
761 
762  // Register singleton columns (that were not already counted as singleton rows)
763  // Check to see if any columns disappeared because all associated rows were eliminated
764  int NumMyColSingletonstmp = 0;
765  for (j=0; j<NumMyCols; j++) {
766  int i = RowIDs[j];
767  if ( ColProfiles[j]==1 && rowMapColors[i]!=1) {
768  ColSingletonRowLIDs_[NumMyColSingletonstmp] = i;
769  ColSingletonColLIDs_[NumMyColSingletonstmp] = j;
770  NumMyColSingletonstmp++;
771  }
772  // Also check for columns that were eliminated implicitly by
773  // having all associated row eliminated
774  else if (NewColProfiles[j]==0 && ColHasRowWithSingleton[j]!=1 && rowMapColors[i]==0) {
775  colMapColors[j] = 1;
776  }
777  }
778 
779  assert(NumMyColSingletonstmp==NumMyColSingletons_); //Sanity check
780  Epetra_Util sorter;
781  sorter.Sort(true, NumMyColSingletons_, ColSingletonRowLIDs_, 0, 0, 1, &ColSingletonColLIDs_, 0, 0);
782 
783  return(0);
784 }
785 
786 #endif // EPETRA_NO_32BIT_GLOBAL_INDICES
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int MaxValue()
Find maximum value.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
Epetra_MapColoring: A class for coloring Epetra_Map and Epetra_BlockMap objects.
Epetra_CrsSingletonFilter()
Epetra_CrsSingletonFilter default constructor.
int Size(int Length_in)
Set length of a Epetra_SerialDenseVector object; init values to zero.
virtual const Epetra_Map & RowMatrixRowMap() const =0
Returns the Epetra_Map object associated with the rows of this matrix.
Epetra_LongLongVector: A class for constructing and using dense integer vectors on a parallel compute...
long long NumGlobalRows64() const
Epetra_MultiVector * GetLHS() const
Get a pointer to the left-hand-side X.
Epetra_MultiVector * GetRHS() const
Get a pointer to the right-hand-side B.
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
bool SameAs(const Epetra_BlockMap &Map) const
Returns true if this and Map are identical maps.
Epetra_SerialDenseVector Values_
Epetra_IntVector: A class for constructing and using dense integer vectors on a parallel computer...
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
long long NumGlobalElements64() const
int UpdateReducedProblem(Epetra_LinearProblem *Problem)
Update a reduced linear problem using new values.
Epetra_LinearProblem * FullProblem_
virtual ~Epetra_CrsSingletonFilter()
Epetra_CrsSingletonFilter Destructor.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Insert a list of elements in a given global row of the matrix.
int ConstructRedistributeExporter(Epetra_Map *SourceMap, Epetra_Map *TargetMap, Epetra_Export *&RedistributeExporter, Epetra_Map *&RedistributeMap)
int Analyze(Epetra_RowMatrix *FullMatrix)
Analyze the input matrix, removing row/column pairs that have singletons.
#define EPETRA_CHK_ERR(a)
const Epetra_Map & FullMatrixColMap() const
int ComputeFullSolution()
Compute a solution for the full problem using the solution of the reduced problem, put in LHS of FullProblem().
Epetra_Export: This class builds an export object for efficient exporting of off-processor elements...
Definition: Epetra_Export.h:62
int PutValue(int Value)
Set all elements of the vector to Value.
Epetra_LinearProblem * ReducedProblem_
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
Epetra_Import: This class builds an import object for efficient importing of off-processor elements...
Definition: Epetra_Import.h:63
int * Values() const
Returns a pointer to an array containing the values of this vector.
int NumVectors() const
Returns the number of vectors in the multi-vector.
virtual int NumMyCols() const =0
Returns the number of matrix columns owned by the calling processor.
int IndexBase() const
Index base for this map.
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
Returns a view of the specified local row values via pointers to internal data.
long long NumGlobalNonzeros64() const
const Epetra_Map & FullMatrixRangeMap() const
Epetra_Util: The Epetra Util Wrapper Class.
Definition: Epetra_Util.h:79
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
int NumMyElements() const
Number of elements on the calling processor.
virtual int MaxNumEntries() const =0
Returns the maximum of NumMyRowEntries() over all rows.
long long GID64(int LID) const
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Exports an Epetra_DistObject using the Epetra_Import object.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
virtual const Epetra_Import * RowMatrixImporter() const =0
Returns the Epetra_Import object that contains the import operations for distributed operations...
Epetra_Map * ReducedMatrixDomainMap() const
Returns pointer to Epetra_Map describing the domain map for the reduced system.
int CreatePostSolveArrays(const Epetra_IntVector &RowIDs, const Epetra_MapColoring &RowMapColors, const Epetra_IntVector &ColProfiles, const Epetra_IntVector &NewColProfiles, const Epetra_IntVector &ColHasRowWithSingleton)
Epetra_Map * ReducedMatrixRowMap() const
Returns pointer to Epetra_Map describing the reduced system row distribution.
virtual int NumMyRows() const =0
Returns the number of matrix rows owned by the calling processor.
const Epetra_Map & FullMatrixDomainMap() const
int GetRow(int Row, int &NumIndices, int *&Indices)
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
Epetra_Map * ReducedMatrixRangeMap() const
Returns pointer to Epetra_Map describing the range map for the reduced system.
Epetra_RowMatrix * FullMatrix() const
Returns pointer to Epetra_CrsMatrix from full problem.
Epetra_CrsMatrix * FullCrsMatrix() const
virtual long long NumGlobalNonzeros64() const =0
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
Epetra_Map * GenerateMap(int Color) const
Generates an Epetra_Map of the GIDs associated with the specified color.
Epetra_RowMatrix * GetMatrix() const
Get a pointer to the matrix A.
Epetra_Map * ReducedMatrixColMap() const
Returns pointer to Epetra_Map describing the reduced system column distribution.
int ConstructReducedProblem(Epetra_LinearProblem *Problem)
Return a reduced linear problem based on results of Analyze().
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
Epetra_LinearProblem * FullProblem() const
Returns pointer to the original unreduced Epetra_LinearProblem.
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_Util Sort Routine (Shell sort)
Epetra_CrsMatrix * ReducedMatrix() const
Returns pointer to Epetra_CrsMatrix from full problem.
Epetra_LinearProblem: The Epetra Linear Problem Class.
double * Values() const
Returns pointer to the values in vector.
long long * MyGlobalElements64() const
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices...
int GetRowGCIDs(int Row, int &NumIndices, double *&Values, int *&GlobalIndices)
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Replace specified existing values with this list of entries for a given global row of the matrix...
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Imports an Epetra_DistObject using the Epetra_Import object.
const Epetra_Map & FullMatrixRowMap() const
virtual long long NumGlobalRows64() const =0
long long * Values() const
Returns a pointer to an array containing the values of this vector.
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.