yat/utility/Matrix.cc

Code
Comments
Other
Rev Date Author Line
42 26 Feb 04 jari 1 // $Id$
12 19 Jun 03 daniel 2
570 05 Apr 06 jari 3 /*
570 05 Apr 06 jari 4   Copyright (C) 2003 Daniel Dalevi, Peter Johansson
2119 12 Dec 09 peter 5   Copyright (C) 2004 Jari Häkkinen, Peter Johansson
2119 12 Dec 09 peter 6   Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
2119 12 Dec 09 peter 7   Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 8   Copyright (C) 2009, 2010, 2011, 2012, 2017, 2021, 2022, 2023 Peter Johansson
570 05 Apr 06 jari 9
1437 25 Aug 08 peter 10   This file is part of the yat library, http://dev.thep.lu.se/yat
570 05 Apr 06 jari 11
675 10 Oct 06 jari 12   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 13   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 14   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 15   License, or (at your option) any later version.
570 05 Apr 06 jari 16
675 10 Oct 06 jari 17   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 18   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 19   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
570 05 Apr 06 jari 20   General Public License for more details.
570 05 Apr 06 jari 21
570 05 Apr 06 jari 22   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 23   along with yat. If not, see <http://www.gnu.org/licenses/>.
570 05 Apr 06 jari 24 */
570 05 Apr 06 jari 25
2881 18 Nov 12 peter 26 #include <config.h>
2881 18 Nov 12 peter 27
1121 22 Feb 08 peter 28 #include "Matrix.h"
2210 05 Mar 10 peter 29
3657 13 Jul 17 peter 30 #include "DiagonalMatrix.h"
2210 05 Mar 10 peter 31 #include "Exception.h"
3657 13 Jul 17 peter 32 #include "SVD.h"
1121 22 Feb 08 peter 33 #include "Vector.h"
1017 01 Feb 08 peter 34 #include "VectorBase.h"
1027 02 Feb 08 peter 35 #include "VectorConstView.h"
1015 01 Feb 08 peter 36 #include "VectorView.h"
680 11 Oct 06 jari 37 #include "utility.h"
273 14 Apr 05 peter 38
1792 11 Feb 09 peter 39 #include <gsl/gsl_blas.h>
1792 11 Feb 09 peter 40
2325 22 Sep 10 peter 41 #include <algorithm>
792 11 Mar 07 jari 42 #include <cassert>
1783 08 Feb 09 peter 43 #include <climits>
293 26 Apr 05 peter 44 #include <cmath>
1792 11 Feb 09 peter 45 #include <iostream>
1920 24 Apr 09 peter 46 #include <limits>
42 26 Feb 04 jari 47 #include <sstream>
17 08 Jul 03 peter 48 #include <vector>
17 08 Jul 03 peter 49
42 26 Feb 04 jari 50 namespace theplu {
680 11 Oct 06 jari 51 namespace yat {
616 31 Aug 06 jari 52 namespace utility {
12 19 Jun 03 daniel 53
12 19 Jun 03 daniel 54
1121 22 Feb 08 peter 55   Matrix::Matrix(void)
4129 19 Jan 22 peter 56     : m_(nullptr)
703 18 Dec 06 jari 57   {
703 18 Dec 06 jari 58   }
12 19 Jun 03 daniel 59
703 18 Dec 06 jari 60
1121 22 Feb 08 peter 61   Matrix::Matrix(const size_t& r, const size_t& c, double init_value)
4129 19 Jan 22 peter 62     : m_(nullptr)
703 18 Dec 06 jari 63   {
4124 13 Jan 22 peter 64     assert(static_cast<bool>(r) == static_cast<bool>(c));
1598 24 Oct 08 peter 65     resize(r,c,init_value);
4124 13 Jan 22 peter 66     assert(rows() == r);
4124 13 Jan 22 peter 67     assert(columns() == c);
703 18 Dec 06 jari 68   }
703 18 Dec 06 jari 69
703 18 Dec 06 jari 70
1121 22 Feb 08 peter 71   Matrix::Matrix(const Matrix& o)
4129 19 Jan 22 peter 72     : m_(nullptr)
703 18 Dec 06 jari 73   {
4124 13 Jan 22 peter 74     m_ = detail::create_copy(o.gsl_matrix_p());
703 18 Dec 06 jari 75   }
703 18 Dec 06 jari 76
703 18 Dec 06 jari 77
4124 13 Jan 22 peter 78   Matrix::Matrix(const MatrixBase& o)
4129 19 Jan 22 peter 79     : m_(nullptr)
4124 13 Jan 22 peter 80   {
4124 13 Jan 22 peter 81     m_ = detail::create_copy(o.gsl_matrix_p());
4124 13 Jan 22 peter 82   }
4124 13 Jan 22 peter 83
4124 13 Jan 22 peter 84
3691 14 Sep 17 peter 85   Matrix::Matrix(Matrix&& other) noexcept
4129 19 Jan 22 peter 86     : m_(nullptr)
3585 19 Jan 17 peter 87   {
3585 19 Jan 17 peter 88     std::swap(m_, other.m_);
3585 19 Jan 17 peter 89   }
3585 19 Jan 17 peter 90
3585 19 Jan 17 peter 91
4124 13 Jan 22 peter 92   Matrix::Matrix(MatrixMutable&& other)
4129 19 Jan 22 peter 93     : m_(nullptr)
4124 13 Jan 22 peter 94   {
4378 11 Oct 23 peter 95     move_assign(std::move(other));
4124 13 Jan 22 peter 96   }
4124 13 Jan 22 peter 97
4124 13 Jan 22 peter 98
42 26 Feb 04 jari 99   // Constructor that gets data from istream
2689 27 Feb 12 peter 100   Matrix::Matrix(std::istream& is, char sep)
4129 19 Jan 22 peter 101     : m_(nullptr)
42 26 Feb 04 jari 102   {
1147 25 Feb 08 peter 103     if (!is.good())
1147 25 Feb 08 peter 104       throw utility::IO_error("Matrix: istream is not good");
1147 25 Feb 08 peter 105
42 26 Feb 04 jari 106     // read the data file and store in stl vectors (dynamically
42 26 Feb 04 jari 107     // expandable)
42 26 Feb 04 jari 108     std::vector<std::vector<double> > data_matrix;
1927 30 Apr 09 peter 109     try {
1928 30 Apr 09 peter 110       load(is, data_matrix, sep, '\n', true);
1927 30 Apr 09 peter 111     }
1928 30 Apr 09 peter 112     catch (utility::IO_error& e) {
4179 02 Jun 22 peter 113       std::stringstream ss;
4179 02 Jun 22 peter 114       ss << "Matrix(std::istream&): invalid dimensions";
4060 10 May 21 peter 115       std::throw_with_nested(IO_error(ss.str()));
1928 30 Apr 09 peter 116     }
2210 05 Mar 10 peter 117     catch (runtime_error& e) {
4179 02 Jun 22 peter 118       std::stringstream ss;
4179 02 Jun 22 peter 119       ss << "Matrix(std::istream&): invalid matrix element";
4060 10 May 21 peter 120       std::throw_with_nested(IO_error(ss.str()));
1927 30 Apr 09 peter 121     }
1927 30 Apr 09 peter 122
1928 30 Apr 09 peter 123     unsigned int nof_rows = data_matrix.size();
1147 25 Feb 08 peter 124     // if stream was empty, create nothing
1928 30 Apr 09 peter 125     if (!nof_rows)
1147 25 Feb 08 peter 126       return;
1147 25 Feb 08 peter 127
1928 30 Apr 09 peter 128     unsigned int nof_columns=data_matrix[0].size();
1928 30 Apr 09 peter 129
42 26 Feb 04 jari 130     // convert the data to a gsl matrix
3651 28 Jun 17 peter 131     m_ = detail::allocate_matrix(nof_rows, nof_columns);
4124 13 Jan 22 peter 132     assert(columns() == nof_columns);
1927 30 Apr 09 peter 133     for (size_t i=0; i<data_matrix.size(); ++i) {
1927 30 Apr 09 peter 134       assert(data_matrix[i].size()==columns());
1928 30 Apr 09 peter 135       assert(i<rows());
1928 30 Apr 09 peter 136       std::copy(data_matrix[i].begin(), data_matrix[i].end(), begin_row(i));
1927 30 Apr 09 peter 137     }
42 26 Feb 04 jari 138   }
12 19 Jun 03 daniel 139
12 19 Jun 03 daniel 140
4174 18 May 22 peter 141   Matrix::~Matrix(void)
4174 18 May 22 peter 142   {
4174 18 May 22 peter 143     detail::deallocate(m_);
4174 18 May 22 peter 144   }
4174 18 May 22 peter 145
4174 18 May 22 peter 146
4129 19 Jan 22 peter 147   void Matrix::copy_assign(const gsl_matrix* rhs)
4129 19 Jan 22 peter 148   {
4129 19 Jan 22 peter 149     // copy via blas_result and then swap result to m_
4129 19 Jan 22 peter 150     detail::copy(blas_result_, rhs);
4129 19 Jan 22 peter 151     std::swap(blas_result_, m_);
4129 19 Jan 22 peter 152   }
4129 19 Jan 22 peter 153
4129 19 Jan 22 peter 154
4129 19 Jan 22 peter 155   gsl_matrix* Matrix::gsl_matrix_p(void)
4129 19 Jan 22 peter 156   {
4129 19 Jan 22 peter 157     return m_;
4129 19 Jan 22 peter 158   }
4129 19 Jan 22 peter 159
4129 19 Jan 22 peter 160
4129 19 Jan 22 peter 161   const gsl_matrix* Matrix::gsl_matrix_p(void) const
4129 19 Jan 22 peter 162   {
4129 19 Jan 22 peter 163     return m_;
4129 19 Jan 22 peter 164   }
4129 19 Jan 22 peter 165
4129 19 Jan 22 peter 166
4129 19 Jan 22 peter 167   void Matrix::move_assign(MatrixMutable&& rhs)
4129 19 Jan 22 peter 168   {
4129 19 Jan 22 peter 169     detail::Mover mover(*this);
4129 19 Jan 22 peter 170     mover(rhs);
4129 19 Jan 22 peter 171   }
4129 19 Jan 22 peter 172
4129 19 Jan 22 peter 173
4129 19 Jan 22 peter 174   void Matrix::move_assign(gsl_matrix*&& rhs)
4129 19 Jan 22 peter 175   {
4129 19 Jan 22 peter 176     // if there is a use case for when rhs is a view, we could allow
4129 19 Jan 22 peter 177     // that with an if statement, but ATM it seems a waste.
4129 19 Jan 22 peter 178     assert(!rhs || rhs->owner);
4129 19 Jan 22 peter 179     std::swap(m_, rhs);
4129 19 Jan 22 peter 180   }
4129 19 Jan 22 peter 181
4129 19 Jan 22 peter 182
1121 22 Feb 08 peter 183   void Matrix::resize(size_t r, size_t c, double init_value)
808 15 Mar 07 peter 184   {
3651 28 Jun 17 peter 185     detail::reallocate(m_, r, c);
808 15 Mar 07 peter 186
3651 28 Jun 17 peter 187     if (m_)
1172 27 Feb 08 peter 188       all(init_value);
808 15 Mar 07 peter 189
810 15 Mar 07 jari 190     // no need to delete blas_result_ if the number of rows fit, it
810 15 Mar 07 jari 191     // may be useful later.
4124 13 Jan 22 peter 192
4124 13 Jan 22 peter 193     // Peter, actually it might be useful regardless of dimension
4124 13 Jan 22 peter 194     // since if we assign from MatrixExpression the rhs might have any
4124 13 Jan 22 peter 195     // dimension. The condition above was based on the fact that
4124 13 Jan 22 peter 196     // blas_result_ was mostly used in operator*= where number of rows
4124 13 Jan 22 peter 197     // is constant.
810 15 Mar 07 jari 198     if (blas_result_ && (blas_result_->size1!=rows())) {
3651 28 Jun 17 peter 199       detail::deallocate(blas_result_);
810 15 Mar 07 jari 200     }
810 15 Mar 07 jari 201   }
808 15 Mar 07 peter 202
810 15 Mar 07 jari 203
1121 22 Feb 08 peter 204   void Matrix::transpose(void)
42 26 Feb 04 jari 205   {
792 11 Mar 07 jari 206     assert(m_);
4124 13 Jan 22 peter 207     if (!m_) // allow empty matrix
4124 13 Jan 22 peter 208       return;
420 02 Dec 05 jari 209     if (columns()==rows())
754 17 Feb 07 jari 210       gsl_matrix_transpose(m_); // this never fails
420 02 Dec 05 jari 211     else {
3651 28 Jun 17 peter 212       gsl_matrix* transposed = detail::allocate_matrix(columns(), rows());
754 17 Feb 07 jari 213       // next line never fails if allocation above succeeded.
420 02 Dec 05 jari 214       gsl_matrix_transpose_memcpy(transposed,m_);
4124 13 Jan 22 peter 215       detail::deallocate(m_);
1098 18 Feb 08 jari 216       m_ = transposed;
3651 28 Jun 17 peter 217       detail::deallocate(blas_result_);
420 02 Dec 05 jari 218     }
42 26 Feb 04 jari 219   }
21 11 Jul 03 daniel 220
12 19 Jun 03 daniel 221
4129 19 Jan 22 peter 222   void Matrix::visit(detail::Mover& mover)
12 19 Jun 03 daniel 223   {
4129 19 Jan 22 peter 224     std::swap(mover.lhs().m_, m_);
42 26 Feb 04 jari 225   }
12 19 Jun 03 daniel 226
12 19 Jun 03 daniel 227
4129 19 Jan 22 peter 228   const Matrix& Matrix::operator=( const Matrix& other )
4124 13 Jan 22 peter 229   {
4124 13 Jan 22 peter 230     copy_assign(other);
4124 13 Jan 22 peter 231     return *this;
4124 13 Jan 22 peter 232   }
4124 13 Jan 22 peter 233
4124 13 Jan 22 peter 234
3585 19 Jan 17 peter 235   Matrix& Matrix::operator=(Matrix&& other)
3585 19 Jan 17 peter 236   {
4129 19 Jan 22 peter 237     std::swap(m_, other.m_);
3585 19 Jan 17 peter 238     return *this;
3585 19 Jan 17 peter 239   }
3585 19 Jan 17 peter 240
3585 19 Jan 17 peter 241
1121 22 Feb 08 peter 242   const Matrix& Matrix::operator*=(const Matrix& other)
42 26 Feb 04 jari 243   {
792 11 Mar 07 jari 244     assert(m_);
3651 28 Jun 17 peter 245     assert(other.columns());
4124 13 Jan 22 peter 246     assert(other.rows() == columns());
3651 28 Jun 17 peter 247     detail::reallocate(blas_result_, rows(), other.columns());
4124 13 Jan 22 peter 248     assert(detail::rows(blas_result_) == detail::rows(m_));
4124 13 Jan 22 peter 249     assert(detail::columns(blas_result_) == detail::columns(other.m_));
2661 21 Nov 11 peter 250     gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0,
2325 22 Sep 10 peter 251                    blas_result_);
2325 22 Sep 10 peter 252     std::swap(m_, blas_result_);
420 02 Dec 05 jari 253     return *this;
42 26 Feb 04 jari 254   }
20 11 Jul 03 daniel 255
12 19 Jun 03 daniel 256
4140 29 Jan 22 peter 257   void inverse_svd(const MatrixBase& A, Matrix& result)
3657 13 Jul 17 peter 258   {
3657 13 Jul 17 peter 259     assert(A.rows() == A.columns());
3657 13 Jul 17 peter 260     SVD svd(A);
3657 13 Jul 17 peter 261     svd.decompose();
4286 30 Jan 23 peter 262     svd.inverse(result);
3657 13 Jul 17 peter 263     assert(result.columns() == result.rows());
3657 13 Jul 17 peter 264     assert(result.rows() == A.rows());
3657 13 Jul 17 peter 265   }
3657 13 Jul 17 peter 266
3657 13 Jul 17 peter 267
4138 28 Jan 22 peter 268   bool nan(const MatrixBase& templat, Matrix& flag)
774 01 Mar 07 jari 269   {
1786 09 Feb 09 peter 270     if (templat.rows()!=flag.rows() || templat.columns()!=flag.columns())
1786 09 Feb 09 peter 271       flag.resize(templat.rows(),templat.columns(),1.0);
1807 18 Feb 09 peter 272     return binary_weight(templat.begin(), templat.end(), flag.begin());
774 01 Mar 07 jari 273   }
774 01 Mar 07 jari 274
774 01 Mar 07 jari 275
4129 19 Jan 22 peter 276   namespace detail {
4129 19 Jan 22 peter 277
4129 19 Jan 22 peter 278     Mover::Mover(Matrix& m)
4129 19 Jan 22 peter 279       : lhs_(m) {}
4129 19 Jan 22 peter 280
4129 19 Jan 22 peter 281
4129 19 Jan 22 peter 282     void Mover::operator()(MatrixMutable& rhs)
4129 19 Jan 22 peter 283     {
4129 19 Jan 22 peter 284       rhs.visit(*this);
4129 19 Jan 22 peter 285     }
4129 19 Jan 22 peter 286
4129 19 Jan 22 peter 287
4129 19 Jan 22 peter 288     void Mover::copy_assign(const MatrixMutable& rhs)
4129 19 Jan 22 peter 289     {
4129 19 Jan 22 peter 290       lhs_.copy_assign(rhs);
4129 19 Jan 22 peter 291     }
4129 19 Jan 22 peter 292
4129 19 Jan 22 peter 293
4129 19 Jan 22 peter 294     Matrix& Mover::lhs(void)
4129 19 Jan 22 peter 295     {
4129 19 Jan 22 peter 296       return lhs_;
4129 19 Jan 22 peter 297     }
4129 19 Jan 22 peter 298   }
4129 19 Jan 22 peter 299
680 11 Oct 06 jari 300 }}} // of namespace utility, yat and thep

Generated on Sat Oct 14 02:24:45 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)