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 "LOCA_Extended_MultiVector.H"
00043 #include "LOCA_Extended_Vector.H"
00044 #include "LOCA_GlobalData.H"
00045 #include "LOCA_ErrorCheck.H"
00046
00047 LOCA::Extended::MultiVector::MultiVector(
00048 const LOCA::Extended::MultiVector& source,
00049 NOX::CopyType type) :
00050 globalData(source.globalData),
00051 numColumns(source.numColumns),
00052 numMultiVecRows(source.numMultiVecRows),
00053 numScalarRows(source.numScalarRows),
00054 multiVectorPtrs(numMultiVecRows),
00055 scalarsPtr(),
00056 extendedVectorPtrs(numColumns),
00057 isView(false)
00058 {
00059
00060 for (int i=0; i<numMultiVecRows; i++)
00061 multiVectorPtrs[i] = source.multiVectorPtrs[i]->clone(type);
00062
00063
00064 scalarsPtr =
00065 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(*source.scalarsPtr));
00066
00067 for (int i=0; i<numColumns; i++) {
00068 extendedVectorPtrs[i] = Teuchos::null;
00069 }
00070 }
00071
00072 LOCA::Extended::MultiVector::MultiVector(
00073 const LOCA::Extended::MultiVector& source,
00074 int nColumns) :
00075 globalData(source.globalData),
00076 numColumns(nColumns),
00077 numMultiVecRows(source.numMultiVecRows),
00078 numScalarRows(source.numScalarRows),
00079 multiVectorPtrs(numMultiVecRows),
00080 scalarsPtr(),
00081 extendedVectorPtrs(numColumns),
00082 isView(false)
00083 {
00084
00085 for (int i=0; i<numMultiVecRows; i++)
00086 multiVectorPtrs[i] = source.multiVectorPtrs[i]->clone(numColumns);
00087
00088 for (int i=0; i<numColumns; i++) {
00089 extendedVectorPtrs[i] = Teuchos::null;
00090 }
00091
00092 scalarsPtr =
00093 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(numScalarRows,
00094 numColumns));
00095 }
00096
00097 LOCA::Extended::MultiVector::MultiVector(
00098 const LOCA::Extended::MultiVector& source,
00099 const vector<int>& index, bool view) :
00100 globalData(source.globalData),
00101 numColumns(index.size()),
00102 numMultiVecRows(source.numMultiVecRows),
00103 numScalarRows(source.numScalarRows),
00104 multiVectorPtrs(numMultiVecRows),
00105 scalarsPtr(),
00106 extendedVectorPtrs(numColumns),
00107 isView(view)
00108 {
00109
00110 for (unsigned int j=0; j<index.size(); j++)
00111 source.checkIndex("LOCA::Extended::MultiVector()", index[j]);
00112
00113 for (int i=0; i<numColumns; i++) {
00114 extendedVectorPtrs[i] = Teuchos::null;
00115 }
00116
00117
00118 bool isCont = isContiguous(index);
00119
00120 if (view) {
00121
00122
00123 for (int i=0; i<numMultiVecRows; i++)
00124 multiVectorPtrs[i] = source.multiVectorPtrs[i]->subView(index);
00125
00126
00127 if (isCont) {
00128 double *vals = source.scalarsPtr->values() +
00129 source.scalarsPtr->numRows()*index[0];
00130 scalarsPtr =
00131 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::View,
00132 vals,
00133 numScalarRows,
00134 numScalarRows,
00135 numColumns));
00136 }
00137 else {
00138 globalData->locaErrorCheck->throwError(
00139 "LOCA::Extended::MultiVector()",
00140 "Sub-view with non-contiguous indices is not supported");
00141 }
00142
00143 }
00144 else {
00145
00146
00147 for (int i=0; i<numMultiVecRows; i++)
00148 multiVectorPtrs[i] = source.multiVectorPtrs[i]->subCopy(index);
00149
00150
00151 if (isCont) {
00152 double *vals = source.scalarsPtr->values() +
00153 source.scalarsPtr->numRows()*index[0];
00154 scalarsPtr =
00155 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::Copy,
00156 vals,
00157 numScalarRows,
00158 numScalarRows,
00159 numColumns));
00160 }
00161 else {
00162 scalarsPtr =
00163 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(numScalarRows,
00164 numColumns));
00165 for (int j=0; j<numColumns; j++)
00166 for (int i=0; i<numScalarRows; i++)
00167 (*scalarsPtr)(i,j) = (*source.scalarsPtr)(i,index[j]);
00168 }
00169 }
00170 }
00171
00172 LOCA::Extended::MultiVector::~MultiVector()
00173 {
00174 }
00175
00176 NOX::Abstract::MultiVector&
00177 LOCA::Extended::MultiVector::init(double gamma)
00178 {
00179
00180 for (int i=0; i<numMultiVecRows; i++)
00181 multiVectorPtrs[i]->init(gamma);
00182
00183
00184 scalarsPtr->putScalar(gamma);
00185
00186 return *this;
00187 }
00188
00189 NOX::Abstract::MultiVector&
00190 LOCA::Extended::MultiVector::random(bool useSeed, int seed)
00191 {
00192
00193 multiVectorPtrs[0]->random(useSeed, seed);
00194 for (int i=1; i<numMultiVecRows; i++)
00195 multiVectorPtrs[i]->random();
00196
00197
00198 scalarsPtr->random();
00199
00200 return *this;
00201 }
00202
00203 NOX::Abstract::MultiVector&
00204 LOCA::Extended::MultiVector::operator=(
00205 const NOX::Abstract::MultiVector& source)
00206 {
00207 operator=(dynamic_cast<const LOCA::Extended::MultiVector&>(source));
00208 return *this;
00209 }
00210
00211 LOCA::Extended::MultiVector&
00212 LOCA::Extended::MultiVector::operator=(
00213 const LOCA::Extended::MultiVector& source)
00214 {
00215 if (this != &source) {
00216
00217
00218 checkDimensions("LOCA::Extended::MultiVector::operator=()", source);
00219
00220 globalData = source.globalData;
00221
00222
00223 for (int i=0; i<numMultiVecRows; i++)
00224 *(multiVectorPtrs[i]) = *(source.multiVectorPtrs[i]);
00225
00226
00227 scalarsPtr->assign(*source.scalarsPtr);
00228 }
00229
00230 return *this;
00231 }
00232
00233 NOX::Abstract::MultiVector&
00234 LOCA::Extended::MultiVector::setBlock(const NOX::Abstract::MultiVector& source,
00235 const vector<int>& index)
00236 {
00237 return setBlock(dynamic_cast<const LOCA::Extended::MultiVector&>(source),
00238 index);
00239 }
00240
00241 NOX::Abstract::MultiVector&
00242 LOCA::Extended::MultiVector::setBlock(
00243 const LOCA::Extended::MultiVector& source,
00244 const vector<int>& index)
00245 {
00246
00247 if (source.numMultiVecRows != numMultiVecRows ||
00248 source.numScalarRows != numScalarRows)
00249 globalData->locaErrorCheck->throwError(
00250 "LOCA::Extended::MultiVector::setBlock()",
00251 "Size of supplied multivector is incompatible with this multivector");
00252 if (static_cast<unsigned int>(source.numColumns) != index.size()) {
00253 globalData->locaErrorCheck->throwError(
00254 "LOCA::Extended::MultiVector::setBlock()",
00255 "Size of supplied index vector is incompatible with this multivector");
00256 }
00257
00258
00259 for (int i=0; i<numMultiVecRows; i++)
00260 multiVectorPtrs[i]->setBlock(*(source.multiVectorPtrs[i]),index);
00261
00262
00263 for (unsigned int j=0; j<index.size(); j++) {
00264 checkIndex("LOCA::Extended::MultiVector::augment()", index[j]);
00265 for (int i=0; i<numScalarRows; i++)
00266 (*scalarsPtr)(i,index[j]) = (*source.scalarsPtr)(i,j);
00267 }
00268
00269 return *this;
00270 }
00271
00272 NOX::Abstract::MultiVector&
00273 LOCA::Extended::MultiVector::augment(const NOX::Abstract::MultiVector& source)
00274 {
00275 return augment(dynamic_cast<const LOCA::Extended::MultiVector&>(source));
00276 }
00277
00278 NOX::Abstract::MultiVector&
00279 LOCA::Extended::MultiVector::augment(const LOCA::Extended::MultiVector& source)
00280 {
00281 if (isView) {
00282 globalData->locaErrorCheck->throwError(
00283 "LOCA::Extended::MultiVector::augment()",
00284 "Augmenting a multivector view is not supported");
00285 }
00286
00287
00288 if (source.numMultiVecRows != numMultiVecRows ||
00289 source.numScalarRows != numScalarRows)
00290 globalData->locaErrorCheck->throwError(
00291 "LOCA::Extended::MultiVector::augment()",
00292 "Size of supplied multivector is incompatible with this multivector");
00293
00294
00295 for (int i=0; i<numMultiVecRows; i++)
00296 multiVectorPtrs[i]->augment(*(source.multiVectorPtrs[i]));
00297
00298
00299 scalarsPtr->reshape(numScalarRows, numColumns + source.numColumns);
00300 for (int j=0; j<source.numColumns; j++)
00301 for (int i=0; i<numScalarRows; i++)
00302 (*scalarsPtr)(i,numColumns + j) = (*source.scalarsPtr)(i,j);
00303
00304
00305 extendedVectorPtrs.resize(numColumns + source.numColumns);
00306 for (int i=0; i<source.numColumns; i++)
00307 extendedVectorPtrs[numColumns + i] = Teuchos::null;
00308
00309
00310 numColumns += source.numColumns;
00311
00312 return *this;
00313 }
00314
00315 NOX::Abstract::Vector&
00316 LOCA::Extended::MultiVector::operator [] (int i)
00317 {
00318 return *(getVector(i));
00319 }
00320
00321 const NOX::Abstract::Vector&
00322 LOCA::Extended::MultiVector::operator [] (int i) const
00323 {
00324 return *(getVector(i));
00325 }
00326
00327 NOX::Abstract::MultiVector&
00328 LOCA::Extended::MultiVector::scale(double gamma)
00329 {
00330
00331 for (int i=0; i<numMultiVecRows; i++)
00332 multiVectorPtrs[i]->scale(gamma);
00333
00334
00335 scalarsPtr->scale(gamma);
00336
00337 return *this;
00338 }
00339
00340 NOX::Abstract::MultiVector&
00341 LOCA::Extended::MultiVector::update(double alpha,
00342 const NOX::Abstract::MultiVector& a,
00343 double gamma)
00344 {
00345 return update(alpha, dynamic_cast<const LOCA::Extended::MultiVector&>(a),
00346 gamma);
00347 }
00348
00349 NOX::Abstract::MultiVector&
00350 LOCA::Extended::MultiVector::update(double alpha,
00351 const LOCA::Extended::MultiVector& a,
00352 double gamma)
00353 {
00354
00355 checkDimensions("LOCA::Extended::MultiVector::update()", a);
00356
00357
00358 for (int i=0; i<numMultiVecRows; i++)
00359 multiVectorPtrs[i]->update(alpha, *(a.multiVectorPtrs[i]), gamma);
00360
00361
00362 for (int j=0; j<numColumns; j++)
00363 for (int i=0; i<numScalarRows; i++)
00364 (*scalarsPtr)(i,j) =
00365 gamma * (*scalarsPtr)(i,j) + alpha * (*a.scalarsPtr)(i,j);
00366
00367 return *this;
00368 }
00369
00370 NOX::Abstract::MultiVector&
00371 LOCA::Extended::MultiVector::update(double alpha,
00372 const NOX::Abstract::MultiVector& a,
00373 double beta,
00374 const NOX::Abstract::MultiVector& b,
00375 double gamma)
00376 {
00377 return update(alpha, dynamic_cast<const LOCA::Extended::MultiVector&>(a),
00378 beta, dynamic_cast<const LOCA::Extended::MultiVector&>(b),
00379 gamma);
00380 }
00381
00382 NOX::Abstract::MultiVector&
00383 LOCA::Extended::MultiVector::update(double alpha,
00384 const LOCA::Extended::MultiVector& a,
00385 double beta,
00386 const LOCA::Extended::MultiVector& b,
00387 double gamma)
00388 {
00389
00390 checkDimensions("LOCA::Extended::MultiVector::update()", a);
00391 checkDimensions("LOCA::Extended::MultiVector::update()", b);
00392
00393
00394 for (int i=0; i<numMultiVecRows; i++)
00395 multiVectorPtrs[i]->update(alpha, *(a.multiVectorPtrs[i]),
00396 beta, *(b.multiVectorPtrs[i]),
00397 gamma);
00398
00399
00400 for (int j=0; j<numColumns; j++)
00401 for (int i=0; i<numScalarRows; i++)
00402 (*scalarsPtr)(i,j) =
00403 gamma * (*scalarsPtr)(i,j) + alpha * (*a.scalarsPtr)(i,j) +
00404 beta * (*b.scalarsPtr)(i,j);
00405
00406 return *this;
00407 }
00408
00409 NOX::Abstract::MultiVector&
00410 LOCA::Extended::MultiVector::update(
00411 Teuchos::ETransp transb,
00412 double alpha,
00413 const NOX::Abstract::MultiVector& a,
00414 const NOX::Abstract::MultiVector::DenseMatrix& b,
00415 double gamma)
00416 {
00417 return update(transb, alpha,
00418 dynamic_cast<const LOCA::Extended::MultiVector&>(a),
00419 b, gamma);
00420 }
00421
00422 NOX::Abstract::MultiVector&
00423 LOCA::Extended::MultiVector::update(
00424 Teuchos::ETransp transb,
00425 double alpha,
00426 const LOCA::Extended::MultiVector& a,
00427 const NOX::Abstract::MultiVector::DenseMatrix& b,
00428 double gamma)
00429 {
00430
00431 if (a.numMultiVecRows != numMultiVecRows ||
00432 a.numScalarRows != numScalarRows)
00433 globalData->locaErrorCheck->throwError(
00434 "LOCA::Extended::MultiVector::update()",
00435 "Size of supplied multivector is incompatible with this multivector");
00436
00437 if (transb == Teuchos::NO_TRANS) {
00438 if (a.numColumns != b.numRows() || numColumns != b.numCols())
00439 globalData->locaErrorCheck->throwError(
00440 "LOCA::Extended::MultiVector::update()",
00441 "Size of supplied matrix is incompatible with this multivector");
00442 }
00443 else
00444 if (a.numColumns != b.numCols() || numColumns != b.numRows())
00445 globalData->locaErrorCheck->throwError(
00446 "LOCA::Extended::MultiVector::update()",
00447 "Size of supplied matrix is incompatible with this multivector");
00448
00449
00450
00451 for (int i=0; i<numMultiVecRows; i++)
00452 multiVectorPtrs[i]->update(transb, alpha, *(a.multiVectorPtrs[i]), b,
00453 gamma);
00454
00455
00456 if (numScalarRows > 0)
00457 scalarsPtr->multiply(Teuchos::NO_TRANS, transb, alpha, *(a.scalarsPtr),
00458 b, gamma);
00459
00460 return *this;
00461 }
00462
00463 Teuchos::RCP<NOX::Abstract::MultiVector>
00464 LOCA::Extended::MultiVector::clone(NOX::CopyType type) const
00465 {
00466 return Teuchos::rcp(new LOCA::Extended::MultiVector(*this, type));
00467 }
00468
00469 Teuchos::RCP<NOX::Abstract::MultiVector>
00470 LOCA::Extended::MultiVector::clone(int numvecs) const
00471 {
00472 return Teuchos::rcp(new LOCA::Extended::MultiVector(*this, numvecs));
00473 }
00474
00475 Teuchos::RCP<NOX::Abstract::MultiVector>
00476 LOCA::Extended::MultiVector::subCopy(const vector<int>& index) const
00477 {
00478 return Teuchos::rcp(new LOCA::Extended::MultiVector(*this, index, false));
00479 }
00480
00481 Teuchos::RCP<NOX::Abstract::MultiVector>
00482 LOCA::Extended::MultiVector::subView(const vector<int>& index) const
00483 {
00484 return Teuchos::rcp(new LOCA::Extended::MultiVector(*this, index, true));
00485 }
00486
00487 void
00488 LOCA::Extended::MultiVector::norm(vector<double>& result,
00489 NOX::Abstract::Vector::NormType type) const
00490 {
00491
00492
00493 if (result.size() != static_cast<unsigned int>(numColumns))
00494 result.resize(numColumns);
00495
00496
00497 for (int i=0; i<numColumns; i++)
00498 result[i] = 0.0;
00499
00500
00501 vector<double> vecNorm(result);
00502
00503 switch (type) {
00504
00505 case NOX::Abstract::Vector::MaxNorm:
00506
00507
00508 for (int i=0; i<numMultiVecRows; i++) {
00509 multiVectorPtrs[i]->norm(vecNorm, type);
00510 for (int j=0; j<numColumns; j++) {
00511 if (result[j] < vecNorm[j])
00512 result[j] = vecNorm[j];
00513 }
00514 }
00515
00516
00517 for (int j=0; j<numColumns; j++) {
00518 for (int i=0; i<numScalarRows; i++)
00519 if (result[j] < (*scalarsPtr)(i,j))
00520 result[j] = (*scalarsPtr)(i,j);
00521 }
00522 break;
00523
00524 case NOX::Abstract::Vector::OneNorm:
00525
00526
00527 for (int i=0; i<numMultiVecRows; i++) {
00528 multiVectorPtrs[i]->norm(vecNorm, type);
00529 for (int j=0; j<numColumns; j++) {
00530 result[j] += vecNorm[j];
00531 }
00532 }
00533
00534
00535 for (int j=0; j<numColumns; j++) {
00536 for (int i=0; i<numScalarRows; i++)
00537 result[j] += fabs((*scalarsPtr)(i,j));
00538 }
00539 break;
00540
00541 case NOX::Abstract::Vector::TwoNorm:
00542 default:
00543
00544
00545 for (int i=0; i<numMultiVecRows; i++) {
00546 multiVectorPtrs[i]->norm(vecNorm, type);
00547 for (int j=0; j<numColumns; j++) {
00548 result[j] += vecNorm[j]*vecNorm[j];
00549 }
00550 }
00551
00552
00553 for (int j=0; j<numColumns; j++) {
00554 for (int i=0; i<numScalarRows; i++)
00555 result[j] += (*scalarsPtr)(i,j) * (*scalarsPtr)(i,j);
00556 result[j] = sqrt(result[j]);
00557 }
00558 break;
00559 }
00560 }
00561
00562 void
00563 LOCA::Extended::MultiVector::multiply(
00564 double alpha,
00565 const NOX::Abstract::MultiVector& y,
00566 NOX::Abstract::MultiVector::DenseMatrix& b) const
00567 {
00568 multiply(alpha, dynamic_cast<const LOCA::Extended::MultiVector&>(y), b);
00569 }
00570
00571 void
00572 LOCA::Extended::MultiVector::multiply(
00573 double alpha,
00574 const LOCA::Extended::MultiVector& y,
00575 NOX::Abstract::MultiVector::DenseMatrix& b) const
00576 {
00577
00578 if (y.numMultiVecRows != numMultiVecRows || y.numColumns != b.numRows() ||
00579 y.numScalarRows != numScalarRows || numColumns != b.numCols())
00580 globalData->locaErrorCheck->throwError(
00581 "LOCA::Extended::MultiVector::multiply()",
00582 "Size of supplied multivector/matrix is incompatible with this multivector");
00583
00584
00585 b.putScalar(0.0);
00586
00587
00588 NOX::Abstract::MultiVector::DenseMatrix tmp(b);
00589
00590
00591 for (int i=0; i<numMultiVecRows; i++) {
00592 multiVectorPtrs[i]->multiply(alpha, *(y.multiVectorPtrs[i]), tmp);
00593 b += tmp;
00594 }
00595
00596
00597 if (numScalarRows > 0)
00598 b.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, alpha, *y.scalarsPtr,
00599 *scalarsPtr, 1.0);
00600
00601 }
00602
00603 int
00604 LOCA::Extended::MultiVector::length() const
00605 {
00606 int len = 0;
00607
00608
00609 for (int i=0; i<numMultiVecRows; i++)
00610 len += multiVectorPtrs[i]->length();
00611
00612
00613 len += numScalarRows;
00614
00615 return len;
00616 }
00617
00618 int
00619 LOCA::Extended::MultiVector::numVectors() const
00620 {
00621 return numColumns;
00622 }
00623
00624 void
00625 LOCA::Extended::MultiVector::print(std::ostream& stream) const
00626 {
00627
00628 for (int i=0; i<numMultiVecRows; i++)
00629 multiVectorPtrs[i]->print(stream);
00630
00631
00632 scalarsPtr->print(stream);
00633 }
00634
00635 Teuchos::RCP<const NOX::Abstract::MultiVector>
00636 LOCA::Extended::MultiVector::getMultiVector(int i) const
00637 {
00638 checkVectorRowIndex("LOCA::Extended::MultiVector::getMultiVector()", i);
00639
00640 return multiVectorPtrs[i];
00641 }
00642
00643 Teuchos::RCP<NOX::Abstract::MultiVector>
00644 LOCA::Extended::MultiVector::getMultiVector(int i)
00645 {
00646 checkVectorRowIndex("LOCA::Extended::MultiVector::getMultiVector()", i);
00647
00648 return multiVectorPtrs[i];
00649 }
00650
00651 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix>
00652 LOCA::Extended::MultiVector::getScalars() const
00653 {
00654 return scalarsPtr;
00655 }
00656
00657 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix>
00658 LOCA::Extended::MultiVector::getScalars()
00659 {
00660 return scalarsPtr;
00661 }
00662
00663 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix>
00664 LOCA::Extended::MultiVector::getScalarRows(int num_rows, int row) const
00665 {
00666 return
00667 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::View,
00668 *scalarsPtr,
00669 num_rows,
00670 numColumns,
00671 row,
00672 0));
00673 }
00674
00675 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix>
00676 LOCA::Extended::MultiVector::getScalarRows(int num_rows, int row)
00677 {
00678 return
00679 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::View,
00680 *scalarsPtr,
00681 num_rows,
00682 numColumns,
00683 row,
00684 0));
00685 }
00686
00687 const double&
00688 LOCA::Extended::MultiVector::getScalar(int i, int j) const
00689 {
00690 checkIndex("LOCA::Extended::MultiVector::getScalar()", i, j);
00691
00692 return (*scalarsPtr)(i,j);
00693 }
00694
00695 double&
00696 LOCA::Extended::MultiVector::getScalar(int i, int j)
00697 {
00698 checkIndex("LOCA::Extended::MultiVector::getScalar()", i, j);
00699
00700 return (*scalarsPtr)(i,j);
00701 }
00702
00703 Teuchos::RCP<LOCA::Extended::Vector>
00704 LOCA::Extended::MultiVector::getVector(int i)
00705 {
00706
00707 checkIndex("LOCA::Extended::MultiVector::vector()", i);
00708
00709
00710 if (extendedVectorPtrs[i] == Teuchos::null) {
00711 extendedVectorPtrs[i] = generateVector(numMultiVecRows, numScalarRows);
00712 for (int k=0; k<numMultiVecRows; k++)
00713 extendedVectorPtrs[i]->setVectorView(
00714 k,
00715 Teuchos::rcp(&(*multiVectorPtrs[k])[i],
00716 false));
00717 if (numScalarRows > 0)
00718 extendedVectorPtrs[i]->setScalarArray((*scalarsPtr)[i]);
00719 }
00720
00721 return extendedVectorPtrs[i];
00722 }
00723
00724 Teuchos::RCP<const LOCA::Extended::Vector>
00725 LOCA::Extended::MultiVector::getVector(int i) const
00726 {
00727
00728 checkIndex("LOCA::Extended::MultiVector::vector()", i);
00729
00730
00731 if (extendedVectorPtrs[i] == Teuchos::null) {
00732 extendedVectorPtrs[i] = generateVector(numMultiVecRows, numScalarRows);
00733 for (int k=0; k<numMultiVecRows; k++)
00734 extendedVectorPtrs[i]->setVectorView(
00735 k,
00736 Teuchos::rcp(&(*multiVectorPtrs[k])[i],
00737 false));
00738 if (numScalarRows > 0)
00739 extendedVectorPtrs[i]->setScalarArray((*scalarsPtr)[i]);
00740 }
00741
00742 return extendedVectorPtrs[i];
00743 }
00744
00745 int
00746 LOCA::Extended::MultiVector::getNumScalarRows() const
00747 {
00748 return numScalarRows;
00749 }
00750
00751 int
00752 LOCA::Extended::MultiVector::getNumMultiVectors() const
00753 {
00754 return numMultiVecRows;
00755 }
00756
00757 LOCA::Extended::MultiVector::MultiVector(
00758 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00759 int nColumns, int nVectorRows,
00760 int nScalarRows) :
00761 globalData(global_data),
00762 numColumns(nColumns),
00763 numMultiVecRows(nVectorRows),
00764 numScalarRows(nScalarRows),
00765 multiVectorPtrs(numMultiVecRows),
00766 scalarsPtr(),
00767 extendedVectorPtrs(numColumns),
00768 isView(false)
00769 {
00770 for (int i=0; i<numColumns; i++) {
00771 extendedVectorPtrs[i] = Teuchos::null;
00772 }
00773
00774 scalarsPtr =
00775 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(numScalarRows,
00776 numColumns));
00777 }
00778
00779 Teuchos::RCP<LOCA::Extended::Vector>
00780 LOCA::Extended::MultiVector::generateVector(int nVecs, int nScalarRows) const
00781 {
00782 return Teuchos::rcp(new LOCA::Extended::Vector(globalData, nVecs,
00783 nScalarRows));
00784 }
00785
00786 void
00787 LOCA::Extended::MultiVector::setMultiVectorPtr(
00788 int i,
00789 Teuchos::RCP<NOX::Abstract::MultiVector> v)
00790 {
00791 checkVectorRowIndex("LOCA::Extended::MultiVector::setMultiVectorPtr()",i);
00792
00793 multiVectorPtrs[i] = v;
00794 }
00795
00796 void
00797 LOCA::Extended::MultiVector::checkDimensions(
00798 const string& callingFunction,
00799 const LOCA::Extended::MultiVector& a) const
00800 {
00801 if (a.numMultiVecRows != numMultiVecRows || a.numColumns != numColumns ||
00802 a.numScalarRows != numScalarRows)
00803 globalData->locaErrorCheck->throwError(callingFunction,
00804 "Size of supplied multivector is incompatible with this multivector");
00805 }
00806
00807 void
00808 LOCA::Extended::MultiVector::checkIndex(const string& callingFunction,
00809 int i) const
00810 {
00811 if ( i < 0 || i >= numColumns )
00812 globalData->locaErrorCheck->throwError(callingFunction,
00813 "Invalid column index");
00814 }
00815
00816 void
00817 LOCA::Extended::MultiVector::checkVectorRowIndex(const string& callingFunction,
00818 int i) const
00819 {
00820 if ( i < 0 || i >= numMultiVecRows)
00821 globalData->locaErrorCheck->throwError(callingFunction,
00822 "Invalid vector row index");
00823 }
00824 void
00825 LOCA::Extended::MultiVector::checkIndex(const string& callingFunction,
00826 int i, int j) const
00827 {
00828 if ( i < 0 || i >= numScalarRows )
00829 globalData->locaErrorCheck->throwError(callingFunction,
00830 "Invalid row index");
00831 if ( j < 0 || j >= numColumns )
00832 globalData->locaErrorCheck->throwError(callingFunction,
00833 "Invalid column index");
00834 }
00835
00836 bool
00837 LOCA::Extended::MultiVector::isContiguous(const vector<int>& index) const
00838 {
00839 for (unsigned int i=0; i<index.size(); i++) {
00840 if (static_cast<unsigned int>(index[i]) != index[0] + i)
00841 return false;
00842 }
00843 return true;
00844 }