EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraExt_CrsSingletonFilter_LinearProblem.cpp
Go to the documentation of this file.
1 
2 //@HEADER
3 // ***********************************************************************
4 //
5 // EpetraExt: Epetra Extended - 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 #include "Epetra_ConfigDefs.h"
44 #include "Epetra_Map.h"
45 #include "Epetra_Util.h"
46 #include "Epetra_Export.h"
47 #include "Epetra_Import.h"
48 #include "Epetra_MultiVector.h"
49 #include "Epetra_Vector.h"
50 #include "Epetra_GIDTypeVector.h"
51 #include "Epetra_Comm.h"
52 #include "Epetra_LinearProblem.h"
53 #include "Epetra_MapColoring.h"
55 
56 namespace EpetraExt {
57 
58 //==============================================================================
61 : verbose_(verbose)
62 {
64 }
65 //==============================================================================
68 {
69  if (ReducedLHS_!=0) delete ReducedLHS_;
70  if (ReducedRHS_!=0) delete ReducedRHS_;
80  if (RowMapColors_!=0) delete RowMapColors_;
81  if (ColMapColors_!=0) delete ColMapColors_;
82 
83  if (ColSingletonRowLIDs_ != 0) delete [] ColSingletonRowLIDs_;
84  if (ColSingletonColLIDs_ != 0) delete [] ColSingletonColLIDs_;
86  if (ColSingletonPivots_ != 0) delete [] ColSingletonPivots_;
87  if (tempExportX_ != 0) delete tempExportX_;
88  if (Indices_int_ != 0) delete [] Indices_int_;
89 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
90  if (Indices_LL_ != 0) delete [] Indices_LL_;
91 #endif
92  if (tempX_ != 0) delete tempX_;
93  if (tempB_ != 0) delete tempB_;
94 
95 }
96 
97 //==============================================================================
101 {
102  analyze( orig );
103  return construct();
104 }
105 
106 //==============================================================================
107 bool
110 {
111  origObj_ = &orig;
112 
113  FullMatrix_ = orig.GetMatrix();
114 
115 #ifdef NDEBUG
116  (void) Analyze( FullMatrix_ );
117 #else
118  // assert() statements go away if NDEBUG is defined. Don't declare
119  // the 'flag' variable if it never gets used.
120  int flag = Analyze( FullMatrix_ );
121  assert( flag >= 0 );
122 #endif // NDEBUG
123 
124  if ( verbose_ && FullMatrix_->Comm().MyPID()==0 ) {
125  std::cout << "\nAnalyzed Singleton Problem:\n";
126  std::cout << "---------------------------\n";
127  }
128  if ( SingletonsDetected() ) {
129  if ( verbose_ && FullMatrix_->Comm().MyPID()==0 ) {
130  std::cout << "Singletons Detected!" << std::endl;;
131  std::cout << "Num Singletons: " << NumSingletons() << std::endl;
132  }
133  }
134  else {
135  if ( verbose_ && FullMatrix_->Comm().MyPID()==0 )
136  std::cout << "No Singletons Detected!" << std::endl;
137  }
138 /*
139  std::cout << "List of Row Singletons:\n";
140  int * slist = RowMapColors_->ColorLIDList(1);
141  for( int i = 0; i < RowMapColors_->NumElementsWithColor(1); ++i )
142  std::cout << slist[i] << " ";
143  std::cout << "\n";
144  std::cout << "List of Col Singletons:\n";
145  slist = ColMapColors_->ColorLIDList(1);
146  for( int i = 0; i < ColMapColors_->NumElementsWithColor(1); ++i )
147  std::cout << slist[i] << " ";
148  std::cout << "\n";
149 */
150  if ( verbose_ && FullMatrix_->Comm().MyPID()==0 )
151  std::cout << "---------------------------\n\n";
152 
153  return true;
154 }
155 
156 //==============================================================================
160 {
161  if( !origObj_ ) abort();
162 
163 #ifdef NDEBUG
165 #else
166  // assert() statements go away if NDEBUG is defined. Don't declare
167  // the 'flag' variable if it never gets used.
168  int flag = ConstructReducedProblem( origObj_ );
169  assert( flag >= 0 );
170 #endif // NDEBUG
171 
173 
174  if( verbose_ && SingletonsDetected() && FullMatrix_->Comm().MyPID()==0 )
175  {
176  std::cout << "\nConstructedSingleton Problem:\n";
177  std::cout << "---------------------------\n";
178  std::cout << "RatioOfDimensions: " << RatioOfDimensions() << std::endl;
179  std::cout << "RatioOfNonzeros: " << RatioOfNonzeros() << std::endl;
180  std::cout << "---------------------------\n\n";
181  }
182 
183  return *newObj_;
184 }
185 
186 //==============================================================================
188 {
189  int ierr = UpdateReducedProblem( FullProblem_ );
190  if( ierr ) std::cout << "EDT_LinearProblem_CrsSingletonFilter::UpdateReducedProblem FAILED!\n";
191 
192  return (ierr==0);
193 }
194 
195 //==============================================================================
197 {
198  int ierr = ComputeFullSolution();
199  if( ierr ) std::cout << "EDT_LinearProblem_CrsSingletonFilter::ComputeFullSolution FAILED!\n";
200 
201  return (ierr==0);
202 }
203 
204 //==============================================================================
206 
207 // Initialize all attributes that have trivial default values
208 
209  FullProblem_ = 0;
210  FullMatrix_ = 0;
211  ReducedRHS_ = 0;
212  ReducedLHS_ = 0;
221 
226 
227  AbsoluteThreshold_ = 0;
228  RelativeThreshold_ = 0;
229 
230  NumMyRowSingletons_ = -1;
231  NumMyColSingletons_ = -1;
234  RatioOfDimensions_ = -1.0;
235  RatioOfNonzeros_ = -1.0;
236 
237  HaveReducedProblem_ = false;
239  AnalysisDone_ = false;
240  SymmetricElimination_ = true;
241 
242  tempExportX_ = 0;
243  tempX_ = 0;
244  tempB_ = 0;
245 
246  Indices_int_ = 0;
247 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
248  Indices_LL_ = 0;
249 #endif
250 
251  RowMapColors_ = 0;
252  ColMapColors_ = 0;
253 
254  FullMatrixIsCrsMatrix_ = false;
255  MaxNumMyEntries_ = 0;
256  return;
257 }
258 //==============================================================================
260  FullMatrix_ = fullMatrix;
261 
262  if (AnalysisDone_) EPETRA_CHK_ERR(-1); // Analysis already done once. Cannot do it again
263  if (fullMatrix==0) EPETRA_CHK_ERR(-2); // Input matrix pointer is zero
264  if (fullMatrix->NumGlobalRows64()==0) EPETRA_CHK_ERR(-3); // Full matrix has zero dimension.
265  if (fullMatrix->NumGlobalNonzeros64()==0) EPETRA_CHK_ERR(-4); // Full matrix has no nonzero terms.
266 
267  // First check for columns with single entries and find columns with singleton rows
268  Epetra_IntVector ColProfiles(FullMatrixColMap()); ColProfiles.PutValue(0);
269  Epetra_IntVector ColHasRowWithSingleton(FullMatrixColMap()); ColHasRowWithSingleton.PutValue(0);
270 
271  // RowIDs[j] will contain the local row ID associated with the jth column,
272  // if the jth col has a single entry
273  Epetra_IntVector RowIDs(FullMatrixColMap()); RowIDs.PutValue(-1);
274 
275  // Define MapColoring objects
276  RowMapColors_ = new Epetra_MapColoring(FullMatrixRowMap()); // Initial colors are all 0
278  Epetra_MapColoring & rowMapColors = *RowMapColors_;
279  Epetra_MapColoring & colMapColors = *ColMapColors_;
280 
281  int NumMyRows = fullMatrix->NumMyRows();
282  int NumMyCols = fullMatrix->NumMyCols();
283 
284  // Set up for accessing full matrix. Will do so row-by-row.
286 
287  // Scan matrix for singleton rows, build up column profiles
288  int NumIndices;
289  int * Indices;
291  for (int i=0; i<NumMyRows; i++) {
292  // Get ith row
293  EPETRA_CHK_ERR(GetRow(i, NumIndices, Indices));
294  for (int j=0; j<NumIndices; j++) {
295  int ColumnIndex = Indices[j];
296  ColProfiles[ColumnIndex]++; // Increment column count
297  // Record local row ID for current column
298  // will use to identify row to eliminate if column is a singleton
299  RowIDs[ColumnIndex] = i;
300  }
301  // If row has single entry, color it and associated column with color=1
302  if (NumIndices==1) {
303  int j2 = Indices[0];
304  ColHasRowWithSingleton[j2]++;
305  rowMapColors[i] = 1;
306  colMapColors[j2] = 1;
308  }
309  }
310 
311  // 1) The vector ColProfiles has column nonzero counts for each processor's contribution
312  // Combine these to get total column profile information and then redistribute to processors
313  // so each can determine if it is the owner of the row associated with the singleton column
314  // 2) The vector ColHasRowWithSingleton[i] contain count of singleton rows that are associated with
315  // the ith column on this processor. Must tell other processors that they should also eliminate
316  // these columns.
317 
318  // Make a copy of ColProfiles for later use when detecting columns that disappear locally
319 
320  Epetra_IntVector NewColProfiles(ColProfiles);
321 
322  // If RowMatrixImporter is non-trivial, we need to perform a gather/scatter to accumulate results
323 
324  if (fullMatrix->RowMatrixImporter()!=0) {
325  Epetra_IntVector tmpVec(FullMatrixDomainMap()); // Use for gather/scatter of column vectors
326  EPETRA_CHK_ERR(tmpVec.PutValue(0));
327  EPETRA_CHK_ERR(tmpVec.Export(ColProfiles, *fullMatrix->RowMatrixImporter(), Add));
328  EPETRA_CHK_ERR(ColProfiles.Import(tmpVec, *fullMatrix->RowMatrixImporter(), Insert));
329 
330  EPETRA_CHK_ERR(tmpVec.PutValue(0));
331  EPETRA_CHK_ERR(tmpVec.Export(ColHasRowWithSingleton, *fullMatrix->RowMatrixImporter(), Add));
332  EPETRA_CHK_ERR(ColHasRowWithSingleton.Import(tmpVec, *fullMatrix->RowMatrixImporter(), Insert));
333  }
334  // ColProfiles now contains the nonzero column entry count for all columns that have
335  // an entry on this processor.
336  // ColHasRowWithSingleton now contains a count of singleton rows associated with the corresponding
337  // local column. Next we check to make sure no column is associated with more than one singleton row.
338 
339  if (ColHasRowWithSingleton.MaxValue()>1) {
340  EPETRA_CHK_ERR(-2); // At least one col is associated with two singleton rows, can't handle it.
341  }
342 
343  Epetra_IntVector RowHasColWithSingleton(fullMatrix->RowMatrixRowMap()); // Use to check for errors
344  RowHasColWithSingleton.PutValue(0);
345 
347  // Count singleton columns (that were not already counted as singleton rows)
348  for (int j=0; j<NumMyCols; j++) {
349  int i2 = RowIDs[j];
350  // Check if column is a singleton
351  if (ColProfiles[j]==1 ) {
352  // Check to make sure RowID is not invalid
353 // assert(i!=-1);
354  // Check to see if this column already eliminated by the row check above
355  if (rowMapColors[i2]!=1) {
356  RowHasColWithSingleton[i2]++; // Increment col singleton counter for ith row
357  rowMapColors[i2] = 2; // Use 2 for now, to distinguish between row eliminated directly or via column singletons
358  colMapColors[j] = 1;
360  // If we delete a row, we need to keep track of associated column entries that were also deleted
361  // in case all entries in a column are eventually deleted, in which case the column should
362  // also be deleted.
363  EPETRA_CHK_ERR(GetRow(i2, NumIndices, Indices));
364  for (int jj=0; jj<NumIndices; jj++) {
365  NewColProfiles[Indices[jj]]--;
366  }
367  }
368  }
369  // Check if some other processor eliminated this column
370  else if (ColHasRowWithSingleton[j]==1 && rowMapColors[i2]!=1) {
371  colMapColors[j] = 1;
372  }
373  }
374  if (RowHasColWithSingleton.MaxValue()>1) {
375  EPETRA_CHK_ERR(-3); // At lease one row is associated with two singleton cols, can't handle it.
376  }
377 
378  // Generate arrays that keep track of column singleton row, col and pivot info needed for post-solve phase
379  EPETRA_CHK_ERR(CreatePostSolveArrays(RowIDs, rowMapColors, ColProfiles, NewColProfiles,
380  ColHasRowWithSingleton));
381 
382  for (int i=0; i<NumMyRows; i++) {
383  if (rowMapColors[i]==2) rowMapColors[i] = 1; // Convert all eliminated rows to same color
384  }
385 
388  AnalysisDone_ = true;
389  return(0);
390 }
391 //==============================================================================
392 template<typename int_type>
394  int i, j;
395  if (HaveReducedProblem_) EPETRA_CHK_ERR(-1); // Setup already done once. Cannot do it again
396  if (Problem==0) EPETRA_CHK_ERR(-2); // Null problem pointer
397 
398  FullProblem_ = Problem;
399  FullMatrix_ = dynamic_cast<Epetra_RowMatrix *>(Problem->GetMatrix());
400  if (FullMatrix_==0) EPETRA_CHK_ERR(-3); // Need a RowMatrix
401  if (Problem->GetRHS()==0) EPETRA_CHK_ERR(-4); // Need a RHS
402  if (Problem->GetLHS()==0) EPETRA_CHK_ERR(-5); // Need a LHS
403  // Generate reduced row and column maps
404 
405  if ( SingletonsDetected() ) {
406 
407  Epetra_MapColoring & rowMapColors = *RowMapColors_;
408  Epetra_MapColoring & colMapColors = *ColMapColors_;
409 
410  ReducedMatrixRowMap_ = rowMapColors.GenerateMap(0);
411  ReducedMatrixColMap_ = colMapColors.GenerateMap(0);
412 
413  // Create domain and range map colorings by exporting map coloring of column and row maps
414 
415  if (FullMatrix()->RowMatrixImporter()!=0) {
416  Epetra_MapColoring DomainMapColors(FullMatrixDomainMap());
417  EPETRA_CHK_ERR(DomainMapColors.Export(*ColMapColors_, *FullMatrix()->RowMatrixImporter(), AbsMax));
418  OrigReducedMatrixDomainMap_ = DomainMapColors.GenerateMap(0);
419  }
420  else
422 
424  if (FullCrsMatrix()->Exporter()!=0) { // Non-trivial exporter
425  Epetra_MapColoring RangeMapColors(FullMatrixRangeMap());
426  EPETRA_CHK_ERR(RangeMapColors.Export(*RowMapColors_, *FullCrsMatrix()->Exporter(),
427  AbsMax));
428  ReducedMatrixRangeMap_ = RangeMapColors.GenerateMap(0);
429  }
430  else
432  }
433  else
435 
436  // Check to see if the reduced system domain and range maps are the same.
437  // If not, we need to remap entries of the LHS multivector so that they are distributed
438  // conformally with the rows of the reduced matrix and the RHS multivector
443  else {
447  }
448 
449  // Create pointer to Full RHS, LHS
450  Epetra_MultiVector * FullRHS = FullProblem()->GetRHS();
451  Epetra_MultiVector * FullLHS = FullProblem()->GetLHS();
452  int NumVectors = FullLHS->NumVectors();
453 
454  // Create importers
457 
458  // Construct Reduced Matrix
460 
461  // Create storage for temporary X values due to explicit elimination of rows
462  tempExportX_ = new Epetra_MultiVector(FullMatrixColMap(), NumVectors);
463 
464  int NumEntries;
465  double * Values;
466  int NumMyRows = FullMatrix()->NumMyRows();
467  int ColSingletonCounter = 0;
468  for (i=0; i<NumMyRows; i++) {
469  int_type curGRID = (int_type) FullMatrixRowMap().GID64(i);
470  if (ReducedMatrixRowMap()->MyGID(curGRID)) { // Check if this row should go into reduced matrix
471  int_type * Indices;
472  EPETRA_CHK_ERR(GetRowGCIDs(i, NumEntries, Values, Indices)); // Get current row (Indices are global)
473 
474  int ierr = ReducedMatrix()->InsertGlobalValues(curGRID, NumEntries,
475  Values, Indices); // Insert into reduce matrix
476  // Positive errors will occur because we are submitting col entries that are not part of
477  // reduced system. However, because we specified a column map to the ReducedMatrix constructor
478  // these extra column entries will be ignored and we will be politely reminded by a positive
479  // error code
480  if (ierr<0) EPETRA_CHK_ERR(ierr);
481  }
482  else {
483  int * Indices;
484  EPETRA_CHK_ERR(GetRow(i, NumEntries, Values, Indices)); // Get current row
485  if (NumEntries==1) {
486  double pivot = Values[0];
487  if (pivot==0.0) EPETRA_CHK_ERR(-1); // Encountered zero row, unable to continue
488  int indX = Indices[0];
489  for (j=0; j<NumVectors; j++)
490  (*tempExportX_)[j][indX] = (*FullRHS)[j][i]/pivot;
491  }
492  // Otherwise, this is a singleton column and we will scan for the pivot element needed
493  // for post-solve equations
494  else {
495  int targetCol = ColSingletonColLIDs_[ColSingletonCounter];
496  for (j=0; j<NumEntries; j++) {
497  if (Indices[j]==targetCol) {
498  double pivot = Values[j];
499  if (pivot==0.0) EPETRA_CHK_ERR(-2); // Encountered zero column, unable to continue
500  ColSingletonPivotLIDs_[ColSingletonCounter] = j; // Save for later use
501  ColSingletonPivots_[ColSingletonCounter] = pivot;
502  ColSingletonCounter++;
503  break;
504  }
505  }
506  }
507  }
508  }
509 
510  // Now convert to local indexing. We have constructed things so that the domain and range of the
511  // matrix will have the same map. If the reduced matrix domain and range maps were not the same, the
512  // differences were addressed in the ConstructRedistributeExporter() method
514 
515  // 1) The vector ColProfiles has column nonzero counts for each processor's contribution
516  // Construct Reduced LHS (Puts any initial guess values into reduced system)
517 
520  FullLHS->PutScalar(0.0); // zero out Full LHS since we will inject values as we get them
521 
522  // Construct Reduced RHS
523 
524  // First compute influence of already-known values of X on RHS
525  tempX_ = new Epetra_MultiVector(FullMatrixDomainMap(), NumVectors);
526  tempB_ = new Epetra_MultiVector(FullRHS->Map(), NumVectors);
527 
528  //Inject known X values into tempX for purpose of computing tempB = FullMatrix*tempX
529  // Also inject into full X since we already know the solution
530 
531  if (FullMatrix()->RowMatrixImporter()!=0) {
532  EPETRA_CHK_ERR(tempX_->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
533  EPETRA_CHK_ERR(FullLHS->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
534  }
535  else {
536  tempX_->Update(1.0, *tempExportX_, 0.0);
537  FullLHS->Update(1.0, *tempExportX_, 0.0);
538  }
539 
540  EPETRA_CHK_ERR(FullMatrix()->Multiply(false, *tempX_, *tempB_));
541 
542  EPETRA_CHK_ERR(tempB_->Update(1.0, *FullRHS, -1.0)); // tempB now has influence of already-known X values
543 
544  ReducedRHS_ = new Epetra_MultiVector(*ReducedMatrixRowMap(), FullRHS->NumVectors());
545  EPETRA_CHK_ERR(ReducedRHS_->Import(*tempB_, *Full2ReducedRHSImporter_, Insert));
546 
547  // Finally construct Reduced Linear Problem
549  }
550  else {
551 
552  // There are no singletons, so don't bother building a reduced problem.
553  ReducedProblem_ = Teuchos::rcp( Problem, false );
554  ReducedMatrix_ = Teuchos::rcp( dynamic_cast<Epetra_CrsMatrix *>(Problem->GetMatrix()), false );
555  }
556 
557  double fn = (double) FullMatrix()->NumGlobalRows64();
558  double fnnz = (double) FullMatrix()->NumGlobalNonzeros64();
559  double rn = (double) ReducedMatrix()->NumGlobalRows64();
560  double rnnz = (double) ReducedMatrix()->NumGlobalNonzeros64();
561 
562  RatioOfDimensions_ = rn/fn;
563  RatioOfNonzeros_ = rnnz/fnnz;
564  HaveReducedProblem_ = true;
565 
566  return(0);
567 }
568 
570 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
571  if(Problem->GetMatrix()->RowMatrixRowMap().GlobalIndicesInt()) {
572  return TConstructReducedProblem<int>(Problem);
573  }
574  else
575 #endif
576 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
577  if(Problem->GetMatrix()->RowMatrixRowMap().GlobalIndicesLongLong()) {
578  return TConstructReducedProblem<long long>(Problem);
579  }
580  else
581 #endif
582  throw "LinearProblem_CrsSingletonFilter::ConstructReducedProblem: ERROR, GlobalIndices type unknown.";
583 }
584 
585 //==============================================================================
586 template<typename int_type>
588 
589  int i, j;
590 
591  if (Problem==0) EPETRA_CHK_ERR(-1); // Null problem pointer
592 
593  FullProblem_ = Problem;
594  FullMatrix_ = dynamic_cast<Epetra_RowMatrix *>(Problem->GetMatrix());
595  if (FullMatrix_==0) EPETRA_CHK_ERR(-2); // Need a RowMatrix
596  if (Problem->GetRHS()==0) EPETRA_CHK_ERR(-3); // Need a RHS
597  if (Problem->GetLHS()==0) EPETRA_CHK_ERR(-4); // Need a LHS
598  if (!HaveReducedProblem_) EPETRA_CHK_ERR(-5); // Must have set up reduced problem
599 
600  if ( SingletonsDetected() ) {
601 
602  // Create pointer to Full RHS, LHS
603  Epetra_MultiVector * FullRHS = FullProblem()->GetRHS();
604  Epetra_MultiVector * FullLHS = FullProblem()->GetLHS();
605 
606  int NumVectors = FullLHS->NumVectors();
607  tempExportX_->PutScalar(0.0);
608 
609  int NumEntries;
610  double * Values;
611  int NumMyRows = FullMatrix()->NumMyRows();
612  int ColSingletonCounter = 0;
613  for (i=0; i<NumMyRows; i++) {
614  int_type curGRID = (int_type) FullMatrixRowMap().GID64(i);
615  if (ReducedMatrixRowMap()->MyGID(curGRID)) { // Check if this row should go into reduced matrix
616  int_type * Indices;
617  EPETRA_CHK_ERR(GetRowGCIDs(i, NumEntries, Values, Indices)); // Get current row (indices global)
618  int ierr = ReducedMatrix()->ReplaceGlobalValues(curGRID, NumEntries,
619  Values, Indices);
620  // Positive errors will occur because we are submitting col entries that are not part of
621  // reduced system. However, because we specified a column map to the ReducedMatrix constructor
622  // these extra column entries will be ignored and we will be politely reminded by a positive
623  // error code
624  if (ierr<0) EPETRA_CHK_ERR(ierr);
625  }
626  // Otherwise if singleton row we explicitly eliminate this row and solve for corresponding X value
627  else {
628  int * Indices;
629  EPETRA_CHK_ERR(GetRow(i, NumEntries, Values, Indices)); // Get current row
630  if (NumEntries==1) {
631  double pivot = Values[0];
632  if (pivot==0.0) EPETRA_CHK_ERR(-1); // Encountered zero row, unable to continue
633  int indX = Indices[0];
634  for (j=0; j<NumVectors; j++)
635  (*tempExportX_)[j][indX] = (*FullRHS)[j][i]/pivot;
636  }
637  // Otherwise, this is a singleton column and we will scan for the pivot element needed
638  // for post-solve equations
639  else {
640  j = ColSingletonPivotLIDs_[ColSingletonCounter];
641  double pivot = Values[j];
642  if (pivot==0.0) EPETRA_CHK_ERR(-2); // Encountered zero column, unable to continue
643  ColSingletonPivots_[ColSingletonCounter] = pivot;
644  ColSingletonCounter++;
645  }
646  }
647  }
648 
649  assert(ColSingletonCounter==NumMyColSingletons_); // Sanity test
650 
651  // Update Reduced LHS (Puts any initial guess values into reduced system)
652 
653  ReducedLHS_->PutScalar(0.0); // zero out Reduced LHS
655  FullLHS->PutScalar(0.0); // zero out Full LHS since we will inject values as we get them
656 
657  // Construct Reduced RHS
658 
659  // Zero out temp space
660  tempX_->PutScalar(0.0);
661  tempB_->PutScalar(0.0);
662 
663  //Inject known X values into tempX for purpose of computing tempB = FullMatrix*tempX
664  // Also inject into full X since we already know the solution
665 
666  if (FullMatrix()->RowMatrixImporter()!=0) {
667  EPETRA_CHK_ERR(tempX_->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
668  EPETRA_CHK_ERR(FullLHS->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
669  }
670  else {
671  tempX_->Update(1.0, *tempExportX_, 0.0);
672  FullLHS->Update(1.0, *tempExportX_, 0.0);
673  }
674 
675  EPETRA_CHK_ERR(FullMatrix()->Multiply(false, *tempX_, *tempB_));
676 
677  EPETRA_CHK_ERR(tempB_->Update(1.0, *FullRHS, -1.0)); // tempB now has influence of already-known X values
678 
679  ReducedRHS_->PutScalar(0.0);
681  }
682  else {
683 
684  // There are no singletons, so don't bother building a reduced problem.
685  ReducedProblem_ = Teuchos::rcp( Problem, false );
686  ReducedMatrix_ = Teuchos::rcp( dynamic_cast<Epetra_CrsMatrix *>(Problem->GetMatrix()), false );
687  }
688 
689  return(0);
690 }
691 
693 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
694  if(Problem->GetMatrix()->RowMatrixRowMap().GlobalIndicesInt()) {
695  return TUpdateReducedProblem<int>(Problem);
696  }
697  else
698 #endif
699 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
700  if(Problem->GetMatrix()->RowMatrixRowMap().GlobalIndicesLongLong()) {
701  return TUpdateReducedProblem<long long>(Problem);
702  }
703  else
704 #endif
705  throw "LinearProblem_CrsSingletonFilter::UpdateReducedProblem: ERROR, GlobalIndices type unknown.";
706 }
707 //==============================================================================
708 template<typename int_type>
710  Epetra_Export * & RedistributeExporter,
711  Epetra_Map * & RedistributeMap) {
712 
713  int_type IndexBase = (int_type) SourceMap->IndexBase64();
714  if (IndexBase!=(int_type) TargetMap->IndexBase64()) EPETRA_CHK_ERR(-1);
715 
716  const Epetra_Comm & Comm = TargetMap->Comm();
717 
718  int TargetNumMyElements = TargetMap->NumMyElements();
719  int SourceNumMyElements = SourceMap->NumMyElements();
720 
721  // ContiguousTargetMap has same number of elements per PE as TargetMap, but uses contigious indexing
722  Epetra_Map ContiguousTargetMap((int_type) -1, TargetNumMyElements, IndexBase,Comm);
723 
724  // Same for ContiguousSourceMap
725  Epetra_Map ContiguousSourceMap((int_type) -1, SourceNumMyElements, IndexBase, Comm);
726 
727  assert(ContiguousSourceMap.NumGlobalElements64()==ContiguousTargetMap.NumGlobalElements64());
728 
729  // Now create a vector that contains the global indices of the Source Epetra_MultiVector
730  int_type* SourceMapMyGlobalElements = 0;
731  SourceMap->MyGlobalElementsPtr(SourceMapMyGlobalElements);
732  typename Epetra_GIDTypeVector<int_type>::impl SourceIndices(View, ContiguousSourceMap, SourceMapMyGlobalElements);
733 
734  // Create an exporter to send the SourceMap global IDs to the target distribution
735  Epetra_Export Exporter(ContiguousSourceMap, ContiguousTargetMap);
736 
737  // Create a vector to catch the global IDs in the target distribution
738  typename Epetra_GIDTypeVector<int_type>::impl TargetIndices(ContiguousTargetMap);
739  TargetIndices.Export(SourceIndices, Exporter, Insert);
740 
741  // Create a new map that describes how the Source MultiVector should be laid out so that it has
742  // the same number of elements on each processor as the TargetMap
743  RedistributeMap = new Epetra_Map((int_type) -1, TargetNumMyElements, TargetIndices.Values(), IndexBase, Comm);
744 
745  // This exporter will finally redistribute the Source MultiVector to the same layout as the TargetMap
746  RedistributeExporter = new Epetra_Export(*SourceMap, *RedistributeMap);
747  return(0);
748 }
749 
751  Epetra_Export * & RedistributeExporter,
752  Epetra_Map * & RedistributeMap) {
753 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
754  if(SourceMap->GlobalIndicesInt() && TargetMap->GlobalIndicesInt()) {
755  return TConstructRedistributeExporter<int>(SourceMap, TargetMap, RedistributeExporter, RedistributeMap);
756  }
757  else
758 #endif
759 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
760  if(SourceMap->GlobalIndicesLongLong() && TargetMap->GlobalIndicesLongLong()) {
761  return TConstructRedistributeExporter<long long>(SourceMap, TargetMap, RedistributeExporter, RedistributeMap);
762  }
763  else
764 #endif
765  throw "LinearProblem_CrsSingletonFilter::ConstructRedistributeExporter: ERROR, GlobalIndices type unknown.";
766 }
767 //==============================================================================
769 
770  if ( SingletonsDetected() ) {
771  int jj, k;
772 
773  Epetra_MultiVector * FullLHS = FullProblem()->GetLHS();
774  Epetra_MultiVector * FullRHS = FullProblem()->GetRHS();
775 
776  tempX_->PutScalar(0.0); tempExportX_->PutScalar(0.0);
777  // Inject values that the user computed for the reduced problem into the full solution vector
779 
780  FullLHS->Update(1.0, *tempX_, 1.0);
781 
782  // Next we will use our full solution vector which is populated with pre-filter solution
783  // values and reduced system solution values to compute the sum of the row contributions
784  // that must be subtracted to get the post-filter solution values
785 
786  EPETRA_CHK_ERR(FullMatrix()->Multiply(false, *FullLHS, *tempB_));
787 
788  // Finally we loop through the local rows that were associated with column singletons and compute the
789  // solution for these equations.
790 
791  int NumVectors = tempB_->NumVectors();
792  for (k=0; k<NumMyColSingletons_; k++) {
793  int i = ColSingletonRowLIDs_[k];
794  int j = ColSingletonColLIDs_[k];
795  double pivot = ColSingletonPivots_[k];
796  for (jj=0; jj<NumVectors; jj++)
797  (*tempExportX_)[jj][j]= ((*FullRHS)[jj][i] - (*tempB_)[jj][i])/pivot;
798  }
799 
800  // Finally, insert values from post-solve step and we are done!!!!
801 
802  if (FullMatrix()->RowMatrixImporter()!=0) {
803  EPETRA_CHK_ERR(tempX_->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));
804  }
805  else {
806  tempX_->Update(1.0, *tempExportX_, 0.0);
807  }
808 
809  FullLHS->Update(1.0, *tempX_, 1.0);
810  }
811 
812  return(0);
813 }
814 //==============================================================================
816 
818 
819  // Cast to CrsMatrix, if possible. Can save some work.
820  FullCrsMatrix_ = dynamic_cast<Epetra_CrsMatrix *>(FullMatrix());
821  FullMatrixIsCrsMatrix_ = (FullCrsMatrix_!=0); // Pointer is non-zero if cast worked
822  Indices_int_ = new int[MaxNumMyEntries_];
823 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
825  Indices_LL_ = new long long[MaxNumMyEntries_];
826 #endif
828 
829  return(0);
830 }
831 //==============================================================================
832 int LinearProblem_CrsSingletonFilter::GetRow(int Row, int & NumIndices, int * & Indices) {
833 
834  if (FullMatrixIsCrsMatrix_) { // View of current row
835  EPETRA_CHK_ERR(FullCrsMatrix()->Graph().ExtractMyRowView(Row, NumIndices, Indices));
836  }
837  else { // Copy of current row (we must get the values, but we ignore them)
838  EPETRA_CHK_ERR(FullMatrix()->ExtractMyRowCopy(Row, MaxNumMyEntries_, NumIndices,
840  Indices = Indices_int_;
841  }
842  return(0);
843 }
844 //==============================================================================
845 int LinearProblem_CrsSingletonFilter::GetRow(int Row, int & NumIndices,
846  double * & Values, int * & Indices) {
847 
848  if (FullMatrixIsCrsMatrix_) { // View of current row
849  EPETRA_CHK_ERR(FullCrsMatrix_->ExtractMyRowView(Row, NumIndices, Values, Indices));
850  }
851  else { // Copy of current row (we must get the values, but we ignore them)
852  EPETRA_CHK_ERR(FullMatrix()->ExtractMyRowCopy(Row, MaxNumMyEntries_, NumIndices,
854  Values = Values_.Values();
855  Indices = Indices_int_;
856  }
857  return(0);
858 }
859 //==============================================================================
860 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
861 int LinearProblem_CrsSingletonFilter::GetRowGCIDs(int Row, int & NumIndices,
862  double * & Values, int * & GlobalIndices) {
863 
864  EPETRA_CHK_ERR(FullMatrix()->ExtractMyRowCopy(Row, MaxNumMyEntries_, NumIndices,
866  for (int j=0; j<NumIndices; j++) Indices_int_[j] = FullMatrixColMap().GID(Indices_int_[j]);
867  Values = Values_.Values();
868  GlobalIndices = Indices_int_;
869  return(0);
870 }
871 #endif
872 
873 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
874 int LinearProblem_CrsSingletonFilter::GetRowGCIDs(int Row, int & NumIndices,
875  double * & Values, long long * & GlobalIndices) {
876  EPETRA_CHK_ERR(FullMatrix()->ExtractMyRowCopy(Row, MaxNumMyEntries_, NumIndices,
878  for (int j=0; j<NumIndices; j++) Indices_LL_[j] = FullMatrixColMap().GID64(Indices_int_[j]);
879  Values = Values_.Values();
880  GlobalIndices = Indices_LL_;
881  return(0);
882 }
883 #endif
884 //==============================================================================
886  const Epetra_MapColoring & rowMapColors,
887  const Epetra_IntVector & ColProfiles,
888  const Epetra_IntVector & NewColProfiles,
889  const Epetra_IntVector & ColHasRowWithSingleton) {
890 
891  int j;
892 
893  if (NumMyColSingletons_==0) return(0); // Nothing to do
894 
895  Epetra_MapColoring & colMapColors = *ColMapColors_;
896 
897  int NumMyCols = FullMatrix()->NumMyCols();
898 
899  // We will need these arrays for the post-solve phase
904 
905  // Register singleton columns (that were not already counted as singleton rows)
906  // Check to see if any columns disappeared because all associated rows were eliminated
907  int NumMyColSingletonstmp = 0;
908  for (j=0; j<NumMyCols; j++) {
909  int i = RowIDs[j];
910  if ( ColProfiles[j]==1 && rowMapColors[i]!=1 ) {
911  ColSingletonRowLIDs_[NumMyColSingletonstmp] = i;
912  ColSingletonColLIDs_[NumMyColSingletonstmp] = j;
913  NumMyColSingletonstmp++;
914  }
915  // Also check for columns that were eliminated implicitly by
916  // having all associated row eliminated
917  else if (NewColProfiles[j]==0 && ColHasRowWithSingleton[j]!=1 && rowMapColors[i]==0) {
918  colMapColors[j] = 1;
919  }
920  }
921 
922  assert(NumMyColSingletonstmp==NumMyColSingletons_); //Sanity check
923  Epetra_Util sorter;
925 
926  return(0);
927 }
928 
929 } //namespace EpetraExt
930 
int GetRowGCIDs(int Row, int &NumIndices, double *&Values, int *&GlobalIndices)
Epetra_Map * ReducedMatrixDomainMap() const
Returns pointer to Epetra_Map describing the domain map for the reduced system.
int Size(int Length_in)
Epetra_Map * ReducedMatrixRangeMap() const
Returns pointer to Epetra_Map describing the range map for the reduced system.
double RatioOfDimensions() const
Returns ratio of reduced system to full system dimensions, returns -1.0 if reduced problem not constr...
virtual const Epetra_Map & RowMatrixRowMap() const =0
long long NumGlobalRows64() const
Epetra_MultiVector * GetLHS() const
Epetra_MultiVector * GetRHS() const
bool GlobalIndicesLongLong() const
NewTypeRef operator()(OriginalTypeRef orig)
Analysis of transform operation on original object and construction of new object.
bool fwd()
Forward transfer of data from orig object input in the operator() method call to the new object creat...
bool SameAs(const Epetra_BlockMap &Map) const
long long IndexBase64() const
long long NumGlobalElements64() const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
NewTypeRef construct()
Construction of new object as a result of the transform.
virtual ~LinearProblem_CrsSingletonFilter()
Epetra_CrsSingletonFilter Destructor.
#define EPETRA_CHK_ERR(a)
bool rvs()
Reverse transfer of data from new object created in the operator() method call to the orig object inp...
int PutValue(int Value)
bool GlobalIndicesInt() const
virtual int MyPID() const =0
bool SingletonsDetected() const
Returns true if singletons were detected in this matrix (must be called after Analyze() to be effecti...
int Analyze(Epetra_RowMatrix *FullMatrix)
Analyze the input matrix, removing row/column pairs that have singletons.
virtual int NumMyCols() const =0
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
long long NumGlobalNonzeros64() const
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
int NumMyElements() const
virtual int MaxNumEntries() const =0
long long GID64(int LID) const
virtual const Epetra_Comm & Comm() const =0
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual const Epetra_Import * RowMatrixImporter() const =0
Epetra_Map * ReducedMatrixRowMap() const
Returns pointer to Epetra_Map describing the reduced system row distribution.
Epetra_CrsMatrix * ReducedMatrix() const
Returns pointer to Epetra_CrsMatrix from full problem.
int ComputeFullSolution()
Compute a solution for the full problem using the solution of the reduced problem, put in LHS of FullProblem().
int GID(int LID) const
virtual const Epetra_BlockMap & Map() const =0
LinearProblem_CrsSingletonFilter(bool verbose=false)
Epetra_CrsSingletonFilter default constructor.
virtual int NumMyRows() const =0
Epetra_LinearProblem * ReducedProblem() const
Returns pointer to the derived reduced Epetra_LinearProblem.
int MyGlobalElementsPtr(int *&MyGlobalElementList) const
bool analyze(OriginalTypeRef orig)
Initial analysis phase of transform.
const Epetra_Comm & Comm() const
virtual long long NumGlobalNonzeros64() const =0
double RatioOfNonzeros() const
Returns ratio of reduced system to full system nonzero count, returns -1.0 if reduced problem not con...
int ConstructRedistributeExporter(Epetra_Map *SourceMap, Epetra_Map *TargetMap, Epetra_Export *&RedistributeExporter, Epetra_Map *&RedistributeMap)
Epetra_Map * GenerateMap(int Color) const
Epetra_RowMatrix * GetMatrix() const
int CreatePostSolveArrays(const Epetra_IntVector &RowIDs, const Epetra_MapColoring &RowMapColors, const Epetra_IntVector &ColProfiles, const Epetra_IntVector &NewColProfiles, const Epetra_IntVector &ColHasRowWithSingleton)
void Sort(bool SortAscending, int NumKeys, T *Keys, int NumDoubleCompanions, double **DoubleCompanions, int NumIntCompanions, int **IntCompanions, int NumLongLongCompanions, long long **LongLongCompanions)
Epetra_RowMatrix * FullMatrix() const
Returns pointer to Epetra_CrsMatrix from full problem.
double * Values() const
int TConstructRedistributeExporter(Epetra_Map *SourceMap, Epetra_Map *TargetMap, Epetra_Export *&RedistributeExporter, Epetra_Map *&RedistributeMap)
int UpdateReducedProblem(Epetra_LinearProblem *Problem)
Update a reduced linear problem using new values.
Epetra_Map * ReducedMatrixColMap() const
Returns pointer to Epetra_Map describing the reduced system column distribution.
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int ConstructReducedProblem(Epetra_LinearProblem *Problem)
Return a reduced linear problem based on results of Analyze().
Epetra_LinearProblem * FullProblem() const
Returns pointer to the original unreduced Epetra_LinearProblem.
virtual long long NumGlobalRows64() const =0
int NumSingletons() const
Return total number of singletons detected, returns -1 if Analysis not performed yet.