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_Thyra_MultiVector.H"
00043 #include "NOX_Thyra_Vector.H"
00044 #include "Thyra_VectorSpaceBase.hpp"
00045 #include "Thyra_MultiVectorStdOps.hpp"
00046 #include "Thyra_DetachedMultiVectorView.hpp"
00047
00048 NOX::Thyra::MultiVector::
00049 MultiVector(
00050 const Teuchos::RCP< ::Thyra::MultiVectorBase<double> >& source)
00051 : thyraMultiVec(source),
00052 noxThyraVectors(thyraMultiVec->domain()->dim())
00053 {
00054 }
00055
00056 NOX::Thyra::MultiVector::
00057 MultiVector(const ::Thyra::MultiVectorBase<double>& source)
00058 : thyraMultiVec(source.clone_mv()),
00059 noxThyraVectors(thyraMultiVec->domain()->dim())
00060 {
00061 }
00062
00063 NOX::Thyra::MultiVector::
00064 MultiVector(const NOX::Thyra::MultiVector& source, NOX::CopyType type)
00065 : thyraMultiVec(source.thyraMultiVec->clone_mv()),
00066 noxThyraVectors(thyraMultiVec->domain()->dim())
00067 {
00068 }
00069
00070 NOX::Thyra::MultiVector::
00071 ~MultiVector()
00072 {
00073 }
00074
00075 NOX::Abstract::MultiVector&
00076 NOX::Thyra::MultiVector::
00077 operator=(const NOX::Abstract::MultiVector& src)
00078 {
00079 const NOX::Thyra::MultiVector& source =
00080 dynamic_cast<const NOX::Thyra::MultiVector&>(src);
00081 ::Thyra::assign(thyraMultiVec.get(), *source.thyraMultiVec);
00082 return *this;
00083 }
00084
00085 Teuchos::RCP< ::Thyra::MultiVectorBase<double> >
00086 NOX::Thyra::MultiVector::
00087 getThyraMultiVector()
00088 {
00089 return thyraMultiVec;
00090 }
00091
00092 Teuchos::RCP< const ::Thyra::MultiVectorBase<double> >
00093 NOX::Thyra::MultiVector::
00094 getThyraMultiVector() const
00095 {
00096 return thyraMultiVec;
00097 }
00098
00099 NOX::Abstract::MultiVector&
00100 NOX::Thyra::MultiVector::
00101 init(double value)
00102 {
00103 ::Thyra::assign(thyraMultiVec.get(), value);
00104 return *this;
00105 }
00106
00107 NOX::Abstract::MultiVector&
00108 NOX::Thyra::MultiVector::
00109 random(bool useSeed, int seed)
00110 {
00111 if (useSeed)
00112 ::Thyra::seed_randomize<double>(seed);
00113 ::Thyra::randomize(-1.0, 1.0, thyraMultiVec.get());
00114 return *this;
00115 }
00116
00117 NOX::Abstract::MultiVector&
00118 NOX::Thyra::MultiVector::
00119 setBlock(const NOX::Abstract::MultiVector& src, const vector<int>& index)
00120 {
00121 const NOX::Thyra::MultiVector& source =
00122 dynamic_cast<const NOX::Thyra::MultiVector&>(src);
00123
00124
00125 Teuchos::RCP< ::Thyra::MultiVectorBase<double> > v;
00126 if (isContiguous(index))
00127 v = thyraMultiVec->subView(::Thyra::Range1D(index[0],
00128 index[index.size()-1]));
00129 else
00130 v = thyraMultiVec->subView(index.size(), &index[0]);
00131
00132
00133 ::Thyra::assign(v.get(), *source.thyraMultiVec);
00134
00135 return *this;
00136 }
00137
00138 NOX::Abstract::MultiVector&
00139 NOX::Thyra::MultiVector::
00140 augment(const NOX::Abstract::MultiVector& src)
00141 {
00142 const NOX::Thyra::MultiVector& source =
00143 dynamic_cast<const NOX::Thyra::MultiVector&>(src);
00144
00145
00146 int num_cols = thyraMultiVec->domain()->dim() +
00147 source.thyraMultiVec->domain()->dim();
00148
00149
00150 Teuchos::RCP< ::Thyra::MultiVectorBase<double> > new_mv =
00151 ::Thyra::createMembers(thyraMultiVec->range(), num_cols);
00152
00153
00154 Teuchos::RCP< ::Thyra::MultiVectorBase<double> > v1 =
00155 new_mv->subView(::Thyra::Range1D(0,thyraMultiVec->domain()->dim()-1));
00156 Teuchos::RCP< ::Thyra::MultiVectorBase<double> > v2 =
00157 new_mv->subView(::Thyra::Range1D(thyraMultiVec->domain()->dim(),
00158 num_cols-1));
00159
00160
00161 ::Thyra::assign(v1.get(), *thyraMultiVec);
00162 ::Thyra::assign(v2.get(), *source.thyraMultiVec);
00163
00164
00165 thyraMultiVec = new_mv;
00166
00167 return *this;
00168 }
00169
00170 NOX::Abstract::Vector&
00171 NOX::Thyra::MultiVector::
00172 operator [] (int i)
00173 {
00174 if (noxThyraVectors[i] == Teuchos::null)
00175 noxThyraVectors[i] =
00176 Teuchos::rcp(new NOX::Thyra::Vector(thyraMultiVec->col(i)));
00177 return *noxThyraVectors[i];
00178 }
00179
00180 const NOX::Abstract::Vector&
00181 NOX::Thyra::MultiVector::
00182 operator [] (int i) const
00183 {
00184 if (noxThyraVectors[i] == Teuchos::null)
00185 noxThyraVectors[i] =
00186 Teuchos::rcp(new NOX::Thyra::Vector(thyraMultiVec->col(i)));
00187 return *noxThyraVectors[i];
00188 }
00189
00190 NOX::Abstract::MultiVector&
00191 NOX::Thyra::MultiVector::
00192 scale(double gamma)
00193 {
00194 ::Thyra::scale(gamma, thyraMultiVec.get());
00195 return *this;
00196 }
00197
00198 NOX::Abstract::MultiVector&
00199 NOX::Thyra::MultiVector::
00200 update(double alpha, const NOX::Abstract::MultiVector& a, double gamma)
00201 {
00202 const NOX::Thyra::MultiVector& aa =
00203 dynamic_cast<const NOX::Thyra::MultiVector&>(a);
00204 const ::Thyra::MultiVectorBase<double> *tv = aa.thyraMultiVec.get();
00205 ::Thyra::linear_combination(1, &alpha, &tv, gamma, thyraMultiVec.get());
00206 return *this;
00207 }
00208
00209 NOX::Abstract::MultiVector&
00210 NOX::Thyra::MultiVector::
00211 update(double alpha,
00212 const NOX::Abstract::MultiVector& a,
00213 double beta,
00214 const NOX::Abstract::MultiVector& b,
00215 double gamma)
00216 {
00217 const NOX::Thyra::MultiVector& aa =
00218 dynamic_cast<const NOX::Thyra::MultiVector&>(a);
00219 const NOX::Thyra::MultiVector& bb =
00220 dynamic_cast<const NOX::Thyra::MultiVector&>(b);
00221 double c[] = {alpha, beta};
00222 const ::Thyra::MultiVectorBase<double>* z[] = {aa.thyraMultiVec.get(),
00223 bb.thyraMultiVec.get()};
00224 ::Thyra::linear_combination(2, c, z, gamma, thyraMultiVec.get());
00225 return *this;
00226 }
00227
00228 NOX::Abstract::MultiVector&
00229 NOX::Thyra::MultiVector::
00230 update(Teuchos::ETransp transb, double alpha,
00231 const NOX::Abstract::MultiVector& a,
00232 const NOX::Abstract::MultiVector::DenseMatrix& b, double gamma)
00233 {
00234 const NOX::Thyra::MultiVector& aa =
00235 dynamic_cast<const NOX::Thyra::MultiVector&>(a);
00236
00237 int m = b.numRows();
00238 int n = b.numCols();
00239 Teuchos::RCP<const ::Thyra::MultiVectorBase<double> > bb;
00240 if (transb == Teuchos::NO_TRANS) {
00241 bb = ::Thyra::createMembersView(
00242 aa.thyraMultiVec->domain(),
00243 RTOpPack::ConstSubMultiVectorView<double>(0,m,0,n,&b(0,0),b.stride()));
00244 }
00245 else {
00246
00247 NOX::Abstract::MultiVector::DenseMatrix btrans(n, m);
00248 for (int i=0; i<m; i++)
00249 for (int j=0; j<n; j++)
00250 btrans(j,i) = b(i,j);
00251 bb = ::Thyra::createMembersView(
00252 aa.thyraMultiVec->domain(),
00253 RTOpPack::ConstSubMultiVectorView<double>(0,n,0,n,&btrans(0,0),
00254 btrans.stride()));
00255 }
00256 aa.thyraMultiVec->apply(::Thyra::NONCONJ_ELE, *bb, thyraMultiVec.get(),
00257 alpha, gamma);
00258
00259 return *this;
00260 }
00261
00262 Teuchos::RCP<NOX::Abstract::MultiVector>
00263 NOX::Thyra::MultiVector::
00264 clone(CopyType type) const
00265 {
00266 return Teuchos::rcp(new NOX::Thyra::MultiVector(*this, type));
00267 }
00268
00269 Teuchos::RCP<NOX::Abstract::MultiVector>
00270 NOX::Thyra::MultiVector::
00271 clone(int numvecs) const
00272 {
00273 Teuchos::RCP< ::Thyra::MultiVectorBase<double> > mv =
00274 ::Thyra::createMembers(thyraMultiVec->range(), numvecs);
00275 return Teuchos::rcp(new NOX::Thyra::MultiVector(mv));
00276 }
00277
00278 Teuchos::RCP<NOX::Abstract::MultiVector>
00279 NOX::Thyra::MultiVector::
00280 subCopy(const vector<int>& index) const
00281 {
00282 Teuchos::RCP< ::Thyra::MultiVectorBase<double> > mv;
00283 if (isContiguous(index))
00284 mv = thyraMultiVec->subView(::Thyra::Range1D(index[0],
00285 index[index.size()-1]));
00286 else
00287 mv = thyraMultiVec->subView(index.size(), &index[0]);
00288
00289 return Teuchos::rcp(new NOX::Thyra::MultiVector(*mv));
00290 }
00291
00292 Teuchos::RCP<NOX::Abstract::MultiVector>
00293 NOX::Thyra::MultiVector::
00294 subView(const vector<int>& index) const
00295 {
00296 Teuchos::RCP< ::Thyra::MultiVectorBase<double> > mv;
00297 if (isContiguous(index))
00298 mv = thyraMultiVec->subView(::Thyra::Range1D(index[0],
00299 index[index.size()-1]));
00300 else
00301 mv = thyraMultiVec->subView(index.size(), &index[0]);
00302
00303 return Teuchos::rcp(new NOX::Thyra::MultiVector(mv));
00304 }
00305
00306 void
00307 NOX::Thyra::MultiVector::
00308 norm(vector<double>& result, NOX::Abstract::Vector::NormType type) const
00309 {
00310 if (type == NOX::Abstract::Vector::TwoNorm)
00311 ::Thyra::norms_2(*thyraMultiVec, &result[0]);
00312 else if (type == NOX::Abstract::Vector::OneNorm)
00313 ::Thyra::norms_1(*thyraMultiVec, &result[0]);
00314 else
00315 ::Thyra::norms_inf(*thyraMultiVec, &result[0]);
00316 }
00317
00318 void
00319 NOX::Thyra::MultiVector::
00320 multiply(double alpha,
00321 const NOX::Abstract::MultiVector& y,
00322 NOX::Abstract::MultiVector::DenseMatrix& b) const
00323 {
00324 const NOX::Thyra::MultiVector& yy =
00325 dynamic_cast<const NOX::Thyra::MultiVector&>(y);
00326
00327 int m = b.numRows();
00328 int n = b.numCols();
00329 Teuchos::RCP< ::Thyra::MultiVectorBase<double> > bb =
00330 ::Thyra::createMembersView(
00331 yy.thyraMultiVec->domain(),
00332 RTOpPack::SubMultiVectorView<double>(0,m,0,n,&b(0,0),b.stride()));
00333
00334 yy.thyraMultiVec->applyTranspose(::Thyra::NONCONJ_ELE, *thyraMultiVec,
00335 bb.get(), alpha, 0.0);
00336 }
00337
00338 int
00339 NOX::Thyra::MultiVector::
00340 length() const
00341 {
00342 return thyraMultiVec->range()->dim();
00343 }
00344
00345 int
00346 NOX::Thyra::MultiVector::
00347 numVectors() const
00348 {
00349 return thyraMultiVec->domain()->dim();
00350 }
00351
00352 void
00353 NOX::Thyra::MultiVector::
00354 print(std::ostream& stream) const
00355 {
00356
00357 thyraMultiVec->describe(*Teuchos::fancyOStream<char>(Teuchos::rcp(&stream,false)), Teuchos::VERB_EXTREME);
00358 return;
00359 }
00360
00361 bool
00362 NOX::Thyra::MultiVector::
00363 isContiguous(const vector<int>& index) const
00364 {
00365 if (index.size()==0)
00366 return true;
00367 int i0 = index[0];
00368 for (std::size_t i=1; i<index.size(); i++)
00369 if (index[i] != static_cast<int>(i0+i))
00370 return false;
00371 return true;
00372 }