source: trunk/yat/utility/Vector.cc @ 1136

Last change on this file since 1136 was 1136, checked in by Peter, 14 years ago

avoid reallocation when sizes match in Vector assignment

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.2 KB
Line 
1// $Id: Vector.cc 1136 2008-02-23 23:46:31Z peter $
2
3/*
4  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
5  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
6  Copyright (C) 2005, 2006 Jari Häkkinen, Markus Ringnér, Peter Johansson
7  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
8
9  This file is part of the yat library, http://trac.thep.lu.se/yat
10
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 2 of the
14  License, or (at your option) any later version.
15
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with this program; if not, write to the Free Software
23  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24  02111-1307, USA.
25*/
26
27#include "Vector.h"
28#include "utility.h"
29#include "yat/random/random.h"
30
31#include <algorithm>
32#include <cassert>
33#include <cmath>
34#include <iostream>
35#include <sstream>
36#include <utility>
37#include <vector>
38
39namespace theplu {
40namespace yat {
41namespace utility {
42
43
44  Vector::Vector(void)
45    : VectorMutable()
46  {
47  }
48
49
50  Vector::Vector(size_t n, double init_value)
51    : VectorMutable(gsl_vector_alloc(n))
52  {
53    if (!vec_)
54      throw utility::GSL_error("Vector::Vector failed to allocate memory");
55    assert(const_vec_);
56    all(init_value);
57  }
58
59
60  Vector::Vector(const Vector& other)
61    : VectorMutable(create_gsl_vector_copy(other))
62  {
63  }
64
65
66  Vector::Vector(const VectorBase& other)
67    : VectorMutable(create_gsl_vector_copy(other))
68  {
69  }
70
71
72  Vector::Vector(std::istream& is, char sep)
73    throw (utility::IO_error, std::exception)
74    : VectorMutable()
75  {
76    // read the data file and store in stl vectors (dynamically
77    // expandable)
78    std::vector<std::vector<double> > data_matrix;
79    u_int nof_columns=0;
80    u_int nof_rows=0;
81    std::string line;
82    while(getline(is, line, '\n')) {
83      // Empty lines
84      if (!line.size())
85          continue;
86      nof_rows++;
87
88      std::vector<double> v;
89      std::string element;
90      std::stringstream ss(line);
91      bool ok=true;
92      while(ok) {
93        if(sep=='\0')
94          ok=(ss>>element);
95        else 
96          ok=getline(ss, element, sep);
97        if(!ok)
98          break;       
99
100        if(utility::is_double(element)) {
101          v.push_back(atof(element.c_str()));
102        }
103        else if (!element.size() || utility::is_nan(element)) {
104          v.push_back(std::numeric_limits<double>::quiet_NaN());
105        }
106        else {
107          std::stringstream ss("Warning: '");
108          ss << element << "' is not a double.";
109          throw IO_error(ss.str());
110        }
111      } 
112      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
113          v.push_back(std::numeric_limits<double>::quiet_NaN());
114      if (!nof_columns)
115        nof_columns=v.size();
116      else if (nof_rows && (nof_columns>1)) {
117        std::ostringstream s;
118        s << "Vector::Vector(std::istream&) data file error:\n"
119          << "    File has inconsistent number of rows (" << nof_rows
120          << ") and columns (" << nof_columns
121          << ").\n    Expected a row or a column Vector.";
122        throw utility::IO_error(s.str());
123      }
124      else if (v.size()!=nof_columns) {
125        std::ostringstream s;
126        s << "Vector::Vector(std::istream&) data file error:\n"
127          << "    Line " << nof_rows << " has " << v.size()
128          << " columns; expected " << nof_columns << " column.";
129        throw utility::IO_error(s.str());
130      }
131      data_matrix.push_back(v);
132    }
133
134    // manipulate the state of the stream to be good
135    is.clear(std::ios::goodbit);
136    // convert the data to a gsl vector
137    vec_ = gsl_vector_alloc(nof_rows*nof_columns);
138    if (!vec_)
139      throw utility::GSL_error("Vector::Vector failed to allocate memory");
140    size_t n=0;
141    // if gsl error handler disabled, out of bounds index will not
142    // abort the program.
143    for (size_t i=0; i<nof_rows; i++)
144      for (size_t j=0; j<nof_columns; j++) 
145        gsl_vector_set( vec_, n++, data_matrix[i][j] );
146    const_vec_ = vec_;
147  }
148
149
150  Vector::~Vector(void)
151  {
152    delete_allocated_memory();
153  }
154
155
156  const Vector& Vector::assign(const VectorBase& other)
157  {
158    // avoid self assignment
159    if (this == &other)
160      return *this;
161    if (!other.size()){
162      delete_allocated_memory();
163      return *this;
164    }
165    // indirect self assignment if begin <= other.begin() < end
166    if (begin() <= other.begin() && other.begin()<end()) {
167      if (size() != other.size()) {
168        gsl_vector* tmp = create_gsl_vector_copy(other);
169        delete_allocated_memory();
170        const_vec_ = vec_ = tmp;
171      }
172      return *this;
173    } 
174    if (size()==other.size()){
175      if (gsl_vector_memcpy(vec_, other.gsl_vector_p()))
176        throw utility::GSL_error("Vector::assign memcpy failed");
177    }
178    else {
179      delete_allocated_memory();
180      vec_ = create_gsl_vector_copy(other);
181    }
182    const_vec_ = vec_;
183    return *this;
184  } 
185
186
187  gsl_vector* Vector::create_gsl_vector_copy(const VectorBase& other) const
188  {
189    gsl_vector* vec = gsl_vector_alloc(other.size());
190    if (!vec)
191      throw utility::GSL_error("Vector::create_gsl_vector_copy failed to allocate memory");
192    if (gsl_vector_memcpy(vec, other.gsl_vector_p()))
193      throw utility::GSL_error("Vector::create_gsl_vector_copy memcpy failed");
194    return vec;
195  }
196
197
198  void Vector::delete_allocated_memory(void) 
199  {
200    if (vec_)
201      gsl_vector_free(vec_);
202    const_vec_ = vec_ = NULL;
203  }
204
205
206  bool Vector::isview(void) const
207  {
208    return false;
209  }
210
211
212  void Vector::resize(size_t n, double init_value)
213  {
214    delete_allocated_memory();
215    const_vec_ = vec_ = gsl_vector_alloc(n);
216    if (!vec_)
217      throw utility::GSL_error("Vector::resize failed to allocate memory");
218    all(init_value);
219  }
220
221
222  void set_basis(Vector& v, size_t i)
223  {
224    assert(v.gsl_vector_p());
225    gsl_vector_set_basis(v.gsl_vector_p(),i);
226  }
227
228
229  void sort(Vector& v)
230  {
231    assert(v.gsl_vector_p());
232    std::sort(v.begin(), v.end());
233  }
234
235
236  void swap(Vector& v, Vector& w)
237  {
238    assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
239    int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
240    if (status)
241      throw utility::GSL_error(std::string("swap(Vector&,Vector&)",status));
242  }
243
244
245  const Vector& Vector::operator=( const VectorBase& other )
246  {
247    return assign(other);
248  } 
249
250
251  const Vector& Vector::operator=( const Vector& other )
252  {
253    return assign(other);
254  }
255
256
257}}} // of namespace utility, yat, and thep
Note: See TracBrowser for help on using the repository browser.