EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
EpetraExt_Scale_LinearProblem.cpp
Go to the documentation of this file.
1//@HEADER
2// ***********************************************************************
3//
4// EpetraExt: Epetra Extended - Linear Algebra Services Package
5// Copyright (2011) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41
43
45#include <Epetra_CrsMatrix.h>
46#include <Epetra_Vector.h>
47
48namespace EpetraExt {
49
52{
53 int lsize = (int) lScaleVecs_.size();
54 for( int i = 0; i < lsize; ++i )
55 delete lScaleVecs_[i];
56 int rsize = (int) rScaleVecs_.size();
57 for( int i = 0; i < rsize; ++i )
58 delete rScaleVecs_[i];
59}
60
61bool
63fwd()
64{
65 Epetra_CrsMatrix & Matrix = *(dynamic_cast<Epetra_CrsMatrix*>(origObj_->GetMatrix()));
66
67 const Epetra_BlockMap & RHSMap = origObj_->GetRHS()->Map();
68 const Epetra_BlockMap & LHSMap = origObj_->GetLHS()->Map();
69
70 if( iters_ > 0 )
71 {
72 if( lScale_ != None && !lScaleVecs_.size() )
73 {
74 lScaleVecs_.resize(iters_);
75 for( int i = 0; i < iters_; ++i )
76 lScaleVecs_[i] = new Epetra_Vector( RHSMap );
77 }
78 if( rScale_ != None && !rScaleVecs_.size() )
79 {
80 rScaleVecs_.resize(iters_);
81 for( int i = 0; i < iters_; ++i )
82 rScaleVecs_[i] = new Epetra_Vector( LHSMap );
83 }
84
85 for( int i = 0; i < iters_; ++i )
86 {
87 if( lScale_ != None )
88 {
89 switch( lScale_ )
90 {
91 case Max: Matrix.InvRowMaxs( *(lScaleVecs_[i]) );
92 break;
93 case Sum: Matrix.InvRowSums( *(lScaleVecs_[i]) );
94 break;
95 case Diag: Matrix.ExtractDiagonalCopy( *(lScaleVecs_[i]) );
96 lScaleVecs_[i]->Reciprocal( *(lScaleVecs_[i]) );
97 break;
98 default: break;
99 };
100 if( expFac_ != 1.0 )
101 {
102 int numVals = RHSMap.NumMyPoints();
103 for( int j = 0; j < numVals; ++j ) (*(lScaleVecs_[i]))[j] = pow( (*(lScaleVecs_[i]))[j], expFac_ );
104 }
105 newObj_->LeftScale( *lScaleVecs_[i] );
106 }
107 if( rScale_ != None )
108 {
109 switch( rScale_ )
110 {
111 case Max: Matrix.InvColMaxs( *(rScaleVecs_[i]) );
112 break;
113 case Sum: Matrix.InvColSums( *(rScaleVecs_[i]) );
114 break;
115 case Diag: Matrix.ExtractDiagonalCopy( *(rScaleVecs_[i]) );
116 rScaleVecs_[i]->Reciprocal( *(rScaleVecs_[i]) );
117 break;
118 default: break;
119 };
120 if( expFac_ != 1.0 )
121 {
122 int numVals = LHSMap.NumMyPoints();
123 for( int j = 0; j < numVals; ++j ) (*(rScaleVecs_[i]))[j] = pow( (*(rScaleVecs_[i]))[j], expFac_ );
124 }
125 newObj_->RightScale( *rScaleVecs_[i] );
126 }
127 }
128 }
129
130 scaled_ = true;
131
132 return true;
133}
134
135bool
137rvs()
138{
139 if( !scaled_ ) std::cout << "EpetraExt::LinearProblem_Scale::rvs() : Problem Not Scaled!\n";
140
141 if( iters_ > 0 )
142 {
143 for( int i = 0; i < iters_; ++i )
144 {
145 int loc = iters_-i-1;
146 if( rScale_ != None )
147 {
148 rScaleVecs_[loc]->Reciprocal( *(rScaleVecs_[loc]) );
149 newObj_->RightScale( *(rScaleVecs_[loc]) );
150 }
151 if( lScale_ != None )
152 {
153 lScaleVecs_[loc]->Reciprocal( *(lScaleVecs_[loc]) );
154 newObj_->LeftScale( *(lScaleVecs_[loc]) );
155 }
156 }
157 }
158 return true;
159}
160
161} //namespace EpetraExt
162
std::vector< Epetra_Vector * > rScaleVecs_
std::vector< Epetra_Vector * > lScaleVecs_
int NumMyPoints() const
int InvRowSums(Epetra_Vector &x) const
int InvColMaxs(Epetra_Vector &x) const
int InvRowMaxs(Epetra_Vector &x) const
int InvColSums(Epetra_Vector &x) const
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.