intern/iksolver/intern/TNT/vec.h

   1 /**
   2  */
   3
   4 /*
   5
   6 *
   7 * Template Numerical Toolkit (TNT): Linear Algebra Module
   8 *
   9 * Mathematical and Computational Sciences Division
  10 * National Institute of Technology,
  11 * Gaithersburg, MD USA
  12 *
  13 *
  14 * This software was developed at the National Institute of Standards and
  15 * Technology (NIST) by employees of the Federal Government in the course
  16 * of their official duties. Pursuant to title 17 Section 105 of the
  17 * United States Code, this software is not subject to copyright protection
  18 * and is in the public domain.  The Template Numerical Toolkit (TNT) is
  19 * an experimental system.  NIST assumes no responsibility whatsoever for
  20 * its use by other parties, and makes no guarantees, expressed or implied,
  21 * about its quality, reliability, or any other characteristic.
  22 *
  23 * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
  24 * see http://math.nist.gov/tnt for latest updates.
  25 *
  26 */
  27
  28
  29 // Basic TNT  numerical vector (0-based [i] AND 1-based (i) indexing )
  30 //
  31
  32 #ifndef VEC_H
  33 #define VEC_H
  34
  35 #include "subscript.h"
  36 #include <stdlib.h>
  37 #include <assert.h>
  38 #include <iostream>
  39
  40 namespace TNT
  41 {
  42
  43 template <class T>
  44 class Vector
  45 {
  46
  47
  48   public:
  49
  50     typedef Subscript   size_type;
  51     typedef         T   value_type;
  52     typedef         T   element_type;
  53     typedef         T*  pointer;
  54     typedef         T*  iterator;
  55     typedef         T&  reference;
  56     typedef const   T*  const_iterator;
  57     typedef const   T&  const_reference;
  58
  59     Subscript lbound() const { return 1;}
  60
  61   protected:
  62     T* v_;
  63     T* vm1_;        // pointer adjustment for optimzied 1-offset indexing
  64     Subscript n_;
  65
  66     // internal helper function to create the array
  67     // of row pointers
  68
  69     void initialize(Subscript N)
  70     {
  71         // adjust pointers so that they are 1-offset:
  72         // v_[] is the internal contiguous array, it is still 0-offset
  73         //
  74         assert(v_ == NULL);
  75         v_ = new T[N];
  76         assert(v_  != NULL);
  77         vm1_ = v_-1;
  78         n_ = N;
  79     }
  80
  81     void copy(const T*  v)
  82     {
  83         Subscript N = n_;
  84         Subscript i;
  85
  86 #ifdef TNT_UNROLL_LOOPS
  87         Subscript Nmod4 = N & 3;
  88         Subscript N4 = N - Nmod4;
  89
  90         for (i=0; i<N4; i+=4)
  91         {
  92             v_[i] = v[i];
  93             v_[i+1] = v[i+1];
  94             v_[i+2] = v[i+2];
  95             v_[i+3] = v[i+3];
  96         }
  97
  98         for (i=N4; i< N; i++)
  99             v_[i] = v[i];
 100 #else
 101
 102         for (i=0; i< N; i++)
 103             v_[i] = v[i];
 104 #endif
 105     }
 106
 107     void set(const T& val)
 108     {
 109         Subscript N = n_;
 110         Subscript i;
 111
 112 #ifdef TNT_UNROLL_LOOPS
 113         Subscript Nmod4 = N & 3;
 114         Subscript N4 = N - Nmod4;
 115
 116         for (i=0; i<N4; i+=4)
 117         {
 118             v_[i] = val;
 119             v_[i+1] = val;
 120             v_[i+2] = val;
 121             v_[i+3] = val;
 122         }
 123
 124         for (i=N4; i< N; i++)
 125             v_[i] = val;
 126 #else
 127
 128         for (i=0; i< N; i++)
 129             v_[i] = val;
 130
 131 #endif
 132     }
 133
 134
 135
 136     void destroy()
 137     {
 138         /* do nothing, if no memory has been previously allocated */
 139         if (v_ == NULL) return ;
 140
 141         /* if we are here, then matrix was previously allocated */
 142         delete [] (v_);
 143
 144         v_ = NULL;
 145         vm1_ = NULL;
 146     }
 147
 148
 149   public:
 150
 151     // access
 152
 153     iterator begin() { return v_;}
 154     iterator end()   { return v_ + n_; }
 155     const iterator begin() const { return v_;}
 156     const iterator end() const  { return v_ + n_; }
 157
 158     // destructor
 159
 160     ~Vector()
 161     {
 162         destroy();
 163     }
 164
 165     // constructors
 166
 167     Vector() : v_(0), vm1_(0), n_(0)  {}
 168
 169     Vector(const Vector<T> &A) : v_(0), vm1_(0), n_(0)
 170     {
 171         initialize(A.n_);
 172         copy(A.v_);
 173     }
 174
 175     Vector(Subscript N, const T& value = T()) :  v_(0), vm1_(0), n_(0)
 176     {
 177         initialize(N);
 178         set(value);
 179     }
 180
 181     Vector(Subscript N, const T* v) :  v_(0), vm1_(0), n_(0)
 182     {
 183         initialize(N);
 184         copy(v);
 185     }
 186
 187
 188     // methods
 189     //
 190     Vector<T>& newsize(Subscript N)
 191     {
 192         if (n_ == N) return *this;
 193
 194         destroy();
 195         initialize(N);
 196
 197         return *this;
 198     }
 199
 200
 201     // assignments
 202     //
 203     Vector<T>& operator=(const Vector<T> &A)
 204     {
 205         if (v_ == A.v_)
 206             return *this;
 207
 208         if (n_ == A.n_)         // no need to re-alloc
 209             copy(A.v_);
 210
 211         else
 212         {
 213             destroy();
 214             initialize(A.n_);
 215             copy(A.v_);
 216         }
 217
 218         return *this;
 219     }
 220
 221     Vector<T>& operator=(const T& scalar)
 222     {
 223         set(scalar);
 224         return *this;
 225     }
 226
 227     inline Subscript dim() const
 228     {
 229         return  n_;
 230     }
 231
 232     inline Subscript size() const
 233     {
 234         return  n_;
 235     }
 236
 237
 238     inline reference operator()(Subscript i)
 239     {
 240 #ifdef TNT_BOUNDS_CHECK
 241         assert(1<=i);
 242         assert(i <= n_) ;
 243 #endif
 244         return vm1_[i];
 245     }
 246
 247     inline const_reference operator() (Subscript i) const
 248     {
 249 #ifdef TNT_BOUNDS_CHECK
 250         assert(1<=i);
 251         assert(i <= n_) ;
 252 #endif
 253         return vm1_[i];
 254     }
 255
 256     inline reference operator[](Subscript i)
 257     {
 258 #ifdef TNT_BOUNDS_CHECK
 259         assert(0<=i);
 260         assert(i < n_) ;
 261 #endif
 262         return v_[i];
 263     }
 264
 265     inline const_reference operator[](Subscript i) const
 266     {
 267 #ifdef TNT_BOUNDS_CHECK
 268         assert(0<=i) ;
 269         assert(i < n_) ;
 270 #endif
 271         return v_[i];
 272     }
 273
 274
 275
 276 };
 277
 278
 279 /* ***************************  I/O  ********************************/
 280
 281 template <class T>
 282 std::ostream& operator<<(std::ostream &s, const Vector<T> &A)
 283 {
 284     Subscript N=A.dim();
 285
 286     s <<  N << std::endl;
 287
 288     for (Subscript i=0; i<N; i++)
 289         s   << A[i] << " " << std::endl;
 290     s << std::endl;
 291
 292     return s;
 293 }
 294
 295 template <class T>
 296 std::istream & operator>>(std::istream &s, Vector<T> &A)
 297 {
 298
 299     Subscript N;
 300
 301     s >> N;
 302
 303     if ( !(N == A.size() ))
 304     {
 305         A.newsize(N);
 306     }
 307
 308
 309     for (Subscript i=0; i<N; i++)
 310             s >>  A[i];
 311
 312
 313     return s;
 314 }
 315
 316 // *******************[ basic matrix algorithms ]***************************
 317
 318
 319 template <class T>
 320 Vector<T> operator+(const Vector<T> &A,
 321     const Vector<T> &B)
 322 {
 323     Subscript N = A.dim();
 324
 325     assert(N==B.dim());
 326
 327     Vector<T> tmp(N);
 328     Subscript i;
 329
 330     for (i=0; i<N; i++)
 331             tmp[i] = A[i] + B[i];
 332
 333     return tmp;
 334 }
 335
 336 template <class T>
 337 Vector<T> operator-(const Vector<T> &A,
 338     const Vector<T> &B)
 339 {
 340     Subscript N = A.dim();
 341
 342     assert(N==B.dim());
 343
 344     Vector<T> tmp(N);
 345     Subscript i;
 346
 347     for (i=0; i<N; i++)
 348             tmp[i] = A[i] - B[i];
 349
 350     return tmp;
 351 }
 352
 353 template <class T>
 354 Vector<T> operator*(const Vector<T> &A,
 355     const Vector<T> &B)
 356 {
 357     Subscript N = A.dim();
 358
 359     assert(N==B.dim());
 360
 361     Vector<T> tmp(N);
 362     Subscript i;
 363
 364     for (i=0; i<N; i++)
 365             tmp[i] = A[i] * B[i];
 366
 367     return tmp;
 368 }
 369
 370 template <class T>
 371 Vector<T> operator*(const Vector<T> &A,
 372     const T &B)
 373 {
 374     Subscript N = A.dim();
 375
 376     Vector<T> tmp(N);
 377     Subscript i;
 378
 379     for (i=0; i<N; i++)
 380             tmp[i] = A[i] * B;
 381
 382     return tmp;
 383 }
 384
 385
 386 template <class T>
 387 T dot_prod(const Vector<T> &A, const Vector<T> &B)
 388 {
 389     Subscript N = A.dim();
 390     assert(N == B.dim());
 391
 392     Subscript i;
 393     T sum = 0;
 394
 395     for (i=0; i<N; i++)
 396         sum += A[i] * B[i];
 397
 398     return sum;
 399 }
 400
 401 // inplace versions of the above template functions
 402
 403 // A = A + B
 404
 405 template <class T>
 406 void vectoradd(
 407         Vector<T> &A,
 408     const Vector<T> &B)
 409 {
 410     const Subscript N = A.dim();
 411     assert(N==B.dim());
 412     Subscript i;
 413
 414     for (i=0; i<N; i++)
 415             A[i] += B[i];
 416 }
 417
 418 // same with separate output vector
 419
 420 template <class T>
 421 void vectoradd(
 422         Vector<T> &C,
 423         const Vector<T> &A,
 424     const Vector<T> &B)
 425 {
 426     const Subscript N = A.dim();
 427     assert(N==B.dim());
 428     Subscript i;
 429
 430     for (i=0; i<N; i++)
 431             C[i] = A[i] + B[i];
 432 }
 433
 434 // A = A - B
 435
 436 template <class T>
 437 void vectorsub(
 438         Vector<T> &A,
 439     const Vector<T> &B)
 440 {
 441     const Subscript N = A.dim();
 442     assert(N==B.dim());
 443     Subscript i;
 444
 445     for (i=0; i<N; i++)
 446             A[i] -= B[i];
 447 }
 448
 449 template <class T>
 450 void vectorsub(
 451         Vector<T> &C,
 452         const Vector<T> &A,
 453     const Vector<T> &B)
 454 {
 455     const Subscript N = A.dim();
 456     assert(N==B.dim());
 457     Subscript i;
 458
 459     for (i=0; i<N; i++)
 460             C[i] = A[i] - B[i];
 461 }
 462
 463 template <class T>
 464 void vectorscale(
 465         Vector<T> &C,
 466         const Vector<T> &A,
 467     const T &B)
 468 {
 469     const Subscript N = A.dim();
 470     Subscript i;
 471
 472     for (i=0; i<N; i++)
 473             C[i] = A[i] * B;
 474 }
 475
 476 template <class T>
 477 void vectorscale(
 478         Vector<T> &A,
 479     const T &B)
 480 {
 481     const Subscript N = A.dim();
 482     Subscript i;
 483
 484     for (i=0; i<N; i++)
 485             A[i] *= B;
 486 }
 487
 488 }   /* namespace TNT */
 489
 490 #endif // VEC_H
 491