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_Random.H"
00043 #include "LOCA_Extended_Vector.H"
00044 #include "LOCA_Extended_MultiVector.H"
00045 #include "LOCA_GlobalData.H"
00046 #include "LOCA_ErrorCheck.H"
00047
00048 LOCA::Extended::Vector::Vector(
00049 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00050 int nvecs, int nscalars) :
00051 globalData(global_data),
00052 vectorPtrs(nvecs),
00053 isView(nvecs),
00054 numScalars(nscalars),
00055 scalarsPtr()
00056 {
00057 scalarsPtr =
00058 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(numScalars, 1));
00059 }
00060
00061 LOCA::Extended::Vector::Vector(const LOCA::Extended::Vector& source,
00062 NOX::CopyType type) :
00063 globalData(source.globalData),
00064 vectorPtrs(source.vectorPtrs.size()),
00065 isView(source.vectorPtrs.size()),
00066 numScalars(source.numScalars),
00067 scalarsPtr()
00068 {
00069 for (unsigned int i=0; i<vectorPtrs.size(); i++) {
00070 vectorPtrs[i] = source.vectorPtrs[i]->clone(type);
00071 isView[i] = false;
00072 }
00073
00074 scalarsPtr =
00075 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(*source.scalarsPtr));
00076
00077 if (type != NOX::DeepCopy)
00078 init(0.0);
00079 }
00080
00081
00082 LOCA::Extended::Vector::~Vector()
00083 {
00084 }
00085
00086 NOX::Abstract::Vector&
00087 LOCA::Extended::Vector::operator=(const NOX::Abstract::Vector& y)
00088 {
00089 operator=(dynamic_cast<const LOCA::Extended::Vector&>(y));
00090 return *this;
00091 }
00092
00093 LOCA::Extended::Vector&
00094 LOCA::Extended::Vector::operator=(const LOCA::Extended::Vector& y)
00095 {
00096 if (this != &y) {
00097
00098 if (numScalars != y.numScalars)
00099 globalData->locaErrorCheck->throwError(
00100 "LOCA::Extended::Vector::operator=()",
00101 "Number of scalars must match in assignment");
00102 if (vectorPtrs.size() != y.vectorPtrs.size())
00103 globalData->locaErrorCheck->throwError(
00104 "LOCA::Extended::Vector::operator=()",
00105 "Number of vectors must match in assignment");
00106
00107 globalData = y.globalData;
00108
00109 for (unsigned int i=0; i<vectorPtrs.size(); i++)
00110 *(vectorPtrs[i]) = *(y.vectorPtrs[i]);
00111
00112 numScalars = y.numScalars;
00113
00114
00115 scalarsPtr->assign(*y.scalarsPtr);
00116
00117 }
00118 return *this;
00119 }
00120
00121 Teuchos::RCP<NOX::Abstract::Vector>
00122 LOCA::Extended::Vector::clone(NOX::CopyType type) const
00123 {
00124 return Teuchos::rcp(new LOCA::Extended::Vector(*this, type));
00125 }
00126
00127 Teuchos::RCP<NOX::Abstract::MultiVector>
00128 LOCA::Extended::Vector::createMultiVector(
00129 const NOX::Abstract::Vector* const* vecs,
00130 int numVecs,
00131 NOX::CopyType type) const
00132 {
00133
00134 const NOX::Abstract::Vector** subvecs =
00135 new const NOX::Abstract::Vector*[numVecs+1];
00136
00137
00138 Teuchos::RCP<NOX::Abstract::MultiVector> submvec;
00139
00140
00141 Teuchos::RCP<LOCA::Extended::MultiVector> mvec =
00142 generateMultiVector(numVecs+1, vectorPtrs.size(), numScalars);
00143
00144 const LOCA::Extended::Vector *evec;
00145
00146
00147 for (unsigned int i=0; i<vectorPtrs.size(); i++) {
00148
00149
00150 subvecs[0] = vectorPtrs[i].get();
00151 for (int j=0; j<numVecs; j++) {
00152 evec = dynamic_cast<const LOCA::Extended::Vector*>(vecs[j]);
00153 subvecs[j+1] = evec->vectorPtrs[i].get();
00154 }
00155
00156
00157 submvec = vectorPtrs[i]->createMultiVector(subvecs, numVecs+1, type);
00158
00159
00160 mvec->setMultiVectorPtr(i, submvec);
00161
00162 }
00163
00164
00165 for (int i=0; i<numScalars; i++)
00166 mvec->getScalar(i,0) = (*scalarsPtr)(i,0);
00167 for (int j=0; j<numVecs; j++) {
00168 evec = dynamic_cast<const LOCA::Extended::Vector*>(vecs[j]);
00169 for (int i=0; i<numScalars; i++)
00170 mvec->getScalar(i,j+1) = (*evec->scalarsPtr)(i,0);
00171 }
00172
00173 delete [] subvecs;
00174
00175 return mvec;
00176
00177 }
00178
00179 Teuchos::RCP<NOX::Abstract::MultiVector>
00180 LOCA::Extended::Vector::createMultiVector(int numVecs,
00181 NOX::CopyType type) const
00182 {
00183
00184 Teuchos::RCP<NOX::Abstract::MultiVector> submvec;
00185
00186
00187 Teuchos::RCP<LOCA::Extended::MultiVector> mvec =
00188 generateMultiVector(numVecs, vectorPtrs.size(), numScalars);
00189
00190
00191 for (unsigned int i=0; i<vectorPtrs.size(); i++) {
00192
00193
00194 submvec = vectorPtrs[i]->createMultiVector(numVecs, type);
00195
00196
00197 mvec->setMultiVectorPtr(i, submvec);
00198
00199 }
00200
00201
00202 if (type == NOX::DeepCopy) {
00203 for (int j=0; j<numVecs; j++) {
00204 for (int i=0; i<numScalars; i++)
00205 mvec->getScalar(i,j) = (*scalarsPtr)(i,0);
00206 }
00207 }
00208
00209 return mvec;
00210 }
00211
00212 NOX::Abstract::Vector&
00213 LOCA::Extended::Vector::init(double gamma)
00214 {
00215 for (unsigned int i=0; i<vectorPtrs.size(); i++)
00216 vectorPtrs[i]->init(gamma);
00217 scalarsPtr->putScalar(gamma);
00218 return *this;
00219 }
00220
00221 NOX::Abstract::Vector&
00222 LOCA::Extended::Vector::random(bool useSeed, int seed) {
00223 if (useSeed)
00224 NOX::Random::setSeed(seed);
00225 if (vectorPtrs.size() > 0)
00226 vectorPtrs[0]->random(useSeed, seed);
00227 for (unsigned int i=1; i<vectorPtrs.size(); i++)
00228 vectorPtrs[i]->random();
00229 scalarsPtr->random();
00230 return *this;
00231 }
00232
00233 NOX::Abstract::Vector&
00234 LOCA::Extended::Vector::abs(const NOX::Abstract::Vector& y)
00235 {
00236 const LOCA::Extended::Vector& Y =
00237 dynamic_cast<const LOCA::Extended::Vector&>(y);
00238
00239 for (unsigned int i=0; i<vectorPtrs.size(); i++)
00240 vectorPtrs[i]->abs(*(Y.vectorPtrs[i]));
00241 for (int i=0; i<numScalars; i++)
00242 (*scalarsPtr)(i,0) = fabs((*Y.scalarsPtr)(i,0));
00243
00244 return *this;
00245 }
00246
00247 NOX::Abstract::Vector&
00248 LOCA::Extended::Vector::reciprocal(const NOX::Abstract::Vector& y)
00249 {
00250 const LOCA::Extended::Vector& Y =
00251 dynamic_cast<const LOCA::Extended::Vector&>(y);
00252
00253 for (unsigned int i=0; i<vectorPtrs.size(); i++)
00254 vectorPtrs[i]->reciprocal(*(Y.vectorPtrs[i]));
00255 for (int i=0; i<numScalars; i++)
00256 (*scalarsPtr)(i,0) = 1.0 / (*Y.scalarsPtr)(i,0);
00257
00258 return *this;
00259 }
00260
00261 NOX::Abstract::Vector&
00262 LOCA::Extended::Vector::scale(double gamma)
00263 {
00264 for (unsigned int i=0; i<vectorPtrs.size(); i++)
00265 vectorPtrs[i]->scale(gamma);
00266 scalarsPtr->scale(gamma);
00267
00268 return *this;
00269 }
00270
00271 NOX::Abstract::Vector&
00272 LOCA::Extended::Vector::scale(const NOX::Abstract::Vector& y)
00273 {
00274 const LOCA::Extended::Vector& Y =
00275 dynamic_cast<const LOCA::Extended::Vector&>(y);
00276
00277 for (unsigned int i=0; i<vectorPtrs.size(); i++)
00278 vectorPtrs[i]->scale(*(Y.vectorPtrs[i]));
00279 scalarsPtr->scale(*Y.scalarsPtr);
00280
00281 return *this;
00282 }
00283
00284 NOX::Abstract::Vector&
00285 LOCA::Extended::Vector::update(double alpha, const NOX::Abstract::Vector& y,
00286 double gamma)
00287 {
00288 const LOCA::Extended::Vector& Y =
00289 dynamic_cast<const LOCA::Extended::Vector&>(y);
00290
00291 for (unsigned int i=0; i<vectorPtrs.size(); i++)
00292 vectorPtrs[i]->update(alpha, *(Y.vectorPtrs[i]), gamma);
00293 for (int i=0; i<numScalars; i++)
00294 (*scalarsPtr)(i,0) = alpha*((*Y.scalarsPtr)(i,0)) +
00295 gamma*((*scalarsPtr)(i,0));
00296
00297 return *this;
00298 }
00299
00300 NOX::Abstract::Vector&
00301 LOCA::Extended::Vector::update(double alpha,
00302 const NOX::Abstract::Vector& y,
00303 double beta, const NOX::Abstract::Vector& b,
00304 double gamma)
00305 {
00306 const LOCA::Extended::Vector& Y =
00307 dynamic_cast<const LOCA::Extended::Vector&>(y);
00308 const LOCA::Extended::Vector& B =
00309 dynamic_cast<const LOCA::Extended::Vector&>(b);
00310
00311 for (unsigned int i=0; i<vectorPtrs.size(); i++)
00312 vectorPtrs[i]->update(alpha, *(Y.vectorPtrs[i]), beta, *(B.vectorPtrs[i]),
00313 gamma);
00314 for (int i=0; i<numScalars; i++)
00315 (*scalarsPtr)(i,0) = alpha*((*Y.scalarsPtr)(i,0)) +
00316 beta*((*B.scalarsPtr)(i,0)) + gamma*((*scalarsPtr)(i,0));
00317
00318 return *this;
00319 }
00320
00321 double
00322 LOCA::Extended::Vector::norm(NormType type) const
00323 {
00324 double n = 0.0;
00325 double nv;
00326
00327 switch (type) {
00328
00329 case NOX::Abstract::Vector::MaxNorm:
00330 for (unsigned int i=0; i<vectorPtrs.size(); i++)
00331 n = NOX_MAX(n, vectorPtrs[i]->norm(type));
00332 n = NOX_MAX(n, scalarsPtr->normInf());
00333 break;
00334
00335 case NOX::Abstract::Vector::OneNorm:
00336 for (unsigned int i=0; i<vectorPtrs.size(); i++)
00337 n += vectorPtrs[i]->norm(type);
00338 n += scalarsPtr->normOne();
00339 break;
00340
00341 case NOX::Abstract::Vector::TwoNorm:
00342 default:
00343 for (unsigned int i=0; i<vectorPtrs.size(); i++) {
00344 nv = vectorPtrs[i]->norm(type);
00345 n += nv*nv;
00346 }
00347 nv = scalarsPtr->normFrobenius();
00348 n += nv*nv;
00349 n = sqrt(n);
00350 break;
00351
00352 }
00353
00354 return n;
00355 }
00356
00357 double
00358 LOCA::Extended::Vector::norm(const NOX::Abstract::Vector& weights) const
00359 {
00360 const LOCA::Extended::Vector& W =
00361 dynamic_cast<const LOCA::Extended::Vector&>(weights);
00362 double n = 0.0;
00363 double nv;
00364
00365 for (unsigned int i=0; i<vectorPtrs.size(); i++) {
00366 nv = vectorPtrs[i]->norm(*(W.vectorPtrs[i]));
00367 n += nv*nv;
00368 }
00369 for (int i=0; i<numScalars; i++)
00370 n += ((*W.scalarsPtr)(i,0))*((*scalarsPtr)(i,0))*((*scalarsPtr)(i,0));
00371 n = sqrt(n);
00372
00373 return n;
00374 }
00375
00376 double
00377 LOCA::Extended::Vector::innerProduct(const NOX::Abstract::Vector& y) const
00378 {
00379 const LOCA::Extended::Vector& Y =
00380 dynamic_cast<const LOCA::Extended::Vector&>(y);
00381 double d = 0.0;
00382
00383 for (unsigned int i=0; i<vectorPtrs.size(); i++)
00384 d += vectorPtrs[i]->innerProduct(*(Y.vectorPtrs[i]));
00385 for (int i=0; i<numScalars; i++)
00386 d += ((*scalarsPtr)(i,0))*((*Y.scalarsPtr)(i,0));
00387
00388 return d;
00389 }
00390
00391 int
00392 LOCA::Extended::Vector::length() const
00393 {
00394 int l = 0;
00395 for (unsigned int i=0; i<vectorPtrs.size(); i++)
00396 l += vectorPtrs[i]->length();
00397 l += numScalars;
00398
00399 return l;
00400 }
00401
00402 void
00403 LOCA::Extended::Vector::print(std::ostream& stream) const
00404 {
00405 for (unsigned int i=0; i<vectorPtrs.size(); i++)
00406 vectorPtrs[i]->print(stream);
00407 scalarsPtr->print(stream);
00408 stream << std::endl;
00409
00410 }
00411
00412 void
00413 LOCA::Extended::Vector::setVector(int i, const NOX::Abstract::Vector& v)
00414 {
00415 if (vectorPtrs[i] != Teuchos::null)
00416 *(vectorPtrs[i]) = v;
00417 else
00418 vectorPtrs[i] = v.clone(NOX::DeepCopy);
00419 isView[i] = false;
00420 }
00421
00422 void
00423 LOCA::Extended::Vector::setVectorView(
00424 int i,
00425 const Teuchos::RCP<NOX::Abstract::Vector>& v)
00426 {
00427 vectorPtrs[i] = v;
00428 isView[i] = true;
00429 }
00430
00431 void
00432 LOCA::Extended::Vector::setScalar(int i, double s)
00433 {
00434 (*scalarsPtr)(i,0) = s;
00435 }
00436
00437 void
00438 LOCA::Extended::Vector::setScalarArray(double *sv)
00439 {
00440 scalarsPtr =
00441 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::View,
00442 sv,
00443 numScalars,
00444 numScalars,
00445 1));
00446 }
00447
00448 Teuchos::RCP<const NOX::Abstract::Vector>
00449 LOCA::Extended::Vector::getVector(int i) const
00450 {
00451 return vectorPtrs[i];
00452 }
00453
00454 Teuchos::RCP<NOX::Abstract::Vector>
00455 LOCA::Extended::Vector::getVector(int i)
00456 {
00457 return vectorPtrs[i];
00458 }
00459
00460 double
00461 LOCA::Extended::Vector::getScalar(int i) const
00462 {
00463 return (*scalarsPtr)(i,0);
00464 }
00465
00466 double&
00467 LOCA::Extended::Vector::getScalar(int i)
00468 {
00469 return (*scalarsPtr)(i,0);
00470 }
00471
00472 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix>
00473 LOCA::Extended::Vector::getScalars() const
00474 {
00475 return scalarsPtr;
00476 }
00477
00478 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix>
00479 LOCA::Extended::Vector::getScalars()
00480 {
00481 return scalarsPtr;
00482 }
00483
00484 int
00485 LOCA::Extended::Vector::getNumScalars() const
00486 {
00487 return numScalars;
00488 }
00489
00490 int
00491 LOCA::Extended::Vector::getNumVectors() const
00492 {
00493 return vectorPtrs.size();
00494 }
00495
00496 Teuchos::RCP<LOCA::Extended::MultiVector>
00497 LOCA::Extended::Vector::generateMultiVector(int nColumns, int nVectorRows,
00498 int nScalarRows) const
00499 {
00500 return Teuchos::rcp(new LOCA::Extended::MultiVector(globalData,
00501 nColumns,
00502 nVectorRows,
00503 nScalarRows));
00504 }
00505