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 "Epetra_Comm.h"
00040 #include "Epetra_Map.h"
00041 #include "Epetra_Vector.h"
00042 #include "Epetra_RowMatrix.h"
00043 #include "NOX_Abstract_Group.H"
00044 #include "NOX_Epetra_Vector.H"
00045 #include "NOX_Epetra_VectorSpace.H"
00046 #include "NOX_Epetra_Interface_Required.H"
00047 #include "NOX_Utils.H"
00048
00049 #include "NOX_Epetra_MatrixFree.H"
00050
00051 using namespace NOX;
00052 using namespace NOX::Epetra;
00053
00054 MatrixFree::MatrixFree(Teuchos::ParameterList& printParams,
00055 const Teuchos::RCP<NOX::Epetra::Interface::Required>& i,
00056 const NOX::Epetra::Vector& x, bool p) :
00057 label("NOX::Matrix-Free"),
00058 interface(i),
00059 currentX(x),
00060 perturbX(x),
00061 fo(x),
00062 fp(x),
00063 diffType(Forward),
00064 lambda(1.0e-6),
00065 eta(0.0),
00066 userEta(1.0e-6),
00067 computeEta(true),
00068 useGroupForComputeF(false),
00069 useNewPerturbation(p),
00070 utils(printParams)
00071 {
00072
00073 perturbX.init(0.0);
00074 fo.init(0.0);
00075 fp.init(0.0);
00076
00077
00078
00079
00080 const Epetra_Map* testMap = 0;
00081 testMap = dynamic_cast<const Epetra_Map*>(¤tX.getEpetraVector().Map());
00082 if (testMap != 0) {
00083 epetraMap = Teuchos::rcp(new Epetra_Map(*testMap));
00084 }
00085 else {
00086 int size = currentX.getEpetraVector().Map().NumGlobalPoints();
00087 int mySize = currentX.getEpetraVector().Map().NumMyPoints();
00088 int indexBase = currentX.getEpetraVector().Map().IndexBase();
00089 const Epetra_Comm& comm = currentX.getEpetraVector().Map().Comm();
00090 epetraMap = Teuchos::rcp(new Epetra_Map(size, mySize, indexBase, comm));
00091 }
00092
00093 }
00094
00095 MatrixFree::~MatrixFree()
00096 {
00097
00098 }
00099
00100 int MatrixFree::SetUseTranspose(bool UseTranspose)
00101 {
00102 if (UseTranspose == true) {
00103 utils.out() << "ERROR: NOX::Epetra::MatrixFree::SetUseTranspose() - Transpose is "
00104 << "unavailable in Matrix-Free mode!" << endl;
00105 throw "NOX Error";
00106 }
00107 return (-1);
00108 }
00109
00110 int MatrixFree::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00111 {
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 Teuchos::RCP<Epetra_Vector> wrappedX =
00129 Teuchos::rcp(new Epetra_Vector(View, X, 0));
00130 Teuchos::RCP<Epetra_Vector> wrappedY =
00131 Teuchos::rcp(new Epetra_Vector(View, Y, 0));
00132 NOX::Epetra::Vector nevX(wrappedX, NOX::Epetra::Vector::CreateView);
00133 NOX::Epetra::Vector nevY(wrappedY, NOX::Epetra::Vector::CreateView);
00134
00135
00136
00137
00138 double solutionNorm = 1.0;
00139 double vectorNorm = 1.0;
00140
00141 solutionNorm = currentX.norm();
00142 vectorNorm = currentX.getVectorSpace()->norm(*wrappedX);
00143
00144
00145 if (vectorNorm == 0.0) {
00146
00147 vectorNorm = 1.0;
00148 }
00149
00150
00151 if ( diffType == Centered )
00152 if ( Teuchos::is_null(fmPtr) )
00153 fmPtr = Teuchos::rcp(new NOX::Epetra::Vector(fo));
00154
00155 double scaleFactor = 1.0;
00156 if ( diffType == Backward )
00157 scaleFactor = -1.0;
00158
00159 if (computeEta) {
00160 if (useNewPerturbation) {
00161 double dotprod = currentX.getVectorSpace()->
00162 innerProduct(currentX.getEpetraVector(), *wrappedX);
00163 if (dotprod==0.0)
00164 dotprod = 1.0e-12;
00165 eta = lambda*(1.0e-12/lambda + fabs(dotprod)/(vectorNorm * vectorNorm))
00166 * dotprod/fabs(dotprod);
00167 }
00168 else
00169 eta = lambda*(lambda + solutionNorm/vectorNorm);
00170 }
00171 else
00172 eta = userEta;
00173
00174
00175 perturbX = currentX;
00176 Y = X;
00177 Y.Scale(eta);
00178 perturbX.update(1.0,nevY,1.0);
00179
00180 if (!useGroupForComputeF)
00181 interface->computeF(perturbX.getEpetraVector(), fp.getEpetraVector(),
00182 NOX::Epetra::Interface::Required::MF_Res);
00183 else{
00184 groupPtr->setX(perturbX);
00185 groupPtr->computeF();
00186 fp = dynamic_cast<const NOX::Epetra::Vector&>
00187 (groupPtr->getF());
00188 }
00189
00190 if ( diffType == Centered ) {
00191 Y.Scale(-2.0);
00192 perturbX.update(scaleFactor,nevY,1.0);
00193 if (!useGroupForComputeF)
00194 interface->computeF(perturbX.getEpetraVector(), fmPtr->getEpetraVector(),
00195 NOX::Epetra::Interface::Required::MF_Res);
00196 else{
00197 groupPtr->setX(perturbX);
00198 groupPtr->computeF();
00199 *fmPtr = dynamic_cast<const NOX::Epetra::Vector&>
00200 (groupPtr->getF());
00201 }
00202 }
00203
00204
00205 if ( diffType != Centered ) {
00206 nevY.update(1.0, fp, -1.0, fo, 0.0);
00207 nevY.scale( 1.0/(scaleFactor * eta) );
00208 }
00209 else {
00210 nevY.update(1.0, fp, -1.0, *fmPtr, 0.0);
00211 nevY.scale( 1.0/(2.0 * eta) );
00212 }
00213
00214 return 0;
00215 }
00216
00217 int MatrixFree::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00218 {
00219 utils.out() << "ERROR: NOX::MatrixFree::ApplyInverse - Not available for Matrix Free!"
00220 << endl;
00221 throw "NOX Error";
00222 return (-1);
00223 }
00224
00225 double MatrixFree::NormInf() const
00226 {
00227 utils.out() << "ERROR: NOX::Epetra::MatrixFree::NormInf() - Not Available for "
00228 << "Matrix-Free mode!" << endl;
00229 throw "NOX Error";
00230 return 1.0;
00231 }
00232
00233
00234 const char* MatrixFree::Label () const
00235 {
00236 return label.c_str();
00237 }
00238
00239 bool MatrixFree::UseTranspose() const
00240 {
00241 return false;
00242 }
00243
00244 bool MatrixFree::HasNormInf() const
00245 {
00246 return false;
00247 }
00248
00249 const Epetra_Comm & MatrixFree::Comm() const
00250 {
00251 return currentX.getEpetraVector().Map().Comm();
00252 }
00253 const Epetra_Map& MatrixFree::OperatorDomainMap() const
00254 {
00255 return *epetraMap;
00256 }
00257
00258 const Epetra_Map& MatrixFree::OperatorRangeMap() const
00259 {
00260 return *epetraMap;
00261 }
00262
00263 bool MatrixFree::computeJacobian(const Epetra_Vector& x, Epetra_Operator& Jac)
00264 {
00265
00266
00267
00268
00269 currentX = x;
00270
00271 bool ok = false;
00272 if (!useGroupForComputeF)
00273 ok = interface->computeF(x, fo.getEpetraVector(),
00274 NOX::Epetra::Interface::Required::MF_Jac);
00275 else {
00276 groupPtr->setX(currentX);
00277 groupPtr->computeF();
00278 fo = dynamic_cast<const NOX::Epetra::Vector&>
00279 (groupPtr->getF());
00280 ok = true;
00281 }
00282 return ok;
00283 }
00284
00285 void MatrixFree::setDifferenceMethod(DifferenceType diffType_)
00286 {
00287 diffType = diffType_;
00288 }
00289
00290 void MatrixFree::setLambda(double lambda_)
00291 {
00292 lambda = lambda_;
00293 }
00294
00295 void MatrixFree::setComputePerturbation(bool bVal)
00296 {
00297 computeEta = bVal;
00298 }
00299
00300 void MatrixFree::setPerturbation(double eta_)
00301 {
00302 userEta = eta_;
00303 computeEta = false;
00304 }
00305
00306 double MatrixFree::getPerturbation() const
00307 {
00308 return eta;
00309 }
00310
00311 void MatrixFree::setGroupForComputeF(const NOX::Abstract::Group& group)
00312 {
00313 useGroupForComputeF = true;
00314 groupPtr = group.clone();
00315 return;
00316 }