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

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

Made output from matrix and vector coherrent.

  • 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 418 2005-12-01 23:52:12Z 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::dec);
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.