EpetraExt Development
Loading...
Searching...
No Matches
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
56namespace 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
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//==============================================================================
107bool
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//==============================================================================
159construct()
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{
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
229
234 RatioOfDimensions_ = -1.0;
235 RatioOfNonzeros_ = -1.0;
236
237 HaveReducedProblem_ = false;
239 AnalysisDone_ = false;
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
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.
285 EPETRA_CHK_ERR(InitFullMatrixAccess());
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//==============================================================================
392template<typename int_type>
393int LinearProblem_CrsSingletonFilter::TConstructReducedProblem(Epetra_LinearProblem * Problem) {
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
459 ReducedMatrix_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *ReducedMatrixRowMap(), *ReducedMatrixColMap(), 0) );
460
461 // Create storage for temporary X values due to explicit elimination of rows
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
513 EPETRA_CHK_ERR(ReducedMatrix()->FillComplete(*ReducedMatrixDomainMap(), *ReducedMatrixRangeMap()));
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
519 EPETRA_CHK_ERR(ReducedLHS_->Import(*FullLHS, *Full2ReducedLHSImporter_, Insert));
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
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
578 return TConstructReducedProblem<long long>(Problem);
579 }
580 else
581#endif
582 throw "LinearProblem_CrsSingletonFilter::ConstructReducedProblem: ERROR, GlobalIndices type unknown.";
583}
584
585//==============================================================================
586template<typename int_type>
587int LinearProblem_CrsSingletonFilter::TUpdateReducedProblem(Epetra_LinearProblem * Problem) {
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();
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
654 EPETRA_CHK_ERR(ReducedLHS_->Import(*FullLHS, *Full2ReducedLHSImporter_, Insert));
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
680 EPETRA_CHK_ERR(ReducedRHS_->Import(*tempB_, *Full2ReducedRHSImporter_, Insert));
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
701 return TUpdateReducedProblem<long long>(Problem);
702 }
703 else
704#endif
705 throw "LinearProblem_CrsSingletonFilter::UpdateReducedProblem: ERROR, GlobalIndices type unknown.";
706}
707//==============================================================================
708template<typename int_type>
709int LinearProblem_CrsSingletonFilter::TConstructRedistributeExporter(Epetra_Map * SourceMap, Epetra_Map * TargetMap,
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
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
823#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
825 Indices_LL_ = new long long[MaxNumMyEntries_];
826#endif
828
829 return(0);
830}
831//==============================================================================
832int 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//==============================================================================
845int 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
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
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
Insert
Add
bool analyze(OriginalTypeRef orig)
Initial analysis phase of transform.
Epetra_RowMatrix * FullMatrix() const
Returns pointer to Epetra_CrsMatrix from full problem.
virtual ~LinearProblem_CrsSingletonFilter()
Epetra_CrsSingletonFilter Destructor.
Epetra_LinearProblem * ReducedProblem() const
Returns pointer to the derived reduced Epetra_LinearProblem.
int GetRowGCIDs(int Row, int &NumIndices, double *&Values, int *&GlobalIndices)
bool rvs()
Reverse transfer of data from new object created in the operator() method call to the orig object inp...
Epetra_LinearProblem * FullProblem() const
Returns pointer to the original unreduced Epetra_LinearProblem.
int ConstructReducedProblem(Epetra_LinearProblem *Problem)
Return a reduced linear problem based on results of Analyze().
int CreatePostSolveArrays(const Epetra_IntVector &RowIDs, const Epetra_MapColoring &RowMapColors, const Epetra_IntVector &ColProfiles, const Epetra_IntVector &NewColProfiles, const Epetra_IntVector &ColHasRowWithSingleton)
NewTypeRef construct()
Construction of new object as a result of the transform.
Epetra_Map * ReducedMatrixRowMap() const
Returns pointer to Epetra_Map describing the reduced system row distribution.
NewTypeRef operator()(OriginalTypeRef orig)
Analysis of transform operation on original object and construction of new object.
double RatioOfNonzeros() const
Returns ratio of reduced system to full system nonzero count, returns -1.0 if reduced problem not con...
Epetra_Map * ReducedMatrixRangeMap() const
Returns pointer to Epetra_Map describing the range map for the reduced system.
int Analyze(Epetra_RowMatrix *FullMatrix)
Analyze the input matrix, removing row/column pairs that have singletons.
double RatioOfDimensions() const
Returns ratio of reduced system to full system dimensions, returns -1.0 if reduced problem not constr...
bool fwd()
Forward transfer of data from orig object input in the operator() method call to the new object creat...
Epetra_CrsMatrix * ReducedMatrix() const
Returns pointer to Epetra_CrsMatrix from full problem.
LinearProblem_CrsSingletonFilter(bool verbose=false)
Epetra_CrsSingletonFilter default constructor.
int ComputeFullSolution()
Compute a solution for the full problem using the solution of the reduced problem,...
int NumSingletons() const
Return total number of singletons detected, returns -1 if Analysis not performed yet.
Epetra_Map * ReducedMatrixDomainMap() const
Returns pointer to Epetra_Map describing the domain map for the reduced system.
int UpdateReducedProblem(Epetra_LinearProblem *Problem)
Update a reduced linear problem using new values.
bool SingletonsDetected() const
Returns true if singletons were detected in this matrix (must be called after Analyze() to be effecti...
int ConstructRedistributeExporter(Epetra_Map *SourceMap, Epetra_Map *TargetMap, Epetra_Export *&RedistributeExporter, Epetra_Map *&RedistributeMap)
Epetra_Map * ReducedMatrixColMap() const
Returns pointer to Epetra_Map describing the reduced system column distribution.
int GID(int LID) const
bool SameAs(const Epetra_BlockMap &Map) const
bool GlobalIndicesInt() const
const Epetra_Comm & Comm() const
int NumMyElements() const
bool MyGID(int GID_in) const
bool GlobalIndicesLongLong() const
virtual int MyPID() const=0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
const Epetra_BlockMap & Map() const
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int PutValue(int Value)
Epetra_MultiVector * GetRHS() const
Epetra_RowMatrix * GetMatrix() const
Epetra_MultiVector * GetLHS() const
Epetra_Map * GenerateMap(int Color) const
int NumVectors() const
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int PutScalar(double ScalarConstant)
virtual const Epetra_Comm & Comm() const=0
virtual int NumMyRows() const=0
virtual int NumMyCols() const=0
virtual const Epetra_Map & RowMatrixRowMap() const=0
virtual int MaxNumEntries() const=0
virtual const Epetra_Import * RowMatrixImporter() const=0
int Size(int Length_in)
double * Values() 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)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.