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_Common.H"
00043 #include "NOX_LAPACK_Vector.H"
00044 #include "Teuchos_BLAS_wrappers.hpp"
00045 #include "NOX_Random.H"
00046
00047 NOX::LAPACK::Vector::Vector() :
00048 n(0),
00049 x()
00050 {
00051 }
00052
00053 NOX::LAPACK::Vector::Vector(int N) :
00054 n(N),
00055 x(N, 0.0)
00056 {
00057 }
00058
00059 NOX::LAPACK::Vector::Vector(int N, double *v) :
00060 n(N),
00061 x(v,v+N)
00062 {
00063 }
00064
00065 NOX::LAPACK::Vector::Vector(const NOX::LAPACK::Vector& source,
00066 NOX::CopyType type) :
00067 n(source.n),
00068 x(source.x)
00069 {
00070 }
00071
00072 NOX::LAPACK::Vector::~Vector()
00073 {
00074 }
00075
00076 NOX::Abstract::Vector& NOX::LAPACK::Vector::operator=(
00077 const vector<double>& source)
00078 {
00079 x = source;
00080 return *this;
00081 }
00082
00083 NOX::Abstract::Vector& NOX::LAPACK::Vector::operator=(
00084 const NOX::Abstract::Vector& source)
00085 {
00086 return operator=(dynamic_cast<const NOX::LAPACK::Vector&>(source));
00087 }
00088
00089 NOX::Abstract::Vector& NOX::LAPACK::Vector::operator=(
00090 const NOX::LAPACK::Vector& source)
00091 {
00092 x = source.x;
00093 return *this;
00094 }
00095
00096 NOX::Abstract::Vector& NOX::LAPACK::Vector::init(double value)
00097 {
00098 for (int i = 0; i < n; i ++)
00099 x[i] = value;
00100 return *this;
00101 }
00102
00103 NOX::Abstract::Vector& NOX::LAPACK::Vector::random(bool useSeed, int seed)
00104 {
00105 if (useSeed)
00106 NOX::Random::setSeed(seed);
00107
00108 for (int i = 0; i < n; i ++)
00109 x[i] = NOX::Random::number();
00110
00111 return *this;
00112 }
00113
00114 NOX::Abstract::Vector& NOX::LAPACK::Vector::abs(
00115 const NOX::Abstract::Vector& base)
00116 {
00117 return abs(dynamic_cast<const NOX::LAPACK::Vector&>(base));
00118 }
00119
00120 NOX::Abstract::Vector& NOX::LAPACK::Vector::abs(
00121 const NOX::LAPACK::Vector& base)
00122 {
00123 for (int i = 0; i < n; i ++)
00124 x[i] = fabs(base[i]);
00125 return *this;
00126 }
00127
00128 NOX::Abstract::Vector& NOX::LAPACK::Vector::reciprocal(
00129 const NOX::Abstract::Vector& base)
00130 {
00131 return reciprocal(dynamic_cast<const NOX::LAPACK::Vector&>(base));
00132 }
00133
00134 NOX::Abstract::Vector& NOX::LAPACK::Vector::reciprocal(
00135 const NOX::LAPACK::Vector& base)
00136 {
00137 for (int i = 0; i < n; i ++)
00138 x[i] = 1.0 / base[i];
00139 return *this;
00140 }
00141
00142 NOX::Abstract::Vector& NOX::LAPACK::Vector::scale(double alpha)
00143 {
00144 for (int i = 0; i <n; i ++)
00145 x[i] *= alpha;
00146 return *this;
00147 }
00148
00149 NOX::Abstract::Vector& NOX::LAPACK::Vector::update(
00150 double alpha,
00151 const NOX::Abstract::Vector& a,
00152 double gamma)
00153 {
00154 return update(alpha, dynamic_cast<const NOX::LAPACK::Vector&>(a), gamma);
00155 }
00156
00157 NOX::Abstract::Vector& NOX::LAPACK::Vector::update(
00158 double alpha,
00159 const NOX::LAPACK::Vector& a,
00160 double gamma)
00161 {
00162 for (int i = 0; i < n; i ++)
00163 x[i] = alpha * a[i] + gamma * x[i];
00164 return *this;
00165 }
00166
00167 NOX::Abstract::Vector& NOX::LAPACK::Vector::update(
00168 double alpha,
00169 const NOX::Abstract::Vector& a,
00170 double beta,
00171 const NOX::Abstract::Vector& b,
00172 double gamma)
00173 {
00174 return update(alpha, dynamic_cast<const NOX::LAPACK::Vector&>(a),
00175 beta, dynamic_cast<const NOX::LAPACK::Vector&>(b), gamma);
00176 }
00177
00178 NOX::Abstract::Vector& NOX::LAPACK::Vector::update(
00179 double alpha,
00180 const NOX::LAPACK::Vector& a,
00181 double beta,
00182 const NOX::LAPACK::Vector& b,
00183 double gamma)
00184 {
00185 for (int i = 0; i < n; i ++)
00186 x[i] = alpha * a[i] + beta * b[i] + gamma * x[i];
00187 return *this;
00188 }
00189
00190 NOX::Abstract::Vector& NOX::LAPACK::Vector::scale(
00191 const NOX::Abstract::Vector& a)
00192 {
00193 return scale(dynamic_cast<const NOX::LAPACK::Vector&>(a));
00194 }
00195
00196 NOX::Abstract::Vector& NOX::LAPACK::Vector::scale(const NOX::LAPACK::Vector& a)
00197 {
00198 for (int i = 0; i < n; i ++)
00199 x[i] = a[i] * x[i];
00200 return *this;
00201 }
00202
00203 Teuchos::RCP<NOX::Abstract::Vector> NOX::LAPACK::Vector::
00204 clone(CopyType type) const
00205 {
00206 Teuchos::RCP<NOX::Abstract::Vector> tmp;
00207 tmp = Teuchos::rcp(new NOX::LAPACK::Vector(*this, type));
00208 return tmp;
00209 }
00210
00211 double NOX::LAPACK::Vector::norm(NOX::Abstract::Vector::NormType type) const
00212 {
00213
00214 if (n == 0)
00215 return 0.0;
00216
00217 int i;
00218 double value;
00219
00220 switch (type) {
00221 case MaxNorm:
00222 value = fabs(x[0]);
00223 for (i = 1; i < n; i ++)
00224 if (value < fabs(x[i]))
00225 value = fabs(x[i]);
00226 break;
00227 case OneNorm:
00228 value = DASUM_F77(&n, &x[0], &i_one);
00229 break;
00230 case TwoNorm:
00231 default:
00232 value = DNRM2_F77(&n, &x[0], &i_one);
00233 break;
00234 }
00235
00236 return value;
00237 }
00238
00239 double NOX::LAPACK::Vector::norm(const NOX::Abstract::Vector& weights) const
00240 {
00241 return norm(dynamic_cast<const NOX::LAPACK::Vector&>(weights));
00242 }
00243
00244 double NOX::LAPACK::Vector::norm(const NOX::LAPACK::Vector& weights) const
00245 {
00246 if (weights.length() != n) {
00247 cerr << "NOX::LAPACK::Vector::norm - size mismatch for weights vector" << endl;
00248 throw "NOX::LAPACK Error";
00249 }
00250
00251 double value = 0;
00252
00253 for (int i = 0; i < n; i ++)
00254 value += weights[i] * x[i] * x[i];
00255
00256 value = sqrt(value);
00257
00258 return value;
00259 }
00260
00261 double NOX::LAPACK::Vector::innerProduct(const NOX::Abstract::Vector& y) const
00262 {
00263 return innerProduct(dynamic_cast<const NOX::LAPACK::Vector&>(y));
00264 }
00265
00266 double NOX::LAPACK::Vector::innerProduct(const NOX::LAPACK::Vector& y) const
00267 {
00268 if (y.length() != n) {
00269 cerr << "NOX::LAPACK::Vector::innerProduct - size mismatch for y vector"
00270 << endl;
00271 throw "NOX::LAPACK Error";
00272 }
00273
00274 return DDOT_F77(&n, &x[0], &i_one, &y[0], &i_one);
00275 }
00276
00277 int NOX::LAPACK::Vector::length() const
00278 {
00279 return n;
00280 }
00281
00282 double& NOX::LAPACK::Vector::operator[] (int i)
00283 {
00284 return x[i];
00285 }
00286
00287 const double& NOX::LAPACK::Vector::operator[] (int i) const
00288 {
00289 return x[i];
00290 }
00291
00292 double& NOX::LAPACK::Vector::operator() (int i)
00293 {
00294 return x[i];
00295 }
00296
00297 const double& NOX::LAPACK::Vector::operator() (int i) const
00298 {
00299 return x[i];
00300 }
00301
00302 ostream& NOX::LAPACK::Vector::leftshift(ostream& stream) const
00303 {
00304 stream << "[ ";
00305 for (int i = 0; i < n; i ++)
00306 stream << x[i] << " ";
00307 stream << "]";
00308 return stream;
00309 }
00310
00311 ostream& operator<<(ostream& stream, const NOX::LAPACK::Vector& v)
00312 {
00313 return v.leftshift(stream);
00314 }
00315
00316 void NOX::LAPACK::Vector::print(std::ostream& stream) const
00317 {
00318 stream << *this << endl;
00319 }
00320