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