Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
55#endif
56
57#include "Epetra_Comm.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
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
126
131 RatioOfDimensions_ = -1.0;
132 RatioOfNonzeros_ = -1.0;
133
134 HaveReducedProblem_ = false;
136 AnalysisDone_ = false;
138
139 tempExportX_ = 0;
140 tempX_ = 0;
141 tempB_ = 0;
142
143 Indices_ = 0;
144
145 RowMapColors_ = 0;
146 ColMapColors_ = 0;
147
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
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
440
441 // Finally construct Reduced Linear Problem
442
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
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
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//==============================================================================
703int 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,
711 Indices = Indices_;
712 }
713 return(0);
714}
715//==============================================================================
716int 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,
725 Values = Values_.Values();
726 Indices = Indices_;
727 }
728 return(0);
729}
730//==============================================================================
731int Epetra_CrsSingletonFilter::GetRowGCIDs(int Row, int & NumIndices,
732 double * & Values, int * & GlobalIndices) {
733
734 EPETRA_CHK_ERR(FullMatrix()->ExtractMyRowCopy(Row, MaxNumMyEntries_, NumIndices,
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;
782
783 return(0);
784}
785
786#endif // EPETRA_NO_32BIT_GLOBAL_INDICES
@ Insert
@ AbsMax
#define EPETRA_CHK_ERR(a)
@ View
@ Copy
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
long long * MyGlobalElements64() const
int IndexBase() const
Index base for this map.
long long GID64(int LID) const
long long NumGlobalElements64() const
bool SameAs(const Epetra_BlockMap &Map) const
Returns true if this and Map are identical maps.
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int NumMyElements() const
Number of elements on the calling processor.
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
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.
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.
long long NumGlobalRows64() const
long long NumGlobalNonzeros64() const
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 ComputeFullSolution()
Compute a solution for the full problem using the solution of the reduced problem,...
int Analyze(Epetra_RowMatrix *FullMatrix)
Analyze the input matrix, removing row/column pairs that have singletons.
int GetRowGCIDs(int Row, int &NumIndices, double *&Values, int *&GlobalIndices)
const Epetra_Map & FullMatrixRowMap() const
int ConstructRedistributeExporter(Epetra_Map *SourceMap, Epetra_Map *TargetMap, Epetra_Export *&RedistributeExporter, Epetra_Map *&RedistributeMap)
int ConstructReducedProblem(Epetra_LinearProblem *Problem)
Return a reduced linear problem based on results of Analyze().
int GetRow(int Row, int &NumIndices, int *&Indices)
Epetra_CrsSingletonFilter()
Epetra_CrsSingletonFilter default constructor.
Epetra_Map * ReducedMatrixDomainMap() const
Returns pointer to Epetra_Map describing the domain map for the reduced system.
Epetra_Map * ReducedMatrixColMap() const
Returns pointer to Epetra_Map describing the reduced system column distribution.
const Epetra_Map & FullMatrixDomainMap() const
const Epetra_Map & FullMatrixColMap() const
Epetra_LinearProblem * FullProblem_
Epetra_CrsMatrix * ReducedMatrix() const
Returns pointer to Epetra_CrsMatrix from full problem.
Epetra_LinearProblem * FullProblem() const
Returns pointer to the original unreduced Epetra_LinearProblem.
virtual ~Epetra_CrsSingletonFilter()
Epetra_CrsSingletonFilter Destructor.
const Epetra_Map & FullMatrixRangeMap() const
int UpdateReducedProblem(Epetra_LinearProblem *Problem)
Update a reduced linear problem using new values.
int CreatePostSolveArrays(const Epetra_IntVector &RowIDs, const Epetra_MapColoring &RowMapColors, const Epetra_IntVector &ColProfiles, const Epetra_IntVector &NewColProfiles, const Epetra_IntVector &ColHasRowWithSingleton)
Epetra_RowMatrix * FullMatrix() const
Returns pointer to Epetra_CrsMatrix from full problem.
Epetra_Map * ReducedMatrixRowMap() const
Returns pointer to Epetra_Map describing the reduced system row distribution.
Epetra_LinearProblem * ReducedProblem_
Epetra_SerialDenseVector Values_
Epetra_Map * ReducedMatrixRangeMap() const
Returns pointer to Epetra_Map describing the range map for the reduced system.
Epetra_CrsMatrix * FullCrsMatrix() const
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
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.
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_Export: This class builds an export object for efficient exporting of off-processor elements.
Definition: Epetra_Export.h:62
Epetra_Import: This class builds an import object for efficient importing of off-processor elements.
Definition: Epetra_Import.h:63
Epetra_IntVector: A class for constructing and using dense integer vectors on a parallel computer.
int MaxValue()
Find maximum value.
int PutValue(int Value)
Set all elements of the vector to Value.
int * Values() const
Returns a pointer to an array containing the values of this vector.
Epetra_LinearProblem: The Epetra Linear Problem Class.
Epetra_MultiVector * GetRHS() const
Get a pointer to the right-hand-side B.
Epetra_RowMatrix * GetMatrix() const
Get a pointer to the matrix A.
Epetra_MultiVector * GetLHS() const
Get a pointer to the left-hand-side X.
Epetra_LongLongVector: A class for constructing and using dense integer vectors on a parallel compute...
long long * Values() const
Returns a pointer to an array containing the values of this vector.
Epetra_MapColoring: A class for coloring Epetra_Map and Epetra_BlockMap objects.
Epetra_Map * GenerateMap(int Color) const
Generates an Epetra_Map of the GIDs associated with the specified color.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int NumVectors() const
Returns the number of vectors in the multi-vector.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices.
virtual int NumMyRows() const =0
Returns the number of matrix rows owned by the calling processor.
virtual int NumMyCols() const =0
Returns the number of matrix columns owned by the calling processor.
virtual long long NumGlobalNonzeros64() const =0
virtual long long NumGlobalRows64() const =0
virtual const Epetra_Map & RowMatrixRowMap() const =0
Returns the Epetra_Map object associated with the rows of this matrix.
virtual int MaxNumEntries() const =0
Returns the maximum of NumMyRowEntries() over all rows.
virtual const Epetra_Import * RowMatrixImporter() const =0
Returns the Epetra_Import object that contains the import operations for distributed operations.
int Size(int Length_in)
Set length of a Epetra_SerialDenseVector object; init values to zero.
double * Values() const
Returns pointer to the values in vector.
Epetra_Util: The Epetra Util Wrapper Class.
Definition: Epetra_Util.h:79
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)