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 "petscvec.h"
00043
00044 #include "NOX_Common.H"
00045 #include "NOX_Petsc_Vector.H"
00046
00047 using namespace NOX;
00048 using namespace NOX::Petsc;
00049
00050 Vector::Vector(const Vec& source, CopyType type)
00051 {
00052 allocate(source, type);
00053 }
00054
00055 Vector::Vector(const Vector& source, CopyType type)
00056 {
00057 allocate(source.getPetscVector(), type);
00058 }
00059
00060 Vector::Vector(const Vec& source, string Name, CopyType type) :
00061 name(Name)
00062 {
00063 allocate(source, type);
00064 }
00065
00066 Vector::~Vector()
00067 {
00068 if(isAlloc)
00069 {
00070
00071
00072
00073 isAlloc = false;
00074 }
00075 }
00076
00077 Abstract::Vector&
00078 Vector::operator=(const Vec& source)
00079 {
00080 VecCopy(source, petscVec);
00081 name = "Unnamed";
00082 isAlloc = true;
00083 return *this;
00084 }
00085
00086 Abstract::Vector&
00087 Vector::operator=(const Abstract::Vector& source)
00088 {
00089 return operator=(dynamic_cast<const Vector&>(source));
00090 }
00091
00092 Abstract::Vector&
00093 Vector::operator=(const Vector& source)
00094 {
00095 VecCopy(source.getPetscVector(), petscVec);
00096 isAlloc = source.isAlloc;
00097 name = source.name;
00098 return *this;
00099 }
00100
00101 int
00102 Vector::allocate(const Vec& source, CopyType type)
00103 {
00104
00105 int ierr = VecDuplicate(source, &petscVec);
00106 isAlloc = true;
00107
00108 switch (type) {
00109
00110 case DeepCopy:
00111
00112 ierr += VecCopy(source, petscVec);
00113 break;
00114
00115 case ShapeCopy:
00116
00117 break;
00118 }
00119
00120 if(ierr)
00121 cout << "ERROR: value " << ierr << " returned during "
00122 << "NOX::Petsc::Vector allocation !!" << endl;
00123
00124 return ierr;
00125 }
00126
00127 Vec&
00128 Vector::getPetscVector()
00129 {
00130 return petscVec;
00131 }
00132
00133 const Vec&
00134 Vector::getPetscVector() const
00135 {
00136 return petscVec;
00137 }
00138
00139 Abstract::Vector&
00140 Vector::init(double value)
00141 {
00142 VecSet(petscVec, value);
00143 return *this;
00144 }
00145
00146 Abstract::Vector&
00147 Vector::abs(const Abstract::Vector& base)
00148 {
00149 return abs(dynamic_cast<const Vector&>(base));
00150 }
00151
00152 Abstract::Vector&
00153 Vector::abs(const Vector& base)
00154 {
00155 VecCopy(base.getPetscVector(), petscVec);
00156 VecAbs(petscVec);
00157 return *this;
00158 }
00159
00160 Abstract::Vector&
00161 Vector::reciprocal(const Abstract::Vector& base)
00162 {
00163 return reciprocal(dynamic_cast<const Vector&>(base));
00164 }
00165
00166 Abstract::Vector&
00167 Vector::reciprocal(const Vector& base)
00168 {
00169 VecCopy(base.getPetscVector(), petscVec);
00170 VecReciprocal(petscVec);
00171 return *this;
00172 }
00173
00174 Abstract::Vector&
00175 Vector::scale(double alpha)
00176 {
00177 VecScale(petscVec, alpha);
00178 return *this;
00179 }
00180
00181 Abstract::Vector&
00182 Vector::scale(const Abstract::Vector& a)
00183 {
00184 return scale(dynamic_cast<const Petsc::Vector&>(a));
00185 }
00186
00187 Abstract::Vector&
00188 Vector::scale(const Vector& a)
00189 {
00190 VecPointwiseMult(petscVec, a.getPetscVector(), petscVec);
00191 return *this;
00192 }
00193
00194 Abstract::Vector&
00195 Vector::update(double alpha, const Abstract::Vector& a,
00196 double gammaval)
00197 {
00198 return update(alpha, dynamic_cast<const Vector&>(a), gammaval);
00199 }
00200
00201 Abstract::Vector&
00202 Vector::update(double alpha, const Vector& a,
00203 double gammaval)
00204 {
00205 VecAXPBY(petscVec, alpha, gammaval, a.getPetscVector());
00206 return *this;
00207 }
00208
00209 Abstract::Vector&
00210 Vector::update(double alpha, const Abstract::Vector& a,
00211 double beta, const Abstract::Vector& b,
00212 double gammaval)
00213 {
00214 return update(alpha, dynamic_cast<const Vector&>(a),
00215 beta, dynamic_cast<const Vector&>(b), gammaval);
00216 }
00217
00218 Abstract::Vector&
00219 Vector::update(double alpha, const Vector& a,
00220 double beta, const Vector& b,
00221 double gammaval)
00222 {
00223 VecAXPBY(petscVec, alpha, gammaval, a.getPetscVector());
00224 VecAXPY(petscVec, beta, b.getPetscVector());
00225 return *this;
00226 }
00227
00228
00229 Teuchos::RCP<NOX::Abstract::Vector>
00230 Vector::clone(CopyType type) const
00231 {
00232 Teuchos::RCP<NOX::Abstract::Vector> newVec =
00233 Teuchos::rcp(new NOX::Petsc::Vector(petscVec, type));
00234 return newVec;
00235 }
00236
00237 double
00238 Vector::norm(Abstract::Vector::NormType type) const
00239 {
00240 double n;
00241 switch (type) {
00242 case MaxNorm:
00243 VecNorm(petscVec, NORM_INFINITY, &n);
00244 break;
00245 case OneNorm:
00246 VecNorm(petscVec, NORM_1, &n);
00247 break;
00248 case TwoNorm:
00249 default:
00250 VecNorm(petscVec, NORM_2, &n);
00251 break;
00252 }
00253 return n;
00254 }
00255
00256 double
00257 Vector::norm(const Abstract::Vector& weights) const
00258 {
00259 return norm(dynamic_cast<const Vector&>(weights));
00260 }
00261
00262 double
00263 Vector::norm(const Vector& weights) const
00264 {
00265 double n = 0.0;
00266 cerr << "Norm type not supported for weighted norm" << endl;
00267 throw;
00268 return n;
00269 }
00270
00271 double
00272 Vector::innerProduct(const Abstract::Vector& y) const
00273 {
00274 return innerProduct(dynamic_cast<const Vector&>(y));
00275 }
00276
00277 double
00278 Vector::innerProduct(const Vector& y) const
00279 {
00280 double dotprod;
00281 VecDot(y.getPetscVector(), petscVec, &dotprod);
00282 return dotprod;
00283 }
00284
00285 int
00286 Vector::length() const
00287 {
00288 int size;
00289 VecGetSize(petscVec, &size);
00290 return size;
00291 }