source: branches/better_matrix_class/lib/gslapi/vector.cc @ 413

Last change on this file since 413 was 413, checked in by Jari Häkkinen, 17 years ago

All check passes, but more matrix checks needed.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.0 KB
Line 
1// $Id: vector.cc 413 2005-11-30 21:35:19Z jari $
2
3#include <c++_tools/gslapi/vector.h>
4#include <c++_tools/gslapi/matrix.h>
5#include <c++_tools/utility/stl_utility.h>
6
7#include <sstream>
8#include <vector>
9#include <utility>
10
11
12namespace theplu {
13namespace gslapi {
14
15
16
17  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
18    : const_view_(NULL)
19  {
20    // Jari, exception handling needed here. Failure in setting up a
21    // proper gsl_vector_view is communicated by NULL pointer in the
22    // view structure (cf. GSL manual). How about GSL error state?
23    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
24                                                                 stride,n));
25    v_ = &(view_->vector);
26  }
27
28
29
30  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
31    : view_(NULL)
32  {
33    // Jari, exception handling needed here. Failure in setting up a
34    // proper gsl_vector_view is communicated by NULL pointer in the
35    // view structure (cf. GSL manual). How about GSL error state?
36    const_view_ = new gsl_vector_const_view(
37                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
38    v_ = const_cast<gsl_vector*>(&(const_view_->vector));
39  }
40
41
42
43  vector::vector(matrix& m, size_t i, bool row)
44    : const_view_(NULL)
45  {
46    view_ = new gsl_vector_view(row ?
47                                gsl_matrix_row   (m.gsl_matrix_pointer(),i) :
48                                gsl_matrix_column(m.gsl_matrix_pointer(),i) );
49    v_ = &(view_->vector);
50  }
51
52
53
54  vector::vector(const matrix& m, size_t i, bool row)
55    : view_(NULL)
56  {
57    const_view_ = new gsl_vector_const_view(row ?
58                             gsl_matrix_const_row   (m.gsl_matrix_pointer(),i) :
59                             gsl_matrix_const_column(m.gsl_matrix_pointer(),i) );
60    v_ = const_cast<gsl_vector*>(&(const_view_->vector));
61  }
62
63
64
65  vector::vector(std::istream& is) throw (utility::IO_error,std::exception)
66    : view_(NULL), const_view_(NULL)
67  {
68    // read the data file and store in stl vectors (dynamically
69    // expandable)
70    std::vector<std::vector<double> > data_matrix;
71    u_int nof_columns=0;
72    u_int nof_rows=0;
73    std::vector<double> v;
74    for (nof_rows = 0; utility::read_to_double(is, v); nof_rows++) {
75      // Ignoring empty lines
76      if (!v.size()) {
77        nof_rows--;
78        continue;
79      }
80      if (!nof_columns)
81        nof_columns=v.size();
82      else if (nof_rows && (nof_columns>1)) {
83        std::ostringstream s;
84        s << "vector::vector(std::istream&) data file error:\n"
85          << "    File has inconsistent number of rows (" << nof_rows+1
86          << ") and columns (" << nof_columns
87          << ").\n    Expected a row or a column vector.";
88        throw utility::IO_error(s.str());
89      }
90      else if (v.size()!=nof_columns) {
91        std::ostringstream s;
92        s << "vector::vector(std::istream&) data file error:\n"
93          << "    Line " << (nof_rows+1) << " has " << v.size()
94          << " columns; expected " << nof_columns << " column.";
95        throw utility::IO_error(s.str());
96      }
97      data_matrix.push_back(v);
98    }
99
100    // manipulate the state of the stream to be good
101    is.clear(std::ios::goodbit);
102    // convert the data to a gsl vector
103    v_ = gsl_vector_alloc(nof_rows*nof_columns);
104    size_t n=0;
105    for (size_t i=0; i<nof_rows; i++)
106      for (size_t j=0; j<nof_columns; j++) 
107        gsl_vector_set( v_, n++, data_matrix[i][j] );
108  }
109
110
111
112  vector::~vector(void)
113  {
114    if (view_)
115      delete view_;
116    else if (const_view_)
117      delete const_view_;
118    else if (v_)
119      gsl_vector_free(v_);
120    v_=NULL;
121  }
122
123
124
125  gsl_vector* vector::gsl_vector_copy(void) const
126  {
127    gsl_vector* vec = gsl_vector_alloc(size());
128    gsl_vector_memcpy(vec,v_);
129    return vec;
130  }
131
132
133
134  std::pair<double,double> vector::minmax(void) const
135  {
136    double min, max;
137    gsl_vector_minmax(v_,&min,&max);
138    return std::pair<double,double>(min,max);
139  }
140
141
142
143  std::pair<size_t,size_t> vector::minmax_index(void) const
144  {
145    size_t min_index, max_index;
146    gsl_vector_minmax_index(v_,&min_index,&max_index);
147    return std::pair<size_t,size_t>(min_index,max_index);
148  }
149
150
151
152  double vector::sum(void) const
153  {
154    double sum = 0;
155    for (size_t i=0; i<size(); i++)
156      sum += gsl_vector_get( v_, i );
157    return( sum );
158  } 
159
160
161
162  double vector::operator*( const vector &other ) const
163  {
164    // Jari, check for gsl_support
165    double res = 0.0;;
166    for ( size_t i = 0; i < size(); ++i ) 
167      res += other[i] * (*this)[i];
168    return res;
169  }
170
171
172
173  vector vector::operator+(const vector &other) const
174  {
175    vector res(*this);
176    gsl_vector_add(res.v_,other.v_);
177    return res;
178  }
179
180
181
182  vector vector::operator-( const vector &other ) const
183  {
184    vector res( *this );
185    gsl_vector_sub(res.v_,other.v_);
186    return res;
187  }
188
189
190
191  bool vector::operator==(const vector& a) const
192  {
193    if (size()!=a.size())
194      return false;
195    for (size_t i=0; i<size(); ++i)
196      if (gsl_vector_get(v_,i)!=a(i))
197        return false;
198    return true;
199  }
200
201
202
203  const vector& vector::operator=( const vector& other )
204  {
205    if( this != &other ) {
206      if ( v_ )
207        gsl_vector_free( v_ );
208      v_ = other.gsl_vector_copy();
209    }
210    return *this;
211  } 
212
213
214
215  std::ostream& operator<<(std::ostream& s, const vector& a)
216  {
217    s.setf(std::ios::fixed);
218    s.precision(12);
219    for (size_t j = 0; j < a.size(); ++j) {
220      s << a[j];
221      if ( (j+1)<a.size() )
222        s << " ";
223    }
224
225    return s;
226  }
227
228
229}} // of namespace gslapi and namespace thep
Note: See TracBrowser for help on using the repository browser.