Jaesung Oh / Mbed 2 deprecated MPC_current_control

quadprog++ and Eigen library test

Dependencies:   mbed Eigen FastPWM

Array.h@2:e843c1b0b25c, 2019-09-24 (annotated)

Committer:
jsoh91
Date:
Tue Sep 24 00:20:12 2019 +0000
Revision:
2:e843c1b0b25c
QP and Eigen library test project;

Who changed what in which revision?

UserRevisionLine numberNew contents of line
jsoh91 2:e843c1b0b25c 1 // $Id: Array.hh 249 2008-11-20 09:58:23Z schaerf $
jsoh91 2:e843c1b0b25c 2 // This file is part of EasyLocalpp: a C++ Object-Oriented framework
jsoh91 2:e843c1b0b25c 3 // aimed at easing the development of Local Search algorithms.
jsoh91 2:e843c1b0b25c 4 // Copyright (C) 2001--2008 Andrea Schaerf, Luca Di Gaspero.
jsoh91 2:e843c1b0b25c 5 //
jsoh91 2:e843c1b0b25c 6 // This software may be modified and distributed under the terms
jsoh91 2:e843c1b0b25c 7 // of the MIT license. See the LICENSE file for details.
jsoh91 2:e843c1b0b25c 8
jsoh91 2:e843c1b0b25c 9 #if !defined(_ARRAY_HH)
jsoh91 2:e843c1b0b25c 10 #define _ARRAY_HH
jsoh91 2:e843c1b0b25c 11
jsoh91 2:e843c1b0b25c 12 #include <set>
jsoh91 2:e843c1b0b25c 13 #include <stdexcept>
jsoh91 2:e843c1b0b25c 14 #include <iostream>
jsoh91 2:e843c1b0b25c 15 #include <iomanip>
jsoh91 2:e843c1b0b25c 16 #include <cmath>
jsoh91 2:e843c1b0b25c 17 #include <cstdlib>
jsoh91 2:e843c1b0b25c 18
jsoh91 2:e843c1b0b25c 19 namespace quadprogpp {
jsoh91 2:e843c1b0b25c 20
jsoh91 2:e843c1b0b25c 21 enum MType { DIAG };
jsoh91 2:e843c1b0b25c 22
jsoh91 2:e843c1b0b25c 23 template <typename T>
jsoh91 2:e843c1b0b25c 24 class Vector
jsoh91 2:e843c1b0b25c 25 {
jsoh91 2:e843c1b0b25c 26 public:
jsoh91 2:e843c1b0b25c 27 Vector();
jsoh91 2:e843c1b0b25c 28 Vector(const unsigned int n);
jsoh91 2:e843c1b0b25c 29 Vector(const T& a, const unsigned int n); //initialize to constant value
jsoh91 2:e843c1b0b25c 30 Vector(const T* a, const unsigned int n); // Initialize to array
jsoh91 2:e843c1b0b25c 31 Vector(const Vector &rhs); // copy constructor
jsoh91 2:e843c1b0b25c 32 ~Vector(); // destructor
jsoh91 2:e843c1b0b25c 33
jsoh91 2:e843c1b0b25c 34 inline void set(const T* a, const unsigned int n);
jsoh91 2:e843c1b0b25c 35 Vector<T> extract(const std::set<unsigned int>& indexes) const;
jsoh91 2:e843c1b0b25c 36 inline T& operator[](const unsigned int& i); //i-th element
jsoh91 2:e843c1b0b25c 37 inline const T& operator[](const unsigned int& i) const;
jsoh91 2:e843c1b0b25c 38
jsoh91 2:e843c1b0b25c 39 inline unsigned int size() const;
jsoh91 2:e843c1b0b25c 40 inline void resize(const unsigned int n);
jsoh91 2:e843c1b0b25c 41 inline void resize(const T& a, const unsigned int n);
jsoh91 2:e843c1b0b25c 42
jsoh91 2:e843c1b0b25c 43 Vector<T>& operator=(const Vector<T>& rhs); //assignment
jsoh91 2:e843c1b0b25c 44 Vector<T>& operator=(const T& a); //assign a to every element
jsoh91 2:e843c1b0b25c 45 inline Vector<T>& operator+=(const Vector<T>& rhs);
jsoh91 2:e843c1b0b25c 46 inline Vector<T>& operator-=(const Vector<T>& rhs);
jsoh91 2:e843c1b0b25c 47 inline Vector<T>& operator*=(const Vector<T>& rhs);
jsoh91 2:e843c1b0b25c 48 inline Vector<T>& operator/=(const Vector<T>& rhs);
jsoh91 2:e843c1b0b25c 49 inline Vector<T>& operator^=(const Vector<T>& rhs);
jsoh91 2:e843c1b0b25c 50 inline Vector<T>& operator+=(const T& a);
jsoh91 2:e843c1b0b25c 51 inline Vector<T>& operator-=(const T& a);
jsoh91 2:e843c1b0b25c 52 inline Vector<T>& operator*=(const T& a);
jsoh91 2:e843c1b0b25c 53 inline Vector<T>& operator/=(const T& a);
jsoh91 2:e843c1b0b25c 54 inline Vector<T>& operator^=(const T& a);
jsoh91 2:e843c1b0b25c 55 private:
jsoh91 2:e843c1b0b25c 56 unsigned int n; // size of array. upper index is n-1
jsoh91 2:e843c1b0b25c 57 T* v; // storage for data
jsoh91 2:e843c1b0b25c 58 };
jsoh91 2:e843c1b0b25c 59
jsoh91 2:e843c1b0b25c 60 template <typename T>
jsoh91 2:e843c1b0b25c 61 Vector<T>::Vector()
jsoh91 2:e843c1b0b25c 62 : n(0), v(0)
jsoh91 2:e843c1b0b25c 63 {}
jsoh91 2:e843c1b0b25c 64
jsoh91 2:e843c1b0b25c 65 template <typename T>
jsoh91 2:e843c1b0b25c 66 Vector<T>::Vector(const unsigned int n)
jsoh91 2:e843c1b0b25c 67 : v(new T[n])
jsoh91 2:e843c1b0b25c 68 {
jsoh91 2:e843c1b0b25c 69 this->n = n;
jsoh91 2:e843c1b0b25c 70 }
jsoh91 2:e843c1b0b25c 71
jsoh91 2:e843c1b0b25c 72 template <typename T>
jsoh91 2:e843c1b0b25c 73 Vector<T>::Vector(const T& a, const unsigned int n)
jsoh91 2:e843c1b0b25c 74 : v(new T[n])
jsoh91 2:e843c1b0b25c 75 {
jsoh91 2:e843c1b0b25c 76 this->n = n;
jsoh91 2:e843c1b0b25c 77 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 78 v[i] = a;
jsoh91 2:e843c1b0b25c 79 }
jsoh91 2:e843c1b0b25c 80
jsoh91 2:e843c1b0b25c 81 template <typename T>
jsoh91 2:e843c1b0b25c 82 Vector<T>::Vector(const T* a, const unsigned int n)
jsoh91 2:e843c1b0b25c 83 : v(new T[n])
jsoh91 2:e843c1b0b25c 84 {
jsoh91 2:e843c1b0b25c 85 this->n = n;
jsoh91 2:e843c1b0b25c 86 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 87 v[i] = *a++;
jsoh91 2:e843c1b0b25c 88 }
jsoh91 2:e843c1b0b25c 89
jsoh91 2:e843c1b0b25c 90 template <typename T>
jsoh91 2:e843c1b0b25c 91 Vector<T>::Vector(const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 92 : v(new T[rhs.n])
jsoh91 2:e843c1b0b25c 93 {
jsoh91 2:e843c1b0b25c 94 this->n = rhs.n;
jsoh91 2:e843c1b0b25c 95 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 96 v[i] = rhs[i];
jsoh91 2:e843c1b0b25c 97 }
jsoh91 2:e843c1b0b25c 98
jsoh91 2:e843c1b0b25c 99 template <typename T>
jsoh91 2:e843c1b0b25c 100 Vector<T>::~Vector()
jsoh91 2:e843c1b0b25c 101 {
jsoh91 2:e843c1b0b25c 102 if (v != 0)
jsoh91 2:e843c1b0b25c 103 delete[] (v);
jsoh91 2:e843c1b0b25c 104 }
jsoh91 2:e843c1b0b25c 105
jsoh91 2:e843c1b0b25c 106 template <typename T>
jsoh91 2:e843c1b0b25c 107 void Vector<T>::resize(const unsigned int n)
jsoh91 2:e843c1b0b25c 108 {
jsoh91 2:e843c1b0b25c 109 if (n == this->n)
jsoh91 2:e843c1b0b25c 110 return;
jsoh91 2:e843c1b0b25c 111 if (v != 0)
jsoh91 2:e843c1b0b25c 112 delete[] (v);
jsoh91 2:e843c1b0b25c 113 v = new T[n];
jsoh91 2:e843c1b0b25c 114 this->n = n;
jsoh91 2:e843c1b0b25c 115 }
jsoh91 2:e843c1b0b25c 116
jsoh91 2:e843c1b0b25c 117 template <typename T>
jsoh91 2:e843c1b0b25c 118 void Vector<T>::resize(const T& a, const unsigned int n)
jsoh91 2:e843c1b0b25c 119 {
jsoh91 2:e843c1b0b25c 120 resize(n);
jsoh91 2:e843c1b0b25c 121 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 122 v[i] = a;
jsoh91 2:e843c1b0b25c 123 }
jsoh91 2:e843c1b0b25c 124
jsoh91 2:e843c1b0b25c 125
jsoh91 2:e843c1b0b25c 126 template <typename T>
jsoh91 2:e843c1b0b25c 127 inline Vector<T>& Vector<T>::operator=(const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 128 // postcondition: normal assignment via copying has been performed;
jsoh91 2:e843c1b0b25c 129 // if vector and rhs were different sizes, vector
jsoh91 2:e843c1b0b25c 130 // has been resized to match the size of rhs
jsoh91 2:e843c1b0b25c 131 {
jsoh91 2:e843c1b0b25c 132 if (this != &rhs)
jsoh91 2:e843c1b0b25c 133 {
jsoh91 2:e843c1b0b25c 134 resize(rhs.n);
jsoh91 2:e843c1b0b25c 135 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 136 v[i] = rhs[i];
jsoh91 2:e843c1b0b25c 137 }
jsoh91 2:e843c1b0b25c 138 return *this;
jsoh91 2:e843c1b0b25c 139 }
jsoh91 2:e843c1b0b25c 140
jsoh91 2:e843c1b0b25c 141 template <typename T>
jsoh91 2:e843c1b0b25c 142 inline Vector<T> & Vector<T>::operator=(const T& a) //assign a to every element
jsoh91 2:e843c1b0b25c 143 {
jsoh91 2:e843c1b0b25c 144 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 145 v[i] = a;
jsoh91 2:e843c1b0b25c 146 return *this;
jsoh91 2:e843c1b0b25c 147 }
jsoh91 2:e843c1b0b25c 148
jsoh91 2:e843c1b0b25c 149 template <typename T>
jsoh91 2:e843c1b0b25c 150 inline T & Vector<T>::operator[](const unsigned int& i) //subscripting
jsoh91 2:e843c1b0b25c 151 {
jsoh91 2:e843c1b0b25c 152 return v[i];
jsoh91 2:e843c1b0b25c 153 }
jsoh91 2:e843c1b0b25c 154
jsoh91 2:e843c1b0b25c 155 template <typename T>
jsoh91 2:e843c1b0b25c 156 inline const T& Vector<T>::operator[](const unsigned int& i) const //subscripting
jsoh91 2:e843c1b0b25c 157 {
jsoh91 2:e843c1b0b25c 158 return v[i];
jsoh91 2:e843c1b0b25c 159 }
jsoh91 2:e843c1b0b25c 160
jsoh91 2:e843c1b0b25c 161 template <typename T>
jsoh91 2:e843c1b0b25c 162 inline unsigned int Vector<T>::size() const
jsoh91 2:e843c1b0b25c 163 {
jsoh91 2:e843c1b0b25c 164 return n;
jsoh91 2:e843c1b0b25c 165 }
jsoh91 2:e843c1b0b25c 166
jsoh91 2:e843c1b0b25c 167 template <typename T>
jsoh91 2:e843c1b0b25c 168 inline void Vector<T>::set(const T* a, unsigned int n)
jsoh91 2:e843c1b0b25c 169 {
jsoh91 2:e843c1b0b25c 170 resize(n);
jsoh91 2:e843c1b0b25c 171 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 172 v[i] = a[i];
jsoh91 2:e843c1b0b25c 173 }
jsoh91 2:e843c1b0b25c 174
jsoh91 2:e843c1b0b25c 175 template <typename T>
jsoh91 2:e843c1b0b25c 176 inline Vector<T> Vector<T>::extract(const std::set<unsigned int>& indexes) const
jsoh91 2:e843c1b0b25c 177 {
jsoh91 2:e843c1b0b25c 178 Vector<T> tmp(indexes.size());
jsoh91 2:e843c1b0b25c 179 unsigned int i = 0;
jsoh91 2:e843c1b0b25c 180
jsoh91 2:e843c1b0b25c 181 for (std::set<unsigned int>::const_iterator el = indexes.begin(); el != indexes.end(); el++)
jsoh91 2:e843c1b0b25c 182 {
jsoh91 2:e843c1b0b25c 183 if (*el >= n)
jsoh91 2:e843c1b0b25c 184 throw std::logic_error("Error extracting subvector: the indexes are out of vector bounds");
jsoh91 2:e843c1b0b25c 185 tmp[i++] = v[*el];
jsoh91 2:e843c1b0b25c 186 }
jsoh91 2:e843c1b0b25c 187
jsoh91 2:e843c1b0b25c 188 return tmp;
jsoh91 2:e843c1b0b25c 189 }
jsoh91 2:e843c1b0b25c 190
jsoh91 2:e843c1b0b25c 191 template <typename T>
jsoh91 2:e843c1b0b25c 192 inline Vector<T>& Vector<T>::operator+=(const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 193 {
jsoh91 2:e843c1b0b25c 194 if (this->size() != rhs.size())
jsoh91 2:e843c1b0b25c 195 throw std::logic_error("Operator+=: vectors have different sizes");
jsoh91 2:e843c1b0b25c 196 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 197 v[i] += rhs[i];
jsoh91 2:e843c1b0b25c 198
jsoh91 2:e843c1b0b25c 199 return *this;
jsoh91 2:e843c1b0b25c 200 }
jsoh91 2:e843c1b0b25c 201
jsoh91 2:e843c1b0b25c 202
jsoh91 2:e843c1b0b25c 203 template <typename T>
jsoh91 2:e843c1b0b25c 204 inline Vector<T>& Vector<T>::operator+=(const T& a)
jsoh91 2:e843c1b0b25c 205 {
jsoh91 2:e843c1b0b25c 206 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 207 v[i] += a;
jsoh91 2:e843c1b0b25c 208
jsoh91 2:e843c1b0b25c 209 return *this;
jsoh91 2:e843c1b0b25c 210 }
jsoh91 2:e843c1b0b25c 211
jsoh91 2:e843c1b0b25c 212 template <typename T>
jsoh91 2:e843c1b0b25c 213 inline Vector<T> operator+(const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 214 {
jsoh91 2:e843c1b0b25c 215 return rhs;
jsoh91 2:e843c1b0b25c 216 }
jsoh91 2:e843c1b0b25c 217
jsoh91 2:e843c1b0b25c 218 template <typename T>
jsoh91 2:e843c1b0b25c 219 inline Vector<T> operator+(const Vector<T>& lhs, const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 220 {
jsoh91 2:e843c1b0b25c 221 if (lhs.size() != rhs.size())
jsoh91 2:e843c1b0b25c 222 throw std::logic_error("Operator+: vectors have different sizes");
jsoh91 2:e843c1b0b25c 223 Vector<T> tmp(lhs.size());
jsoh91 2:e843c1b0b25c 224 for (unsigned int i = 0; i < lhs.size(); i++)
jsoh91 2:e843c1b0b25c 225 tmp[i] = lhs[i] + rhs[i];
jsoh91 2:e843c1b0b25c 226
jsoh91 2:e843c1b0b25c 227 return tmp;
jsoh91 2:e843c1b0b25c 228 }
jsoh91 2:e843c1b0b25c 229
jsoh91 2:e843c1b0b25c 230 template <typename T>
jsoh91 2:e843c1b0b25c 231 inline Vector<T> operator+(const Vector<T>& lhs, const T& a)
jsoh91 2:e843c1b0b25c 232 {
jsoh91 2:e843c1b0b25c 233 Vector<T> tmp(lhs.size());
jsoh91 2:e843c1b0b25c 234 for (unsigned int i = 0; i < lhs.size(); i++)
jsoh91 2:e843c1b0b25c 235 tmp[i] = lhs[i] + a;
jsoh91 2:e843c1b0b25c 236
jsoh91 2:e843c1b0b25c 237 return tmp;
jsoh91 2:e843c1b0b25c 238 }
jsoh91 2:e843c1b0b25c 239
jsoh91 2:e843c1b0b25c 240 template <typename T>
jsoh91 2:e843c1b0b25c 241 inline Vector<T> operator+(const T& a, const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 242 {
jsoh91 2:e843c1b0b25c 243 Vector<T> tmp(rhs.size());
jsoh91 2:e843c1b0b25c 244 for (unsigned int i = 0; i < rhs.size(); i++)
jsoh91 2:e843c1b0b25c 245 tmp[i] = a + rhs[i];
jsoh91 2:e843c1b0b25c 246
jsoh91 2:e843c1b0b25c 247 return tmp;
jsoh91 2:e843c1b0b25c 248 }
jsoh91 2:e843c1b0b25c 249
jsoh91 2:e843c1b0b25c 250 template <typename T>
jsoh91 2:e843c1b0b25c 251 inline Vector<T>& Vector<T>::operator-=(const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 252 {
jsoh91 2:e843c1b0b25c 253 if (this->size() != rhs.size())
jsoh91 2:e843c1b0b25c 254 throw std::logic_error("Operator-=: vectors have different sizes");
jsoh91 2:e843c1b0b25c 255 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 256 v[i] -= rhs[i];
jsoh91 2:e843c1b0b25c 257
jsoh91 2:e843c1b0b25c 258 return *this;
jsoh91 2:e843c1b0b25c 259 }
jsoh91 2:e843c1b0b25c 260
jsoh91 2:e843c1b0b25c 261
jsoh91 2:e843c1b0b25c 262 template <typename T>
jsoh91 2:e843c1b0b25c 263 inline Vector<T>& Vector<T>::operator-=(const T& a)
jsoh91 2:e843c1b0b25c 264 {
jsoh91 2:e843c1b0b25c 265 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 266 v[i] -= a;
jsoh91 2:e843c1b0b25c 267
jsoh91 2:e843c1b0b25c 268 return *this;
jsoh91 2:e843c1b0b25c 269 }
jsoh91 2:e843c1b0b25c 270
jsoh91 2:e843c1b0b25c 271 template <typename T>
jsoh91 2:e843c1b0b25c 272 inline Vector<T> operator-(const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 273 {
jsoh91 2:e843c1b0b25c 274 return (T)(-1) * rhs;
jsoh91 2:e843c1b0b25c 275 }
jsoh91 2:e843c1b0b25c 276
jsoh91 2:e843c1b0b25c 277 template <typename T>
jsoh91 2:e843c1b0b25c 278 inline Vector<T> operator-(const Vector<T>& lhs, const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 279 {
jsoh91 2:e843c1b0b25c 280 if (lhs.size() != rhs.size())
jsoh91 2:e843c1b0b25c 281 throw std::logic_error("Operator-: vectors have different sizes");
jsoh91 2:e843c1b0b25c 282 Vector<T> tmp(lhs.size());
jsoh91 2:e843c1b0b25c 283 for (unsigned int i = 0; i < lhs.size(); i++)
jsoh91 2:e843c1b0b25c 284 tmp[i] = lhs[i] - rhs[i];
jsoh91 2:e843c1b0b25c 285
jsoh91 2:e843c1b0b25c 286 return tmp;
jsoh91 2:e843c1b0b25c 287 }
jsoh91 2:e843c1b0b25c 288
jsoh91 2:e843c1b0b25c 289 template <typename T>
jsoh91 2:e843c1b0b25c 290 inline Vector<T> operator-(const Vector<T>& lhs, const T& a)
jsoh91 2:e843c1b0b25c 291 {
jsoh91 2:e843c1b0b25c 292 Vector<T> tmp(lhs.size());
jsoh91 2:e843c1b0b25c 293 for (unsigned int i = 0; i < lhs.size(); i++)
jsoh91 2:e843c1b0b25c 294 tmp[i] = lhs[i] - a;
jsoh91 2:e843c1b0b25c 295
jsoh91 2:e843c1b0b25c 296 return tmp;
jsoh91 2:e843c1b0b25c 297 }
jsoh91 2:e843c1b0b25c 298
jsoh91 2:e843c1b0b25c 299 template <typename T>
jsoh91 2:e843c1b0b25c 300 inline Vector<T> operator-(const T& a, const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 301 {
jsoh91 2:e843c1b0b25c 302 Vector<T> tmp(rhs.size());
jsoh91 2:e843c1b0b25c 303 for (unsigned int i = 0; i < rhs.size(); i++)
jsoh91 2:e843c1b0b25c 304 tmp[i] = a - rhs[i];
jsoh91 2:e843c1b0b25c 305
jsoh91 2:e843c1b0b25c 306 return tmp;
jsoh91 2:e843c1b0b25c 307 }
jsoh91 2:e843c1b0b25c 308
jsoh91 2:e843c1b0b25c 309 template <typename T>
jsoh91 2:e843c1b0b25c 310 inline Vector<T>& Vector<T>::operator*=(const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 311 {
jsoh91 2:e843c1b0b25c 312 if (this->size() != rhs.size())
jsoh91 2:e843c1b0b25c 313 throw std::logic_error("Operator*=: vectors have different sizes");
jsoh91 2:e843c1b0b25c 314 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 315 v[i] *= rhs[i];
jsoh91 2:e843c1b0b25c 316
jsoh91 2:e843c1b0b25c 317 return *this;
jsoh91 2:e843c1b0b25c 318 }
jsoh91 2:e843c1b0b25c 319
jsoh91 2:e843c1b0b25c 320
jsoh91 2:e843c1b0b25c 321 template <typename T>
jsoh91 2:e843c1b0b25c 322 inline Vector<T>& Vector<T>::operator*=(const T& a)
jsoh91 2:e843c1b0b25c 323 {
jsoh91 2:e843c1b0b25c 324 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 325 v[i] *= a;
jsoh91 2:e843c1b0b25c 326
jsoh91 2:e843c1b0b25c 327 return *this;
jsoh91 2:e843c1b0b25c 328 }
jsoh91 2:e843c1b0b25c 329
jsoh91 2:e843c1b0b25c 330 template <typename T>
jsoh91 2:e843c1b0b25c 331 inline Vector<T> operator*(const Vector<T>& lhs, const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 332 {
jsoh91 2:e843c1b0b25c 333 if (lhs.size() != rhs.size())
jsoh91 2:e843c1b0b25c 334 throw std::logic_error("Operator*: vectors have different sizes");
jsoh91 2:e843c1b0b25c 335 Vector<T> tmp(lhs.size());
jsoh91 2:e843c1b0b25c 336 for (unsigned int i = 0; i < lhs.size(); i++)
jsoh91 2:e843c1b0b25c 337 tmp[i] = lhs[i] * rhs[i];
jsoh91 2:e843c1b0b25c 338
jsoh91 2:e843c1b0b25c 339 return tmp;
jsoh91 2:e843c1b0b25c 340 }
jsoh91 2:e843c1b0b25c 341
jsoh91 2:e843c1b0b25c 342 template <typename T>
jsoh91 2:e843c1b0b25c 343 inline Vector<T> operator*(const Vector<T>& lhs, const T& a)
jsoh91 2:e843c1b0b25c 344 {
jsoh91 2:e843c1b0b25c 345 Vector<T> tmp(lhs.size());
jsoh91 2:e843c1b0b25c 346 for (unsigned int i = 0; i < lhs.size(); i++)
jsoh91 2:e843c1b0b25c 347 tmp[i] = lhs[i] * a;
jsoh91 2:e843c1b0b25c 348
jsoh91 2:e843c1b0b25c 349 return tmp;
jsoh91 2:e843c1b0b25c 350 }
jsoh91 2:e843c1b0b25c 351
jsoh91 2:e843c1b0b25c 352 template <typename T>
jsoh91 2:e843c1b0b25c 353 inline Vector<T> operator*(const T& a, const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 354 {
jsoh91 2:e843c1b0b25c 355 Vector<T> tmp(rhs.size());
jsoh91 2:e843c1b0b25c 356 for (unsigned int i = 0; i < rhs.size(); i++)
jsoh91 2:e843c1b0b25c 357 tmp[i] = a * rhs[i];
jsoh91 2:e843c1b0b25c 358
jsoh91 2:e843c1b0b25c 359 return tmp;
jsoh91 2:e843c1b0b25c 360 }
jsoh91 2:e843c1b0b25c 361
jsoh91 2:e843c1b0b25c 362 template <typename T>
jsoh91 2:e843c1b0b25c 363 inline Vector<T>& Vector<T>::operator/=(const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 364 {
jsoh91 2:e843c1b0b25c 365 if (this->size() != rhs.size())
jsoh91 2:e843c1b0b25c 366 throw std::logic_error("Operator/=: vectors have different sizes");
jsoh91 2:e843c1b0b25c 367 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 368 v[i] /= rhs[i];
jsoh91 2:e843c1b0b25c 369
jsoh91 2:e843c1b0b25c 370 return *this;
jsoh91 2:e843c1b0b25c 371 }
jsoh91 2:e843c1b0b25c 372
jsoh91 2:e843c1b0b25c 373
jsoh91 2:e843c1b0b25c 374 template <typename T>
jsoh91 2:e843c1b0b25c 375 inline Vector<T>& Vector<T>::operator/=(const T& a)
jsoh91 2:e843c1b0b25c 376 {
jsoh91 2:e843c1b0b25c 377 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 378 v[i] /= a;
jsoh91 2:e843c1b0b25c 379
jsoh91 2:e843c1b0b25c 380 return *this;
jsoh91 2:e843c1b0b25c 381 }
jsoh91 2:e843c1b0b25c 382
jsoh91 2:e843c1b0b25c 383 template <typename T>
jsoh91 2:e843c1b0b25c 384 inline Vector<T> operator/(const Vector<T>& lhs, const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 385 {
jsoh91 2:e843c1b0b25c 386 if (lhs.size() != rhs.size())
jsoh91 2:e843c1b0b25c 387 throw std::logic_error("Operator/: vectors have different sizes");
jsoh91 2:e843c1b0b25c 388 Vector<T> tmp(lhs.size());
jsoh91 2:e843c1b0b25c 389 for (unsigned int i = 0; i < lhs.size(); i++)
jsoh91 2:e843c1b0b25c 390 tmp[i] = lhs[i] / rhs[i];
jsoh91 2:e843c1b0b25c 391
jsoh91 2:e843c1b0b25c 392 return tmp;
jsoh91 2:e843c1b0b25c 393 }
jsoh91 2:e843c1b0b25c 394
jsoh91 2:e843c1b0b25c 395 template <typename T>
jsoh91 2:e843c1b0b25c 396 inline Vector<T> operator/(const Vector<T>& lhs, const T& a)
jsoh91 2:e843c1b0b25c 397 {
jsoh91 2:e843c1b0b25c 398 Vector<T> tmp(lhs.size());
jsoh91 2:e843c1b0b25c 399 for (unsigned int i = 0; i < lhs.size(); i++)
jsoh91 2:e843c1b0b25c 400 tmp[i] = lhs[i] / a;
jsoh91 2:e843c1b0b25c 401
jsoh91 2:e843c1b0b25c 402 return tmp;
jsoh91 2:e843c1b0b25c 403 }
jsoh91 2:e843c1b0b25c 404
jsoh91 2:e843c1b0b25c 405 template <typename T>
jsoh91 2:e843c1b0b25c 406 inline Vector<T> operator/(const T& a, const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 407 {
jsoh91 2:e843c1b0b25c 408 Vector<T> tmp(rhs.size());
jsoh91 2:e843c1b0b25c 409 for (unsigned int i = 0; i < rhs.size(); i++)
jsoh91 2:e843c1b0b25c 410 tmp[i] = a / rhs[i];
jsoh91 2:e843c1b0b25c 411
jsoh91 2:e843c1b0b25c 412 return tmp;
jsoh91 2:e843c1b0b25c 413 }
jsoh91 2:e843c1b0b25c 414
jsoh91 2:e843c1b0b25c 415 template <typename T>
jsoh91 2:e843c1b0b25c 416 inline Vector<T> operator^(const Vector<T>& lhs, const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 417 {
jsoh91 2:e843c1b0b25c 418 if (lhs.size() != rhs.size())
jsoh91 2:e843c1b0b25c 419 throw std::logic_error("Operator^: vectors have different sizes");
jsoh91 2:e843c1b0b25c 420 Vector<T> tmp(lhs.size());
jsoh91 2:e843c1b0b25c 421 for (unsigned int i = 0; i < lhs.size(); i++)
jsoh91 2:e843c1b0b25c 422 tmp[i] = pow(lhs[i], rhs[i]);
jsoh91 2:e843c1b0b25c 423
jsoh91 2:e843c1b0b25c 424 return tmp;
jsoh91 2:e843c1b0b25c 425 }
jsoh91 2:e843c1b0b25c 426
jsoh91 2:e843c1b0b25c 427 template <typename T>
jsoh91 2:e843c1b0b25c 428 inline Vector<T> operator^(const Vector<T>& lhs, const T& a)
jsoh91 2:e843c1b0b25c 429 {
jsoh91 2:e843c1b0b25c 430 Vector<T> tmp(lhs.size());
jsoh91 2:e843c1b0b25c 431 for (unsigned int i = 0; i < lhs.size(); i++)
jsoh91 2:e843c1b0b25c 432 tmp[i] = pow(lhs[i], a);
jsoh91 2:e843c1b0b25c 433
jsoh91 2:e843c1b0b25c 434 return tmp;
jsoh91 2:e843c1b0b25c 435 }
jsoh91 2:e843c1b0b25c 436
jsoh91 2:e843c1b0b25c 437 template <typename T>
jsoh91 2:e843c1b0b25c 438 inline Vector<T> operator^(const T& a, const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 439 {
jsoh91 2:e843c1b0b25c 440 Vector<T> tmp(rhs.size());
jsoh91 2:e843c1b0b25c 441 for (unsigned int i = 0; i < rhs.size(); i++)
jsoh91 2:e843c1b0b25c 442 tmp[i] = pow(a, rhs[i]);
jsoh91 2:e843c1b0b25c 443
jsoh91 2:e843c1b0b25c 444 return tmp;
jsoh91 2:e843c1b0b25c 445 }
jsoh91 2:e843c1b0b25c 446
jsoh91 2:e843c1b0b25c 447 template <typename T>
jsoh91 2:e843c1b0b25c 448 inline Vector<T>& Vector<T>::operator^=(const Vector<T>& rhs)
jsoh91 2:e843c1b0b25c 449 {
jsoh91 2:e843c1b0b25c 450 if (this->size() != rhs.size())
jsoh91 2:e843c1b0b25c 451 throw std::logic_error("Operator^=: vectors have different sizes");
jsoh91 2:e843c1b0b25c 452 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 453 v[i] = pow(v[i], rhs[i]);
jsoh91 2:e843c1b0b25c 454
jsoh91 2:e843c1b0b25c 455 return *this;
jsoh91 2:e843c1b0b25c 456 }
jsoh91 2:e843c1b0b25c 457
jsoh91 2:e843c1b0b25c 458 template <typename T>
jsoh91 2:e843c1b0b25c 459 inline Vector<T>& Vector<T>::operator^=(const T& a)
jsoh91 2:e843c1b0b25c 460 {
jsoh91 2:e843c1b0b25c 461 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 462 v[i] = pow(v[i], a);
jsoh91 2:e843c1b0b25c 463
jsoh91 2:e843c1b0b25c 464 return *this;
jsoh91 2:e843c1b0b25c 465 }
jsoh91 2:e843c1b0b25c 466
jsoh91 2:e843c1b0b25c 467 template <typename T>
jsoh91 2:e843c1b0b25c 468 inline bool operator==(const Vector<T>& v, const Vector<T>& w)
jsoh91 2:e843c1b0b25c 469 {
jsoh91 2:e843c1b0b25c 470 if (v.size() != w.size())
jsoh91 2:e843c1b0b25c 471 throw std::logic_error("Vectors of different size are not confrontable");
jsoh91 2:e843c1b0b25c 472 for (unsigned i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 473 if (v[i] != w[i])
jsoh91 2:e843c1b0b25c 474 return false;
jsoh91 2:e843c1b0b25c 475 return true;
jsoh91 2:e843c1b0b25c 476 }
jsoh91 2:e843c1b0b25c 477
jsoh91 2:e843c1b0b25c 478 template <typename T>
jsoh91 2:e843c1b0b25c 479 inline bool operator!=(const Vector<T>& v, const Vector<T>& w)
jsoh91 2:e843c1b0b25c 480 {
jsoh91 2:e843c1b0b25c 481 if (v.size() != w.size())
jsoh91 2:e843c1b0b25c 482 throw std::logic_error("Vectors of different size are not confrontable");
jsoh91 2:e843c1b0b25c 483 for (unsigned i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 484 if (v[i] != w[i])
jsoh91 2:e843c1b0b25c 485 return true;
jsoh91 2:e843c1b0b25c 486 return false;
jsoh91 2:e843c1b0b25c 487 }
jsoh91 2:e843c1b0b25c 488
jsoh91 2:e843c1b0b25c 489 template <typename T>
jsoh91 2:e843c1b0b25c 490 inline bool operator<(const Vector<T>& v, const Vector<T>& w)
jsoh91 2:e843c1b0b25c 491 {
jsoh91 2:e843c1b0b25c 492 if (v.size() != w.size())
jsoh91 2:e843c1b0b25c 493 throw std::logic_error("Vectors of different size are not confrontable");
jsoh91 2:e843c1b0b25c 494 for (unsigned i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 495 if (v[i] >= w[i])
jsoh91 2:e843c1b0b25c 496 return false;
jsoh91 2:e843c1b0b25c 497 return true;
jsoh91 2:e843c1b0b25c 498 }
jsoh91 2:e843c1b0b25c 499
jsoh91 2:e843c1b0b25c 500 template <typename T>
jsoh91 2:e843c1b0b25c 501 inline bool operator<=(const Vector<T>& v, const Vector<T>& w)
jsoh91 2:e843c1b0b25c 502 {
jsoh91 2:e843c1b0b25c 503 if (v.size() != w.size())
jsoh91 2:e843c1b0b25c 504 throw std::logic_error("Vectors of different size are not confrontable");
jsoh91 2:e843c1b0b25c 505 for (unsigned i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 506 if (v[i] > w[i])
jsoh91 2:e843c1b0b25c 507 return false;
jsoh91 2:e843c1b0b25c 508 return true;
jsoh91 2:e843c1b0b25c 509 }
jsoh91 2:e843c1b0b25c 510
jsoh91 2:e843c1b0b25c 511 template <typename T>
jsoh91 2:e843c1b0b25c 512 inline bool operator>(const Vector<T>& v, const Vector<T>& w)
jsoh91 2:e843c1b0b25c 513 {
jsoh91 2:e843c1b0b25c 514 if (v.size() != w.size())
jsoh91 2:e843c1b0b25c 515 throw std::logic_error("Vectors of different size are not confrontable");
jsoh91 2:e843c1b0b25c 516 for (unsigned i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 517 if (v[i] <= w[i])
jsoh91 2:e843c1b0b25c 518 return false;
jsoh91 2:e843c1b0b25c 519 return true;
jsoh91 2:e843c1b0b25c 520 }
jsoh91 2:e843c1b0b25c 521
jsoh91 2:e843c1b0b25c 522 template <typename T>
jsoh91 2:e843c1b0b25c 523 inline bool operator>=(const Vector<T>& v, const Vector<T>& w)
jsoh91 2:e843c1b0b25c 524 {
jsoh91 2:e843c1b0b25c 525 if (v.size() != w.size())
jsoh91 2:e843c1b0b25c 526 throw std::logic_error("Vectors of different size are not confrontable");
jsoh91 2:e843c1b0b25c 527 for (unsigned i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 528 if (v[i] < w[i])
jsoh91 2:e843c1b0b25c 529 return false;
jsoh91 2:e843c1b0b25c 530 return true;
jsoh91 2:e843c1b0b25c 531 }
jsoh91 2:e843c1b0b25c 532
jsoh91 2:e843c1b0b25c 533 /**
jsoh91 2:e843c1b0b25c 534 Input/Output
jsoh91 2:e843c1b0b25c 535 */
jsoh91 2:e843c1b0b25c 536 template <typename T>
jsoh91 2:e843c1b0b25c 537 inline std::ostream& operator<<(std::ostream& os, const Vector<T>& v)
jsoh91 2:e843c1b0b25c 538 {
jsoh91 2:e843c1b0b25c 539 os << std::endl << v.size() << std::endl;
jsoh91 2:e843c1b0b25c 540 for (unsigned int i = 0; i < v.size() - 1; i++)
jsoh91 2:e843c1b0b25c 541 os << std::setw(20) << std::setprecision(16) << v[i] << ", ";
jsoh91 2:e843c1b0b25c 542 os << std::setw(20) << std::setprecision(16) << v[v.size() - 1] << std::endl;
jsoh91 2:e843c1b0b25c 543
jsoh91 2:e843c1b0b25c 544 return os;
jsoh91 2:e843c1b0b25c 545 }
jsoh91 2:e843c1b0b25c 546
jsoh91 2:e843c1b0b25c 547 template <typename T>
jsoh91 2:e843c1b0b25c 548 std::istream& operator>>(std::istream& is, Vector<T>& v)
jsoh91 2:e843c1b0b25c 549 {
jsoh91 2:e843c1b0b25c 550 int elements;
jsoh91 2:e843c1b0b25c 551 char comma;
jsoh91 2:e843c1b0b25c 552 is >> elements;
jsoh91 2:e843c1b0b25c 553 v.resize(elements);
jsoh91 2:e843c1b0b25c 554 for (unsigned int i = 0; i < elements; i++)
jsoh91 2:e843c1b0b25c 555 is >> v[i] >> comma;
jsoh91 2:e843c1b0b25c 556
jsoh91 2:e843c1b0b25c 557 return is;
jsoh91 2:e843c1b0b25c 558 }
jsoh91 2:e843c1b0b25c 559
jsoh91 2:e843c1b0b25c 560 /**
jsoh91 2:e843c1b0b25c 561 Index utilities
jsoh91 2:e843c1b0b25c 562 */
jsoh91 2:e843c1b0b25c 563
jsoh91 2:e843c1b0b25c 564 std::set<unsigned int> seq(unsigned int s, unsigned int e);
jsoh91 2:e843c1b0b25c 565
jsoh91 2:e843c1b0b25c 566 std::set<unsigned int> singleton(unsigned int i);
jsoh91 2:e843c1b0b25c 567
jsoh91 2:e843c1b0b25c 568 template <typename T>
jsoh91 2:e843c1b0b25c 569 class CanonicalBaseVector : public Vector<T>
jsoh91 2:e843c1b0b25c 570 {
jsoh91 2:e843c1b0b25c 571 public:
jsoh91 2:e843c1b0b25c 572 CanonicalBaseVector(unsigned int i, unsigned int n);
jsoh91 2:e843c1b0b25c 573 inline void reset(unsigned int i);
jsoh91 2:e843c1b0b25c 574 private:
jsoh91 2:e843c1b0b25c 575 unsigned int e;
jsoh91 2:e843c1b0b25c 576 };
jsoh91 2:e843c1b0b25c 577
jsoh91 2:e843c1b0b25c 578 template <typename T>
jsoh91 2:e843c1b0b25c 579 CanonicalBaseVector<T>::CanonicalBaseVector(unsigned int i, unsigned int n)
jsoh91 2:e843c1b0b25c 580 : Vector<T>((T)0, n), e(i)
jsoh91 2:e843c1b0b25c 581 { (*this)[e] = (T)1; }
jsoh91 2:e843c1b0b25c 582
jsoh91 2:e843c1b0b25c 583 template <typename T>
jsoh91 2:e843c1b0b25c 584 inline void CanonicalBaseVector<T>::reset(unsigned int i)
jsoh91 2:e843c1b0b25c 585 {
jsoh91 2:e843c1b0b25c 586 (*this)[e] = (T)0;
jsoh91 2:e843c1b0b25c 587 e = i;
jsoh91 2:e843c1b0b25c 588 (*this)[e] = (T)1;
jsoh91 2:e843c1b0b25c 589 }
jsoh91 2:e843c1b0b25c 590
jsoh91 2:e843c1b0b25c 591 #include <stdexcept>
jsoh91 2:e843c1b0b25c 592
jsoh91 2:e843c1b0b25c 593 template <typename T>
jsoh91 2:e843c1b0b25c 594 inline T sum(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 595 {
jsoh91 2:e843c1b0b25c 596 T tmp = (T)0;
jsoh91 2:e843c1b0b25c 597 for (unsigned int i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 598 tmp += v[i];
jsoh91 2:e843c1b0b25c 599
jsoh91 2:e843c1b0b25c 600 return tmp;
jsoh91 2:e843c1b0b25c 601 }
jsoh91 2:e843c1b0b25c 602
jsoh91 2:e843c1b0b25c 603 template <typename T>
jsoh91 2:e843c1b0b25c 604 inline T prod(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 605 {
jsoh91 2:e843c1b0b25c 606 T tmp = (T)1;
jsoh91 2:e843c1b0b25c 607 for (unsigned int i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 608 tmp *= v[i];
jsoh91 2:e843c1b0b25c 609
jsoh91 2:e843c1b0b25c 610 return tmp;
jsoh91 2:e843c1b0b25c 611 }
jsoh91 2:e843c1b0b25c 612
jsoh91 2:e843c1b0b25c 613 template <typename T>
jsoh91 2:e843c1b0b25c 614 inline T mean(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 615 {
jsoh91 2:e843c1b0b25c 616 T sum = (T)0;
jsoh91 2:e843c1b0b25c 617 for (unsigned int i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 618 sum += v[i];
jsoh91 2:e843c1b0b25c 619 return sum / v.size();
jsoh91 2:e843c1b0b25c 620 }
jsoh91 2:e843c1b0b25c 621
jsoh91 2:e843c1b0b25c 622 template <typename T>
jsoh91 2:e843c1b0b25c 623 inline T median(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 624 {
jsoh91 2:e843c1b0b25c 625 Vector<T> tmp = sort(v);
jsoh91 2:e843c1b0b25c 626 if (v.size() % 2 == 1) // it is an odd-sized vector
jsoh91 2:e843c1b0b25c 627 return tmp[v.size() / 2];
jsoh91 2:e843c1b0b25c 628 else
jsoh91 2:e843c1b0b25c 629 return 0.5 * (tmp[v.size() / 2 - 1] + tmp[v.size() / 2]);
jsoh91 2:e843c1b0b25c 630 }
jsoh91 2:e843c1b0b25c 631
jsoh91 2:e843c1b0b25c 632 template <typename T>
jsoh91 2:e843c1b0b25c 633 inline T stdev(const Vector<T>& v, bool sample_correction = false)
jsoh91 2:e843c1b0b25c 634 {
jsoh91 2:e843c1b0b25c 635 return sqrt(var(v, sample_correction));
jsoh91 2:e843c1b0b25c 636 }
jsoh91 2:e843c1b0b25c 637
jsoh91 2:e843c1b0b25c 638 template <typename T>
jsoh91 2:e843c1b0b25c 639 inline T var(const Vector<T>& v, bool sample_correction = false)
jsoh91 2:e843c1b0b25c 640 {
jsoh91 2:e843c1b0b25c 641 T sum = (T)0, ssum = (T)0;
jsoh91 2:e843c1b0b25c 642 unsigned int n = v.size();
jsoh91 2:e843c1b0b25c 643 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 644 {
jsoh91 2:e843c1b0b25c 645 sum += v[i];
jsoh91 2:e843c1b0b25c 646 ssum += (v[i] * v[i]);
jsoh91 2:e843c1b0b25c 647 }
jsoh91 2:e843c1b0b25c 648 if (!sample_correction)
jsoh91 2:e843c1b0b25c 649 return (ssum / n) - (sum / n) * (sum / n);
jsoh91 2:e843c1b0b25c 650 else
jsoh91 2:e843c1b0b25c 651 return n * ((ssum / n) - (sum / n) * (sum / n)) / (n - 1);
jsoh91 2:e843c1b0b25c 652 }
jsoh91 2:e843c1b0b25c 653
jsoh91 2:e843c1b0b25c 654 template <typename T>
jsoh91 2:e843c1b0b25c 655 inline T max(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 656 {
jsoh91 2:e843c1b0b25c 657 T value = v[0];
jsoh91 2:e843c1b0b25c 658 for (unsigned int i = 1; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 659 value = std::max(v[i], value);
jsoh91 2:e843c1b0b25c 660
jsoh91 2:e843c1b0b25c 661 return value;
jsoh91 2:e843c1b0b25c 662 }
jsoh91 2:e843c1b0b25c 663
jsoh91 2:e843c1b0b25c 664 template <typename T>
jsoh91 2:e843c1b0b25c 665 inline T min(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 666 {
jsoh91 2:e843c1b0b25c 667 T value = v[0];
jsoh91 2:e843c1b0b25c 668 for (unsigned int i = 1; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 669 value = std::min(v[i], value);
jsoh91 2:e843c1b0b25c 670
jsoh91 2:e843c1b0b25c 671 return value;
jsoh91 2:e843c1b0b25c 672 }
jsoh91 2:e843c1b0b25c 673
jsoh91 2:e843c1b0b25c 674 template <typename T>
jsoh91 2:e843c1b0b25c 675 inline unsigned int index_max(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 676 {
jsoh91 2:e843c1b0b25c 677 unsigned int max = 0;
jsoh91 2:e843c1b0b25c 678 for (unsigned int i = 1; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 679 if (v[i] > v[max])
jsoh91 2:e843c1b0b25c 680 max = i;
jsoh91 2:e843c1b0b25c 681
jsoh91 2:e843c1b0b25c 682 return max;
jsoh91 2:e843c1b0b25c 683 }
jsoh91 2:e843c1b0b25c 684
jsoh91 2:e843c1b0b25c 685 template <typename T>
jsoh91 2:e843c1b0b25c 686 inline unsigned int index_min(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 687 {
jsoh91 2:e843c1b0b25c 688 unsigned int min = 0;
jsoh91 2:e843c1b0b25c 689 for (unsigned int i = 1; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 690 if (v[i] < v[min])
jsoh91 2:e843c1b0b25c 691 min = i;
jsoh91 2:e843c1b0b25c 692
jsoh91 2:e843c1b0b25c 693 return min;
jsoh91 2:e843c1b0b25c 694 }
jsoh91 2:e843c1b0b25c 695
jsoh91 2:e843c1b0b25c 696
jsoh91 2:e843c1b0b25c 697 template <typename T>
jsoh91 2:e843c1b0b25c 698 inline T dot_prod(const Vector<T>& a, const Vector<T>& b)
jsoh91 2:e843c1b0b25c 699 {
jsoh91 2:e843c1b0b25c 700 T sum = (T)0;
jsoh91 2:e843c1b0b25c 701 if (a.size() != b.size())
jsoh91 2:e843c1b0b25c 702 throw std::logic_error("Dotprod error: the vectors are not the same size");
jsoh91 2:e843c1b0b25c 703 for (unsigned int i = 0; i < a.size(); i++)
jsoh91 2:e843c1b0b25c 704 sum += a[i] * b[i];
jsoh91 2:e843c1b0b25c 705
jsoh91 2:e843c1b0b25c 706 return sum;
jsoh91 2:e843c1b0b25c 707 }
jsoh91 2:e843c1b0b25c 708
jsoh91 2:e843c1b0b25c 709 /**
jsoh91 2:e843c1b0b25c 710 Single element mathematical functions
jsoh91 2:e843c1b0b25c 711 */
jsoh91 2:e843c1b0b25c 712
jsoh91 2:e843c1b0b25c 713 template <typename T>
jsoh91 2:e843c1b0b25c 714 inline Vector<T> exp(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 715 {
jsoh91 2:e843c1b0b25c 716 Vector<T> tmp(v.size());
jsoh91 2:e843c1b0b25c 717 for (unsigned int i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 718 tmp[i] = exp(v[i]);
jsoh91 2:e843c1b0b25c 719
jsoh91 2:e843c1b0b25c 720 return tmp;
jsoh91 2:e843c1b0b25c 721 }
jsoh91 2:e843c1b0b25c 722
jsoh91 2:e843c1b0b25c 723 template <typename T>
jsoh91 2:e843c1b0b25c 724 inline Vector<T> log(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 725 {
jsoh91 2:e843c1b0b25c 726 Vector<T> tmp(v.size());
jsoh91 2:e843c1b0b25c 727 for (unsigned int i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 728 tmp[i] = log(v[i]);
jsoh91 2:e843c1b0b25c 729
jsoh91 2:e843c1b0b25c 730 return tmp;
jsoh91 2:e843c1b0b25c 731 }
jsoh91 2:e843c1b0b25c 732
jsoh91 2:e843c1b0b25c 733 template <typename T>
jsoh91 2:e843c1b0b25c 734 inline Vector<T> vec_sqrt(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 735 {
jsoh91 2:e843c1b0b25c 736 Vector<T> tmp(v.size());
jsoh91 2:e843c1b0b25c 737 for (unsigned int i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 738 tmp[i] = sqrt(v[i]);
jsoh91 2:e843c1b0b25c 739
jsoh91 2:e843c1b0b25c 740 return tmp;
jsoh91 2:e843c1b0b25c 741 }
jsoh91 2:e843c1b0b25c 742
jsoh91 2:e843c1b0b25c 743 template <typename T>
jsoh91 2:e843c1b0b25c 744 inline Vector<T> pow(const Vector<T>& v, double a)
jsoh91 2:e843c1b0b25c 745 {
jsoh91 2:e843c1b0b25c 746 Vector<T> tmp(v.size());
jsoh91 2:e843c1b0b25c 747 for (unsigned int i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 748 tmp[i] = pow(v[i], a);
jsoh91 2:e843c1b0b25c 749
jsoh91 2:e843c1b0b25c 750 return tmp;
jsoh91 2:e843c1b0b25c 751 }
jsoh91 2:e843c1b0b25c 752
jsoh91 2:e843c1b0b25c 753 template <typename T>
jsoh91 2:e843c1b0b25c 754 inline Vector<T> abs(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 755 {
jsoh91 2:e843c1b0b25c 756 Vector<T> tmp(v.size());
jsoh91 2:e843c1b0b25c 757 for (unsigned int i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 758 tmp[i] = (T)fabs(v[i]);
jsoh91 2:e843c1b0b25c 759
jsoh91 2:e843c1b0b25c 760 return tmp;
jsoh91 2:e843c1b0b25c 761 }
jsoh91 2:e843c1b0b25c 762
jsoh91 2:e843c1b0b25c 763 template <typename T>
jsoh91 2:e843c1b0b25c 764 inline Vector<T> sign(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 765 {
jsoh91 2:e843c1b0b25c 766 Vector<T> tmp(v.size());
jsoh91 2:e843c1b0b25c 767 for (unsigned int i = 0; i < v.size(); i++)
jsoh91 2:e843c1b0b25c 768 tmp[i] = v[i] > 0 ? +1 : v[i] == 0 ? 0 : -1;
jsoh91 2:e843c1b0b25c 769
jsoh91 2:e843c1b0b25c 770 return tmp;
jsoh91 2:e843c1b0b25c 771 }
jsoh91 2:e843c1b0b25c 772
jsoh91 2:e843c1b0b25c 773 template <typename T>
jsoh91 2:e843c1b0b25c 774 inline unsigned int partition(Vector<T>& v, unsigned int begin, unsigned int end)
jsoh91 2:e843c1b0b25c 775 {
jsoh91 2:e843c1b0b25c 776 unsigned int i = begin + 1, j = begin + 1;
jsoh91 2:e843c1b0b25c 777 T pivot = v[begin];
jsoh91 2:e843c1b0b25c 778 while (j <= end)
jsoh91 2:e843c1b0b25c 779 {
jsoh91 2:e843c1b0b25c 780 if (v[j] < pivot) {
jsoh91 2:e843c1b0b25c 781 std::swap(v[i], v[j]);
jsoh91 2:e843c1b0b25c 782 i++;
jsoh91 2:e843c1b0b25c 783 }
jsoh91 2:e843c1b0b25c 784 j++;
jsoh91 2:e843c1b0b25c 785 }
jsoh91 2:e843c1b0b25c 786 v[begin] = v[i - 1];
jsoh91 2:e843c1b0b25c 787 v[i - 1] = pivot;
jsoh91 2:e843c1b0b25c 788 return i - 2;
jsoh91 2:e843c1b0b25c 789 }
jsoh91 2:e843c1b0b25c 790
jsoh91 2:e843c1b0b25c 791
jsoh91 2:e843c1b0b25c 792 template <typename T>
jsoh91 2:e843c1b0b25c 793 inline void quicksort(Vector<T>& v, unsigned int begin, unsigned int end)
jsoh91 2:e843c1b0b25c 794 {
jsoh91 2:e843c1b0b25c 795 if (end > begin)
jsoh91 2:e843c1b0b25c 796 {
jsoh91 2:e843c1b0b25c 797 unsigned int index = partition(v, begin, end);
jsoh91 2:e843c1b0b25c 798 quicksort(v, begin, index);
jsoh91 2:e843c1b0b25c 799 quicksort(v, index + 2, end);
jsoh91 2:e843c1b0b25c 800 }
jsoh91 2:e843c1b0b25c 801 }
jsoh91 2:e843c1b0b25c 802
jsoh91 2:e843c1b0b25c 803 template <typename T>
jsoh91 2:e843c1b0b25c 804 inline Vector<T> sort(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 805 {
jsoh91 2:e843c1b0b25c 806 Vector<T> tmp(v);
jsoh91 2:e843c1b0b25c 807
jsoh91 2:e843c1b0b25c 808 quicksort<T>(tmp, 0, tmp.size() - 1);
jsoh91 2:e843c1b0b25c 809
jsoh91 2:e843c1b0b25c 810 return tmp;
jsoh91 2:e843c1b0b25c 811 }
jsoh91 2:e843c1b0b25c 812
jsoh91 2:e843c1b0b25c 813 template <typename T>
jsoh91 2:e843c1b0b25c 814 inline Vector<double> rank(const Vector<T>& v)
jsoh91 2:e843c1b0b25c 815 {
jsoh91 2:e843c1b0b25c 816 Vector<T> tmp(v);
jsoh91 2:e843c1b0b25c 817 Vector<double> tmp_rank(0.0, v.size());
jsoh91 2:e843c1b0b25c 818
jsoh91 2:e843c1b0b25c 819 for (unsigned int i = 0; i < tmp.size(); i++)
jsoh91 2:e843c1b0b25c 820 {
jsoh91 2:e843c1b0b25c 821 unsigned int smaller = 0, equal = 0;
jsoh91 2:e843c1b0b25c 822 for (unsigned int j = 0; j < tmp.size(); j++)
jsoh91 2:e843c1b0b25c 823 if (i == j)
jsoh91 2:e843c1b0b25c 824 continue;
jsoh91 2:e843c1b0b25c 825 else
jsoh91 2:e843c1b0b25c 826 if (tmp[j] < tmp[i])
jsoh91 2:e843c1b0b25c 827 smaller++;
jsoh91 2:e843c1b0b25c 828 else if (tmp[j] == tmp[i])
jsoh91 2:e843c1b0b25c 829 equal++;
jsoh91 2:e843c1b0b25c 830 tmp_rank[i] = smaller + 1;
jsoh91 2:e843c1b0b25c 831 if (equal > 0)
jsoh91 2:e843c1b0b25c 832 {
jsoh91 2:e843c1b0b25c 833 for (unsigned int j = 1; j <= equal; j++)
jsoh91 2:e843c1b0b25c 834 tmp_rank[i] += smaller + 1 + j;
jsoh91 2:e843c1b0b25c 835 tmp_rank[i] /= (double)(equal + 1);
jsoh91 2:e843c1b0b25c 836 }
jsoh91 2:e843c1b0b25c 837 }
jsoh91 2:e843c1b0b25c 838
jsoh91 2:e843c1b0b25c 839 return tmp_rank;
jsoh91 2:e843c1b0b25c 840 }
jsoh91 2:e843c1b0b25c 841
jsoh91 2:e843c1b0b25c 842 //enum MType { DIAG };
jsoh91 2:e843c1b0b25c 843
jsoh91 2:e843c1b0b25c 844 template <typename T>
jsoh91 2:e843c1b0b25c 845 class Matrix
jsoh91 2:e843c1b0b25c 846 {
jsoh91 2:e843c1b0b25c 847 public:
jsoh91 2:e843c1b0b25c 848 Matrix(); // Default constructor
jsoh91 2:e843c1b0b25c 849 Matrix(const unsigned int n, const unsigned int m); // Construct a n x m matrix
jsoh91 2:e843c1b0b25c 850 Matrix(const T& a, const unsigned int n, const unsigned int m); // Initialize the content to constant a
jsoh91 2:e843c1b0b25c 851 Matrix(MType t, const T& a, const T& o, const unsigned int n, const unsigned int m);
jsoh91 2:e843c1b0b25c 852 Matrix(MType t, const Vector<T>& v, const T& o, const unsigned int n, const unsigned int m);
jsoh91 2:e843c1b0b25c 853 Matrix(const T* a, const unsigned int n, const unsigned int m); // Initialize to array
jsoh91 2:e843c1b0b25c 854 Matrix(const Matrix<T>& rhs); // Copy constructor
jsoh91 2:e843c1b0b25c 855 ~Matrix(); // destructor
jsoh91 2:e843c1b0b25c 856
jsoh91 2:e843c1b0b25c 857 inline T* operator[](const unsigned int& i) { return v[i]; } // Subscripting: row i
jsoh91 2:e843c1b0b25c 858 inline const T* operator[](const unsigned int& i) const { return v[i]; }; // const subsctipting
jsoh91 2:e843c1b0b25c 859
jsoh91 2:e843c1b0b25c 860 inline void resize(const unsigned int n, const unsigned int m);
jsoh91 2:e843c1b0b25c 861 inline void resize(const T& a, const unsigned int n, const unsigned int m);
jsoh91 2:e843c1b0b25c 862
jsoh91 2:e843c1b0b25c 863
jsoh91 2:e843c1b0b25c 864 inline Vector<T> extractRow(const unsigned int i) const;
jsoh91 2:e843c1b0b25c 865 inline Vector<T> extractColumn(const unsigned int j) const;
jsoh91 2:e843c1b0b25c 866 inline Vector<T> extractDiag() const;
jsoh91 2:e843c1b0b25c 867 inline Matrix<T> extractRows(const std::set<unsigned int>& indexes) const;
jsoh91 2:e843c1b0b25c 868 inline Matrix<T> extractColumns(const std::set<unsigned int>& indexes) const;
jsoh91 2:e843c1b0b25c 869 inline Matrix<T> extract(const std::set<unsigned int>& r_indexes, const std::set<unsigned int>& c_indexes) const;
jsoh91 2:e843c1b0b25c 870
jsoh91 2:e843c1b0b25c 871 inline void set(const T* a, unsigned int n, unsigned int m);
jsoh91 2:e843c1b0b25c 872 inline void set(const std::set<unsigned int>& r_indexes, const std::set<unsigned int>& c_indexes, const Matrix<T>& m);
jsoh91 2:e843c1b0b25c 873 inline void setRow(const unsigned int index, const Vector<T>& v);
jsoh91 2:e843c1b0b25c 874 inline void setRow(const unsigned int index, const Matrix<T>& v);
jsoh91 2:e843c1b0b25c 875 inline void setRows(const std::set<unsigned int>& indexes, const Matrix<T>& m);
jsoh91 2:e843c1b0b25c 876 inline void setColumn(const unsigned int index, const Vector<T>& v);
jsoh91 2:e843c1b0b25c 877 inline void setColumn(const unsigned int index, const Matrix<T>& v);
jsoh91 2:e843c1b0b25c 878 inline void setColumns(const std::set<unsigned int>& indexes, const Matrix<T>& m);
jsoh91 2:e843c1b0b25c 879
jsoh91 2:e843c1b0b25c 880
jsoh91 2:e843c1b0b25c 881 inline unsigned int nrows() const { return n; } // number of rows
jsoh91 2:e843c1b0b25c 882 inline unsigned int ncols() const { return m; } // number of columns
jsoh91 2:e843c1b0b25c 883
jsoh91 2:e843c1b0b25c 884 inline Matrix<T>& operator=(const Matrix<T>& rhs); // Assignment operator
jsoh91 2:e843c1b0b25c 885 inline Matrix<T>& operator=(const T& a); // Assign to every element value a
jsoh91 2:e843c1b0b25c 886 inline Matrix<T>& operator+=(const Matrix<T>& rhs);
jsoh91 2:e843c1b0b25c 887 inline Matrix<T>& operator-=(const Matrix<T>& rhs);
jsoh91 2:e843c1b0b25c 888 inline Matrix<T>& operator*=(const Matrix<T>& rhs);
jsoh91 2:e843c1b0b25c 889 inline Matrix<T>& operator/=(const Matrix<T>& rhs);
jsoh91 2:e843c1b0b25c 890 inline Matrix<T>& operator^=(const Matrix<T>& rhs);
jsoh91 2:e843c1b0b25c 891 inline Matrix<T>& operator+=(const T& a);
jsoh91 2:e843c1b0b25c 892 inline Matrix<T>& operator-=(const T& a);
jsoh91 2:e843c1b0b25c 893 inline Matrix<T>& operator*=(const T& a);
jsoh91 2:e843c1b0b25c 894 inline Matrix<T>& operator/=(const T& a);
jsoh91 2:e843c1b0b25c 895 inline Matrix<T>& operator^=(const T& a);
jsoh91 2:e843c1b0b25c 896 inline operator Vector<T>();
jsoh91 2:e843c1b0b25c 897 private:
jsoh91 2:e843c1b0b25c 898 unsigned int n; // number of rows
jsoh91 2:e843c1b0b25c 899 unsigned int m; // number of columns
jsoh91 2:e843c1b0b25c 900 T **v; // storage for data
jsoh91 2:e843c1b0b25c 901 };
jsoh91 2:e843c1b0b25c 902
jsoh91 2:e843c1b0b25c 903 template <typename T>
jsoh91 2:e843c1b0b25c 904 Matrix<T>::Matrix()
jsoh91 2:e843c1b0b25c 905 : n(0), m(0), v(0)
jsoh91 2:e843c1b0b25c 906 {}
jsoh91 2:e843c1b0b25c 907
jsoh91 2:e843c1b0b25c 908 template <typename T>
jsoh91 2:e843c1b0b25c 909 Matrix<T>::Matrix(unsigned int n, unsigned int m)
jsoh91 2:e843c1b0b25c 910 : v(new T*[n])
jsoh91 2:e843c1b0b25c 911 {
jsoh91 2:e843c1b0b25c 912 register unsigned int i;
jsoh91 2:e843c1b0b25c 913 this->n = n; this->m = m;
jsoh91 2:e843c1b0b25c 914 v[0] = new T[m * n];
jsoh91 2:e843c1b0b25c 915 for (i = 1; i < n; i++)
jsoh91 2:e843c1b0b25c 916 v[i] = v[i - 1] + m;
jsoh91 2:e843c1b0b25c 917 }
jsoh91 2:e843c1b0b25c 918
jsoh91 2:e843c1b0b25c 919 template <typename T>
jsoh91 2:e843c1b0b25c 920 Matrix<T>::Matrix(const T& a, unsigned int n, unsigned int m)
jsoh91 2:e843c1b0b25c 921 : v(new T*[n])
jsoh91 2:e843c1b0b25c 922 {
jsoh91 2:e843c1b0b25c 923 register unsigned int i, j;
jsoh91 2:e843c1b0b25c 924 this->n = n; this->m = m;
jsoh91 2:e843c1b0b25c 925 v[0] = new T[m * n];
jsoh91 2:e843c1b0b25c 926 for (i = 1; i < n; i++)
jsoh91 2:e843c1b0b25c 927 v[i] = v[i - 1] + m;
jsoh91 2:e843c1b0b25c 928 for (i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 929 for (j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 930 v[i][j] = a;
jsoh91 2:e843c1b0b25c 931 }
jsoh91 2:e843c1b0b25c 932
jsoh91 2:e843c1b0b25c 933 template <class T>
jsoh91 2:e843c1b0b25c 934 Matrix<T>::Matrix(const T* a, unsigned int n, unsigned int m)
jsoh91 2:e843c1b0b25c 935 : v(new T*[n])
jsoh91 2:e843c1b0b25c 936 {
jsoh91 2:e843c1b0b25c 937 register unsigned int i, j;
jsoh91 2:e843c1b0b25c 938 this->n = n; this->m = m;
jsoh91 2:e843c1b0b25c 939 v[0] = new T[m * n];
jsoh91 2:e843c1b0b25c 940 for (i = 1; i < n; i++)
jsoh91 2:e843c1b0b25c 941 v[i] = v[i - 1] + m;
jsoh91 2:e843c1b0b25c 942 for (i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 943 for (j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 944 v[i][j] = *a++;
jsoh91 2:e843c1b0b25c 945 }
jsoh91 2:e843c1b0b25c 946
jsoh91 2:e843c1b0b25c 947 template <class T>
jsoh91 2:e843c1b0b25c 948 Matrix<T>::Matrix(MType t, const T& a, const T& o, unsigned int n, unsigned int m)
jsoh91 2:e843c1b0b25c 949 : v(new T*[n])
jsoh91 2:e843c1b0b25c 950 {
jsoh91 2:e843c1b0b25c 951 register unsigned int i, j;
jsoh91 2:e843c1b0b25c 952 this->n = n; this->m = m;
jsoh91 2:e843c1b0b25c 953 v[0] = new T[m * n];
jsoh91 2:e843c1b0b25c 954 for (i = 1; i < n; i++)
jsoh91 2:e843c1b0b25c 955 v[i] = v[i - 1] + m;
jsoh91 2:e843c1b0b25c 956 switch (t)
jsoh91 2:e843c1b0b25c 957 {
jsoh91 2:e843c1b0b25c 958 case DIAG:
jsoh91 2:e843c1b0b25c 959 for (i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 960 for (j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 961 if (i != j)
jsoh91 2:e843c1b0b25c 962 v[i][j] = o;
jsoh91 2:e843c1b0b25c 963 else
jsoh91 2:e843c1b0b25c 964 v[i][j] = a;
jsoh91 2:e843c1b0b25c 965 break;
jsoh91 2:e843c1b0b25c 966 default:
jsoh91 2:e843c1b0b25c 967 throw std::logic_error("Matrix type not supported");
jsoh91 2:e843c1b0b25c 968 }
jsoh91 2:e843c1b0b25c 969 }
jsoh91 2:e843c1b0b25c 970
jsoh91 2:e843c1b0b25c 971 template <class T>
jsoh91 2:e843c1b0b25c 972 Matrix<T>::Matrix(MType t, const Vector<T>& a, const T& o, unsigned int n, unsigned int m)
jsoh91 2:e843c1b0b25c 973 : v(new T*[n])
jsoh91 2:e843c1b0b25c 974 {
jsoh91 2:e843c1b0b25c 975 register unsigned int i, j;
jsoh91 2:e843c1b0b25c 976 this->n = n; this->m = m;
jsoh91 2:e843c1b0b25c 977 v[0] = new T[m * n];
jsoh91 2:e843c1b0b25c 978 for (i = 1; i < n; i++)
jsoh91 2:e843c1b0b25c 979 v[i] = v[i - 1] + m;
jsoh91 2:e843c1b0b25c 980 switch (t)
jsoh91 2:e843c1b0b25c 981 {
jsoh91 2:e843c1b0b25c 982 case DIAG:
jsoh91 2:e843c1b0b25c 983 for (i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 984 for (j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 985 if (i != j)
jsoh91 2:e843c1b0b25c 986 v[i][j] = o;
jsoh91 2:e843c1b0b25c 987 else
jsoh91 2:e843c1b0b25c 988 v[i][j] = a[i];
jsoh91 2:e843c1b0b25c 989 break;
jsoh91 2:e843c1b0b25c 990 default:
jsoh91 2:e843c1b0b25c 991 throw std::logic_error("Matrix type not supported");
jsoh91 2:e843c1b0b25c 992 }
jsoh91 2:e843c1b0b25c 993 }
jsoh91 2:e843c1b0b25c 994
jsoh91 2:e843c1b0b25c 995 template <typename T>
jsoh91 2:e843c1b0b25c 996 Matrix<T>::Matrix(const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 997 : v(new T*[rhs.n])
jsoh91 2:e843c1b0b25c 998 {
jsoh91 2:e843c1b0b25c 999 register unsigned int i, j;
jsoh91 2:e843c1b0b25c 1000 n = rhs.n; m = rhs.m;
jsoh91 2:e843c1b0b25c 1001 v[0] = new T[m * n];
jsoh91 2:e843c1b0b25c 1002 for (i = 1; i < n; i++)
jsoh91 2:e843c1b0b25c 1003 v[i] = v[i - 1] + m;
jsoh91 2:e843c1b0b25c 1004 for (i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1005 for (j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1006 v[i][j] = rhs[i][j];
jsoh91 2:e843c1b0b25c 1007 }
jsoh91 2:e843c1b0b25c 1008
jsoh91 2:e843c1b0b25c 1009 template <typename T>
jsoh91 2:e843c1b0b25c 1010 Matrix<T>::~Matrix()
jsoh91 2:e843c1b0b25c 1011 {
jsoh91 2:e843c1b0b25c 1012 if (v != 0) {
jsoh91 2:e843c1b0b25c 1013 delete[] (v[0]);
jsoh91 2:e843c1b0b25c 1014 delete[] (v);
jsoh91 2:e843c1b0b25c 1015 }
jsoh91 2:e843c1b0b25c 1016 }
jsoh91 2:e843c1b0b25c 1017
jsoh91 2:e843c1b0b25c 1018 template <typename T>
jsoh91 2:e843c1b0b25c 1019 inline Matrix<T>& Matrix<T>::operator=(const Matrix<T> &rhs)
jsoh91 2:e843c1b0b25c 1020 // postcondition: normal assignment via copying has been performed;
jsoh91 2:e843c1b0b25c 1021 // if matrix and rhs were different sizes, matrix
jsoh91 2:e843c1b0b25c 1022 // has been resized to match the size of rhs
jsoh91 2:e843c1b0b25c 1023 {
jsoh91 2:e843c1b0b25c 1024 register unsigned int i, j;
jsoh91 2:e843c1b0b25c 1025 if (this != &rhs)
jsoh91 2:e843c1b0b25c 1026 {
jsoh91 2:e843c1b0b25c 1027 resize(rhs.n, rhs.m);
jsoh91 2:e843c1b0b25c 1028 for (i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1029 for (j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1030 v[i][j] = rhs[i][j];
jsoh91 2:e843c1b0b25c 1031 }
jsoh91 2:e843c1b0b25c 1032 return *this;
jsoh91 2:e843c1b0b25c 1033 }
jsoh91 2:e843c1b0b25c 1034
jsoh91 2:e843c1b0b25c 1035 template <typename T>
jsoh91 2:e843c1b0b25c 1036 inline Matrix<T>& Matrix<T>::operator=(const T& a) // assign a to every element
jsoh91 2:e843c1b0b25c 1037 {
jsoh91 2:e843c1b0b25c 1038 register unsigned int i, j;
jsoh91 2:e843c1b0b25c 1039 for (i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1040 for (j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1041 v[i][j] = a;
jsoh91 2:e843c1b0b25c 1042 return *this;
jsoh91 2:e843c1b0b25c 1043 }
jsoh91 2:e843c1b0b25c 1044
jsoh91 2:e843c1b0b25c 1045
jsoh91 2:e843c1b0b25c 1046 template <typename T>
jsoh91 2:e843c1b0b25c 1047 inline void Matrix<T>::resize(const unsigned int n, const unsigned int m)
jsoh91 2:e843c1b0b25c 1048 {
jsoh91 2:e843c1b0b25c 1049 register unsigned int i;
jsoh91 2:e843c1b0b25c 1050 if (n == this->n && m == this->m)
jsoh91 2:e843c1b0b25c 1051 return;
jsoh91 2:e843c1b0b25c 1052 if (v != 0)
jsoh91 2:e843c1b0b25c 1053 {
jsoh91 2:e843c1b0b25c 1054 delete[] (v[0]);
jsoh91 2:e843c1b0b25c 1055 delete[] (v);
jsoh91 2:e843c1b0b25c 1056 }
jsoh91 2:e843c1b0b25c 1057 this->n = n; this->m = m;
jsoh91 2:e843c1b0b25c 1058 v = new T*[n];
jsoh91 2:e843c1b0b25c 1059 v[0] = new T[m * n];
jsoh91 2:e843c1b0b25c 1060 for (i = 1; i < n; i++)
jsoh91 2:e843c1b0b25c 1061 v[i] = v[i - 1] + m;
jsoh91 2:e843c1b0b25c 1062 }
jsoh91 2:e843c1b0b25c 1063
jsoh91 2:e843c1b0b25c 1064 template <typename T>
jsoh91 2:e843c1b0b25c 1065 inline void Matrix<T>::resize(const T& a, const unsigned int n, const unsigned int m)
jsoh91 2:e843c1b0b25c 1066 {
jsoh91 2:e843c1b0b25c 1067 register unsigned int i, j;
jsoh91 2:e843c1b0b25c 1068 resize(n, m);
jsoh91 2:e843c1b0b25c 1069 for (i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1070 for (j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1071 v[i][j] = a;
jsoh91 2:e843c1b0b25c 1072 }
jsoh91 2:e843c1b0b25c 1073
jsoh91 2:e843c1b0b25c 1074
jsoh91 2:e843c1b0b25c 1075
jsoh91 2:e843c1b0b25c 1076 template <typename T>
jsoh91 2:e843c1b0b25c 1077 inline Vector<T> Matrix<T>::extractRow(const unsigned int i) const
jsoh91 2:e843c1b0b25c 1078 {
jsoh91 2:e843c1b0b25c 1079 if (i >= n)
jsoh91 2:e843c1b0b25c 1080 throw std::logic_error("Error in extractRow: trying to extract a row out of matrix bounds");
jsoh91 2:e843c1b0b25c 1081 Vector<T> tmp(v[i], m);
jsoh91 2:e843c1b0b25c 1082
jsoh91 2:e843c1b0b25c 1083 return tmp;
jsoh91 2:e843c1b0b25c 1084 }
jsoh91 2:e843c1b0b25c 1085
jsoh91 2:e843c1b0b25c 1086 template <typename T>
jsoh91 2:e843c1b0b25c 1087 inline Vector<T> Matrix<T>::extractColumn(const unsigned int j) const
jsoh91 2:e843c1b0b25c 1088 {
jsoh91 2:e843c1b0b25c 1089 register unsigned int i;
jsoh91 2:e843c1b0b25c 1090 if (j >= m)
jsoh91 2:e843c1b0b25c 1091 throw std::logic_error("Error in extractRow: trying to extract a row out of matrix bounds");
jsoh91 2:e843c1b0b25c 1092 Vector<T> tmp(n);
jsoh91 2:e843c1b0b25c 1093
jsoh91 2:e843c1b0b25c 1094 for (i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1095 tmp[i] = v[i][j];
jsoh91 2:e843c1b0b25c 1096
jsoh91 2:e843c1b0b25c 1097 return tmp;
jsoh91 2:e843c1b0b25c 1098 }
jsoh91 2:e843c1b0b25c 1099
jsoh91 2:e843c1b0b25c 1100 template <typename T>
jsoh91 2:e843c1b0b25c 1101 inline Vector<T> Matrix<T>::extractDiag() const
jsoh91 2:e843c1b0b25c 1102 {
jsoh91 2:e843c1b0b25c 1103 register unsigned int d = std::min(n, m), i;
jsoh91 2:e843c1b0b25c 1104
jsoh91 2:e843c1b0b25c 1105 Vector<T> tmp(d);
jsoh91 2:e843c1b0b25c 1106
jsoh91 2:e843c1b0b25c 1107 for (i = 0; i < d; i++)
jsoh91 2:e843c1b0b25c 1108 tmp[i] = v[i][i];
jsoh91 2:e843c1b0b25c 1109
jsoh91 2:e843c1b0b25c 1110 return tmp;
jsoh91 2:e843c1b0b25c 1111
jsoh91 2:e843c1b0b25c 1112 }
jsoh91 2:e843c1b0b25c 1113
jsoh91 2:e843c1b0b25c 1114 template <typename T>
jsoh91 2:e843c1b0b25c 1115 inline Matrix<T> Matrix<T>::extractRows(const std::set<unsigned int>& indexes) const
jsoh91 2:e843c1b0b25c 1116 {
jsoh91 2:e843c1b0b25c 1117 Matrix<T> tmp(indexes.size(), m);
jsoh91 2:e843c1b0b25c 1118 register unsigned int i = 0, j;
jsoh91 2:e843c1b0b25c 1119
jsoh91 2:e843c1b0b25c 1120 for (std::set<unsigned int>::const_iterator el = indexes.begin(); el != indexes.end(); el++)
jsoh91 2:e843c1b0b25c 1121 {
jsoh91 2:e843c1b0b25c 1122 for (j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1123 {
jsoh91 2:e843c1b0b25c 1124 if (*el >= n)
jsoh91 2:e843c1b0b25c 1125 throw std::logic_error("Error extracting rows: the indexes are out of matrix bounds");
jsoh91 2:e843c1b0b25c 1126 tmp[i][j] = v[*el][j];
jsoh91 2:e843c1b0b25c 1127 }
jsoh91 2:e843c1b0b25c 1128 i++;
jsoh91 2:e843c1b0b25c 1129 }
jsoh91 2:e843c1b0b25c 1130
jsoh91 2:e843c1b0b25c 1131 return tmp;
jsoh91 2:e843c1b0b25c 1132 }
jsoh91 2:e843c1b0b25c 1133
jsoh91 2:e843c1b0b25c 1134 template <typename T>
jsoh91 2:e843c1b0b25c 1135 inline Matrix<T> Matrix<T>::extractColumns(const std::set<unsigned int>& indexes) const
jsoh91 2:e843c1b0b25c 1136 {
jsoh91 2:e843c1b0b25c 1137 Matrix<T> tmp(n, indexes.size());
jsoh91 2:e843c1b0b25c 1138 register unsigned int i, j = 0;
jsoh91 2:e843c1b0b25c 1139
jsoh91 2:e843c1b0b25c 1140 for (std::set<unsigned int>::const_iterator el = indexes.begin(); el != indexes.end(); el++)
jsoh91 2:e843c1b0b25c 1141 {
jsoh91 2:e843c1b0b25c 1142 for (i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1143 {
jsoh91 2:e843c1b0b25c 1144 if (*el >= m)
jsoh91 2:e843c1b0b25c 1145 throw std::logic_error("Error extracting columns: the indexes are out of matrix bounds");
jsoh91 2:e843c1b0b25c 1146 tmp[i][j] = v[i][*el];
jsoh91 2:e843c1b0b25c 1147 }
jsoh91 2:e843c1b0b25c 1148 j++;
jsoh91 2:e843c1b0b25c 1149 }
jsoh91 2:e843c1b0b25c 1150
jsoh91 2:e843c1b0b25c 1151 return tmp;
jsoh91 2:e843c1b0b25c 1152 }
jsoh91 2:e843c1b0b25c 1153
jsoh91 2:e843c1b0b25c 1154 template <typename T>
jsoh91 2:e843c1b0b25c 1155 inline Matrix<T> Matrix<T>::extract(const std::set<unsigned int>& r_indexes, const std::set<unsigned int>& c_indexes) const
jsoh91 2:e843c1b0b25c 1156 {
jsoh91 2:e843c1b0b25c 1157 Matrix<T> tmp(r_indexes.size(), c_indexes.size());
jsoh91 2:e843c1b0b25c 1158 register unsigned int i = 0, j;
jsoh91 2:e843c1b0b25c 1159
jsoh91 2:e843c1b0b25c 1160 for (std::set<unsigned int>::const_iterator r_el = r_indexes.begin(); r_el != r_indexes.end(); r_el++)
jsoh91 2:e843c1b0b25c 1161 {
jsoh91 2:e843c1b0b25c 1162 if (*r_el >= n)
jsoh91 2:e843c1b0b25c 1163 throw std::logic_error("Error extracting submatrix: the indexes are out of matrix bounds");
jsoh91 2:e843c1b0b25c 1164 j = 0;
jsoh91 2:e843c1b0b25c 1165 for (std::set<unsigned int>::const_iterator c_el = c_indexes.begin(); c_el != c_indexes.end(); c_el++)
jsoh91 2:e843c1b0b25c 1166 {
jsoh91 2:e843c1b0b25c 1167 if (*c_el >= m)
jsoh91 2:e843c1b0b25c 1168 throw std::logic_error("Error extracting rows: the indexes are out of matrix bounds");
jsoh91 2:e843c1b0b25c 1169 tmp[i][j] = v[*r_el][*c_el];
jsoh91 2:e843c1b0b25c 1170 j++;
jsoh91 2:e843c1b0b25c 1171 }
jsoh91 2:e843c1b0b25c 1172 i++;
jsoh91 2:e843c1b0b25c 1173 }
jsoh91 2:e843c1b0b25c 1174
jsoh91 2:e843c1b0b25c 1175 return tmp;
jsoh91 2:e843c1b0b25c 1176 }
jsoh91 2:e843c1b0b25c 1177
jsoh91 2:e843c1b0b25c 1178 template <typename T>
jsoh91 2:e843c1b0b25c 1179 inline void Matrix<T>::setRow(unsigned int i, const Vector<T>& a)
jsoh91 2:e843c1b0b25c 1180 {
jsoh91 2:e843c1b0b25c 1181 if (i >= n)
jsoh91 2:e843c1b0b25c 1182 throw std::logic_error("Error in setRow: trying to set a row out of matrix bounds");
jsoh91 2:e843c1b0b25c 1183 if (this->m != a.size())
jsoh91 2:e843c1b0b25c 1184 throw std::logic_error("Error setting matrix row: ranges are not compatible");
jsoh91 2:e843c1b0b25c 1185 for (unsigned int j = 0; j < ncols(); j++)
jsoh91 2:e843c1b0b25c 1186 v[i][j] = a[j];
jsoh91 2:e843c1b0b25c 1187 }
jsoh91 2:e843c1b0b25c 1188
jsoh91 2:e843c1b0b25c 1189 template <typename T>
jsoh91 2:e843c1b0b25c 1190 inline void Matrix<T>::setRow(unsigned int i, const Matrix<T>& a)
jsoh91 2:e843c1b0b25c 1191 {
jsoh91 2:e843c1b0b25c 1192 if (i >= n)
jsoh91 2:e843c1b0b25c 1193 throw std::logic_error("Error in setRow: trying to set a row out of matrix bounds");
jsoh91 2:e843c1b0b25c 1194 if (this->m != a.ncols())
jsoh91 2:e843c1b0b25c 1195 throw std::logic_error("Error setting matrix column: ranges are not compatible");
jsoh91 2:e843c1b0b25c 1196 if (a.nrows() != 1)
jsoh91 2:e843c1b0b25c 1197 throw std::logic_error("Error setting matrix column with a non-row matrix");
jsoh91 2:e843c1b0b25c 1198 for (unsigned int j = 0; j < ncols(); j++)
jsoh91 2:e843c1b0b25c 1199 v[i][j] = a[0][j];
jsoh91 2:e843c1b0b25c 1200 }
jsoh91 2:e843c1b0b25c 1201
jsoh91 2:e843c1b0b25c 1202 template <typename T>
jsoh91 2:e843c1b0b25c 1203 inline void Matrix<T>::setRows(const std::set<unsigned int>& indexes, const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 1204 {
jsoh91 2:e843c1b0b25c 1205 unsigned int i = 0;
jsoh91 2:e843c1b0b25c 1206
jsoh91 2:e843c1b0b25c 1207 if (indexes.size() != m.nrows() || this->m != m.ncols())
jsoh91 2:e843c1b0b25c 1208 throw std::logic_error("Error setting matrix rows: ranges are not compatible");
jsoh91 2:e843c1b0b25c 1209 for (std::set<unsigned int>::const_iterator el = indexes.begin(); el != indexes.end(); el++)
jsoh91 2:e843c1b0b25c 1210 {
jsoh91 2:e843c1b0b25c 1211 for (unsigned int j = 0; j < ncols(); j++)
jsoh91 2:e843c1b0b25c 1212 {
jsoh91 2:e843c1b0b25c 1213 if (*el >= n)
jsoh91 2:e843c1b0b25c 1214 throw std::logic_error("Error in setRows: trying to set a row out of matrix bounds");
jsoh91 2:e843c1b0b25c 1215 v[*el][j] = m[i][j];
jsoh91 2:e843c1b0b25c 1216 }
jsoh91 2:e843c1b0b25c 1217 i++;
jsoh91 2:e843c1b0b25c 1218 }
jsoh91 2:e843c1b0b25c 1219 }
jsoh91 2:e843c1b0b25c 1220
jsoh91 2:e843c1b0b25c 1221 template <typename T>
jsoh91 2:e843c1b0b25c 1222 inline void Matrix<T>::setColumn(unsigned int j, const Vector<T>& a)
jsoh91 2:e843c1b0b25c 1223 {
jsoh91 2:e843c1b0b25c 1224 if (j >= m)
jsoh91 2:e843c1b0b25c 1225 throw std::logic_error("Error in setColumn: trying to set a column out of matrix bounds");
jsoh91 2:e843c1b0b25c 1226 if (this->n != a.size())
jsoh91 2:e843c1b0b25c 1227 throw std::logic_error("Error setting matrix column: ranges are not compatible");
jsoh91 2:e843c1b0b25c 1228 for (unsigned int i = 0; i < nrows(); i++)
jsoh91 2:e843c1b0b25c 1229 v[i][j] = a[i];
jsoh91 2:e843c1b0b25c 1230 }
jsoh91 2:e843c1b0b25c 1231
jsoh91 2:e843c1b0b25c 1232 template <typename T>
jsoh91 2:e843c1b0b25c 1233 inline void Matrix<T>::setColumn(unsigned int j, const Matrix<T>& a)
jsoh91 2:e843c1b0b25c 1234 {
jsoh91 2:e843c1b0b25c 1235 if (j >= m)
jsoh91 2:e843c1b0b25c 1236 throw std::logic_error("Error in setColumn: trying to set a column out of matrix bounds");
jsoh91 2:e843c1b0b25c 1237 if (this->n != a.nrows())
jsoh91 2:e843c1b0b25c 1238 throw std::logic_error("Error setting matrix column: ranges are not compatible");
jsoh91 2:e843c1b0b25c 1239 if (a.ncols() != 1)
jsoh91 2:e843c1b0b25c 1240 throw std::logic_error("Error setting matrix column with a non-column matrix");
jsoh91 2:e843c1b0b25c 1241 for (unsigned int i = 0; i < nrows(); i++)
jsoh91 2:e843c1b0b25c 1242 v[i][j] = a[i][0];
jsoh91 2:e843c1b0b25c 1243 }
jsoh91 2:e843c1b0b25c 1244
jsoh91 2:e843c1b0b25c 1245
jsoh91 2:e843c1b0b25c 1246 template <typename T>
jsoh91 2:e843c1b0b25c 1247 inline void Matrix<T>::setColumns(const std::set<unsigned int>& indexes, const Matrix<T>& a)
jsoh91 2:e843c1b0b25c 1248 {
jsoh91 2:e843c1b0b25c 1249 unsigned int j = 0;
jsoh91 2:e843c1b0b25c 1250
jsoh91 2:e843c1b0b25c 1251 if (indexes.size() != a.ncols() || this->n != a.nrows())
jsoh91 2:e843c1b0b25c 1252 throw std::logic_error("Error setting matrix columns: ranges are not compatible");
jsoh91 2:e843c1b0b25c 1253 for (std::set<unsigned int>::const_iterator el = indexes.begin(); el != indexes.end(); el++)
jsoh91 2:e843c1b0b25c 1254 {
jsoh91 2:e843c1b0b25c 1255 for (unsigned int i = 0; i < nrows(); i++)
jsoh91 2:e843c1b0b25c 1256 {
jsoh91 2:e843c1b0b25c 1257 if (*el >= m)
jsoh91 2:e843c1b0b25c 1258 throw std::logic_error("Error in setColumns: trying to set a column out of matrix bounds");
jsoh91 2:e843c1b0b25c 1259 v[i][*el] = a[i][j];
jsoh91 2:e843c1b0b25c 1260 }
jsoh91 2:e843c1b0b25c 1261 j++;
jsoh91 2:e843c1b0b25c 1262 }
jsoh91 2:e843c1b0b25c 1263 }
jsoh91 2:e843c1b0b25c 1264
jsoh91 2:e843c1b0b25c 1265 template <typename T>
jsoh91 2:e843c1b0b25c 1266 inline void Matrix<T>::set(const std::set<unsigned int>& r_indexes, const std::set<unsigned int>& c_indexes, const Matrix<T>& a)
jsoh91 2:e843c1b0b25c 1267 {
jsoh91 2:e843c1b0b25c 1268 unsigned int i = 0, j;
jsoh91 2:e843c1b0b25c 1269 if (c_indexes.size() != a.ncols() || r_indexes.size() != a.nrows())
jsoh91 2:e843c1b0b25c 1270 throw std::logic_error("Error setting matrix elements: ranges are not compatible");
jsoh91 2:e843c1b0b25c 1271
jsoh91 2:e843c1b0b25c 1272 for (std::set<unsigned int>::const_iterator r_el = r_indexes.begin(); r_el != r_indexes.end(); r_el++)
jsoh91 2:e843c1b0b25c 1273 {
jsoh91 2:e843c1b0b25c 1274 if (*r_el >= n)
jsoh91 2:e843c1b0b25c 1275 throw std::logic_error("Error in set: trying to set a row out of matrix bounds");
jsoh91 2:e843c1b0b25c 1276 j = 0;
jsoh91 2:e843c1b0b25c 1277 for (std::set<unsigned int>::const_iterator c_el = c_indexes.begin(); c_el != c_indexes.end(); c_el++)
jsoh91 2:e843c1b0b25c 1278 {
jsoh91 2:e843c1b0b25c 1279 if (*c_el >= m)
jsoh91 2:e843c1b0b25c 1280 throw std::logic_error("Error in set: trying to set a column out of matrix bounds");
jsoh91 2:e843c1b0b25c 1281 v[*r_el][*c_el] = a[i][j];
jsoh91 2:e843c1b0b25c 1282 j++;
jsoh91 2:e843c1b0b25c 1283 }
jsoh91 2:e843c1b0b25c 1284 i++;
jsoh91 2:e843c1b0b25c 1285 }
jsoh91 2:e843c1b0b25c 1286 }
jsoh91 2:e843c1b0b25c 1287
jsoh91 2:e843c1b0b25c 1288 template <typename T>
jsoh91 2:e843c1b0b25c 1289 inline void Matrix<T>::set(const T* a, unsigned int n, unsigned int m)
jsoh91 2:e843c1b0b25c 1290 {
jsoh91 2:e843c1b0b25c 1291 if (this->n != n || this->m != m)
jsoh91 2:e843c1b0b25c 1292 resize(n, m);
jsoh91 2:e843c1b0b25c 1293 unsigned int k = 0;
jsoh91 2:e843c1b0b25c 1294 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1295 for (unsigned int j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1296 v[i][j] = a[k++];
jsoh91 2:e843c1b0b25c 1297 }
jsoh91 2:e843c1b0b25c 1298
jsoh91 2:e843c1b0b25c 1299
jsoh91 2:e843c1b0b25c 1300 template <typename T>
jsoh91 2:e843c1b0b25c 1301 Matrix<T> operator+(const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1302 {
jsoh91 2:e843c1b0b25c 1303 return rhs;
jsoh91 2:e843c1b0b25c 1304 }
jsoh91 2:e843c1b0b25c 1305
jsoh91 2:e843c1b0b25c 1306 template <typename T>
jsoh91 2:e843c1b0b25c 1307 Matrix<T> operator+(const Matrix<T>& lhs, const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1308 {
jsoh91 2:e843c1b0b25c 1309 if (lhs.ncols() != rhs.ncols() || lhs.nrows() != rhs.nrows())
jsoh91 2:e843c1b0b25c 1310 throw std::logic_error("Operator+: matrices have different sizes");
jsoh91 2:e843c1b0b25c 1311 Matrix<T> tmp(lhs.nrows(), lhs.ncols());
jsoh91 2:e843c1b0b25c 1312 for (unsigned int i = 0; i < lhs.nrows(); i++)
jsoh91 2:e843c1b0b25c 1313 for (unsigned int j = 0; j < lhs.ncols(); j++)
jsoh91 2:e843c1b0b25c 1314 tmp[i][j] = lhs[i][j] + rhs[i][j];
jsoh91 2:e843c1b0b25c 1315
jsoh91 2:e843c1b0b25c 1316 return tmp;
jsoh91 2:e843c1b0b25c 1317 }
jsoh91 2:e843c1b0b25c 1318
jsoh91 2:e843c1b0b25c 1319 template <typename T>
jsoh91 2:e843c1b0b25c 1320 Matrix<T> operator+(const Matrix<T>& lhs, const T& a)
jsoh91 2:e843c1b0b25c 1321 {
jsoh91 2:e843c1b0b25c 1322 Matrix<T> tmp(lhs.nrows(), lhs.ncols());
jsoh91 2:e843c1b0b25c 1323 for (unsigned int i = 0; i < lhs.nrows(); i++)
jsoh91 2:e843c1b0b25c 1324 for (unsigned int j = 0; j < lhs.ncols(); j++)
jsoh91 2:e843c1b0b25c 1325 tmp[i][j] = lhs[i][j] + a;
jsoh91 2:e843c1b0b25c 1326
jsoh91 2:e843c1b0b25c 1327 return tmp;
jsoh91 2:e843c1b0b25c 1328 }
jsoh91 2:e843c1b0b25c 1329
jsoh91 2:e843c1b0b25c 1330 template <typename T>
jsoh91 2:e843c1b0b25c 1331 Matrix<T> operator+(const T& a, const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1332 {
jsoh91 2:e843c1b0b25c 1333 Matrix<T> tmp(rhs.nrows(), rhs.ncols());
jsoh91 2:e843c1b0b25c 1334 for (unsigned int i = 0; i < rhs.nrows(); i++)
jsoh91 2:e843c1b0b25c 1335 for (unsigned int j = 0; j < rhs.ncols(); j++)
jsoh91 2:e843c1b0b25c 1336 tmp[i][j] = a + rhs[i][j];
jsoh91 2:e843c1b0b25c 1337
jsoh91 2:e843c1b0b25c 1338 return tmp;
jsoh91 2:e843c1b0b25c 1339 }
jsoh91 2:e843c1b0b25c 1340
jsoh91 2:e843c1b0b25c 1341 template <typename T>
jsoh91 2:e843c1b0b25c 1342 inline Matrix<T>& Matrix<T>::operator+=(const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1343 {
jsoh91 2:e843c1b0b25c 1344 if (m != rhs.ncols() || n != rhs.nrows())
jsoh91 2:e843c1b0b25c 1345 throw std::logic_error("Operator+=: matrices have different sizes");
jsoh91 2:e843c1b0b25c 1346 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1347 for (unsigned int j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1348 v[i][j] += rhs[i][j];
jsoh91 2:e843c1b0b25c 1349
jsoh91 2:e843c1b0b25c 1350 return *this;
jsoh91 2:e843c1b0b25c 1351 }
jsoh91 2:e843c1b0b25c 1352
jsoh91 2:e843c1b0b25c 1353 template <typename T>
jsoh91 2:e843c1b0b25c 1354 inline Matrix<T>& Matrix<T>::operator+=(const T& a)
jsoh91 2:e843c1b0b25c 1355 {
jsoh91 2:e843c1b0b25c 1356 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1357 for (unsigned int j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1358 v[i][j] += a;
jsoh91 2:e843c1b0b25c 1359
jsoh91 2:e843c1b0b25c 1360 return *this;
jsoh91 2:e843c1b0b25c 1361 }
jsoh91 2:e843c1b0b25c 1362
jsoh91 2:e843c1b0b25c 1363 template <typename T>
jsoh91 2:e843c1b0b25c 1364 Matrix<T> operator-(const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1365 {
jsoh91 2:e843c1b0b25c 1366 return (T)(-1) * rhs;
jsoh91 2:e843c1b0b25c 1367 }
jsoh91 2:e843c1b0b25c 1368
jsoh91 2:e843c1b0b25c 1369 template <typename T>
jsoh91 2:e843c1b0b25c 1370 Matrix<T> operator-(const Matrix<T>& lhs, const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1371 {
jsoh91 2:e843c1b0b25c 1372 if (lhs.ncols() != rhs.ncols() || lhs.nrows() != rhs.nrows())
jsoh91 2:e843c1b0b25c 1373 throw std::logic_error("Operator-: matrices have different sizes");
jsoh91 2:e843c1b0b25c 1374 Matrix<T> tmp(lhs.nrows(), lhs.ncols());
jsoh91 2:e843c1b0b25c 1375 for (unsigned int i = 0; i < lhs.nrows(); i++)
jsoh91 2:e843c1b0b25c 1376 for (unsigned int j = 0; j < lhs.ncols(); j++)
jsoh91 2:e843c1b0b25c 1377 tmp[i][j] = lhs[i][j] - rhs[i][j];
jsoh91 2:e843c1b0b25c 1378
jsoh91 2:e843c1b0b25c 1379 return tmp;
jsoh91 2:e843c1b0b25c 1380 }
jsoh91 2:e843c1b0b25c 1381
jsoh91 2:e843c1b0b25c 1382 template <typename T>
jsoh91 2:e843c1b0b25c 1383 Matrix<T> operator-(const Matrix<T>& lhs, const T& a)
jsoh91 2:e843c1b0b25c 1384 {
jsoh91 2:e843c1b0b25c 1385 Matrix<T> tmp(lhs.nrows(), lhs.ncols());
jsoh91 2:e843c1b0b25c 1386 for (unsigned int i = 0; i < lhs.nrows(); i++)
jsoh91 2:e843c1b0b25c 1387 for (unsigned int j = 0; j < lhs.ncols(); j++)
jsoh91 2:e843c1b0b25c 1388 tmp[i][j] = lhs[i][j] - a;
jsoh91 2:e843c1b0b25c 1389
jsoh91 2:e843c1b0b25c 1390 return tmp;
jsoh91 2:e843c1b0b25c 1391 }
jsoh91 2:e843c1b0b25c 1392
jsoh91 2:e843c1b0b25c 1393 template <typename T>
jsoh91 2:e843c1b0b25c 1394 Matrix<T> operator-(const T& a, const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1395 {
jsoh91 2:e843c1b0b25c 1396 Matrix<T> tmp(rhs.nrows(), rhs.ncols());
jsoh91 2:e843c1b0b25c 1397 for (unsigned int i = 0; i < rhs.nrows(); i++)
jsoh91 2:e843c1b0b25c 1398 for (unsigned int j = 0; j < rhs.ncols(); j++)
jsoh91 2:e843c1b0b25c 1399 tmp[i][j] = a - rhs[i][j];
jsoh91 2:e843c1b0b25c 1400
jsoh91 2:e843c1b0b25c 1401 return tmp;
jsoh91 2:e843c1b0b25c 1402 }
jsoh91 2:e843c1b0b25c 1403
jsoh91 2:e843c1b0b25c 1404 template <typename T>
jsoh91 2:e843c1b0b25c 1405 inline Matrix<T>& Matrix<T>::operator-=(const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1406 {
jsoh91 2:e843c1b0b25c 1407 if (m != rhs.ncols() || n != rhs.nrows())
jsoh91 2:e843c1b0b25c 1408 throw std::logic_error("Operator-=: matrices have different sizes");
jsoh91 2:e843c1b0b25c 1409 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1410 for (unsigned int j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1411 v[i][j] -= rhs[i][j];
jsoh91 2:e843c1b0b25c 1412
jsoh91 2:e843c1b0b25c 1413 return *this;
jsoh91 2:e843c1b0b25c 1414 }
jsoh91 2:e843c1b0b25c 1415
jsoh91 2:e843c1b0b25c 1416 template <typename T>
jsoh91 2:e843c1b0b25c 1417 inline Matrix<T>& Matrix<T>::operator-=(const T& a)
jsoh91 2:e843c1b0b25c 1418 {
jsoh91 2:e843c1b0b25c 1419 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1420 for (unsigned int j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1421 v[i][j] -= a;
jsoh91 2:e843c1b0b25c 1422
jsoh91 2:e843c1b0b25c 1423 return *this;
jsoh91 2:e843c1b0b25c 1424 }
jsoh91 2:e843c1b0b25c 1425
jsoh91 2:e843c1b0b25c 1426 template <typename T>
jsoh91 2:e843c1b0b25c 1427 Matrix<T> operator*(const Matrix<T>& lhs, const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1428 {
jsoh91 2:e843c1b0b25c 1429 if (lhs.ncols() != rhs.ncols() || lhs.nrows() != rhs.nrows())
jsoh91 2:e843c1b0b25c 1430 throw std::logic_error("Operator*: matrices have different sizes");
jsoh91 2:e843c1b0b25c 1431 Matrix<T> tmp(lhs.nrows(), lhs.ncols());
jsoh91 2:e843c1b0b25c 1432 for (unsigned int i = 0; i < lhs.nrows(); i++)
jsoh91 2:e843c1b0b25c 1433 for (unsigned int j = 0; j < lhs.ncols(); j++)
jsoh91 2:e843c1b0b25c 1434 tmp[i][j] = lhs[i][j] * rhs[i][j];
jsoh91 2:e843c1b0b25c 1435
jsoh91 2:e843c1b0b25c 1436 return tmp;
jsoh91 2:e843c1b0b25c 1437 }
jsoh91 2:e843c1b0b25c 1438
jsoh91 2:e843c1b0b25c 1439 template <typename T>
jsoh91 2:e843c1b0b25c 1440 Matrix<T> operator*(const Matrix<T>& lhs, const T& a)
jsoh91 2:e843c1b0b25c 1441 {
jsoh91 2:e843c1b0b25c 1442 Matrix<T> tmp(lhs.nrows(), lhs.ncols());
jsoh91 2:e843c1b0b25c 1443 for (unsigned int i = 0; i < lhs.nrows(); i++)
jsoh91 2:e843c1b0b25c 1444 for (unsigned int j = 0; j < lhs.ncols(); j++)
jsoh91 2:e843c1b0b25c 1445 tmp[i][j] = lhs[i][j] * a;
jsoh91 2:e843c1b0b25c 1446
jsoh91 2:e843c1b0b25c 1447 return tmp;
jsoh91 2:e843c1b0b25c 1448 }
jsoh91 2:e843c1b0b25c 1449
jsoh91 2:e843c1b0b25c 1450 template <typename T>
jsoh91 2:e843c1b0b25c 1451 Matrix<T> operator*(const T& a, const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1452 {
jsoh91 2:e843c1b0b25c 1453 Matrix<T> tmp(rhs.nrows(), rhs.ncols());
jsoh91 2:e843c1b0b25c 1454 for (unsigned int i = 0; i < rhs.nrows(); i++)
jsoh91 2:e843c1b0b25c 1455 for (unsigned int j = 0; j < rhs.ncols(); j++)
jsoh91 2:e843c1b0b25c 1456 tmp[i][j] = a * rhs[i][j];
jsoh91 2:e843c1b0b25c 1457
jsoh91 2:e843c1b0b25c 1458 return tmp;
jsoh91 2:e843c1b0b25c 1459 }
jsoh91 2:e843c1b0b25c 1460
jsoh91 2:e843c1b0b25c 1461 template <typename T>
jsoh91 2:e843c1b0b25c 1462 inline Matrix<T>& Matrix<T>::operator*=(const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1463 {
jsoh91 2:e843c1b0b25c 1464 if (m != rhs.ncols() || n != rhs.nrows())
jsoh91 2:e843c1b0b25c 1465 throw std::logic_error("Operator*=: matrices have different sizes");
jsoh91 2:e843c1b0b25c 1466 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1467 for (unsigned int j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1468 v[i][j] *= rhs[i][j];
jsoh91 2:e843c1b0b25c 1469
jsoh91 2:e843c1b0b25c 1470 return *this;
jsoh91 2:e843c1b0b25c 1471 }
jsoh91 2:e843c1b0b25c 1472
jsoh91 2:e843c1b0b25c 1473 template <typename T>
jsoh91 2:e843c1b0b25c 1474 inline Matrix<T>& Matrix<T>::operator*=(const T& a)
jsoh91 2:e843c1b0b25c 1475 {
jsoh91 2:e843c1b0b25c 1476 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1477 for (unsigned int j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1478 v[i][j] *= a;
jsoh91 2:e843c1b0b25c 1479
jsoh91 2:e843c1b0b25c 1480 return *this;
jsoh91 2:e843c1b0b25c 1481 }
jsoh91 2:e843c1b0b25c 1482
jsoh91 2:e843c1b0b25c 1483 template <typename T>
jsoh91 2:e843c1b0b25c 1484 Matrix<T> operator/(const Matrix<T>& lhs, const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1485 {
jsoh91 2:e843c1b0b25c 1486 if (lhs.ncols() != rhs.ncols() || lhs.nrows() != rhs.nrows())
jsoh91 2:e843c1b0b25c 1487 throw std::logic_error("Operator+: matrices have different sizes");
jsoh91 2:e843c1b0b25c 1488 Matrix<T> tmp(lhs.nrows(), lhs.ncols());
jsoh91 2:e843c1b0b25c 1489 for (unsigned int i = 0; i < lhs.nrows(); i++)
jsoh91 2:e843c1b0b25c 1490 for (unsigned int j = 0; j < lhs.ncols(); j++)
jsoh91 2:e843c1b0b25c 1491 tmp[i][j] = lhs[i][j] / rhs[i][j];
jsoh91 2:e843c1b0b25c 1492
jsoh91 2:e843c1b0b25c 1493 return tmp;
jsoh91 2:e843c1b0b25c 1494 }
jsoh91 2:e843c1b0b25c 1495
jsoh91 2:e843c1b0b25c 1496 template <typename T>
jsoh91 2:e843c1b0b25c 1497 Matrix<T> operator/(const Matrix<T>& lhs, const T& a)
jsoh91 2:e843c1b0b25c 1498 {
jsoh91 2:e843c1b0b25c 1499 Matrix<T> tmp(lhs.nrows(), lhs.ncols());
jsoh91 2:e843c1b0b25c 1500 for (unsigned int i = 0; i < lhs.nrows(); i++)
jsoh91 2:e843c1b0b25c 1501 for (unsigned int j = 0; j < lhs.ncols(); j++)
jsoh91 2:e843c1b0b25c 1502 tmp[i][j] = lhs[i][j] / a;
jsoh91 2:e843c1b0b25c 1503
jsoh91 2:e843c1b0b25c 1504 return tmp;
jsoh91 2:e843c1b0b25c 1505 }
jsoh91 2:e843c1b0b25c 1506
jsoh91 2:e843c1b0b25c 1507 template <typename T>
jsoh91 2:e843c1b0b25c 1508 Matrix<T> operator/(const T& a, const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1509 {
jsoh91 2:e843c1b0b25c 1510 Matrix<T> tmp(rhs.nrows(), rhs.ncols());
jsoh91 2:e843c1b0b25c 1511 for (unsigned int i = 0; i < rhs.nrows(); i++)
jsoh91 2:e843c1b0b25c 1512 for (unsigned int j = 0; j < rhs.ncols(); j++)
jsoh91 2:e843c1b0b25c 1513 tmp[i][j] = a / rhs[i][j];
jsoh91 2:e843c1b0b25c 1514
jsoh91 2:e843c1b0b25c 1515 return tmp;
jsoh91 2:e843c1b0b25c 1516 }
jsoh91 2:e843c1b0b25c 1517
jsoh91 2:e843c1b0b25c 1518 template <typename T>
jsoh91 2:e843c1b0b25c 1519 inline Matrix<T>& Matrix<T>::operator/=(const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1520 {
jsoh91 2:e843c1b0b25c 1521 if (m != rhs.ncols() || n != rhs.nrows())
jsoh91 2:e843c1b0b25c 1522 throw std::logic_error("Operator+=: matrices have different sizes");
jsoh91 2:e843c1b0b25c 1523 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1524 for (unsigned int j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1525 v[i][j] /= rhs[i][j];
jsoh91 2:e843c1b0b25c 1526
jsoh91 2:e843c1b0b25c 1527 return *this;
jsoh91 2:e843c1b0b25c 1528 }
jsoh91 2:e843c1b0b25c 1529
jsoh91 2:e843c1b0b25c 1530 template <typename T>
jsoh91 2:e843c1b0b25c 1531 inline Matrix<T>& Matrix<T>::operator/=(const T& a)
jsoh91 2:e843c1b0b25c 1532 {
jsoh91 2:e843c1b0b25c 1533 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1534 for (unsigned int j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1535 v[i][j] /= a;
jsoh91 2:e843c1b0b25c 1536
jsoh91 2:e843c1b0b25c 1537 return *this;
jsoh91 2:e843c1b0b25c 1538 }
jsoh91 2:e843c1b0b25c 1539
jsoh91 2:e843c1b0b25c 1540 template <typename T>
jsoh91 2:e843c1b0b25c 1541 Matrix<T> operator^(const Matrix<T>& lhs, const T& a)
jsoh91 2:e843c1b0b25c 1542 {
jsoh91 2:e843c1b0b25c 1543 Matrix<T> tmp(lhs.nrows(), lhs.ncols());
jsoh91 2:e843c1b0b25c 1544 for (unsigned int i = 0; i < lhs.nrows(); i++)
jsoh91 2:e843c1b0b25c 1545 for (unsigned int j = 0; j < lhs.ncols(); j++)
jsoh91 2:e843c1b0b25c 1546 tmp[i][j] = pow(lhs[i][j], a);
jsoh91 2:e843c1b0b25c 1547
jsoh91 2:e843c1b0b25c 1548 return tmp;
jsoh91 2:e843c1b0b25c 1549 }
jsoh91 2:e843c1b0b25c 1550
jsoh91 2:e843c1b0b25c 1551 template <typename T>
jsoh91 2:e843c1b0b25c 1552 inline Matrix<T>& Matrix<T>::operator^=(const Matrix<T>& rhs)
jsoh91 2:e843c1b0b25c 1553 {
jsoh91 2:e843c1b0b25c 1554 if (m != rhs.ncols() || n != rhs.nrows())
jsoh91 2:e843c1b0b25c 1555 throw std::logic_error("Operator^=: matrices have different sizes");
jsoh91 2:e843c1b0b25c 1556 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1557 for (unsigned int j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1558 v[i][j] = pow(v[i][j], rhs[i][j]);
jsoh91 2:e843c1b0b25c 1559
jsoh91 2:e843c1b0b25c 1560 return *this;
jsoh91 2:e843c1b0b25c 1561 }
jsoh91 2:e843c1b0b25c 1562
jsoh91 2:e843c1b0b25c 1563
jsoh91 2:e843c1b0b25c 1564 template <typename T>
jsoh91 2:e843c1b0b25c 1565 inline Matrix<T>& Matrix<T>::operator^=(const T& a)
jsoh91 2:e843c1b0b25c 1566 {
jsoh91 2:e843c1b0b25c 1567 for (unsigned int i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1568 for (unsigned int j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1569 v[i][j] = pow(v[i][j], a);
jsoh91 2:e843c1b0b25c 1570
jsoh91 2:e843c1b0b25c 1571 return *this;
jsoh91 2:e843c1b0b25c 1572 }
jsoh91 2:e843c1b0b25c 1573
jsoh91 2:e843c1b0b25c 1574 template <typename T>
jsoh91 2:e843c1b0b25c 1575 inline Matrix<T>::operator Vector<T>()
jsoh91 2:e843c1b0b25c 1576 {
jsoh91 2:e843c1b0b25c 1577 if (n > 1 && m > 1)
jsoh91 2:e843c1b0b25c 1578 throw std::logic_error("Error matrix cast to vector: trying to cast a multi-dimensional matrix");
jsoh91 2:e843c1b0b25c 1579 if (n == 1)
jsoh91 2:e843c1b0b25c 1580 return extractRow(0);
jsoh91 2:e843c1b0b25c 1581 else
jsoh91 2:e843c1b0b25c 1582 return extractColumn(0);
jsoh91 2:e843c1b0b25c 1583 }
jsoh91 2:e843c1b0b25c 1584
jsoh91 2:e843c1b0b25c 1585 template <typename T>
jsoh91 2:e843c1b0b25c 1586 inline bool operator==(const Matrix<T>& a, const Matrix<T>& b)
jsoh91 2:e843c1b0b25c 1587 {
jsoh91 2:e843c1b0b25c 1588 if (a.nrows() != b.nrows() || a.ncols() != b.ncols())
jsoh91 2:e843c1b0b25c 1589 throw std::logic_error("Matrices of different size are not confrontable");
jsoh91 2:e843c1b0b25c 1590 for (unsigned i = 0; i < a.nrows(); i++)
jsoh91 2:e843c1b0b25c 1591 for (unsigned j = 0; j < a.ncols(); j++)
jsoh91 2:e843c1b0b25c 1592 if (a[i][j] != b[i][j])
jsoh91 2:e843c1b0b25c 1593 return false;
jsoh91 2:e843c1b0b25c 1594 return true;
jsoh91 2:e843c1b0b25c 1595 }
jsoh91 2:e843c1b0b25c 1596
jsoh91 2:e843c1b0b25c 1597 template <typename T>
jsoh91 2:e843c1b0b25c 1598 inline bool operator!=(const Matrix<T>& a, const Matrix<T>& b)
jsoh91 2:e843c1b0b25c 1599 {
jsoh91 2:e843c1b0b25c 1600 if (a.nrows() != b.nrows() || a.ncols() != b.ncols())
jsoh91 2:e843c1b0b25c 1601 throw std::logic_error("Matrices of different size are not confrontable");
jsoh91 2:e843c1b0b25c 1602 for (unsigned i = 0; i < a.nrows(); i++)
jsoh91 2:e843c1b0b25c 1603 for (unsigned j = 0; j < a.ncols(); j++)
jsoh91 2:e843c1b0b25c 1604 if (a[i][j] != b[i][j])
jsoh91 2:e843c1b0b25c 1605 return true;
jsoh91 2:e843c1b0b25c 1606 return false;
jsoh91 2:e843c1b0b25c 1607 }
jsoh91 2:e843c1b0b25c 1608
jsoh91 2:e843c1b0b25c 1609
jsoh91 2:e843c1b0b25c 1610
jsoh91 2:e843c1b0b25c 1611 /**
jsoh91 2:e843c1b0b25c 1612 Input/Output
jsoh91 2:e843c1b0b25c 1613 */
jsoh91 2:e843c1b0b25c 1614 template <typename T>
jsoh91 2:e843c1b0b25c 1615 std::ostream& operator<<(std::ostream& os, const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 1616 {
jsoh91 2:e843c1b0b25c 1617 os << std::endl << m.nrows() << " " << m.ncols() << std::endl;
jsoh91 2:e843c1b0b25c 1618 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 1619 {
jsoh91 2:e843c1b0b25c 1620 for (unsigned int j = 0; j < m.ncols() - 1; j++)
jsoh91 2:e843c1b0b25c 1621 os << std::setw(20) << std::setprecision(16) << m[i][j] << ", ";
jsoh91 2:e843c1b0b25c 1622 os << std::setw(20) << std::setprecision(16) << m[i][m.ncols() - 1] << std::endl;
jsoh91 2:e843c1b0b25c 1623 }
jsoh91 2:e843c1b0b25c 1624
jsoh91 2:e843c1b0b25c 1625 return os;
jsoh91 2:e843c1b0b25c 1626 }
jsoh91 2:e843c1b0b25c 1627
jsoh91 2:e843c1b0b25c 1628 template <typename T>
jsoh91 2:e843c1b0b25c 1629 std::istream& operator>>(std::istream& is, Matrix<T>& m)
jsoh91 2:e843c1b0b25c 1630 {
jsoh91 2:e843c1b0b25c 1631 int rows, cols;
jsoh91 2:e843c1b0b25c 1632 char comma;
jsoh91 2:e843c1b0b25c 1633 is >> rows >> cols;
jsoh91 2:e843c1b0b25c 1634 m.resize(rows, cols);
jsoh91 2:e843c1b0b25c 1635 for (unsigned int i = 0; i < rows; i++)
jsoh91 2:e843c1b0b25c 1636 for (unsigned int j = 0; j < cols; j++)
jsoh91 2:e843c1b0b25c 1637 is >> m[i][j] >> comma;
jsoh91 2:e843c1b0b25c 1638
jsoh91 2:e843c1b0b25c 1639 return is;
jsoh91 2:e843c1b0b25c 1640 }
jsoh91 2:e843c1b0b25c 1641
jsoh91 2:e843c1b0b25c 1642 template <typename T>
jsoh91 2:e843c1b0b25c 1643 T sign(const T& v)
jsoh91 2:e843c1b0b25c 1644 {
jsoh91 2:e843c1b0b25c 1645 if (v >= (T)0.0)
jsoh91 2:e843c1b0b25c 1646 return (T)1.0;
jsoh91 2:e843c1b0b25c 1647 else
jsoh91 2:e843c1b0b25c 1648 return (T)-1.0;
jsoh91 2:e843c1b0b25c 1649 }
jsoh91 2:e843c1b0b25c 1650
jsoh91 2:e843c1b0b25c 1651 template <typename T>
jsoh91 2:e843c1b0b25c 1652 T dist(const T& a, const T& b)
jsoh91 2:e843c1b0b25c 1653 {
jsoh91 2:e843c1b0b25c 1654 T abs_a = (T)fabs(a), abs_b = (T)fabs(b);
jsoh91 2:e843c1b0b25c 1655 if (abs_a > abs_b)
jsoh91 2:e843c1b0b25c 1656 return abs_a * sqrt((T)1.0 + (abs_b / abs_a) * (abs_b / abs_a));
jsoh91 2:e843c1b0b25c 1657 else
jsoh91 2:e843c1b0b25c 1658 return (abs_b == (T)0.0 ? (T)0.0 : abs_b * sqrt((T)1.0 + (abs_a / abs_b) * (abs_a / abs_b)));
jsoh91 2:e843c1b0b25c 1659 }
jsoh91 2:e843c1b0b25c 1660
jsoh91 2:e843c1b0b25c 1661 template <typename T>
jsoh91 2:e843c1b0b25c 1662 void svd(const Matrix<T>& A, Matrix<T>& U, Vector<T>& W, Matrix<T>& V)
jsoh91 2:e843c1b0b25c 1663 {
jsoh91 2:e843c1b0b25c 1664 int m = A.nrows(), n = A.ncols(), i, j, k, l, jj, nm;
jsoh91 2:e843c1b0b25c 1665 const unsigned int max_its = 30;
jsoh91 2:e843c1b0b25c 1666 bool flag;
jsoh91 2:e843c1b0b25c 1667 Vector<T> rv1(n);
jsoh91 2:e843c1b0b25c 1668 U = A;
jsoh91 2:e843c1b0b25c 1669 W.resize(n);
jsoh91 2:e843c1b0b25c 1670 V.resize(n, n);
jsoh91 2:e843c1b0b25c 1671 T anorm, c, f, g, h, s, scale, x, y, z;
jsoh91 2:e843c1b0b25c 1672 g = scale = anorm = (T)0.0;
jsoh91 2:e843c1b0b25c 1673
jsoh91 2:e843c1b0b25c 1674 // Householder reduction to bidiagonal form
jsoh91 2:e843c1b0b25c 1675 for (i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1676 {
jsoh91 2:e843c1b0b25c 1677 l = i + 1;
jsoh91 2:e843c1b0b25c 1678 rv1[i] = scale * g;
jsoh91 2:e843c1b0b25c 1679 g = s = scale = (T)0.0;
jsoh91 2:e843c1b0b25c 1680 if (i < m)
jsoh91 2:e843c1b0b25c 1681 {
jsoh91 2:e843c1b0b25c 1682 for (k = i; k < m; k++)
jsoh91 2:e843c1b0b25c 1683 scale += fabs(U[k][i]);
jsoh91 2:e843c1b0b25c 1684 if (scale != (T)0.0)
jsoh91 2:e843c1b0b25c 1685 {
jsoh91 2:e843c1b0b25c 1686 for (k = i; k < m; k++)
jsoh91 2:e843c1b0b25c 1687 {
jsoh91 2:e843c1b0b25c 1688 U[k][i] /= scale;
jsoh91 2:e843c1b0b25c 1689 s += U[k][i] * U[k][i];
jsoh91 2:e843c1b0b25c 1690 }
jsoh91 2:e843c1b0b25c 1691 f = U[i][i];
jsoh91 2:e843c1b0b25c 1692 g = -sign(f) * sqrt(s);
jsoh91 2:e843c1b0b25c 1693 h = f * g - s;
jsoh91 2:e843c1b0b25c 1694 U[i][i] = f - g;
jsoh91 2:e843c1b0b25c 1695 for (j = l; j < n; j++)
jsoh91 2:e843c1b0b25c 1696 {
jsoh91 2:e843c1b0b25c 1697 s = (T)0.0;
jsoh91 2:e843c1b0b25c 1698 for (k = i; k < m; k++)
jsoh91 2:e843c1b0b25c 1699 s += U[k][i] * U[k][j];
jsoh91 2:e843c1b0b25c 1700 f = s / h;
jsoh91 2:e843c1b0b25c 1701 for (k = i; k < m; k++)
jsoh91 2:e843c1b0b25c 1702 U[k][j] += f * U[k][i];
jsoh91 2:e843c1b0b25c 1703 }
jsoh91 2:e843c1b0b25c 1704 for (k = i; k < m; k++)
jsoh91 2:e843c1b0b25c 1705 U[k][i] *= scale;
jsoh91 2:e843c1b0b25c 1706 }
jsoh91 2:e843c1b0b25c 1707 }
jsoh91 2:e843c1b0b25c 1708 W[i] = scale * g;
jsoh91 2:e843c1b0b25c 1709 g = s = scale = (T)0.0;
jsoh91 2:e843c1b0b25c 1710 if (i < m && i != n - 1)
jsoh91 2:e843c1b0b25c 1711 {
jsoh91 2:e843c1b0b25c 1712 for (k = l; k < n; k++)
jsoh91 2:e843c1b0b25c 1713 scale += fabs(U[i][k]);
jsoh91 2:e843c1b0b25c 1714 if (scale != (T)0.0)
jsoh91 2:e843c1b0b25c 1715 {
jsoh91 2:e843c1b0b25c 1716 for (k = l; k < n; k++)
jsoh91 2:e843c1b0b25c 1717 {
jsoh91 2:e843c1b0b25c 1718 U[i][k] /= scale;
jsoh91 2:e843c1b0b25c 1719 s += U[i][k] * U[i][k];
jsoh91 2:e843c1b0b25c 1720 }
jsoh91 2:e843c1b0b25c 1721 f = U[i][l];
jsoh91 2:e843c1b0b25c 1722 g = -sign(f) * sqrt(s);
jsoh91 2:e843c1b0b25c 1723 h = f * g - s;
jsoh91 2:e843c1b0b25c 1724 U[i][l] = f - g;
jsoh91 2:e843c1b0b25c 1725 for (k = l; k <n; k++)
jsoh91 2:e843c1b0b25c 1726 rv1[k] = U[i][k] / h;
jsoh91 2:e843c1b0b25c 1727 for (j = l; j < m; j++)
jsoh91 2:e843c1b0b25c 1728 {
jsoh91 2:e843c1b0b25c 1729 s = (T)0.0;
jsoh91 2:e843c1b0b25c 1730 for (k = l; k < n; k++)
jsoh91 2:e843c1b0b25c 1731 s += U[j][k] * U[i][k];
jsoh91 2:e843c1b0b25c 1732 for (k = l; k < n; k++)
jsoh91 2:e843c1b0b25c 1733 U[j][k] += s * rv1[k];
jsoh91 2:e843c1b0b25c 1734 }
jsoh91 2:e843c1b0b25c 1735 for (k = l; k < n; k++)
jsoh91 2:e843c1b0b25c 1736 U[i][k] *= scale;
jsoh91 2:e843c1b0b25c 1737 }
jsoh91 2:e843c1b0b25c 1738 }
jsoh91 2:e843c1b0b25c 1739 anorm = std::max(anorm, fabs(W[i]) + fabs(rv1[i]));
jsoh91 2:e843c1b0b25c 1740 }
jsoh91 2:e843c1b0b25c 1741 // Accumulation of right-hand transformations
jsoh91 2:e843c1b0b25c 1742 for (i = n - 1; i >= 0; i--)
jsoh91 2:e843c1b0b25c 1743 {
jsoh91 2:e843c1b0b25c 1744 if (i < n - 1)
jsoh91 2:e843c1b0b25c 1745 {
jsoh91 2:e843c1b0b25c 1746 if (g != (T)0.0)
jsoh91 2:e843c1b0b25c 1747 {
jsoh91 2:e843c1b0b25c 1748 for (j = l; j < n; j++)
jsoh91 2:e843c1b0b25c 1749 V[j][i] = (U[i][j] / U[i][l]) / g;
jsoh91 2:e843c1b0b25c 1750 for (j = l; j < n; j++)
jsoh91 2:e843c1b0b25c 1751 {
jsoh91 2:e843c1b0b25c 1752 s = (T)0.0;
jsoh91 2:e843c1b0b25c 1753 for (k = l; k < n; k++)
jsoh91 2:e843c1b0b25c 1754 s += U[i][k] * V[k][j];
jsoh91 2:e843c1b0b25c 1755 for (k = l; k < n; k++)
jsoh91 2:e843c1b0b25c 1756 V[k][j] += s * V[k][i];
jsoh91 2:e843c1b0b25c 1757 }
jsoh91 2:e843c1b0b25c 1758 }
jsoh91 2:e843c1b0b25c 1759 for (j = l; j < n; j++)
jsoh91 2:e843c1b0b25c 1760 V[i][j] = V[j][i] = (T)0.0;
jsoh91 2:e843c1b0b25c 1761 }
jsoh91 2:e843c1b0b25c 1762 V[i][i] = (T)1.0;
jsoh91 2:e843c1b0b25c 1763 g = rv1[i];
jsoh91 2:e843c1b0b25c 1764 l = i;
jsoh91 2:e843c1b0b25c 1765 }
jsoh91 2:e843c1b0b25c 1766 // Accumulation of left-hand transformations
jsoh91 2:e843c1b0b25c 1767 for (i = std::min(m, n) - 1; i >= 0; i--)
jsoh91 2:e843c1b0b25c 1768 {
jsoh91 2:e843c1b0b25c 1769 l = i + 1;
jsoh91 2:e843c1b0b25c 1770 g = W[i];
jsoh91 2:e843c1b0b25c 1771 for (j = l; j < n; j++)
jsoh91 2:e843c1b0b25c 1772 U[i][j] = (T)0.0;
jsoh91 2:e843c1b0b25c 1773 if (g != (T)0.0)
jsoh91 2:e843c1b0b25c 1774 {
jsoh91 2:e843c1b0b25c 1775 g = (T)1.0 / g;
jsoh91 2:e843c1b0b25c 1776 for (j = l; j < n; j++)
jsoh91 2:e843c1b0b25c 1777 {
jsoh91 2:e843c1b0b25c 1778 s = (T)0.0;
jsoh91 2:e843c1b0b25c 1779 for (k = l; k < m; k++)
jsoh91 2:e843c1b0b25c 1780 s += U[k][i] * U[k][j];
jsoh91 2:e843c1b0b25c 1781 f = (s / U[i][i]) * g;
jsoh91 2:e843c1b0b25c 1782 for (k = i; k < m; k++)
jsoh91 2:e843c1b0b25c 1783 U[k][j] += f * U[k][i];
jsoh91 2:e843c1b0b25c 1784 }
jsoh91 2:e843c1b0b25c 1785 for (j = i; j < m; j++)
jsoh91 2:e843c1b0b25c 1786 U[j][i] *= g;
jsoh91 2:e843c1b0b25c 1787 }
jsoh91 2:e843c1b0b25c 1788 else
jsoh91 2:e843c1b0b25c 1789 for (j = i; j < m; j++)
jsoh91 2:e843c1b0b25c 1790 U[j][i] = (T)0.0;
jsoh91 2:e843c1b0b25c 1791 U[i][i]++;
jsoh91 2:e843c1b0b25c 1792 }
jsoh91 2:e843c1b0b25c 1793 // Diagonalization of the bidiagonal form: loop over singular values, and over allowed iterations.
jsoh91 2:e843c1b0b25c 1794 for (k = n - 1; k >= 0; k--)
jsoh91 2:e843c1b0b25c 1795 {
jsoh91 2:e843c1b0b25c 1796 for (unsigned int its = 0; its < max_its; its++)
jsoh91 2:e843c1b0b25c 1797 {
jsoh91 2:e843c1b0b25c 1798 flag = true;
jsoh91 2:e843c1b0b25c 1799 for (l = k; l >= 0; l--) // FIXME: in NR it was l >= 1 but there subscripts start from one
jsoh91 2:e843c1b0b25c 1800 { // Test for splitting
jsoh91 2:e843c1b0b25c 1801 nm = l - 1; // Note that rV[0] is always zero
jsoh91 2:e843c1b0b25c 1802 if ((T)(fabs(rv1[l]) + anorm) == anorm)
jsoh91 2:e843c1b0b25c 1803 {
jsoh91 2:e843c1b0b25c 1804 flag = false;
jsoh91 2:e843c1b0b25c 1805 break;
jsoh91 2:e843c1b0b25c 1806 }
jsoh91 2:e843c1b0b25c 1807 if ((T)(fabs(W[nm]) + anorm) == anorm)
jsoh91 2:e843c1b0b25c 1808 break;
jsoh91 2:e843c1b0b25c 1809 }
jsoh91 2:e843c1b0b25c 1810 if (flag)
jsoh91 2:e843c1b0b25c 1811 {
jsoh91 2:e843c1b0b25c 1812 // Cancellation of rv1[l], if l > 0 FIXME: it was l > 1 in NR
jsoh91 2:e843c1b0b25c 1813 c = (T)0.0;
jsoh91 2:e843c1b0b25c 1814 s = (T)1.0;
jsoh91 2:e843c1b0b25c 1815 for (i = l; i <= k; i++)
jsoh91 2:e843c1b0b25c 1816 {
jsoh91 2:e843c1b0b25c 1817 f = s * rv1[i];
jsoh91 2:e843c1b0b25c 1818 rv1[i] *= c;
jsoh91 2:e843c1b0b25c 1819 if ((T)(fabs(f) + anorm) == anorm)
jsoh91 2:e843c1b0b25c 1820 break;
jsoh91 2:e843c1b0b25c 1821 g = W[i];
jsoh91 2:e843c1b0b25c 1822 h = dist(f, g);
jsoh91 2:e843c1b0b25c 1823 W[i] = h;
jsoh91 2:e843c1b0b25c 1824 h = (T)1.0 / h;
jsoh91 2:e843c1b0b25c 1825 c = g * h;
jsoh91 2:e843c1b0b25c 1826 s = -f * h;
jsoh91 2:e843c1b0b25c 1827 for (j = 0; j < m; j++)
jsoh91 2:e843c1b0b25c 1828 {
jsoh91 2:e843c1b0b25c 1829 y = U[j][nm];
jsoh91 2:e843c1b0b25c 1830 z = U[j][i];
jsoh91 2:e843c1b0b25c 1831 U[j][nm] = y * c + z * s;
jsoh91 2:e843c1b0b25c 1832 U[j][i] = z * c - y * s;
jsoh91 2:e843c1b0b25c 1833 }
jsoh91 2:e843c1b0b25c 1834 }
jsoh91 2:e843c1b0b25c 1835 }
jsoh91 2:e843c1b0b25c 1836 z = W[k];
jsoh91 2:e843c1b0b25c 1837 if (l == k)
jsoh91 2:e843c1b0b25c 1838 { // Convergence
jsoh91 2:e843c1b0b25c 1839 if (z < (T)0.0)
jsoh91 2:e843c1b0b25c 1840 { // Singular value is made nonnegative
jsoh91 2:e843c1b0b25c 1841 W[k] = -z;
jsoh91 2:e843c1b0b25c 1842 for (j = 0; j < n; j++)
jsoh91 2:e843c1b0b25c 1843 V[j][k] = -V[j][k];
jsoh91 2:e843c1b0b25c 1844 }
jsoh91 2:e843c1b0b25c 1845 break;
jsoh91 2:e843c1b0b25c 1846 }
jsoh91 2:e843c1b0b25c 1847 if (its == max_its)
jsoh91 2:e843c1b0b25c 1848 throw std::logic_error("Error svd: no convergence in the maximum number of iterations");
jsoh91 2:e843c1b0b25c 1849 x = W[l];
jsoh91 2:e843c1b0b25c 1850 nm = k - 1;
jsoh91 2:e843c1b0b25c 1851 y = W[nm];
jsoh91 2:e843c1b0b25c 1852 g = rv1[nm];
jsoh91 2:e843c1b0b25c 1853 h = rv1[k];
jsoh91 2:e843c1b0b25c 1854 f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
jsoh91 2:e843c1b0b25c 1855 g = dist(f, (T)1.0);
jsoh91 2:e843c1b0b25c 1856 f = ((x - z) * (x + z) + h * ((y / (f + sign(f)*fabs(g))) - h)) / x;
jsoh91 2:e843c1b0b25c 1857 c = s = (T)1.0; // Next QR transformation
jsoh91 2:e843c1b0b25c 1858 for (j = l; j <= nm; j++)
jsoh91 2:e843c1b0b25c 1859 {
jsoh91 2:e843c1b0b25c 1860 i = j + 1;
jsoh91 2:e843c1b0b25c 1861 g = rv1[i];
jsoh91 2:e843c1b0b25c 1862 y = W[i];
jsoh91 2:e843c1b0b25c 1863 h = s * g;
jsoh91 2:e843c1b0b25c 1864 g *= c;
jsoh91 2:e843c1b0b25c 1865 z = dist(f, h);
jsoh91 2:e843c1b0b25c 1866 rv1[j] = z;
jsoh91 2:e843c1b0b25c 1867 c = f / z;
jsoh91 2:e843c1b0b25c 1868 s = h / z;
jsoh91 2:e843c1b0b25c 1869 f = x * c + g * s;
jsoh91 2:e843c1b0b25c 1870 g = g * c - x * s;
jsoh91 2:e843c1b0b25c 1871 h = y * s;
jsoh91 2:e843c1b0b25c 1872 y *= c;
jsoh91 2:e843c1b0b25c 1873 for (jj = 0; jj < n; jj++)
jsoh91 2:e843c1b0b25c 1874 {
jsoh91 2:e843c1b0b25c 1875 x = V[jj][j];
jsoh91 2:e843c1b0b25c 1876 z = V[jj][i];
jsoh91 2:e843c1b0b25c 1877 V[jj][j] = x * c + z * s;
jsoh91 2:e843c1b0b25c 1878 V[jj][i] = z * c - x * s;
jsoh91 2:e843c1b0b25c 1879 }
jsoh91 2:e843c1b0b25c 1880 z = dist(f, h);
jsoh91 2:e843c1b0b25c 1881 W[j] = z;
jsoh91 2:e843c1b0b25c 1882 if (z != 0) // Rotation can be arbitrary if z = 0
jsoh91 2:e843c1b0b25c 1883 {
jsoh91 2:e843c1b0b25c 1884 z = (T)1.0 / z;
jsoh91 2:e843c1b0b25c 1885 c = f * z;
jsoh91 2:e843c1b0b25c 1886 s = h * z;
jsoh91 2:e843c1b0b25c 1887 }
jsoh91 2:e843c1b0b25c 1888 f = c * g + s * y;
jsoh91 2:e843c1b0b25c 1889 x = c * y - s * g;
jsoh91 2:e843c1b0b25c 1890 for (jj = 0; jj < m; jj++)
jsoh91 2:e843c1b0b25c 1891 {
jsoh91 2:e843c1b0b25c 1892 y = U[jj][j];
jsoh91 2:e843c1b0b25c 1893 z = U[jj][i];
jsoh91 2:e843c1b0b25c 1894 U[jj][j] = y * c + z * s;
jsoh91 2:e843c1b0b25c 1895 U[jj][i] = z * c - y * s;
jsoh91 2:e843c1b0b25c 1896 }
jsoh91 2:e843c1b0b25c 1897 }
jsoh91 2:e843c1b0b25c 1898 rv1[l] = (T)0.0;
jsoh91 2:e843c1b0b25c 1899 rv1[k] = f;
jsoh91 2:e843c1b0b25c 1900 W[k] = x;
jsoh91 2:e843c1b0b25c 1901 }
jsoh91 2:e843c1b0b25c 1902 }
jsoh91 2:e843c1b0b25c 1903 }
jsoh91 2:e843c1b0b25c 1904
jsoh91 2:e843c1b0b25c 1905 template <typename T>
jsoh91 2:e843c1b0b25c 1906 Matrix<T> pinv(const Matrix<T>& A)
jsoh91 2:e843c1b0b25c 1907 {
jsoh91 2:e843c1b0b25c 1908 Matrix<T> U, V, x, tmp(A.ncols(), A.nrows());
jsoh91 2:e843c1b0b25c 1909 Vector<T> W;
jsoh91 2:e843c1b0b25c 1910 CanonicalBaseVector<T> e(0, A.nrows());
jsoh91 2:e843c1b0b25c 1911 svd(A, U, W, V);
jsoh91 2:e843c1b0b25c 1912 for (unsigned int i = 0; i < A.nrows(); i++)
jsoh91 2:e843c1b0b25c 1913 {
jsoh91 2:e843c1b0b25c 1914 e.reset(i);
jsoh91 2:e843c1b0b25c 1915 tmp.setColumn(i, dot_prod(dot_prod(dot_prod(V, Matrix<double>(DIAG, 1.0 / W, 0.0, W.size(), W.size())), t(U)), e));
jsoh91 2:e843c1b0b25c 1916 }
jsoh91 2:e843c1b0b25c 1917
jsoh91 2:e843c1b0b25c 1918 return tmp;
jsoh91 2:e843c1b0b25c 1919 }
jsoh91 2:e843c1b0b25c 1920
jsoh91 2:e843c1b0b25c 1921 template <typename T>
jsoh91 2:e843c1b0b25c 1922 int lu(const Matrix<T>& A, Matrix<T>& LU, Vector<unsigned int>& index)
jsoh91 2:e843c1b0b25c 1923 {
jsoh91 2:e843c1b0b25c 1924 if (A.ncols() != A.nrows())
jsoh91 2:e843c1b0b25c 1925 throw std::logic_error("Error in LU decomposition: matrix must be squared");
jsoh91 2:e843c1b0b25c 1926 int i, p, j, k, n = A.ncols(), ex;
jsoh91 2:e843c1b0b25c 1927 T val, tmp;
jsoh91 2:e843c1b0b25c 1928 Vector<T> d(n);
jsoh91 2:e843c1b0b25c 1929 LU = A;
jsoh91 2:e843c1b0b25c 1930 index.resize(n);
jsoh91 2:e843c1b0b25c 1931
jsoh91 2:e843c1b0b25c 1932 ex = 1;
jsoh91 2:e843c1b0b25c 1933 for (i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 1934 {
jsoh91 2:e843c1b0b25c 1935 index[i] = i;
jsoh91 2:e843c1b0b25c 1936 val = (T)0.0;
jsoh91 2:e843c1b0b25c 1937 for (j = 0; j < n; j++)
jsoh91 2:e843c1b0b25c 1938 val = std::max(val, (T)fabs(LU[i][j]));
jsoh91 2:e843c1b0b25c 1939 if (val == (T)0.0)
jsoh91 2:e843c1b0b25c 1940 std::logic_error("Error in LU decomposition: matrix was singular");
jsoh91 2:e843c1b0b25c 1941 d[i] = val;
jsoh91 2:e843c1b0b25c 1942 }
jsoh91 2:e843c1b0b25c 1943
jsoh91 2:e843c1b0b25c 1944 for (k = 0; k < n - 1; k++)
jsoh91 2:e843c1b0b25c 1945 {
jsoh91 2:e843c1b0b25c 1946 p = k;
jsoh91 2:e843c1b0b25c 1947 val = fabs(LU[k][k]) / d[k];
jsoh91 2:e843c1b0b25c 1948 for (i = k + 1; i < n; i++)
jsoh91 2:e843c1b0b25c 1949 {
jsoh91 2:e843c1b0b25c 1950 tmp = fabs(LU[i][k]) / d[i];
jsoh91 2:e843c1b0b25c 1951 if (tmp > val)
jsoh91 2:e843c1b0b25c 1952 {
jsoh91 2:e843c1b0b25c 1953 val = tmp;
jsoh91 2:e843c1b0b25c 1954 p = i;
jsoh91 2:e843c1b0b25c 1955 }
jsoh91 2:e843c1b0b25c 1956 }
jsoh91 2:e843c1b0b25c 1957 if (val == (T)0.0)
jsoh91 2:e843c1b0b25c 1958 std::logic_error("Error in LU decomposition: matrix was singular");
jsoh91 2:e843c1b0b25c 1959 if (p > k)
jsoh91 2:e843c1b0b25c 1960 {
jsoh91 2:e843c1b0b25c 1961 ex = -ex;
jsoh91 2:e843c1b0b25c 1962 std::swap(index[k], index[p]);
jsoh91 2:e843c1b0b25c 1963 std::swap(d[k], d[p]);
jsoh91 2:e843c1b0b25c 1964 for (j = 0; j < n; j++)
jsoh91 2:e843c1b0b25c 1965 std::swap(LU[k][j], LU[p][j]);
jsoh91 2:e843c1b0b25c 1966 }
jsoh91 2:e843c1b0b25c 1967
jsoh91 2:e843c1b0b25c 1968 for (i = k + 1; i < n; i++)
jsoh91 2:e843c1b0b25c 1969 {
jsoh91 2:e843c1b0b25c 1970 LU[i][k] /= LU[k][k];
jsoh91 2:e843c1b0b25c 1971 for (j = k + 1; j < n; j++)
jsoh91 2:e843c1b0b25c 1972 LU[i][j] -= LU[i][k] * LU[k][j];
jsoh91 2:e843c1b0b25c 1973 }
jsoh91 2:e843c1b0b25c 1974 }
jsoh91 2:e843c1b0b25c 1975 if (LU[n - 1][n - 1] == (T)0.0)
jsoh91 2:e843c1b0b25c 1976 std::logic_error("Error in LU decomposition: matrix was singular");
jsoh91 2:e843c1b0b25c 1977
jsoh91 2:e843c1b0b25c 1978 return ex;
jsoh91 2:e843c1b0b25c 1979 }
jsoh91 2:e843c1b0b25c 1980
jsoh91 2:e843c1b0b25c 1981 template <typename T>
jsoh91 2:e843c1b0b25c 1982 Vector<T> lu_solve(const Matrix<T>& LU, const Vector<T>& b, Vector<unsigned int>& index)
jsoh91 2:e843c1b0b25c 1983 {
jsoh91 2:e843c1b0b25c 1984 if (LU.ncols() != LU.nrows())
jsoh91 2:e843c1b0b25c 1985 throw std::logic_error("Error in LU solve: LU matrix should be squared");
jsoh91 2:e843c1b0b25c 1986 unsigned int n = LU.ncols();
jsoh91 2:e843c1b0b25c 1987 if (b.size() != n)
jsoh91 2:e843c1b0b25c 1988 throw std::logic_error("Error in LU solve: b vector must be of the same dimensions of LU matrix");
jsoh91 2:e843c1b0b25c 1989 Vector<T> x((T)0.0, n);
jsoh91 2:e843c1b0b25c 1990 int i, j, p;
jsoh91 2:e843c1b0b25c 1991 T sum;
jsoh91 2:e843c1b0b25c 1992
jsoh91 2:e843c1b0b25c 1993 p = index[0];
jsoh91 2:e843c1b0b25c 1994 x[0] = b[p];
jsoh91 2:e843c1b0b25c 1995
jsoh91 2:e843c1b0b25c 1996 for (i = 1; i < n; i++)
jsoh91 2:e843c1b0b25c 1997 {
jsoh91 2:e843c1b0b25c 1998 sum = (T)0.0;
jsoh91 2:e843c1b0b25c 1999 for (j = 0; j < i; j++)
jsoh91 2:e843c1b0b25c 2000 sum += LU[i][j] * x[j];
jsoh91 2:e843c1b0b25c 2001 p = index[i];
jsoh91 2:e843c1b0b25c 2002 x[i] = b[p] - sum;
jsoh91 2:e843c1b0b25c 2003 }
jsoh91 2:e843c1b0b25c 2004 x[n - 1] /= LU[n - 1][n - 1];
jsoh91 2:e843c1b0b25c 2005 for (i = n - 2; i >= 0; i--)
jsoh91 2:e843c1b0b25c 2006 {
jsoh91 2:e843c1b0b25c 2007 sum = (T)0.0;
jsoh91 2:e843c1b0b25c 2008 for (j = i + 1; j < n; j++)
jsoh91 2:e843c1b0b25c 2009 sum += LU[i][j] * x[j];
jsoh91 2:e843c1b0b25c 2010 x[i] = (x[i] - sum) / LU[i][i];
jsoh91 2:e843c1b0b25c 2011 }
jsoh91 2:e843c1b0b25c 2012 return x;
jsoh91 2:e843c1b0b25c 2013 }
jsoh91 2:e843c1b0b25c 2014
jsoh91 2:e843c1b0b25c 2015 template <typename T>
jsoh91 2:e843c1b0b25c 2016 void lu_solve(const Matrix<T>& LU, Vector<T>& x, const Vector<T>& b, Vector<unsigned int>& index)
jsoh91 2:e843c1b0b25c 2017 {
jsoh91 2:e843c1b0b25c 2018 x = lu_solve(LU, b, index);
jsoh91 2:e843c1b0b25c 2019 }
jsoh91 2:e843c1b0b25c 2020
jsoh91 2:e843c1b0b25c 2021 template <typename T>
jsoh91 2:e843c1b0b25c 2022 Matrix<T> lu_inverse(const Matrix<T>& A)
jsoh91 2:e843c1b0b25c 2023 {
jsoh91 2:e843c1b0b25c 2024 if (A.ncols() != A.nrows())
jsoh91 2:e843c1b0b25c 2025 throw std::logic_error("Error in LU invert: matrix must be squared");
jsoh91 2:e843c1b0b25c 2026 unsigned int n = A.ncols();
jsoh91 2:e843c1b0b25c 2027 Matrix<T> A1(n, n), LU;
jsoh91 2:e843c1b0b25c 2028 Vector<unsigned int> index;
jsoh91 2:e843c1b0b25c 2029
jsoh91 2:e843c1b0b25c 2030 lu(A, LU, index);
jsoh91 2:e843c1b0b25c 2031 CanonicalBaseVector<T> e(0, n);
jsoh91 2:e843c1b0b25c 2032 for (unsigned i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 2033 {
jsoh91 2:e843c1b0b25c 2034 e.reset(i);
jsoh91 2:e843c1b0b25c 2035 A1.setColumn(i, lu_solve(LU, e, index));
jsoh91 2:e843c1b0b25c 2036 }
jsoh91 2:e843c1b0b25c 2037
jsoh91 2:e843c1b0b25c 2038 return A1;
jsoh91 2:e843c1b0b25c 2039 }
jsoh91 2:e843c1b0b25c 2040
jsoh91 2:e843c1b0b25c 2041 template <typename T>
jsoh91 2:e843c1b0b25c 2042 T lu_det(const Matrix<T>& A)
jsoh91 2:e843c1b0b25c 2043 {
jsoh91 2:e843c1b0b25c 2044 if (A.ncols() != A.nrows())
jsoh91 2:e843c1b0b25c 2045 throw std::logic_error("Error in LU determinant: matrix must be squared");
jsoh91 2:e843c1b0b25c 2046 unsigned int d;
jsoh91 2:e843c1b0b25c 2047 Matrix<T> LU;
jsoh91 2:e843c1b0b25c 2048 Vector<unsigned int> index;
jsoh91 2:e843c1b0b25c 2049
jsoh91 2:e843c1b0b25c 2050 d = lu(A, LU, index);
jsoh91 2:e843c1b0b25c 2051
jsoh91 2:e843c1b0b25c 2052 return d * prod(LU.extractDiag());
jsoh91 2:e843c1b0b25c 2053 }
jsoh91 2:e843c1b0b25c 2054
jsoh91 2:e843c1b0b25c 2055 template <typename T>
jsoh91 2:e843c1b0b25c 2056 void cholesky(const Matrix<T> A, Matrix<T>& LL)
jsoh91 2:e843c1b0b25c 2057 {
jsoh91 2:e843c1b0b25c 2058 if (A.ncols() != A.nrows())
jsoh91 2:e843c1b0b25c 2059 throw std::logic_error("Error in Cholesky decomposition: matrix must be squared");
jsoh91 2:e843c1b0b25c 2060 register int i, j, k, n = A.ncols();
jsoh91 2:e843c1b0b25c 2061 register double sum;
jsoh91 2:e843c1b0b25c 2062 LL = A;
jsoh91 2:e843c1b0b25c 2063
jsoh91 2:e843c1b0b25c 2064 for (i = 0; i < n; i++)
jsoh91 2:e843c1b0b25c 2065 {
jsoh91 2:e843c1b0b25c 2066 for (j = i; j < n; j++)
jsoh91 2:e843c1b0b25c 2067 {
jsoh91 2:e843c1b0b25c 2068 sum = LL[i][j];
jsoh91 2:e843c1b0b25c 2069 for (k = i - 1; k >= 0; k--)
jsoh91 2:e843c1b0b25c 2070 sum -= LL[i][k] * LL[j][k];
jsoh91 2:e843c1b0b25c 2071 if (i == j)
jsoh91 2:e843c1b0b25c 2072 {
jsoh91 2:e843c1b0b25c 2073 if (sum <= 0.0)
jsoh91 2:e843c1b0b25c 2074 throw std::logic_error("Error in Cholesky decomposition: matrix is not postive definite");
jsoh91 2:e843c1b0b25c 2075 LL[i][i] = sqrt(sum);
jsoh91 2:e843c1b0b25c 2076 }
jsoh91 2:e843c1b0b25c 2077 else
jsoh91 2:e843c1b0b25c 2078 LL[j][i] = sum / LL[i][i];
jsoh91 2:e843c1b0b25c 2079 }
jsoh91 2:e843c1b0b25c 2080 for (k = i + 1; k < n; k++)
jsoh91 2:e843c1b0b25c 2081 LL[i][k] = LL[k][i];
jsoh91 2:e843c1b0b25c 2082 }
jsoh91 2:e843c1b0b25c 2083 }
jsoh91 2:e843c1b0b25c 2084
jsoh91 2:e843c1b0b25c 2085 template <typename T>
jsoh91 2:e843c1b0b25c 2086 Matrix<T> cholesky(const Matrix<T> A)
jsoh91 2:e843c1b0b25c 2087 {
jsoh91 2:e843c1b0b25c 2088 Matrix<T> LL;
jsoh91 2:e843c1b0b25c 2089 cholesky(A, LL);
jsoh91 2:e843c1b0b25c 2090
jsoh91 2:e843c1b0b25c 2091 return LL;
jsoh91 2:e843c1b0b25c 2092 }
jsoh91 2:e843c1b0b25c 2093
jsoh91 2:e843c1b0b25c 2094 template <typename T>
jsoh91 2:e843c1b0b25c 2095 Vector<T> cholesky_solve(const Matrix<T>& LL, const Vector<T>& b)
jsoh91 2:e843c1b0b25c 2096 {
jsoh91 2:e843c1b0b25c 2097 if (LL.ncols() != LL.nrows())
jsoh91 2:e843c1b0b25c 2098 throw std::logic_error("Error in Cholesky solve: matrix must be squared");
jsoh91 2:e843c1b0b25c 2099 unsigned int n = LL.ncols();
jsoh91 2:e843c1b0b25c 2100 if (b.size() != n)
jsoh91 2:e843c1b0b25c 2101 throw std::logic_error("Error in Cholesky decomposition: b vector must be of the same dimensions of LU matrix");
jsoh91 2:e843c1b0b25c 2102 Vector<T> x, y;
jsoh91 2:e843c1b0b25c 2103
jsoh91 2:e843c1b0b25c 2104 /* Solve L * y = b */
jsoh91 2:e843c1b0b25c 2105 forward_elimination(LL, y, b);
jsoh91 2:e843c1b0b25c 2106 /* Solve L^T * x = y */
jsoh91 2:e843c1b0b25c 2107 backward_elimination(LL, x, y);
jsoh91 2:e843c1b0b25c 2108
jsoh91 2:e843c1b0b25c 2109 return x;
jsoh91 2:e843c1b0b25c 2110 }
jsoh91 2:e843c1b0b25c 2111
jsoh91 2:e843c1b0b25c 2112 template <typename T>
jsoh91 2:e843c1b0b25c 2113 void cholesky_solve(const Matrix<T>& LL, Vector<T>& x, const Vector<T>& b)
jsoh91 2:e843c1b0b25c 2114 {
jsoh91 2:e843c1b0b25c 2115 x = cholesky_solve(LL, b);
jsoh91 2:e843c1b0b25c 2116 }
jsoh91 2:e843c1b0b25c 2117
jsoh91 2:e843c1b0b25c 2118 template <typename T>
jsoh91 2:e843c1b0b25c 2119 void forward_elimination(const Matrix<T>& L, Vector<T>& y, const Vector<T> b)
jsoh91 2:e843c1b0b25c 2120 {
jsoh91 2:e843c1b0b25c 2121 if (L.ncols() != L.nrows())
jsoh91 2:e843c1b0b25c 2122 throw std::logic_error("Error in Forward elimination: matrix must be squared (lower triangular)");
jsoh91 2:e843c1b0b25c 2123 if (b.size() != L.nrows())
jsoh91 2:e843c1b0b25c 2124 throw std::logic_error("Error in Forward elimination: b vector must be of the same dimensions of L matrix");
jsoh91 2:e843c1b0b25c 2125 register int i, j, n = b.size();
jsoh91 2:e843c1b0b25c 2126 y.resize(n);
jsoh91 2:e843c1b0b25c 2127
jsoh91 2:e843c1b0b25c 2128 y[0] = b[0] / L[0][0];
jsoh91 2:e843c1b0b25c 2129 for (i = 1; i < n; i++)
jsoh91 2:e843c1b0b25c 2130 {
jsoh91 2:e843c1b0b25c 2131 y[i] = b[i];
jsoh91 2:e843c1b0b25c 2132 for (j = 0; j < i; j++)
jsoh91 2:e843c1b0b25c 2133 y[i] -= L[i][j] * y[j];
jsoh91 2:e843c1b0b25c 2134 y[i] = y[i] / L[i][i];
jsoh91 2:e843c1b0b25c 2135 }
jsoh91 2:e843c1b0b25c 2136 }
jsoh91 2:e843c1b0b25c 2137
jsoh91 2:e843c1b0b25c 2138 template <typename T>
jsoh91 2:e843c1b0b25c 2139 Vector<T> forward_elimination(const Matrix<T>& L, const Vector<T> b)
jsoh91 2:e843c1b0b25c 2140 {
jsoh91 2:e843c1b0b25c 2141 Vector<T> y;
jsoh91 2:e843c1b0b25c 2142 forward_elimination(L, y, b);
jsoh91 2:e843c1b0b25c 2143
jsoh91 2:e843c1b0b25c 2144 return y;
jsoh91 2:e843c1b0b25c 2145 }
jsoh91 2:e843c1b0b25c 2146
jsoh91 2:e843c1b0b25c 2147 template <typename T>
jsoh91 2:e843c1b0b25c 2148 void backward_elimination(const Matrix<T>& U, Vector<T>& x, const Vector<T>& y)
jsoh91 2:e843c1b0b25c 2149 {
jsoh91 2:e843c1b0b25c 2150 if (U.ncols() != U.nrows())
jsoh91 2:e843c1b0b25c 2151 throw std::logic_error("Error in Backward elimination: matrix must be squared (upper triangular)");
jsoh91 2:e843c1b0b25c 2152 if (y.size() != U.nrows())
jsoh91 2:e843c1b0b25c 2153 throw std::logic_error("Error in Backward elimination: b vector must be of the same dimensions of U matrix");
jsoh91 2:e843c1b0b25c 2154 register int i, j, n = y.size();
jsoh91 2:e843c1b0b25c 2155 x.resize(n);
jsoh91 2:e843c1b0b25c 2156
jsoh91 2:e843c1b0b25c 2157 x[n - 1] = y[n - 1] / U[n - 1][n - 1];
jsoh91 2:e843c1b0b25c 2158 for (i = n - 2; i >= 0; i--)
jsoh91 2:e843c1b0b25c 2159 {
jsoh91 2:e843c1b0b25c 2160 x[i] = y[i];
jsoh91 2:e843c1b0b25c 2161 for (j = i + 1; j < n; j++)
jsoh91 2:e843c1b0b25c 2162 x[i] -= U[i][j] * x[j];
jsoh91 2:e843c1b0b25c 2163 x[i] = x[i] / U[i][i];
jsoh91 2:e843c1b0b25c 2164 }
jsoh91 2:e843c1b0b25c 2165 }
jsoh91 2:e843c1b0b25c 2166
jsoh91 2:e843c1b0b25c 2167 template <typename T>
jsoh91 2:e843c1b0b25c 2168 Vector<T> backward_elimination(const Matrix<T>& U, const Vector<T> y)
jsoh91 2:e843c1b0b25c 2169 {
jsoh91 2:e843c1b0b25c 2170 Vector<T> x;
jsoh91 2:e843c1b0b25c 2171 forward_elimination(U, x, y);
jsoh91 2:e843c1b0b25c 2172
jsoh91 2:e843c1b0b25c 2173 return x;
jsoh91 2:e843c1b0b25c 2174 }
jsoh91 2:e843c1b0b25c 2175
jsoh91 2:e843c1b0b25c 2176 /* Setting default linear systems machinery */
jsoh91 2:e843c1b0b25c 2177
jsoh91 2:e843c1b0b25c 2178 #define det lu_det
jsoh91 2:e843c1b0b25c 2179 //#define inverse lu_inverse
jsoh91 2:e843c1b0b25c 2180 #define solve lu_solve
jsoh91 2:e843c1b0b25c 2181
jsoh91 2:e843c1b0b25c 2182 /* Random */
jsoh91 2:e843c1b0b25c 2183
jsoh91 2:e843c1b0b25c 2184 template <typename T>
jsoh91 2:e843c1b0b25c 2185 void random(Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2186 {
jsoh91 2:e843c1b0b25c 2187 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2188 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2189 m[i][j] = (T)(rand() / double(RAND_MAX));
jsoh91 2:e843c1b0b25c 2190 }
jsoh91 2:e843c1b0b25c 2191
jsoh91 2:e843c1b0b25c 2192 /**
jsoh91 2:e843c1b0b25c 2193 Aggregate functions
jsoh91 2:e843c1b0b25c 2194 */
jsoh91 2:e843c1b0b25c 2195
jsoh91 2:e843c1b0b25c 2196 template <typename T>
jsoh91 2:e843c1b0b25c 2197 Vector<T> sum(const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2198 {
jsoh91 2:e843c1b0b25c 2199 Vector<T> tmp((T)0, m.ncols());
jsoh91 2:e843c1b0b25c 2200 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2201 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2202 tmp[j] += m[i][j];
jsoh91 2:e843c1b0b25c 2203 return tmp;
jsoh91 2:e843c1b0b25c 2204 }
jsoh91 2:e843c1b0b25c 2205
jsoh91 2:e843c1b0b25c 2206 template <typename T>
jsoh91 2:e843c1b0b25c 2207 Vector<T> r_sum(const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2208 {
jsoh91 2:e843c1b0b25c 2209 Vector<T> tmp((T)0, m.nrows());
jsoh91 2:e843c1b0b25c 2210 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2211 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2212 tmp[i] += m[i][j];
jsoh91 2:e843c1b0b25c 2213 return tmp;
jsoh91 2:e843c1b0b25c 2214 }
jsoh91 2:e843c1b0b25c 2215
jsoh91 2:e843c1b0b25c 2216 template <typename T>
jsoh91 2:e843c1b0b25c 2217 T all_sum(const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2218 {
jsoh91 2:e843c1b0b25c 2219 T tmp = (T)0;
jsoh91 2:e843c1b0b25c 2220 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2221 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2222 tmp += m[i][j];
jsoh91 2:e843c1b0b25c 2223 return tmp;
jsoh91 2:e843c1b0b25c 2224 }
jsoh91 2:e843c1b0b25c 2225
jsoh91 2:e843c1b0b25c 2226 template <typename T>
jsoh91 2:e843c1b0b25c 2227 Vector<T> prod(const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2228 {
jsoh91 2:e843c1b0b25c 2229 Vector<T> tmp((T)1, m.ncols());
jsoh91 2:e843c1b0b25c 2230 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2231 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2232 tmp[j] *= m[i][j];
jsoh91 2:e843c1b0b25c 2233 return tmp;
jsoh91 2:e843c1b0b25c 2234 }
jsoh91 2:e843c1b0b25c 2235
jsoh91 2:e843c1b0b25c 2236 template <typename T>
jsoh91 2:e843c1b0b25c 2237 Vector<T> r_prod(const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2238 {
jsoh91 2:e843c1b0b25c 2239 Vector<T> tmp((T)1, m.nrows());
jsoh91 2:e843c1b0b25c 2240 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2241 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2242 tmp[i] *= m[i][j];
jsoh91 2:e843c1b0b25c 2243 return tmp;
jsoh91 2:e843c1b0b25c 2244 }
jsoh91 2:e843c1b0b25c 2245
jsoh91 2:e843c1b0b25c 2246 template <typename T>
jsoh91 2:e843c1b0b25c 2247 T all_prod(const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2248 {
jsoh91 2:e843c1b0b25c 2249 T tmp = (T)1;
jsoh91 2:e843c1b0b25c 2250 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2251 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2252 tmp *= m[i][j];
jsoh91 2:e843c1b0b25c 2253 return tmp;
jsoh91 2:e843c1b0b25c 2254 }
jsoh91 2:e843c1b0b25c 2255
jsoh91 2:e843c1b0b25c 2256 template <typename T>
jsoh91 2:e843c1b0b25c 2257 Vector<T> mean(const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2258 {
jsoh91 2:e843c1b0b25c 2259 Vector<T> res((T)0, m.ncols());
jsoh91 2:e843c1b0b25c 2260 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2261 {
jsoh91 2:e843c1b0b25c 2262 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2263 res[j] += m[i][j];
jsoh91 2:e843c1b0b25c 2264 res[j] /= m.nrows();
jsoh91 2:e843c1b0b25c 2265 }
jsoh91 2:e843c1b0b25c 2266
jsoh91 2:e843c1b0b25c 2267 return res;
jsoh91 2:e843c1b0b25c 2268 }
jsoh91 2:e843c1b0b25c 2269
jsoh91 2:e843c1b0b25c 2270 template <typename T>
jsoh91 2:e843c1b0b25c 2271 Vector<T> r_mean(const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2272 {
jsoh91 2:e843c1b0b25c 2273 Vector<T> res((T)0, m.rows());
jsoh91 2:e843c1b0b25c 2274 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2275 {
jsoh91 2:e843c1b0b25c 2276 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2277 res[i] += m[i][j];
jsoh91 2:e843c1b0b25c 2278 res[i] /= m.nrows();
jsoh91 2:e843c1b0b25c 2279 }
jsoh91 2:e843c1b0b25c 2280
jsoh91 2:e843c1b0b25c 2281 return res;
jsoh91 2:e843c1b0b25c 2282 }
jsoh91 2:e843c1b0b25c 2283
jsoh91 2:e843c1b0b25c 2284 template <typename T>
jsoh91 2:e843c1b0b25c 2285 T all_mean(const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2286 {
jsoh91 2:e843c1b0b25c 2287 T tmp = (T)0;
jsoh91 2:e843c1b0b25c 2288 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2289 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2290 tmp += m[i][j];
jsoh91 2:e843c1b0b25c 2291 return tmp / (m.nrows() * m.ncols());
jsoh91 2:e843c1b0b25c 2292 }
jsoh91 2:e843c1b0b25c 2293
jsoh91 2:e843c1b0b25c 2294 template <typename T>
jsoh91 2:e843c1b0b25c 2295 Vector<T> var(const Matrix<T>& m, bool sample_correction = false)
jsoh91 2:e843c1b0b25c 2296 {
jsoh91 2:e843c1b0b25c 2297 Vector<T> res((T)0, m.ncols());
jsoh91 2:e843c1b0b25c 2298 unsigned int n = m.nrows();
jsoh91 2:e843c1b0b25c 2299 double sum, ssum;
jsoh91 2:e843c1b0b25c 2300 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2301 {
jsoh91 2:e843c1b0b25c 2302 sum = (T)0.0; ssum = (T)0.0;
jsoh91 2:e843c1b0b25c 2303 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2304 {
jsoh91 2:e843c1b0b25c 2305 sum += m[i][j];
jsoh91 2:e843c1b0b25c 2306 ssum += (m[i][j] * m[i][j]);
jsoh91 2:e843c1b0b25c 2307 }
jsoh91 2:e843c1b0b25c 2308 if (!sample_correction)
jsoh91 2:e843c1b0b25c 2309 res[j] = (ssum / n) - (sum / n) * (sum / n);
jsoh91 2:e843c1b0b25c 2310 else
jsoh91 2:e843c1b0b25c 2311 res[j] = n * ((ssum / n) - (sum / n) * (sum / n)) / (n - 1);
jsoh91 2:e843c1b0b25c 2312 }
jsoh91 2:e843c1b0b25c 2313
jsoh91 2:e843c1b0b25c 2314 return res;
jsoh91 2:e843c1b0b25c 2315 }
jsoh91 2:e843c1b0b25c 2316
jsoh91 2:e843c1b0b25c 2317 template <typename T>
jsoh91 2:e843c1b0b25c 2318 Vector<T> stdev(const Matrix<T>& m, bool sample_correction = false)
jsoh91 2:e843c1b0b25c 2319 {
jsoh91 2:e843c1b0b25c 2320 return vec_sqrt(var(m, sample_correction));
jsoh91 2:e843c1b0b25c 2321 }
jsoh91 2:e843c1b0b25c 2322
jsoh91 2:e843c1b0b25c 2323 template <typename T>
jsoh91 2:e843c1b0b25c 2324 Vector<T> r_var(const Matrix<T>& m, bool sample_correction = false)
jsoh91 2:e843c1b0b25c 2325 {
jsoh91 2:e843c1b0b25c 2326 Vector<T> res((T)0, m.nrows());
jsoh91 2:e843c1b0b25c 2327 double sum, ssum;
jsoh91 2:e843c1b0b25c 2328 unsigned int n = m.ncols();
jsoh91 2:e843c1b0b25c 2329 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2330 {
jsoh91 2:e843c1b0b25c 2331 sum = 0.0; ssum = 0.0;
jsoh91 2:e843c1b0b25c 2332 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2333 {
jsoh91 2:e843c1b0b25c 2334 sum += m[i][j];
jsoh91 2:e843c1b0b25c 2335 ssum += (m[i][j] * m[i][j]);
jsoh91 2:e843c1b0b25c 2336 }
jsoh91 2:e843c1b0b25c 2337 if (!sample_correction)
jsoh91 2:e843c1b0b25c 2338 res[i] = (ssum / n) - (sum / n) * (sum / n);
jsoh91 2:e843c1b0b25c 2339 else
jsoh91 2:e843c1b0b25c 2340 res[i] = n * ((ssum / n) - (sum / n) * (sum / n)) / (n - 1);
jsoh91 2:e843c1b0b25c 2341 }
jsoh91 2:e843c1b0b25c 2342
jsoh91 2:e843c1b0b25c 2343 return res;
jsoh91 2:e843c1b0b25c 2344 }
jsoh91 2:e843c1b0b25c 2345
jsoh91 2:e843c1b0b25c 2346 template <typename T>
jsoh91 2:e843c1b0b25c 2347 Vector<T> r_stdev(const Matrix<T>& m, bool sample_correction = false)
jsoh91 2:e843c1b0b25c 2348 {
jsoh91 2:e843c1b0b25c 2349 return vec_sqrt(r_var(m, sample_correction));
jsoh91 2:e843c1b0b25c 2350 }
jsoh91 2:e843c1b0b25c 2351
jsoh91 2:e843c1b0b25c 2352 template <typename T>
jsoh91 2:e843c1b0b25c 2353 Vector<T> max(const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2354 {
jsoh91 2:e843c1b0b25c 2355 Vector<T> res(m.ncols());
jsoh91 2:e843c1b0b25c 2356 double value;
jsoh91 2:e843c1b0b25c 2357 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2358 {
jsoh91 2:e843c1b0b25c 2359 value = m[0][j];
jsoh91 2:e843c1b0b25c 2360 for (unsigned int i = 1; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2361 value = std::max(m[i][j], value);
jsoh91 2:e843c1b0b25c 2362 res[j] = value;
jsoh91 2:e843c1b0b25c 2363 }
jsoh91 2:e843c1b0b25c 2364
jsoh91 2:e843c1b0b25c 2365 return res;
jsoh91 2:e843c1b0b25c 2366 }
jsoh91 2:e843c1b0b25c 2367
jsoh91 2:e843c1b0b25c 2368 template <typename T>
jsoh91 2:e843c1b0b25c 2369 Vector<T> r_max(const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2370 {
jsoh91 2:e843c1b0b25c 2371 Vector<T> res(m.nrows());
jsoh91 2:e843c1b0b25c 2372 double value;
jsoh91 2:e843c1b0b25c 2373 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2374 {
jsoh91 2:e843c1b0b25c 2375 value = m[i][0];
jsoh91 2:e843c1b0b25c 2376 for (unsigned int j = 1; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2377 value = std::max(m[i][j], value);
jsoh91 2:e843c1b0b25c 2378 res[i] = value;
jsoh91 2:e843c1b0b25c 2379 }
jsoh91 2:e843c1b0b25c 2380
jsoh91 2:e843c1b0b25c 2381 return res;
jsoh91 2:e843c1b0b25c 2382 }
jsoh91 2:e843c1b0b25c 2383
jsoh91 2:e843c1b0b25c 2384 template <typename T>
jsoh91 2:e843c1b0b25c 2385 Vector<T> min(const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2386 {
jsoh91 2:e843c1b0b25c 2387 Vector<T> res(m.ncols());
jsoh91 2:e843c1b0b25c 2388 double value;
jsoh91 2:e843c1b0b25c 2389 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2390 {
jsoh91 2:e843c1b0b25c 2391 value = m[0][j];
jsoh91 2:e843c1b0b25c 2392 for (unsigned int i = 1; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2393 value = std::min(m[i][j], value);
jsoh91 2:e843c1b0b25c 2394 res[j] = value;
jsoh91 2:e843c1b0b25c 2395 }
jsoh91 2:e843c1b0b25c 2396
jsoh91 2:e843c1b0b25c 2397 return res;
jsoh91 2:e843c1b0b25c 2398 }
jsoh91 2:e843c1b0b25c 2399
jsoh91 2:e843c1b0b25c 2400 template <typename T>
jsoh91 2:e843c1b0b25c 2401 Vector<T> r_min(const Matrix<T>& m)
jsoh91 2:e843c1b0b25c 2402 {
jsoh91 2:e843c1b0b25c 2403 Vector<T> res(m.nrows());
jsoh91 2:e843c1b0b25c 2404 double value;
jsoh91 2:e843c1b0b25c 2405 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2406 {
jsoh91 2:e843c1b0b25c 2407 value = m[i][0];
jsoh91 2:e843c1b0b25c 2408 for (unsigned int j = 1; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2409 value = std::min(m[i][j], value);
jsoh91 2:e843c1b0b25c 2410 res[i] = value;
jsoh91 2:e843c1b0b25c 2411 }
jsoh91 2:e843c1b0b25c 2412
jsoh91 2:e843c1b0b25c 2413 return res;
jsoh91 2:e843c1b0b25c 2414 }
jsoh91 2:e843c1b0b25c 2415
jsoh91 2:e843c1b0b25c 2416
jsoh91 2:e843c1b0b25c 2417
jsoh91 2:e843c1b0b25c 2418 /**
jsoh91 2:e843c1b0b25c 2419 Single element mathematical functions
jsoh91 2:e843c1b0b25c 2420 */
jsoh91 2:e843c1b0b25c 2421
jsoh91 2:e843c1b0b25c 2422 template <typename T>
jsoh91 2:e843c1b0b25c 2423 Matrix<T> exp(const Matrix<T>&m)
jsoh91 2:e843c1b0b25c 2424 {
jsoh91 2:e843c1b0b25c 2425 Matrix<T> tmp(m.nrows(), m.ncols());
jsoh91 2:e843c1b0b25c 2426
jsoh91 2:e843c1b0b25c 2427 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2428 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2429 tmp[i][j] = exp(m[i][j]);
jsoh91 2:e843c1b0b25c 2430
jsoh91 2:e843c1b0b25c 2431 return tmp;
jsoh91 2:e843c1b0b25c 2432 }
jsoh91 2:e843c1b0b25c 2433
jsoh91 2:e843c1b0b25c 2434 template <typename T>
jsoh91 2:e843c1b0b25c 2435 Matrix<T> mat_sqrt(const Matrix<T>&m)
jsoh91 2:e843c1b0b25c 2436 {
jsoh91 2:e843c1b0b25c 2437 Matrix<T> tmp(m.nrows(), m.ncols());
jsoh91 2:e843c1b0b25c 2438
jsoh91 2:e843c1b0b25c 2439 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2440 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2441 tmp[i][j] = sqrt(m[i][j]);
jsoh91 2:e843c1b0b25c 2442
jsoh91 2:e843c1b0b25c 2443 return tmp;
jsoh91 2:e843c1b0b25c 2444 }
jsoh91 2:e843c1b0b25c 2445
jsoh91 2:e843c1b0b25c 2446 /**
jsoh91 2:e843c1b0b25c 2447 Matrix operators
jsoh91 2:e843c1b0b25c 2448 */
jsoh91 2:e843c1b0b25c 2449
jsoh91 2:e843c1b0b25c 2450 template <typename T>
jsoh91 2:e843c1b0b25c 2451 Matrix<T> kron(const Vector<T>& b, const Vector<T>& a)
jsoh91 2:e843c1b0b25c 2452 {
jsoh91 2:e843c1b0b25c 2453 Matrix<T> tmp(b.size(), a.size());
jsoh91 2:e843c1b0b25c 2454 for (unsigned int i = 0; i < b.size(); i++)
jsoh91 2:e843c1b0b25c 2455 for (unsigned int j = 0; j < a.size(); j++)
jsoh91 2:e843c1b0b25c 2456 tmp[i][j] = a[j] * b[i];
jsoh91 2:e843c1b0b25c 2457
jsoh91 2:e843c1b0b25c 2458 return tmp;
jsoh91 2:e843c1b0b25c 2459 }
jsoh91 2:e843c1b0b25c 2460
jsoh91 2:e843c1b0b25c 2461 template <typename T>
jsoh91 2:e843c1b0b25c 2462 Matrix<T> t(const Matrix<T>& a)
jsoh91 2:e843c1b0b25c 2463 {
jsoh91 2:e843c1b0b25c 2464 Matrix<T> tmp(a.ncols(), a.nrows());
jsoh91 2:e843c1b0b25c 2465 for (unsigned int i = 0; i < a.nrows(); i++)
jsoh91 2:e843c1b0b25c 2466 for (unsigned int j = 0; j < a.ncols(); j++)
jsoh91 2:e843c1b0b25c 2467 tmp[j][i] = a[i][j];
jsoh91 2:e843c1b0b25c 2468
jsoh91 2:e843c1b0b25c 2469 return tmp;
jsoh91 2:e843c1b0b25c 2470 }
jsoh91 2:e843c1b0b25c 2471
jsoh91 2:e843c1b0b25c 2472 template <typename T>
jsoh91 2:e843c1b0b25c 2473 Matrix<T> dot_prod(const Matrix<T>& a, const Matrix<T>& b)
jsoh91 2:e843c1b0b25c 2474 {
jsoh91 2:e843c1b0b25c 2475 if (a.ncols() != b.nrows())
jsoh91 2:e843c1b0b25c 2476 throw std::logic_error("Error matrix dot product: dimensions of the matrices are not compatible");
jsoh91 2:e843c1b0b25c 2477 Matrix<T> tmp(a.nrows(), b.ncols());
jsoh91 2:e843c1b0b25c 2478 for (unsigned int i = 0; i < tmp.nrows(); i++)
jsoh91 2:e843c1b0b25c 2479 for (unsigned int j = 0; j < tmp.ncols(); j++)
jsoh91 2:e843c1b0b25c 2480 {
jsoh91 2:e843c1b0b25c 2481 tmp[i][j] = (T)0;
jsoh91 2:e843c1b0b25c 2482 for (unsigned int k = 0; k < a.ncols(); k++)
jsoh91 2:e843c1b0b25c 2483 tmp[i][j] += a[i][k] * b[k][j];
jsoh91 2:e843c1b0b25c 2484 }
jsoh91 2:e843c1b0b25c 2485
jsoh91 2:e843c1b0b25c 2486 return tmp;
jsoh91 2:e843c1b0b25c 2487 }
jsoh91 2:e843c1b0b25c 2488
jsoh91 2:e843c1b0b25c 2489 template <typename T>
jsoh91 2:e843c1b0b25c 2490 Matrix<T> dot_prod(const Matrix<T>& a, const Vector<T>& b)
jsoh91 2:e843c1b0b25c 2491 {
jsoh91 2:e843c1b0b25c 2492 if (a.ncols() != b.size())
jsoh91 2:e843c1b0b25c 2493 throw std::logic_error("Error matrix dot product: dimensions of the matrix and the vector are not compatible");
jsoh91 2:e843c1b0b25c 2494 Matrix<T> tmp(a.nrows(), 1);
jsoh91 2:e843c1b0b25c 2495 for (unsigned int i = 0; i < tmp.nrows(); i++)
jsoh91 2:e843c1b0b25c 2496 {
jsoh91 2:e843c1b0b25c 2497 tmp[i][0] = (T)0;
jsoh91 2:e843c1b0b25c 2498 for (unsigned int k = 0; k < a.ncols(); k++)
jsoh91 2:e843c1b0b25c 2499 tmp[i][0] += a[i][k] * b[k];
jsoh91 2:e843c1b0b25c 2500 }
jsoh91 2:e843c1b0b25c 2501
jsoh91 2:e843c1b0b25c 2502 return tmp;
jsoh91 2:e843c1b0b25c 2503 }
jsoh91 2:e843c1b0b25c 2504
jsoh91 2:e843c1b0b25c 2505 template <typename T>
jsoh91 2:e843c1b0b25c 2506 Matrix<T> dot_prod(const Vector<T>& a, const Matrix<T>& b)
jsoh91 2:e843c1b0b25c 2507 {
jsoh91 2:e843c1b0b25c 2508 if (a.size() != b.ncols())
jsoh91 2:e843c1b0b25c 2509 throw std::logic_error("Error matrix dot product: dimensions of the vector and matrix are not compatible");
jsoh91 2:e843c1b0b25c 2510 Matrix<T> tmp(1, b.ncols());
jsoh91 2:e843c1b0b25c 2511 for (unsigned int j = 0; j < tmp.ncols(); j++)
jsoh91 2:e843c1b0b25c 2512 {
jsoh91 2:e843c1b0b25c 2513 tmp[0][j] = (T)0;
jsoh91 2:e843c1b0b25c 2514 for (unsigned int k = 0; k < a.size(); k++)
jsoh91 2:e843c1b0b25c 2515 tmp[0][j] += a[k] * b[k][j];
jsoh91 2:e843c1b0b25c 2516 }
jsoh91 2:e843c1b0b25c 2517
jsoh91 2:e843c1b0b25c 2518 return tmp;
jsoh91 2:e843c1b0b25c 2519 }
jsoh91 2:e843c1b0b25c 2520
jsoh91 2:e843c1b0b25c 2521 template <typename T>
jsoh91 2:e843c1b0b25c 2522 inline Matrix<double> rank(const Matrix<T> m)
jsoh91 2:e843c1b0b25c 2523 {
jsoh91 2:e843c1b0b25c 2524 Matrix<double> tmp(m.nrows(), m.ncols());
jsoh91 2:e843c1b0b25c 2525 for (unsigned int j = 0; j < m.ncols(); j++)
jsoh91 2:e843c1b0b25c 2526 tmp.setColumn(j, rank<T>(m.extractColumn(j)));
jsoh91 2:e843c1b0b25c 2527
jsoh91 2:e843c1b0b25c 2528 return tmp;
jsoh91 2:e843c1b0b25c 2529 }
jsoh91 2:e843c1b0b25c 2530
jsoh91 2:e843c1b0b25c 2531 template <typename T>
jsoh91 2:e843c1b0b25c 2532 inline Matrix<double> r_rank(const Matrix<T> m)
jsoh91 2:e843c1b0b25c 2533 {
jsoh91 2:e843c1b0b25c 2534 Matrix<double> tmp(m.nrows(), m.ncols());
jsoh91 2:e843c1b0b25c 2535 for (unsigned int i = 0; i < m.nrows(); i++)
jsoh91 2:e843c1b0b25c 2536 tmp.setRow(i, rank<T>(m.extractRow(i)));
jsoh91 2:e843c1b0b25c 2537
jsoh91 2:e843c1b0b25c 2538 return tmp;
jsoh91 2:e843c1b0b25c 2539 }
jsoh91 2:e843c1b0b25c 2540
jsoh91 2:e843c1b0b25c 2541 } // namespace quadprogpp
jsoh91 2:e843c1b0b25c 2542
jsoh91 2:e843c1b0b25c 2543 #endif // define _ARRAY_HH_