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

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

Removed some gslapi copy intensive operators.
Made gslapi assignment operators ignore views, and change them
into normal vectors if assigned.
Cleaning up of code aswell.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.8 KB
Line 
1// $Id: vector.cc 415 2005-12-01 15:52:36Z 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_p(),i) :
48                              gsl_matrix_column(m.gsl_matrix_p(),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_p(),i) :
59                                   gsl_matrix_const_column(m.gsl_matrix_p(),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::create_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  bool vector::operator==(const vector& a) const
174  {
175    if (size()!=a.size())
176      return false;
177    for (size_t i=0; i<size(); ++i)
178      if (gsl_vector_get(v_,i)!=a(i))
179        return false;
180    return true;
181  }
182
183
184
185  const vector& vector::operator=( const vector& other )
186  {
187    if( this != &other ) {
188      if (view_)
189        delete view_;
190      else if (const_view_)
191        delete const_view_;
192      else if ( v_ )
193        gsl_vector_free( v_ );
194      v_ = other.create_gsl_vector_copy();
195    }
196    return *this;
197  } 
198
199
200
201  std::ostream& operator<<(std::ostream& s, const vector& a)
202  {
203    s.setf(std::ios::fixed);
204    s.precision(12);
205    for (size_t j = 0; j < a.size(); ++j) {
206      s << a[j];
207      if ( (j+1)<a.size() )
208        s << " ";
209    }
210
211    return s;
212  }
213
214
215}} // of namespace gslapi and namespace thep
Note: See TracBrowser for help on using the repository browser.