MFEM v2.0
vector.cpp
Go to the documentation of this file.
00001 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
00002 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
00003 // reserved. See file COPYRIGHT for details.
00004 //
00005 // This file is part of the MFEM library. For more information and source code
00006 // availability see http://mfem.googlecode.com.
00007 //
00008 // MFEM is free software; you can redistribute it and/or modify it under the
00009 // terms of the GNU Lesser General Public License (as published by the Free
00010 // Software Foundation) version 2.1 dated February 1999.
00011 
00012 // Implementation of data type vector
00013 
00014 #include <iostream>
00015 #include <iomanip>
00016 #include <math.h>
00017 #include <stdlib.h>
00018 #include <time.h>
00019 
00020 #include "vector.hpp"
00021 
00022 Vector::Vector(const Vector &v)
00023 {
00024    int s = v.Size();
00025    if (s > 0)
00026    {
00027       allocsize = size = s;
00028       data = new double[s];
00029       for (int i = 0; i < s; i++)
00030          data[i] = v(i);
00031    }
00032    else
00033    {
00034       allocsize = size = 0;
00035       data = NULL;
00036    }
00037 }
00038 
00039 void Vector::Load(istream **in, int np, int *dim)
00040 {
00041    int i, j, s;
00042 
00043    s = 0;
00044    for (i = 0; i < np; i++)
00045       s += dim[i];
00046 
00047    SetSize(s);
00048 
00049    int p = 0;
00050    for (i = 0; i < np; i++)
00051       for (j = 0; j < dim[i]; j++)
00052          *in[i] >> data[p++];
00053 }
00054 
00055 void Vector::Load(istream &in, int Size)
00056 {
00057    SetSize(Size);
00058 
00059    for (int i = 0; i < size; i++)
00060       in >> data[i];
00061 }
00062 
00063 double &Vector::Elem(int i)
00064 {
00065    return operator()(i);
00066 }
00067 
00068 const double &Vector::Elem(int i) const
00069 {
00070    return operator()(i);
00071 }
00072 
00073 double Vector::operator*(const double *v) const
00074 {
00075    int s = size;
00076    const double *d = data;
00077    double prod = 0.0;
00078 #ifdef MFEM_USE_OPENMP
00079 #pragma omp parallel for reduction(+:prod)
00080 #endif
00081    for (int i = 0; i < s; i++)
00082       prod += d[i] * v[i];
00083    return prod;
00084 }
00085 
00086 double Vector::operator*(const Vector &v) const
00087 {
00088 #ifdef MFEM_DEBUG
00089    if (v.size != size)
00090       mfem_error("Vector::operator*(const Vector &) const");
00091 #endif
00092 
00093    return operator*(v.data);
00094 }
00095 
00096 Vector &Vector::operator=(const Vector &v)
00097 {
00098    SetSize(v.Size());
00099    for (int i = 0; i < size; i++)
00100       data[i] = v.data[i];
00101    return *this;
00102 }
00103 
00104 Vector &Vector::operator=(double value)
00105 {
00106    register int i, s = size;
00107    register double *p = data, v = value;
00108    for (i = 0; i < s; i++)
00109       *(p++) = v;
00110    return *this;
00111 }
00112 
00113 Vector &Vector::operator*=(double c)
00114 {
00115    for (int i = 0; i < size; i++)
00116       data[i] *= c;
00117    return *this;
00118 }
00119 
00120 Vector &Vector::operator/=(double c)
00121 {
00122    double m = 1.0/c;
00123    for (int i = 0; i < size; i++)
00124       data[i] *= m;
00125    return *this;
00126 }
00127 
00128 Vector &Vector::operator-=(double c)
00129 {
00130    for (int i = 0; i < size; i++)
00131       data[i] -= c;
00132    return *this;
00133 }
00134 
00135 Vector &Vector::operator-=(const Vector &v)
00136 {
00137 #ifdef MFEM_DEBUG
00138    if (size != v.size)
00139       mfem_error("Vector::operator-=(const Vector &)");
00140 #endif
00141    for (int i = 0; i < size; i++)
00142       data[i] -= v(i);
00143    return *this;
00144 }
00145 
00146 Vector &Vector::operator+=(const Vector &v)
00147 {
00148 #ifdef MFEM_DEBUG
00149    if (size != v.size)
00150       mfem_error("Vector::operator+=(const Vector &)");
00151 #endif
00152    for (int i = 0; i < size; i++)
00153       data[i] += v(i);
00154    return *this;
00155 }
00156 
00157 Vector &Vector::Add(const double a, const Vector &Va)
00158 {
00159 #ifdef MFEM_DEBUG
00160    if (size != Va.size)
00161       mfem_error("Vector::Add(const double, const Vector &)");
00162 #endif
00163    if (a != 0.0)
00164    {
00165       for (int i = 0; i < size; i++)
00166          data[i] += a * Va(i);
00167    }
00168    return *this;
00169 }
00170 
00171 Vector &Vector::Set(const double a, const Vector &Va)
00172 {
00173 #ifdef MFEM_DEBUG
00174    if (size != Va.size)
00175       mfem_error("Vector::Set(const double, const Vector &)");
00176 #endif
00177    for (int i = 0; i < size; i++)
00178       data[i] = a * Va(i);
00179    return *this;
00180 }
00181 
00182 void Vector::SetVector(const Vector &v, int offset)
00183 {
00184    int vs = v.Size();
00185    double *vp = v.data, *p = data + offset;
00186 
00187 #ifdef MFEM_DEBUG
00188    if (offset+vs > size)
00189       mfem_error("Vector::SetVector(const Vector &, int)");
00190 #endif
00191 
00192    for (int i = 0; i < vs; i++)
00193       p[i] = vp[i];
00194 }
00195 
00196 void Vector::Neg()
00197 {
00198    for (int i = 0; i < size; i++)
00199       data[i] = -data[i];
00200 }
00201 
00202 void swap(Vector *v1, Vector *v2)
00203 {
00204    int size = v1->size, allocsize = v1->allocsize;
00205    double *data = v1->data;
00206 
00207    v1->size      = v2->size;
00208    v1->allocsize = v2->allocsize;
00209    v1->data      = v2->data;
00210 
00211    v2->size      = size;
00212    v2->allocsize = allocsize;
00213    v2->data      = data;
00214 }
00215 
00216 void add(const Vector &v1, const Vector &v2, Vector &v)
00217 {
00218 #ifdef MFEM_DEBUG
00219    if (v.size != v1.size || v.size != v2.size)
00220       mfem_error("add(Vector &v1, Vector &v2, Vector &v)");
00221 #endif
00222 
00223 #ifdef MFEM_USE_OPENMP
00224 #pragma omp parallel for
00225 #endif
00226    for (int i = 0; i < v.size; i++)
00227       v.data[i] = v1.data[i] + v2.data[i];
00228 }
00229 
00230 void add(const Vector &v1, double alpha, const Vector &v2, Vector &v)
00231 {
00232 #ifdef MFEM_DEBUG
00233    if (v.size != v1.size || v.size != v2.size)
00234       mfem_error ("add(Vector &v1, double alpha, Vector &v2, Vector &v)");
00235 #endif
00236    if (alpha == 0.0)
00237    {
00238       v = v1;
00239    }
00240    else if (alpha == 1.0)
00241    {
00242       add(v1, v2, v);
00243    }
00244    else
00245    {
00246       const double *v1p = v1.data, *v2p = v2.data;
00247       double *vp = v.data;
00248       int s = v.size;
00249 #ifdef MFEM_USE_OPENMP
00250 #pragma omp parallel for
00251 #endif
00252       for (int i = 0; i < s; i++)
00253          vp[i] = v1p[i] + alpha*v2p[i];
00254    }
00255 }
00256 
00257 void add(const double a, const Vector &x, const Vector &y, Vector &z)
00258 {
00259 #ifdef MFEM_DEBUG
00260    if (x.size != y.size || x.size != z.size)
00261       mfem_error ("add(const double a, const Vector &x, const Vector &y,"
00262                   " Vector &z)");
00263 #endif
00264    if (a == 0.0)
00265    {
00266       z = 0.0;
00267    }
00268    else if (a == 1.0)
00269    {
00270       add(x, y, z);
00271    }
00272    else
00273    {
00274       const double *xp = x.data;
00275       const double *yp = y.data;
00276       double       *zp = z.data;
00277       int            s = x.size;
00278 
00279 #ifdef MFEM_USE_OPENMP
00280 #pragma omp parallel for
00281 #endif
00282       for (int i = 0; i < s; i++)
00283          zp[i] = a * (xp[i] + yp[i]);
00284    }
00285 }
00286 
00287 void add(const double a, const Vector &x,
00288          const double b, const Vector &y, Vector &z)
00289 {
00290 #ifdef MFEM_DEBUG
00291    if (x.size != y.size || x.size != z.size)
00292       mfem_error("add(const double a, const Vector &x,\n"
00293                  "    const double b, const Vector &y, Vector &z)");
00294 #endif
00295    if (a == 0.0)
00296    {
00297       z.Set(b, y);
00298    }
00299    else if (b == 0.0)
00300    {
00301       z.Set(a, x);
00302    }
00303    else if (a == 1.0)
00304    {
00305       add(x, b, y, z);
00306    }
00307    else if (b == 1.0)
00308    {
00309       add(y, a, x, z);
00310    }
00311    else if (a == b)
00312    {
00313       add(a, x, y, z);
00314    }
00315    else
00316    {
00317       const double *xp = x.data;
00318       const double *yp = y.data;
00319       double       *zp = z.data;
00320       int            s = x.size;
00321 
00322 #ifdef MFEM_USE_OPENMP
00323 #pragma omp parallel for
00324 #endif
00325       for (int i = 0; i < s; i++)
00326          zp[i] = a * xp[i] + b * yp[i];
00327    }
00328 }
00329 
00330 void subtract(const Vector &x, const Vector &y, Vector &z)
00331 {
00332 #ifdef MFEM_DEBUG
00333    if (x.size != y.size || x.size != z.size)
00334       mfem_error ("subtract(const Vector &, const Vector &, Vector &)");
00335 #endif
00336    const double *xp = x.data;
00337    const double *yp = y.data;
00338    double       *zp = z.data;
00339    int            s = x.size;
00340 
00341 #ifdef MFEM_USE_OPENMP
00342 #pragma omp parallel for
00343 #endif
00344    for (int i = 0; i < s; i++)
00345       zp[i] = xp[i] - yp[i];
00346 }
00347 
00348 void subtract(const double a, const Vector &x, const Vector &y, Vector &z)
00349 {
00350 #ifdef MFEM_DEBUG
00351    if (x.size != y.size || x.size != z.size)
00352       mfem_error("subtract(const double a, const Vector &x,"
00353                  " const Vector &y, Vector &z)");
00354 #endif
00355 
00356    if (a == 0.)
00357    {
00358       z = 0.;
00359    }
00360    else if (a == 1.)
00361    {
00362       subtract(x, y, z);
00363    }
00364    else
00365    {
00366       const double *xp = x.data;
00367       const double *yp = y.data;
00368       double       *zp = z.data;
00369       int            s = x.size;
00370 
00371 #ifdef MFEM_USE_OPENMP
00372 #pragma omp parallel for
00373 #endif
00374       for (int i = 0; i < s; i++)
00375          zp[i] = a * (xp[i] - yp[i]);
00376    }
00377 }
00378 
00379 void Vector::GetSubVector(const Array<int> &dofs, Vector &elemvect) const
00380 {
00381    int i, j, n = dofs.Size();
00382 
00383    elemvect.SetSize (n);
00384 
00385    for (i = 0; i < n; i++)
00386       if ((j=dofs[i]) >= 0)
00387          elemvect(i) = data[j];
00388       else
00389          elemvect(i) = -data[-1-j];
00390 }
00391 
00392 void Vector::GetSubVector(const Array<int> &dofs, double *elem_data) const
00393 {
00394    int i, j, n = dofs.Size();
00395 
00396    for (i = 0; i < n; i++)
00397       if ((j=dofs[i]) >= 0)
00398          elem_data[i] = data[j];
00399       else
00400          elem_data[i] = -data[-1-j];
00401 }
00402 
00403 void Vector::SetSubVector(const Array<int> &dofs, const Vector &elemvect)
00404 {
00405    int i, j, n = dofs.Size();
00406 
00407    for (i = 0; i < n; i++)
00408       if ((j=dofs[i]) >= 0)
00409          data[j] = elemvect(i);
00410       else
00411          data[-1-j] = -elemvect(i);
00412 }
00413 
00414 void Vector::SetSubVector(const Array<int> &dofs, double *elem_data)
00415 {
00416    int i, j, n = dofs.Size();
00417 
00418    for (i = 0; i < n; i++)
00419       if ((j=dofs[i]) >= 0)
00420          data[j] = elem_data[i];
00421       else
00422          data[-1-j] = -elem_data[i];
00423 }
00424 
00425 void Vector::AddElementVector(const Array<int> &dofs, const Vector &elemvect)
00426 {
00427    int i, j, n = dofs.Size();
00428 
00429    for (i = 0; i < n; i++)
00430       if ((j=dofs[i]) >= 0)
00431          data[j] += elemvect(i);
00432       else
00433          data[-1-j] -= elemvect(i);
00434 }
00435 
00436 void Vector::AddElementVector(const Array<int> &dofs, double *elem_data)
00437 {
00438    int i, j, n = dofs.Size();
00439 
00440    for (i = 0; i < n; i++)
00441       if ((j=dofs[i]) >= 0)
00442          data[j] += elem_data[i];
00443       else
00444          data[-1-j] -= elem_data[i];
00445 }
00446 
00447 void Vector::AddElementVector(const Array<int> &dofs, const double a,
00448                               const Vector &elemvect)
00449 {
00450    int i, j, n = dofs.Size();
00451 
00452    for (i = 0; i < n; i++)
00453       if ((j=dofs[i]) >= 0)
00454          data[j] += a * elemvect(i);
00455       else
00456          data[-1-j] -= a * elemvect(i);
00457 }
00458 
00459 void Vector::Print(ostream &out, int width) const
00460 {
00461    for (int i = 0; 1; )
00462    {
00463       out << data[i];
00464       i++;
00465       if (i == size)
00466          break;
00467       if ( i % width == 0 )
00468          out << '\n';
00469       else
00470          out << ' ';
00471    }
00472    out << '\n';
00473 }
00474 
00475 void Vector::Print_HYPRE(ostream &out) const
00476 {
00477    int i;
00478    ios::fmtflags old_fmt = out.setf(ios::scientific);
00479    int old_prec = out.precision(14);
00480 
00481    out << size << '\n';  // number of rows
00482 
00483    for (i = 0; i < size; i++)
00484       out << data[i] << '\n';
00485 
00486    out.precision(old_prec);
00487    out.setf(old_fmt);
00488 }
00489 
00490 void Vector::Randomize(int seed)
00491 {
00492    // static unsigned int seed = time(0);
00493    const double max = (double)(RAND_MAX) + 1.;
00494 
00495    if (seed == 0)
00496       seed = time(0);
00497 
00498    // srand(seed++);
00499    srand(seed);
00500 
00501    for (int i = 0; i < size; i++)
00502       data[i] = fabs(rand()/max);
00503 }
00504 
00505 double Vector::Norml2()
00506 {
00507    return sqrt((*this)*(*this));
00508 }
00509 
00510 double Vector::Normlinf()
00511 {
00512    double max = fabs(data[0]);
00513 
00514    for (int i = 1; i < size; i++)
00515       if (fabs(data[i]) > max)
00516          max = fabs(data[i]);
00517 
00518    return max;
00519 }
00520 
00521 double Vector::Norml1()
00522 {
00523    double sum = 0.0;
00524 
00525    for (int i = 0; i < size; i++)
00526       sum += fabs (data[i]);
00527 
00528    return sum;
00529 }
00530 
00531 double Vector::Max()
00532 {
00533    double max = data[0];
00534 
00535    for (int i = 1; i < size; i++)
00536       if (data[i] > max)
00537          max = data[i];
00538 
00539    return max;
00540 }
00541 
00542 double Vector::Min()
00543 {
00544    double min = data[0];
00545 
00546    for (int i = 1; i < size; i++)
00547       if (data[i] < min)
00548          min = data[i];
00549 
00550    return min;
00551 }
00552 
00553 double Vector::DistanceTo(const double *p) const
00554 {
00555    return Distance(data, p, size);
00556 }
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines