Changeset 42 for trunk/src/vector.cc


Ignore:
Timestamp:
Feb 26, 2004, 4:06:20 PM (19 years ago)
Author:
Jari Häkkinen
Message:

Made a major revision of matrix and vector classes. Everything compiles
but the binaries have not been tested.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/vector.cc

    r37 r42  
    11// $Id$
    22
    3 // #include <iostream>
     3#include <iostream>
     4#include <sstream>
    45#include <string>
    5 #include <sstream>
    66#include <vector>
     7
    78#include "vector.h"
    89
    910
    10 using namespace thep_gsl_api;
    11 
    12 
    13 // Constructors and Destructors
    14 ///////////////////////////////
    15 vector::vector() : v_( NULL )
    16 {
    17   is_col_ = true;
    18 }
    19 
    20 
    21 vector::vector( const size_t& N, bool is_col_vector,
    22     bool init_to_zero ) :
    23   is_col_( is_col_vector )
    24 {
    25   if( init_to_zero )
    26     {
    27      v_ = gsl_vector_calloc( N );
     11namespace theplu {
     12namespace gslapi {
     13
     14
     15
     16  vector::vector(void)
     17    : v_(NULL)
     18  {
     19  }
     20
     21
     22
     23  vector::vector(const size_t n, const double init_value)
     24  {
     25    v_ = gsl_vector_alloc(n);
     26    set_all(init_value);
     27  }
     28
     29
     30
     31  vector::vector(const vector& other)
     32  {
     33    gsl_vector_free(v_);
     34    v_ = other.gsl_vector_copy();
     35  }
     36
     37
     38
     39  vector::vector(gsl_vector* vec)
     40    : v_(vec)
     41  {
     42  }
     43
     44
     45
     46  vector::vector(std::istream& is)
     47  {
     48    using namespace std;
     49 
     50    // read the data file and store in stl vectors (dynamically expandable)
     51    std::vector<std::vector<double> > data_matrix;
     52    u_int nof_columns=0;
     53    u_int nof_rows = 0;
     54    string s;
     55    for (nof_rows = 0; getline(is, s, '\n'); nof_rows++) {
     56      istringstream line(s);
     57      std::vector<double> v;
     58      string tmp_string;
     59      while (line >> tmp_string) {
     60        if(!is.good()) {
     61          // Jari. change to throw exception, remove unnecessary includes
     62          cerr << "vector::vector(std::istream&): "
     63               << "error reading data file!" << endl;
     64          exit(1);
     65        }
     66        double t=atof(tmp_string.c_str());
     67        // Here we should check that it actually was a double!!!!!
     68        v.push_back(t);
     69      }
     70      if(nof_columns==0)
     71        nof_columns=v.size();
     72      else if(v.size()!=nof_columns) {
     73        // Jari. change to throw exception, remove unnecessary includes
     74        cerr << "vector::vector(std::istream&) data file error: "
     75             << "line" << nof_rows+1 << " has " << v.size()
     76             << " columns; expected " << nof_columns
     77             << " columns"
     78             << endl;
     79        exit(1);
     80      }
     81      data_matrix.push_back(v);
     82    }
     83 
     84    // manipulate the state of the stream to be good
     85    is.clear(std::ios::goodbit);
     86 
     87    // convert the data to a gsl vector and check that data file is a
     88    // column vector or a row vector
     89    if(nof_columns==1) {
     90      v_ = gsl_vector_alloc ( nof_rows );
     91      for(u_int i=0;i<nof_rows;i++)
     92        gsl_vector_set( v_, i, data_matrix[i][0] );
     93    }
     94    else if(nof_rows==1){
     95      v_ = gsl_vector_alloc ( nof_columns );
     96      for(u_int i=0;i<nof_columns;i++)
     97        gsl_vector_set( v_, i, data_matrix[0][i] );
     98    }
     99    else {
     100      // Jari. change to throw exception, remove unnecessary includes
     101      cerr << "vector::vector(std::istream&) data file error: "
     102           << "file has " << nof_rows << " rows and " << nof_columns
     103           << " columns; expected a row vector or a column vector"
     104           << endl;
     105      exit(1);
     106    }
     107  }
     108
     109
     110
     111  vector::~vector(void)
     112  {
     113    if (v_) {
     114      gsl_vector_free(v_);
     115      v_=NULL;
    28116    }
    29   else     
    30     {
    31      v_ = gsl_vector_alloc ( N );
     117  }
     118
     119
     120
     121  gsl_vector* vector::gsl_vector_copy(void) const
     122  {
     123    gsl_vector* vec = gsl_vector_alloc(size());
     124    gsl_vector_memcpy(vec,v_);
     125    return vec;
     126  }
     127
     128
     129
     130  vector vector::mul_elements( const vector& other ) const
     131  {
     132    vector res( *this );
     133    gsl_vector_mul( res.v_, other.v_ );
     134    return res;
     135  }
     136
     137
     138
     139  double vector::sum(void) const
     140  {
     141    double sum = 0;
     142    for (u_int i=0; i<size(); i++)
     143      sum += gsl_vector_get( v_, i );
     144    return( sum );
     145  } 
     146
     147
     148
     149  vector vector::operator-(void) const
     150  {
     151    vector res( *this );
     152    gsl_vector_scale (res.v_,-1.);
     153    return res;
     154  }
     155
     156
     157
     158  double vector::operator*( const vector &other ) const
     159  {
     160    // Jari, check for gsl_support
     161    double res = 0.0;;
     162    for ( size_t i = 0; i < size(); ++i )
     163      res += other[i] * (*this)[i];
     164    return res;
     165  }
     166
     167
     168
     169  vector vector::operator+(const vector &other) const
     170  {
     171    vector res(*this);
     172    gsl_vector_add(res.v_,other.v_);
     173    return res;
     174  }
     175
     176
     177
     178  vector vector::operator-( const vector &other ) const
     179  {
     180    vector res( *this );
     181    gsl_vector_sub(res.v_,other.v_);
     182    return res;
     183  }
     184
     185
     186
     187  vector& vector::operator=( const vector& other )
     188  {
     189    if( this != &other ) {
     190      if ( v_ )
     191        gsl_vector_free( v_ );
     192      v_ = other.gsl_vector_copy();
    32193    }
    33 }
    34 
    35 
    36 // Is this the way to do it? No copy ...
    37 // internal data ...
    38 vector::vector( gsl_vector* v, bool is_col ) : v_( v ), is_col_( is_col )
    39 {
    40 }
    41 
    42 
    43 // Copy constructor
    44 vector::vector( const vector& other )
    45 {
    46   v_ = new_copy( other.get_gsl_vector() );
    47   is_col_ = other.is_column();
    48 }
    49 
    50 
    51 // Constructor that gets data from istream
    52 vector::vector(std::istream& is)
    53 {
    54   using namespace std;
    55  
    56   // read the data file and store in stl vectors (dynamically expandable)
    57   std::vector<std::vector<double> > data_matrix;
    58   u_int nof_columns=0;
    59   u_int nof_rows = 0;
    60   string s;
    61   for (nof_rows = 0; getline(is, s, '\n'); nof_rows++) {
    62    istringstream line(s);
    63    std::vector<double> v;
    64    string tmp_string;
    65    while (line >> tmp_string) {
    66     if(!is.good()) {
    67      cerr << "vector::vector(std::istream&): "
    68     << "error reading data file!" << endl;
    69      exit(1);
    70     }
    71     double t=atof(tmp_string.c_str());
    72     // Here we should check that it actually was a double!!!!!
    73     v.push_back(t);
    74    }
    75    if(nof_columns==0)
    76      nof_columns=v.size();
    77    else if(v.size()!=nof_columns) {
    78     cerr << "vector::vector(std::istream&) data file error: "
    79    << "line" << nof_rows+1 << " has " << v.size()
    80    << " columns; expected " << nof_columns
    81    << " columns"
    82    << endl;
    83     exit(1);
    84    }
    85    data_matrix.push_back(v);
    86   }
    87  
    88   // manipulate the state of the stream to be good
    89   is.clear(std::ios::goodbit);
    90  
    91   // convert the data to a gsl vector and check that data file is a
    92   // column vector or a row vector
    93   if(nof_columns==1) {
    94    v_ = gsl_vector_alloc ( nof_rows );
    95    is_col_ = true;
    96    for(u_int i=0;i<nof_rows;i++)
    97      gsl_vector_set( v_, i, data_matrix[i][0] );
    98   }
    99   else if(nof_rows==1){
    100    v_ = gsl_vector_alloc ( nof_columns );
    101    is_col_ = false;
    102    for(u_int i=0;i<nof_columns;i++)
    103      gsl_vector_set( v_, i, data_matrix[0][i] );
    104   }
    105   else {
    106     cerr << "vector::vector(std::istream&) data file error: "
    107    << "file has " << nof_rows << " rows and " << nof_columns
    108    << " columns; expected a row vector or a column vector"
    109    << endl;
    110     exit(1);
    111   }
    112 }
    113 
    114  
    115 vector::~vector()
    116 {
    117   if( v_ )
    118     {
    119      gsl_vector_free( v_ );
    120      v_ = NULL;
    121     }
    122 }
    123 
    124 
    125 gsl_vector* vector::new_copy( const gsl_vector* p_other )
    126 {
    127   // Get dimenstions
    128   size_t N = p_other->size;
    129    
    130   // Create new empty vector
    131   gsl_vector* p_res =  gsl_vector_alloc( N );
    132 
    133   // Copy p_others elements into p_res
    134   gsl_vector_memcpy( p_res, p_other );
    135 
    136   return p_res;
    137 }
    138 
    139 
    140 // Operators on vectors
    141 ////////////////////////
    142 vector& vector::operator=( const vector& other )
    143 {
    144   if( this != &other )
    145     {
    146      gsl_vector* v_new = new_copy( other.get_gsl_vector() );
    147      if( v_ ) gsl_vector_free( v_ );
    148      v_ = v_new;
    149     }
    150   return *this;
    151 }
    152 
    153 
    154 vector vector::operator+( const vector &other ) const
    155 {
    156   assert( size() == other.size() );
    157   vector res( *this );
    158   gsl_vector_add( res.get_gsl_vector(), other.get_gsl_vector() );
    159   return res;
    160 }
    161 
    162 
    163 vector vector::operator-( const vector &other ) const
    164 {
    165   assert( size() == other.size() );
    166   vector res( *this );
    167   gsl_vector_sub( res.get_gsl_vector(), other.get_gsl_vector() );
    168   return res;
    169 }
    170 
    171 vector vector::operator-() const
    172 {
    173   vector res( *this );
    174   gsl_vector_scale (res.get_gsl_vector() , -1.);
    175   return res;
    176 }
    177 double vector::operator*( const vector &other ) const
    178 {
    179   assert( size() == other.size() );
    180   double res = 0.0;;
    181   for ( size_t i = 0; i < other.size(); ++i )
    182     {
    183      res += other.get( i ) * get( i );
    184     }
    185   return res;
    186 }
    187 
    188 
    189 int vector::operator/=( const vector& other ) const
    190 {
    191   return gsl_vector_div( v_, other.get_gsl_vector() );
    192 }
    193 
    194 // Each element is divided by the "constant"
    195 int vector::operator/=( const double& constant )
    196 {
    197   assert( constant != 0 );
    198   return gsl_vector_scale( v_, 1 / constant );
    199 }
    200 
    201 // This is not the dot product, hence a matrix
    202 // is returned
    203 // matrix vector::operator*( const vector &other ) const
    204 // {
    205 //   assert( rows() == other.cols() );
    206 //   matrix res( rows(), other.cols() );
    207 //   gsl_linalg_matmult( v_, other.get_gsl_vector(),
    208 //          res.get_gsl_vector() );
    209 //   return res;
    210 // }
    211 
    212 
    213 
    214 vector vector::mul_elements( const vector& other ) const
    215 {
    216   assert( size() == other.size() );
    217   vector res( *this );
    218   gsl_vector_mul( res.get_gsl_vector(), other.get_gsl_vector() );
    219   return res;
    220 }
    221 
    222 std::ostream& thep_gsl_api::operator<< ( std::ostream& s_out,
    223            const vector& a )
    224 {
    225   using namespace std;
    226   s_out.setf( ios::fixed ); 
    227 
    228   if( a.is_column() )
    229     {
    230      for ( size_t i = 0; i < a.size(); ++i )
    231        {
    232   s_out << a.get( i  ) << endl;
    233        }
    234     }
    235   else
    236     {
    237      for ( size_t j = 0; j < a.size(); ++j )
    238        {
    239   s_out << a.get( j ) << " ";
    240        }
    241     }
    242 
    243   return s_out;
    244 }
    245 
    246 double vector::sum() const
    247 {
    248   double sum = 0;
    249   for (u_int i=0; i<size(); i++)
    250    sum += gsl_vector_get( v_, i );
    251   return( sum );
    252 } 
    253  
    254  
    255 
    256 // // Public functions (A-Z)
    257 // /////////////////////////
    258 
    259 
     194    return *this;
     195  }
     196
     197
     198
     199  std::ostream& theplu::gslapi::operator<<(std::ostream& s, const vector& a)
     200  {
     201    s.setf(std::ios::fixed);
     202
     203    for (size_t j = 0; j < a.size(); ++j) {
     204      s << a[j];
     205      if ( (j+1)<a.size() )
     206        s << " ";
     207    }
     208
     209    return s;
     210  }
     211
     212
     213}} // of namespace gslapi and namespace thep
Note: See TracChangeset for help on using the changeset viewer.