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_MultiVector.H"
00043
00044 NOX::MultiVector::MultiVector(int numVecs) :
00045 vecs(numVecs)
00046 {
00047 if (numVecs <= 0) {
00048 cerr << "NOX::MultiVector: Error! Multivector"
00049 << " must have postive number of columns!" << endl;
00050 throw "NOX Error";
00051 }
00052 }
00053
00054 NOX::MultiVector::MultiVector(const NOX::Abstract::Vector& v, int numVecs,
00055 NOX::CopyType type) :
00056 vecs(numVecs)
00057 {
00058 if (numVecs <= 0) {
00059 cerr << "NOX::MultiVector: Error! Multivector"
00060 << " must have postive number of columns!" << endl;
00061 throw "NOX Error";
00062 }
00063
00064 for (int i=0; i<numVecs; i++) {
00065 vecs[i] = v.clone(type);
00066 }
00067 }
00068
00069 NOX::MultiVector::MultiVector(const NOX::Abstract::Vector* const* vs,
00070 int numVecs,
00071 NOX::CopyType type) :
00072 vecs(numVecs)
00073 {
00074 if (numVecs <= 0) {
00075 cerr << "NOX::MultiVector: Error! Multivector"
00076 << " must have postive number of columns!" << endl;
00077 throw "NOX Error";
00078 }
00079
00080 for (int i=0; i<numVecs; i++) {
00081 vecs[i] = vs[i]->clone(type);
00082 }
00083 }
00084
00085 NOX::MultiVector::MultiVector(const NOX::MultiVector& source,
00086 NOX::CopyType type) :
00087 vecs(source.vecs.size())
00088 {
00089 for (unsigned int i=0; i<source.vecs.size(); i++) {
00090 vecs[i] = source.vecs[i]->clone(type);
00091 }
00092 }
00093
00094 NOX::MultiVector::~MultiVector()
00095 {
00096
00097 }
00098
00099 NOX::Abstract::MultiVector&
00100 NOX::MultiVector::init(double gamma)
00101 {
00102 for (unsigned int i=0; i<vecs.size(); i++)
00103 vecs[i]->init(gamma);
00104 return *this;
00105 }
00106
00107 NOX::Abstract::MultiVector&
00108 NOX::MultiVector::random(bool useSeed, int seed)
00109 {
00110 if (vecs.size() > 0)
00111 vecs[0]->random(useSeed, seed);
00112 for (unsigned int i=1; i<vecs.size(); i++)
00113 vecs[i]->random();
00114 return *this;
00115 }
00116
00117 NOX::Abstract::MultiVector&
00118 NOX::MultiVector::operator=(const NOX::Abstract::MultiVector& source)
00119 {
00120 return operator=(dynamic_cast<const NOX::MultiVector&>(source));
00121 }
00122
00123 NOX::Abstract::MultiVector&
00124 NOX::MultiVector::operator=(const NOX::MultiVector& source)
00125 {
00126 if (this != &source) {
00127 checkSize(source.vecs.size());
00128 for (unsigned int i=0; i<vecs.size(); i++)
00129 *(vecs[i]) = *(source.vecs[i]);
00130 }
00131
00132 return *this;
00133 }
00134
00135 NOX::Abstract::MultiVector&
00136 NOX::MultiVector::setBlock(const NOX::Abstract::MultiVector& source,
00137 const vector<int>& index)
00138 {
00139 return setBlock(dynamic_cast<const NOX::MultiVector&>(source), index);
00140 }
00141
00142 NOX::Abstract::MultiVector&
00143 NOX::MultiVector::setBlock(const NOX::MultiVector& source,
00144 const vector<int>& index)
00145 {
00146 int ind;
00147
00148 source.checkIndex(index.size()-1);
00149 for (unsigned int i=0; i<index.size(); i++) {
00150 ind = index[i];
00151 checkIndex(ind);
00152 *(vecs[ind]) = *(source.vecs[i]);
00153 }
00154
00155 return *this;
00156 }
00157
00158 NOX::Abstract::MultiVector&
00159 NOX::MultiVector::augment(const NOX::Abstract::MultiVector& source)
00160 {
00161 return augment(dynamic_cast<const NOX::MultiVector&>(source));
00162 }
00163
00164 NOX::Abstract::MultiVector&
00165 NOX::MultiVector::augment(const NOX::MultiVector& source)
00166 {
00167 int sz = vecs.size();
00168 int newsize = sz + source.vecs.size();
00169 vecs.resize(newsize);
00170
00171 for (unsigned int i=0; i<source.vecs.size(); i++) {
00172 vecs[sz+i] = source.vecs[i]->clone(NOX::DeepCopy);
00173 }
00174
00175 return *this;
00176 }
00177
00178 NOX::Abstract::Vector&
00179 NOX::MultiVector::operator [] (int i)
00180 {
00181 checkIndex(i);
00182 return *(vecs[i]);
00183 }
00184
00185 const NOX::Abstract::Vector&
00186 NOX::MultiVector::operator [] (int i) const
00187 {
00188 checkIndex(i);
00189 return *(vecs[i]);
00190 }
00191
00192 NOX::Abstract::MultiVector&
00193 NOX::MultiVector::scale(double gamma)
00194 {
00195 for (unsigned int i=0; i<vecs.size(); i++)
00196 vecs[i]->scale(gamma);
00197 return *this;
00198 }
00199
00200 NOX::Abstract::MultiVector&
00201 NOX::MultiVector::update(double alpha, const NOX::Abstract::MultiVector& a,
00202 double gamma)
00203 {
00204 return update(alpha, dynamic_cast<const NOX::MultiVector&>(a), gamma);
00205 }
00206
00207 NOX::Abstract::MultiVector&
00208 NOX::MultiVector::update(double alpha, const NOX::MultiVector& a,
00209 double gamma)
00210 {
00211 checkSize(a.vecs.size());
00212 for (unsigned int i=0; i<vecs.size(); i++)
00213 vecs[i]->update(alpha, *(a.vecs[i]), gamma);
00214 return *this;
00215 }
00216
00217 NOX::Abstract::MultiVector&
00218 NOX::MultiVector::update(double alpha, const NOX::Abstract::MultiVector& a,
00219 double beta, const NOX::Abstract::MultiVector& b,
00220 double gamma)
00221 {
00222 return update(alpha, dynamic_cast<const NOX::MultiVector&>(a),
00223 beta, dynamic_cast<const NOX::MultiVector&>(b), gamma);
00224 }
00225
00226 NOX::Abstract::MultiVector&
00227 NOX::MultiVector::update(double alpha, const NOX::MultiVector& a,
00228 double beta, const NOX::MultiVector& b,
00229 double gamma)
00230 {
00231 checkSize(a.vecs.size());
00232 checkSize(b.vecs.size());
00233 for (unsigned int i=0; i<vecs.size(); i++)
00234 vecs[i]->update(alpha, *(a.vecs[i]), beta, *(b.vecs[i]), gamma);
00235 return *this;
00236 }
00237
00238 NOX::Abstract::MultiVector&
00239 NOX::MultiVector::update(Teuchos::ETransp transb, double alpha,
00240 const NOX::Abstract::MultiVector& a,
00241 const NOX::Abstract::MultiVector::DenseMatrix& b,
00242 double gamma)
00243 {
00244 return update(transb, alpha, dynamic_cast<const NOX::MultiVector&>(a), b,
00245 gamma);
00246 }
00247
00248 NOX::Abstract::MultiVector&
00249 NOX::MultiVector::update(Teuchos::ETransp transb, double alpha,
00250 const NOX::MultiVector& a,
00251 const NOX::Abstract::MultiVector::DenseMatrix& b,
00252 double gamma)
00253 {
00254 if (transb == Teuchos::NO_TRANS) {
00255 a.checkSize(b.numRows());
00256 checkSize(b.numCols());
00257 }
00258 else {
00259 a.checkSize(b.numCols());
00260 checkSize(b.numRows());
00261 }
00262
00263 int sz_a = a.vecs.size();
00264 int p = sz_a / 2;
00265 int q = sz_a - 2*p;
00266
00267 if (transb == Teuchos::NO_TRANS) {
00268 for (unsigned int i=0; i<vecs.size(); i++) {
00269
00270 if (p == 0)
00271 vecs[i]->update(alpha*b(0,i), *(a.vecs[0]), gamma);
00272 else {
00273 vecs[i]->update(alpha*b(0,i), *(a.vecs[0]),
00274 alpha*b(1,i), *(a.vecs[1]), gamma);
00275
00276 for (int j=1; j<p; j++)
00277 vecs[i]->update(alpha*b(2*j,i), *(a.vecs[2*j]),
00278 alpha*b(2*j+1,i), *(a.vecs[2*j+1]), 1.0);
00279
00280 if (q > 0)
00281 vecs[i]->update(alpha*b(sz_a-1,i), *(a.vecs[sz_a-1]), 1.0);
00282 }
00283
00284 }
00285
00286 }
00287 else {
00288
00289 for (unsigned int i=0; i<vecs.size(); i++) {
00290
00291 if (p == 0)
00292 vecs[i]->update(alpha*b(i,0), *(a.vecs[0]), gamma);
00293 else {
00294 vecs[i]->update(alpha*b(i,0), *(a.vecs[0]),
00295 alpha*b(i,1), *(a.vecs[1]), gamma);
00296
00297 for (int j=1; j<p; j++)
00298 vecs[i]->update(alpha*b(i,2*j), *(a.vecs[2*j]),
00299 alpha*b(i,2*j+1), *(a.vecs[2*j+1]), 1.0);
00300
00301 if (q > 0)
00302 vecs[i]->update(alpha*b(i,sz_a-1), *(a.vecs[sz_a-1]), 1.0);
00303 }
00304
00305 }
00306
00307 }
00308 return *this;
00309 }
00310
00311 Teuchos::RCP<NOX::Abstract::MultiVector>
00312 NOX::MultiVector::clone(NOX::CopyType type) const
00313 {
00314 Teuchos::RCP<NOX::Abstract::MultiVector> tmp =
00315 Teuchos::rcp(new NOX::MultiVector(*this, type));
00316 return tmp;
00317 }
00318
00319 Teuchos::RCP<NOX::Abstract::MultiVector>
00320 NOX::MultiVector::clone(int numvecs) const
00321 {
00322 Teuchos::RCP<NOX::MultiVector> tmp =
00323 Teuchos::rcp(new NOX::MultiVector(numvecs));
00324 for (int i=0; i<numvecs; i++) {
00325 tmp->vecs[i] = vecs[0]->clone(NOX::ShapeCopy);
00326 }
00327 return tmp;
00328 }
00329
00330 Teuchos::RCP<NOX::Abstract::MultiVector>
00331 NOX::MultiVector::subCopy(const vector<int>& index) const
00332 {
00333 Teuchos::RCP<NOX::MultiVector> tmp =
00334 Teuchos::rcp(new NOX::MultiVector(index.size()));
00335 int ind;
00336 for (unsigned int i=0; i<index.size(); i++) {
00337 ind = index[i];
00338 checkIndex(ind);
00339 tmp->vecs[i] = vecs[ind]->clone(NOX::DeepCopy);
00340 }
00341 return tmp;
00342 }
00343
00344 Teuchos::RCP<NOX::Abstract::MultiVector>
00345 NOX::MultiVector::subView(const vector<int>& index) const
00346 {
00347 Teuchos::RCP<NOX::MultiVector> tmp =
00348 Teuchos::rcp(new NOX::MultiVector(index.size()));
00349 int ind;
00350 for (unsigned int i=0; i<index.size(); i++) {
00351 ind = index[i];
00352 checkIndex(ind);
00353 tmp->vecs[i] = vecs[ind];
00354 }
00355 return tmp;
00356 }
00357
00358 void
00359 NOX::MultiVector::norm(vector<double>& result,
00360 NOX::Abstract::Vector::NormType type) const
00361 {
00362 if (result.size() != vecs.size())
00363 result.resize(vecs.size());
00364
00365 for (unsigned int i=0; i<vecs.size(); i++)
00366 result[i] = vecs[i]->norm(type);
00367 }
00368
00369 void
00370 NOX::MultiVector::multiply(double alpha, const NOX::Abstract::MultiVector& y,
00371 NOX::Abstract::MultiVector::DenseMatrix& b) const
00372 {
00373 multiply(alpha, dynamic_cast<const NOX::MultiVector&>(y), b);
00374 }
00375
00376 void
00377 NOX::MultiVector::multiply(double alpha, const NOX::MultiVector& y,
00378 NOX::Abstract::MultiVector::DenseMatrix& b) const
00379 {
00380 for (unsigned int i=0; i<y.vecs.size(); i++) {
00381 for (unsigned int j=0; j<vecs.size(); j++) {
00382 b(i,j) = alpha*(y.vecs[i]->innerProduct(*(vecs[j])));
00383 }
00384 }
00385 }
00386
00387 int
00388 NOX::MultiVector::length() const
00389 {
00390 return vecs[0]->length();
00391 }
00392
00393 int
00394 NOX::MultiVector::numVectors() const
00395 {
00396 return vecs.size();
00397 }
00398
00399 void
00400 NOX::MultiVector::print(std::ostream& stream) const
00401 {
00402 for (unsigned int i=0; i<vecs.size(); i++)
00403 vecs[i]->print(stream);
00404 }
00405
00406 void
00407 NOX::MultiVector::checkIndex(int idx) const
00408 {
00409 if ( idx < 0 || idx >= static_cast<int>(vecs.size()) ) {
00410 cerr << "NOX::MultiVector: Error! Invalid index " << idx << endl;
00411 throw "NOX Error";
00412 }
00413 }
00414
00415 void
00416 NOX::MultiVector::checkSize(int sz) const
00417 {
00418 if (static_cast<int>(vecs.size()) != sz) {
00419 cerr << "NOX::MultiVector: Error! Size of supplied multivector is"
00420 << " incompatible with this multivector" << endl;
00421 throw "NOX Error";
00422 }
00423 }