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
00040
00041
00042 #include "NOX_Epetra_MultiVector.H"
00043 #include "NOX_Epetra_Vector.H"
00044
00045 #include "Epetra_LocalMap.h"
00046 #include "Epetra_MultiVector.h"
00047 #include "Epetra_Vector.h"
00048
00049 NOX::Epetra::MultiVector::MultiVector(int numvecs)
00050 : noxEpetraVectors(numvecs)
00051 {
00052 for (unsigned int i=0; i<noxEpetraVectors.size(); i++)
00053 noxEpetraVectors[i] = NULL;
00054 }
00055
00056 NOX::Epetra::MultiVector::
00057 MultiVector(const Teuchos::RCP<Epetra_MultiVector>& source,
00058 NOX::CopyType type,
00059 NOX::Epetra::MultiVector::MemoryType memoryType)
00060 : noxEpetraVectors(source->NumVectors())
00061 {
00062 if (memoryType == NOX::Epetra::MultiVector::CreateView) {
00063 epetraMultiVec = source;
00064 }
00065 else {
00066
00067 switch (type) {
00068
00069 case DeepCopy:
00070
00071 epetraMultiVec = Teuchos::rcp(new Epetra_MultiVector(*source));
00072 break;
00073
00074 case ShapeCopy:
00075
00076 epetraMultiVec =
00077 Teuchos::rcp(new Epetra_MultiVector(source->Map(),
00078 source->NumVectors()));
00079 break;
00080 }
00081 }
00082
00083 for (unsigned int i=0; i<noxEpetraVectors.size(); i++)
00084 noxEpetraVectors[i] = NULL;
00085 }
00086
00087 NOX::Epetra::MultiVector::MultiVector(const Epetra_MultiVector& source,
00088 NOX::CopyType type)
00089 : noxEpetraVectors(source.NumVectors())
00090 {
00091 switch (type) {
00092
00093 case DeepCopy:
00094
00095 epetraMultiVec = Teuchos::rcp(new Epetra_MultiVector(source));
00096 break;
00097
00098 case ShapeCopy:
00099
00100 epetraMultiVec =
00101 Teuchos::rcp(new Epetra_MultiVector(source.Map(), source.NumVectors()));
00102 break;
00103
00104 }
00105
00106 for (unsigned int i=0; i<noxEpetraVectors.size(); i++)
00107 noxEpetraVectors[i] = NULL;
00108 }
00109
00110 NOX::Epetra::MultiVector::MultiVector(const NOX::Epetra::MultiVector& source,
00111 NOX::CopyType type)
00112 : noxEpetraVectors(source.epetraMultiVec->NumVectors())
00113 {
00114
00115 switch (type) {
00116
00117 case DeepCopy:
00118
00119 epetraMultiVec = Teuchos::rcp(new Epetra_MultiVector(source.getEpetraMultiVector()));
00120 break;
00121
00122 case ShapeCopy:
00123
00124 epetraMultiVec =
00125 Teuchos::rcp(new Epetra_MultiVector(source.getEpetraMultiVector().Map(),
00126 source.numVectors()));
00127 break;
00128
00129 }
00130
00131 for (unsigned int i=0; i<noxEpetraVectors.size(); i++)
00132 noxEpetraVectors[i] = NULL;
00133 }
00134
00135 NOX::Epetra::MultiVector::~MultiVector()
00136 {
00137 for (unsigned int i=0; i<noxEpetraVectors.size(); i++)
00138 if (noxEpetraVectors[i] != NULL)
00139 delete noxEpetraVectors[i];
00140 }
00141
00142 NOX::Abstract::MultiVector&
00143 NOX::Epetra::MultiVector::operator=(const Epetra_MultiVector& source)
00144 {
00145 *epetraMultiVec = source;
00146 return *this;
00147 }
00148
00149 NOX::Abstract::MultiVector&
00150 NOX::Epetra::MultiVector::operator=(const NOX::Abstract::MultiVector& source)
00151 {
00152 return operator=(dynamic_cast<const NOX::Epetra::MultiVector&>(source));
00153 }
00154
00155 NOX::Abstract::MultiVector&
00156 NOX::Epetra::MultiVector::operator=(const NOX::Epetra::MultiVector& source)
00157 {
00158 *epetraMultiVec = *source.epetraMultiVec;
00159 return *this;
00160 }
00161
00162 Epetra_MultiVector&
00163 NOX::Epetra::MultiVector::getEpetraMultiVector()
00164 {
00165 return *epetraMultiVec;
00166 }
00167
00168 const Epetra_MultiVector&
00169 NOX::Epetra::MultiVector::getEpetraMultiVector() const
00170 {
00171 return *epetraMultiVec;
00172 }
00173
00174 NOX::Abstract::MultiVector&
00175 NOX::Epetra::MultiVector::init(double value)
00176 {
00177 epetraMultiVec->PutScalar(value);
00178 return *this;
00179 }
00180
00181 NOX::Abstract::MultiVector&
00182 NOX::Epetra::MultiVector::random(bool useSeed, int seed)
00183 {
00184 if (useSeed)
00185 epetraMultiVec->SetSeed(seed);
00186 epetraMultiVec->Random();
00187 return *this;
00188 }
00189
00190 NOX::Abstract::MultiVector&
00191 NOX::Epetra::MultiVector::setBlock(const NOX::Abstract::MultiVector& source,
00192 const vector<int>& index) {
00193 return setBlock(dynamic_cast<const NOX::Epetra::MultiVector&>(source),
00194 index);
00195 }
00196
00197 NOX::Abstract::MultiVector&
00198 NOX::Epetra::MultiVector::setBlock(const NOX::Epetra::MultiVector& source,
00199 const vector<int>& index) {
00200 double* vecPtr;
00201 double *sourceVecPtr;
00202 int ind;
00203 int vecLength = epetraMultiVec->MyLength();
00204
00205 source.checkIndex(index.size()-1);
00206 for (unsigned int i=0; i<index.size(); i++) {
00207 ind = index[i];
00208 checkIndex(ind);
00209 vecPtr = epetraMultiVec->operator[] (ind);
00210 sourceVecPtr = source.epetraMultiVec->operator[] (i);
00211 for (int j = 0; j<vecLength; j++)
00212 vecPtr[j] = sourceVecPtr[j];
00213 }
00214 return *this;
00215 }
00216
00217 NOX::Abstract::MultiVector&
00218 NOX::Epetra::MultiVector::augment(const NOX::Abstract::MultiVector& source) {
00219 return augment(dynamic_cast<const NOX::Epetra::MultiVector&>(source));
00220 }
00221
00222 NOX::Abstract::MultiVector&
00223 NOX::Epetra::MultiVector::augment(const NOX::Epetra::MultiVector& source) {
00224 int numVecs = numVectors();
00225 int sourceNumVecs = source.numVectors();
00226 int totalNumVecs = numVecs + sourceNumVecs;
00227 int vecLength = epetraMultiVec->MyLength();
00228 double* vecPtr;
00229 double *tmpVecPtr;
00230
00231 Epetra_MultiVector *tmp = new Epetra_MultiVector(epetraMultiVec->Map(),
00232 totalNumVecs);
00233 for (int i=0; i<numVecs; i++) {
00234 vecPtr = epetraMultiVec->operator[] (i);
00235 tmpVecPtr = tmp->operator[] (i);
00236 for (int j = 0; j<vecLength; j++)
00237 tmpVecPtr[j] = vecPtr[j];
00238 }
00239
00240 for (int i=0; i<sourceNumVecs; i++) {
00241 vecPtr = source.epetraMultiVec->operator[] (i);
00242 tmpVecPtr = tmp->operator[] (numVecs + i);
00243 for (int j = 0; j<vecLength; j++)
00244 tmpVecPtr[j] = vecPtr[j];
00245 }
00246
00247 epetraMultiVec = Teuchos::rcp(tmp);
00248
00249 return *this;
00250 }
00251
00252 NOX::Abstract::Vector&
00253 NOX::Epetra::MultiVector::operator [] (int i)
00254 {
00255 if ( i < 0 || i > (int) noxEpetraVectors.size() ) {
00256 cerr << "NOX::Epetra::MultiVector::operator[]: Error! Invalid index "
00257 << i << endl;
00258 throw "NOX::Epetra Error";
00259 }
00260 if (noxEpetraVectors[i] == NULL) {
00261 Teuchos::RCP<Epetra_Vector> epetra_vec =
00262 Teuchos::rcp(epetraMultiVec->operator() (i), false);
00263 noxEpetraVectors[i] =
00264 new NOX::Epetra::Vector(epetra_vec, NOX::Epetra::Vector::CreateView);
00265 }
00266 return *(noxEpetraVectors[i]);
00267 }
00268
00269 const NOX::Abstract::Vector&
00270 NOX::Epetra::MultiVector::operator [] (int i) const
00271 {
00272 if ( i < 0 || i > (int) noxEpetraVectors.size() ) {
00273 cerr << "NOX::Epetra::MultiVector::operator[]: Error! Invalid index "
00274 << i << endl;
00275 throw "NOX::Epetra Error";
00276 }
00277 if (noxEpetraVectors[i] == NULL) {
00278 Teuchos::RCP<Epetra_Vector> epetra_vec =
00279 Teuchos::rcp(epetraMultiVec->operator() (i), false);
00280 noxEpetraVectors[i] =
00281 new NOX::Epetra::Vector(epetra_vec, NOX::Epetra::Vector::CreateView);
00282 }
00283 return *(noxEpetraVectors[i]);
00284 }
00285
00286 NOX::Abstract::MultiVector&
00287 NOX::Epetra::MultiVector::scale(double gamma)
00288 {
00289 epetraMultiVec->Scale(gamma);
00290 return *this;
00291 }
00292
00293 NOX::Abstract::MultiVector&
00294 NOX::Epetra::MultiVector::update(double alpha,
00295 const NOX::Abstract::MultiVector& a,
00296 double gamma)
00297 {
00298 return update(alpha, dynamic_cast<const NOX::Epetra::MultiVector&>(a),
00299 gamma);
00300 }
00301
00302 NOX::Abstract::MultiVector&
00303 NOX::Epetra::MultiVector::update(double alpha,
00304 const NOX::Epetra::MultiVector& a,
00305 double gamma)
00306 {
00307 epetraMultiVec->Update(alpha, a.getEpetraMultiVector(), gamma);
00308 return *this;
00309 }
00310
00311 NOX::Abstract::MultiVector&
00312 NOX::Epetra::MultiVector::update(double alpha,
00313 const NOX::Abstract::MultiVector& a,
00314 double beta,
00315 const NOX::Abstract::MultiVector& b,
00316 double gamma)
00317 {
00318 return update(alpha, dynamic_cast<const NOX::Epetra::MultiVector&>(a),
00319 beta, dynamic_cast<const NOX::Epetra::MultiVector&>(b), gamma);
00320 }
00321
00322 NOX::Abstract::MultiVector&
00323 NOX::Epetra::MultiVector::update(double alpha,
00324 const NOX::Epetra::MultiVector& a,
00325 double beta,
00326 const NOX::Epetra::MultiVector& b,
00327 double gamma)
00328 {
00329 epetraMultiVec->Update(alpha, a.getEpetraMultiVector(),
00330 beta, b.getEpetraMultiVector(), gamma);
00331 return *this;
00332 }
00333
00334 NOX::Abstract::MultiVector&
00335 NOX::Epetra::MultiVector::update(
00336 Teuchos::ETransp transb,
00337 double alpha,
00338 const NOX::Abstract::MultiVector& a,
00339 const NOX::Abstract::MultiVector::DenseMatrix& b,
00340 double gamma)
00341 {
00342 return update(transb, alpha,
00343 dynamic_cast<const NOX::Epetra::MultiVector&>(a),
00344 b, gamma);
00345 }
00346
00347 NOX::Abstract::MultiVector&
00348 NOX::Epetra::MultiVector::update(
00349 Teuchos::ETransp transb,
00350 double alpha,
00351 const NOX::Epetra::MultiVector& a,
00352 const NOX::Abstract::MultiVector::DenseMatrix& b,
00353 double gamma)
00354 {
00355
00356 NOX::Abstract::MultiVector::DenseMatrix& nc_b =
00357 const_cast<NOX::Abstract::MultiVector::DenseMatrix&>(b);
00358
00359
00360 const int izero = 0;
00361 Epetra_LocalMap localMap(b.numRows(), izero, epetraMultiVec->Map().Comm());
00362 Epetra_MultiVector B(View, localMap, nc_b.values(), b.stride(), b.numCols());
00363
00364 char epetra_trans_b;
00365 if (transb == Teuchos::NO_TRANS)
00366 epetra_trans_b = 'N';
00367 else
00368 epetra_trans_b = 'T';
00369
00370 epetraMultiVec->Multiply('N', epetra_trans_b, alpha,
00371 a.getEpetraMultiVector(), B, gamma);
00372 return *this;
00373 }
00374
00375 Teuchos::RCP<NOX::Abstract::MultiVector>
00376 NOX::Epetra::MultiVector::clone(CopyType type) const
00377 {
00378 Teuchos::RCP<NOX::Abstract::MultiVector> newVec =
00379 Teuchos::rcp(new NOX::Epetra::MultiVector(*epetraMultiVec, type));
00380 return newVec;
00381 }
00382
00383 Teuchos::RCP<NOX::Abstract::MultiVector>
00384 NOX::Epetra::MultiVector::clone(int numvecs) const
00385 {
00386 Teuchos::RCP<NOX::Epetra::MultiVector> newVec =
00387 Teuchos::rcp(new NOX::Epetra::MultiVector(numvecs));
00388 newVec->epetraMultiVec =
00389 Teuchos::rcp(new Epetra_MultiVector(epetraMultiVec->Map(), numvecs));
00390 return newVec;
00391 }
00392
00393 Teuchos::RCP<NOX::Abstract::MultiVector>
00394 NOX::Epetra::MultiVector::subCopy(const vector<int>& index) const
00395 {
00396 int numvecs = index.size();
00397 Teuchos::RCP<NOX::Epetra::MultiVector> newVec =
00398 Teuchos::rcp(new NOX::Epetra::MultiVector(numvecs));
00399 newVec->epetraMultiVec =
00400 Teuchos::rcp(new Epetra_MultiVector(Copy, *epetraMultiVec,
00401 const_cast<int*>(&index[0]),
00402 numvecs));
00403 return newVec;
00404 }
00405
00406 Teuchos::RCP<NOX::Abstract::MultiVector>
00407 NOX::Epetra::MultiVector::subView(const vector<int>& index) const
00408 {
00409 int numvecs = index.size();
00410 Teuchos::RCP<NOX::Epetra::MultiVector> newVec =
00411 Teuchos::rcp(new NOX::Epetra::MultiVector(numvecs));
00412 newVec->epetraMultiVec =
00413 Teuchos::rcp(new Epetra_MultiVector(View, *epetraMultiVec,
00414 const_cast<int*>(&index[0]),
00415 numvecs));
00416 return newVec;
00417 }
00418
00419 void
00420 NOX::Epetra::MultiVector::norm(vector<double>& result,
00421 NOX::Abstract::Vector::NormType type) const
00422 {
00423 switch (type) {
00424 case NOX::Abstract::Vector::MaxNorm:
00425 epetraMultiVec->NormInf(&result[0]);
00426 break;
00427 case NOX::Abstract::Vector::OneNorm:
00428 epetraMultiVec->Norm1(&result[0]);
00429 break;
00430 case NOX::Abstract::Vector::TwoNorm:
00431 default:
00432 epetraMultiVec->Norm2(&result[0]);
00433 break;
00434 }
00435 }
00436
00437 void
00438 NOX::Epetra::MultiVector::multiply(
00439 double alpha,
00440 const NOX::Abstract::MultiVector& y,
00441 NOX::Abstract::MultiVector::DenseMatrix& b) const
00442 {
00443 multiply(alpha, dynamic_cast<const NOX::Epetra::MultiVector&>(y), b);
00444 }
00445
00446 void
00447 NOX::Epetra::MultiVector::multiply(
00448 double alpha,
00449 const NOX::Epetra::MultiVector& y,
00450 NOX::Abstract::MultiVector::DenseMatrix& b) const
00451 {
00452
00453 const int izero = 0;
00454 Epetra_LocalMap localMap(b.numRows(), izero, epetraMultiVec->Map().Comm());
00455 Epetra_MultiVector B(View, localMap, b.values(), b.stride(), b.numCols());
00456
00457
00458 B.Multiply('T', 'N', alpha, *(y.epetraMultiVec), *epetraMultiVec, 0.0);
00459 }
00460
00461 int NOX::Epetra::MultiVector::length() const
00462 {
00463 return epetraMultiVec->GlobalLength();
00464 }
00465
00466 int NOX::Epetra::MultiVector::numVectors() const
00467 {
00468 return epetraMultiVec->NumVectors();
00469 }
00470
00471 void NOX::Epetra::MultiVector::print(std::ostream& stream) const
00472 {
00473 epetraMultiVec->Print(stream);
00474 return;
00475 }
00476
00477 void NOX::Epetra::MultiVector::checkIndex(int idx) const
00478 {
00479 if ( idx < 0 || idx >= epetraMultiVec->NumVectors() ) {
00480 cerr << "NOX::Epetra::MultiVector: Error! Invalid index " << idx << endl;
00481 throw "NOX::Epetra Error";
00482 }
00483 }