00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include "NOX_Epetra_Scaling.H"
00040
00041 #include "Epetra_Vector.h"
00042 #include "Epetra_Operator.h"
00043 #include "Epetra_RowMatrix.h"
00044 #include "Epetra_LinearProblem.h"
00045
00046 #include "NOX_Utils.H"
00047
00048 NOX::Epetra::Scaling::Scaling()
00049 {
00050
00051 }
00052
00053 NOX::Epetra::Scaling::~Scaling()
00054 {
00055
00056 }
00057
00058 void NOX::Epetra::Scaling::addUserScaling(ScaleType type, const Teuchos::RCP<Epetra_Vector>& D)
00059 {
00060 if ( Teuchos::is_null(tmpVectorPtr) )
00061 tmpVectorPtr = Teuchos::rcp(new Epetra_Vector(*D));
00062
00063 scaleType.push_back(type);
00064 sourceType.push_back(UserDefined);
00065 scaleVector.push_back(D);
00066 }
00067
00068 void NOX::Epetra::Scaling::addRowSumScaling(ScaleType type, const Teuchos::RCP<Epetra_Vector>& D)
00069 {
00070 if ( Teuchos::is_null(tmpVectorPtr) )
00071 tmpVectorPtr = Teuchos::rcp(new Epetra_Vector(*D));
00072
00073 scaleType.push_back(type);
00074 sourceType.push_back(RowSum);
00075 scaleVector.push_back(D);
00076 }
00077
00078 void NOX::Epetra::Scaling::addColSumScaling(ScaleType type, const Teuchos::RCP<Epetra_Vector>& D)
00079 {
00080 if ( Teuchos::is_null(tmpVectorPtr) )
00081 tmpVectorPtr = Teuchos::rcp(new Epetra_Vector(*D));
00082
00083 scaleType.push_back(type);
00084 sourceType.push_back(ColSum);
00085 scaleVector.push_back(D);
00086 }
00087
00088 void NOX::Epetra::Scaling::computeScaling(const Epetra_LinearProblem& problem)
00089 {
00090
00091 Epetra_Vector* diagonal = 0;
00092 for (unsigned int i = 0; i < scaleVector.size(); i ++) {
00093
00094 if (sourceType[i] == RowSum) {
00095
00096 diagonal = scaleVector[i].get();
00097
00098
00099
00100 const Epetra_RowMatrix* test = 0;
00101 test = dynamic_cast<const Epetra_RowMatrix*>(problem.GetOperator());
00102 if (test == 0) {
00103 cout << "ERROR: NOX::Epetra::Scaling::scaleLinearSystem() - "
00104 << "For \"Row Sum\" scaling, the Matrix must be an "
00105 << "Epetra_RowMatrix derived object!" << endl;
00106 throw "NOX Error";
00107 }
00108
00109 test->InvRowSums(*diagonal);
00110 diagonal->Reciprocal(*diagonal);
00111
00112 }
00113
00114 else if (sourceType[i] == ColSum) {
00115
00116 diagonal = scaleVector[i].get();
00117
00118
00119
00120 const Epetra_RowMatrix* test = 0;
00121 test = dynamic_cast<const Epetra_RowMatrix*>(problem.GetOperator());
00122 if (test == 0) {
00123 cout << "ERROR: NOX::Epetra::Scaling::scaleLinearSystem() - "
00124 << "For \"Column Sum\" scaling, the Matrix must be an "
00125 << "Epetra_RowMatrix derived object!" << endl;
00126 throw "NOX Error";
00127 }
00128
00129 test->InvColSums(*diagonal);
00130 diagonal->Reciprocal(*diagonal);
00131
00132 }
00133
00134 }
00135
00136 }
00137
00138 void NOX::Epetra::Scaling::scaleLinearSystem(Epetra_LinearProblem& problem)
00139 {
00140 Epetra_Vector* diagonal = 0;
00141 for (unsigned int i = 0; i < scaleVector.size(); i ++) {
00142
00143 diagonal = scaleVector[i].get();
00144
00145 if (scaleType[i] == Left) {
00146
00147 tmpVectorPtr->Reciprocal(*diagonal);
00148 problem.LeftScale(*tmpVectorPtr);
00149
00150 }
00151 else if (scaleType[i] == Right) {
00152
00153 tmpVectorPtr->Reciprocal(*diagonal);
00154 problem.RightScale(*tmpVectorPtr);
00155 }
00156
00157 }
00158 }
00159
00160 void NOX::Epetra::Scaling::unscaleLinearSystem(Epetra_LinearProblem& problem)
00161 {
00162 Epetra_Vector* diagonal = 0;
00163 for (unsigned int i = 0; i < scaleVector.size(); i ++) {
00164
00165 diagonal = scaleVector[i].get();
00166
00167 if (scaleType[i] == Left) {
00168 problem.LeftScale(*diagonal);
00169 }
00170 else if (scaleType[i] == Right) {
00171 problem.RightScale(*diagonal);
00172
00173 }
00174 }
00175 }
00176
00177 void NOX::Epetra::Scaling::applyRightScaling(const Epetra_Vector& input,
00178 Epetra_Vector& result)
00179 {
00180 if (scaleVector.size() == 0) {
00181 result = input;
00182 }
00183 else {
00184 Epetra_Vector* diagonal = 0;
00185 for (unsigned int i = 0; i < scaleVector.size(); i ++) {
00186
00187 if (scaleType[i] == Right) {
00188 diagonal = scaleVector[i].get();
00189
00190 tmpVectorPtr->Reciprocal(*diagonal);
00191
00192 result.Multiply(1.0, input, *tmpVectorPtr, 0.0);
00193 }
00194 }
00195 }
00196 }
00197
00198 void NOX::Epetra::Scaling::applyLeftScaling(const Epetra_Vector& input,
00199 Epetra_Vector& result)
00200 {
00201 if (scaleVector.size() == 0) {
00202 result = input;
00203 }
00204 else {
00205 Epetra_Vector* diagonal = 0;
00206 for (unsigned int i = 0; i < scaleVector.size(); i ++) {
00207
00208 if (scaleType[i] == Left) {
00209 diagonal = scaleVector[i].get();
00210
00211 tmpVectorPtr->Reciprocal(*diagonal);
00212
00213 result.Multiply(1.0, input, *tmpVectorPtr, 0.0);
00214 }
00215 }
00216 }
00217 }
00218
00219 void NOX::Epetra::Scaling::print(ostream& os)
00220 {
00221
00222 os << "\n LINEAR SOLVER SCALING:" << endl;
00223
00224 for (unsigned int i = 0; i < scaleVector.size(); i ++) {
00225
00226 string source = " ";
00227 if (sourceType[i] == UserDefined)
00228 source = "User Defined Vector";
00229 else if (sourceType[i] == RowSum)
00230 source = "Row Sum";
00231 else if (sourceType[i] == ColSum)
00232 source = "Col Sum";
00233
00234 if (scaleType[i] == Left) {
00235 os << " " << (i+1) << ". Left Scaled with " << source << endl;
00236
00237 }
00238 else if (scaleType[i] == Right)
00239 os << " " << (i+1) << ". Right Scaled with " << source << endl;
00240 }
00241
00242 return;
00243 }
00244
00245 ostream& operator<<(ostream& os, NOX::Epetra::Scaling& scalingObject)
00246 {
00247 scalingObject.print(os);
00248 return os;
00249 }