yat/utility/MatrixMutable.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 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
4119 07 Nov 21 peter 28 #include "MatrixMutable.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
4119 07 Nov 21 peter 55   MatrixMutable::MatrixMutable(void)
4129 19 Jan 22 peter 56     : blas_result_(nullptr)
42 26 Feb 04 jari 57   {
42 26 Feb 04 jari 58   }
12 19 Jun 03 daniel 59
12 19 Jun 03 daniel 60
4119 07 Nov 21 peter 61   MatrixMutable::~MatrixMutable(void)
42 26 Feb 04 jari 62   {
3651 28 Jun 17 peter 63     detail::deallocate(blas_result_);
42 26 Feb 04 jari 64   }
12 19 Jun 03 daniel 65
12 19 Jun 03 daniel 66
4119 07 Nov 21 peter 67   void MatrixMutable::all(const double value)
1098 18 Feb 08 jari 68   {
4124 13 Jan 22 peter 69     assert(gsl_matrix_p());
4124 13 Jan 22 peter 70     gsl_matrix_set_all(gsl_matrix_p(), value);
1098 18 Feb 08 jari 71   }
1098 18 Feb 08 jari 72
1098 18 Feb 08 jari 73
4119 07 Nov 21 peter 74   MatrixMutable::iterator MatrixMutable::begin(void)
1064 10 Feb 08 peter 75   {
4136 21 Jan 22 peter 76     return iterator(*this, 0, 0);
1064 10 Feb 08 peter 77   }
1064 10 Feb 08 peter 78
1064 10 Feb 08 peter 79
4119 07 Nov 21 peter 80   MatrixMutable::column_iterator MatrixMutable::begin_column(size_t i)
1064 10 Feb 08 peter 81   {
4133 21 Jan 22 peter 82     return column_iterator(&(*this)(0,i), gsl_matrix_p()->tda);
1064 10 Feb 08 peter 83   }
1064 10 Feb 08 peter 84
1064 10 Feb 08 peter 85
4119 07 Nov 21 peter 86   MatrixMutable::row_iterator MatrixMutable::begin_row(size_t i)
1064 10 Feb 08 peter 87   {
1690 31 Dec 08 peter 88     return row_iterator(&(*this)(i,0), 1);
1064 10 Feb 08 peter 89   }
1064 10 Feb 08 peter 90
1064 10 Feb 08 peter 91
4119 07 Nov 21 peter 92   VectorView MatrixMutable::column_view(size_t col)
1064 10 Feb 08 peter 93   {
1098 18 Feb 08 jari 94     VectorView res(*this, col, false);
1098 18 Feb 08 jari 95     return res;
1098 18 Feb 08 jari 96   }
792 11 Mar 07 jari 97
792 11 Mar 07 jari 98
4119 07 Nov 21 peter 99   void MatrixMutable::div(const MatrixBase& other)
1098 18 Feb 08 jari 100   {
4124 13 Jan 22 peter 101     assert(gsl_matrix_p());
4124 13 Jan 22 peter 102     int status=gsl_matrix_div_elements(gsl_matrix_p(), other.gsl_matrix_p());
755 18 Feb 07 jari 103     if (status)
3743 12 Jul 18 peter 104       throw utility::GSL_error("matrix::div_elements",status);
716 25 Dec 06 jari 105   }
716 25 Dec 06 jari 106
716 25 Dec 06 jari 107
4119 07 Nov 21 peter 108   MatrixMutable::iterator MatrixMutable::end(void)
1064 10 Feb 08 peter 109   {
4136 21 Jan 22 peter 110     return iterator(*this, rows(), 0);
1064 10 Feb 08 peter 111   }
1064 10 Feb 08 peter 112
1064 10 Feb 08 peter 113
4119 07 Nov 21 peter 114   MatrixMutable::column_iterator MatrixMutable::end_column(size_t i)
1064 10 Feb 08 peter 115   {
4133 21 Jan 22 peter 116     return begin_column(i) + rows();
1064 10 Feb 08 peter 117   }
1064 10 Feb 08 peter 118
1064 10 Feb 08 peter 119
4119 07 Nov 21 peter 120   MatrixMutable::row_iterator MatrixMutable::end_row(size_t i)
1064 10 Feb 08 peter 121   {
1103 18 Feb 08 peter 122     return row_iterator(&(*this)(i,0)+columns(), 1);
1064 10 Feb 08 peter 123   }
1064 10 Feb 08 peter 124
1064 10 Feb 08 peter 125
4119 07 Nov 21 peter 126   void MatrixMutable::mul(const MatrixBase& other)
42 26 Feb 04 jari 127   {
4124 13 Jan 22 peter 128     assert(gsl_matrix_p());
4124 13 Jan 22 peter 129     int status=gsl_matrix_mul_elements(gsl_matrix_p(), other.gsl_matrix_p());
755 18 Feb 07 jari 130     if (status)
4119 07 Nov 21 peter 131       throw utility::GSL_error("MatrixMutable::mul_elements",status);
716 25 Dec 06 jari 132   }
516 20 Feb 06 peter 133
716 25 Dec 06 jari 134
4119 07 Nov 21 peter 135   VectorView MatrixMutable::row_view(size_t row)
808 15 Mar 07 peter 136   {
4124 13 Jan 22 peter 137     assert(row < rows());
1015 01 Feb 08 peter 138     VectorView res(*this, row, true);
1009 01 Feb 08 peter 139     return res;
1009 01 Feb 08 peter 140   }
1009 01 Feb 08 peter 141
1009 01 Feb 08 peter 142
4119 07 Nov 21 peter 143   void MatrixMutable::swap_columns(const size_t i, const size_t j)
716 25 Dec 06 jari 144   {
4124 13 Jan 22 peter 145     assert(gsl_matrix_p());
4124 13 Jan 22 peter 146     int status=gsl_matrix_swap_columns(gsl_matrix_p(), i, j);
755 18 Feb 07 jari 147     if (status)
4119 07 Nov 21 peter 148       throw utility::GSL_error("MatrixMutable::swap_columns",status);
716 25 Dec 06 jari 149   }
716 25 Dec 06 jari 150
716 25 Dec 06 jari 151
4119 07 Nov 21 peter 152   void MatrixMutable::swap_rowcol(const size_t i, const size_t j)
716 25 Dec 06 jari 153   {
4124 13 Jan 22 peter 154     assert(gsl_matrix_p());
4124 13 Jan 22 peter 155     int status=gsl_matrix_swap_rowcol(gsl_matrix_p(), i, j);
755 18 Feb 07 jari 156     if (status)
4119 07 Nov 21 peter 157       throw utility::GSL_error("MatrixMutable::swap_rowcol",status);
716 25 Dec 06 jari 158   }
716 25 Dec 06 jari 159
716 25 Dec 06 jari 160
4119 07 Nov 21 peter 161   void MatrixMutable::swap_rows(const size_t i, const size_t j)
716 25 Dec 06 jari 162   {
4124 13 Jan 22 peter 163     assert(gsl_matrix_p());
4124 13 Jan 22 peter 164     int status=gsl_matrix_swap_rows(gsl_matrix_p(), i, j);
755 18 Feb 07 jari 165     if (status)
4119 07 Nov 21 peter 166       throw utility::GSL_error("MatrixMutable::swap_rows",status);
716 25 Dec 06 jari 167   }
716 25 Dec 06 jari 168
716 25 Dec 06 jari 169
4124 13 Jan 22 peter 170   MatrixMutable& MatrixMutable::operator=(const MatrixBase& other)
4124 13 Jan 22 peter 171   {
4124 13 Jan 22 peter 172     copy_assign(other);
4124 13 Jan 22 peter 173     return *this;
4124 13 Jan 22 peter 174   }
4124 13 Jan 22 peter 175
4124 13 Jan 22 peter 176
4124 13 Jan 22 peter 177   MatrixMutable& MatrixMutable::operator=(const MatrixMutable& other)
4124 13 Jan 22 peter 178   {
4124 13 Jan 22 peter 179     copy_assign(other);
4124 13 Jan 22 peter 180     return *this;
4124 13 Jan 22 peter 181   }
4124 13 Jan 22 peter 182
4124 13 Jan 22 peter 183
4124 13 Jan 22 peter 184   MatrixMutable& MatrixMutable::operator=(MatrixMutable&& other)
4124 13 Jan 22 peter 185   {
4124 13 Jan 22 peter 186     move_assign(std::move(other));
4124 13 Jan 22 peter 187     return *this;
4124 13 Jan 22 peter 188   }
4124 13 Jan 22 peter 189
4124 13 Jan 22 peter 190
4119 07 Nov 21 peter 191   double& MatrixMutable::operator()(size_t row, size_t column)
42 26 Feb 04 jari 192   {
4124 13 Jan 22 peter 193     assert(gsl_matrix_p());
813 16 Mar 07 peter 194     assert(row<rows());
813 16 Mar 07 peter 195     assert(column<columns());
4124 13 Jan 22 peter 196     double* d=gsl_matrix_ptr(gsl_matrix_p(), row, column);
755 18 Feb 07 jari 197     if (!d)
4119 07 Nov 21 peter 198       throw utility::GSL_error("MatrixMutable::operator()",GSL_EINVAL);
755 18 Feb 07 jari 199     return *d;
716 25 Dec 06 jari 200   }
12 19 Jun 03 daniel 201
716 25 Dec 06 jari 202
4119 07 Nov 21 peter 203   const MatrixMutable& MatrixMutable::operator+=(const MatrixBase& other)
716 25 Dec 06 jari 204   {
4124 13 Jan 22 peter 205     assert(gsl_matrix_p());
4124 13 Jan 22 peter 206     int status=gsl_matrix_add(gsl_matrix_p(), other.gsl_matrix_p());
773 01 Mar 07 jari 207     if (status)
4119 07 Nov 21 peter 208       throw utility::GSL_error("MatrixMutable::operator+=", status);
703 18 Dec 06 jari 209     return *this;
703 18 Dec 06 jari 210   }
20 11 Jul 03 daniel 211
703 18 Dec 06 jari 212
4119 07 Nov 21 peter 213   const MatrixMutable& MatrixMutable::operator+=(const double d)
703 18 Dec 06 jari 214   {
4124 13 Jan 22 peter 215     assert(gsl_matrix_p());
4124 13 Jan 22 peter 216     gsl_matrix_add_constant(gsl_matrix_p(), d);
703 18 Dec 06 jari 217     return *this;
703 18 Dec 06 jari 218   }
703 18 Dec 06 jari 219
703 18 Dec 06 jari 220
4119 07 Nov 21 peter 221   const MatrixMutable& MatrixMutable::operator-=(const MatrixBase& other)
703 18 Dec 06 jari 222   {
4124 13 Jan 22 peter 223     assert(gsl_matrix_p());
4124 13 Jan 22 peter 224     int status=gsl_matrix_sub(gsl_matrix_p(), other.gsl_matrix_p());
773 01 Mar 07 jari 225     if (status)
4119 07 Nov 21 peter 226       throw utility::GSL_error("MatrixMutable::operator-=", status);
703 18 Dec 06 jari 227     return *this;
703 18 Dec 06 jari 228   }
703 18 Dec 06 jari 229
703 18 Dec 06 jari 230
4119 07 Nov 21 peter 231   const MatrixMutable& MatrixMutable::operator-=(const double d)
773 01 Mar 07 jari 232   {
4124 13 Jan 22 peter 233     assert(gsl_matrix_p());
4124 13 Jan 22 peter 234     gsl_matrix_add_constant(gsl_matrix_p(), -d);
773 01 Mar 07 jari 235     return *this;
773 01 Mar 07 jari 236   }
773 01 Mar 07 jari 237
773 01 Mar 07 jari 238
4119 07 Nov 21 peter 239   const MatrixMutable& MatrixMutable::operator*=(const double d)
42 26 Feb 04 jari 240   {
4124 13 Jan 22 peter 241     assert(gsl_matrix_p());
4124 13 Jan 22 peter 242     gsl_matrix_scale(gsl_matrix_p(), d);
703 18 Dec 06 jari 243     return *this;
703 18 Dec 06 jari 244   }
12 19 Jun 03 daniel 245
703 18 Dec 06 jari 246
4129 19 Jan 22 peter 247   void MatrixMutable::copy_assign(const MatrixBase& other)
4129 19 Jan 22 peter 248   {
4129 19 Jan 22 peter 249     if (this != &other)
4129 19 Jan 22 peter 250       copy_assign(other.gsl_matrix_p());
4129 19 Jan 22 peter 251   }
4129 19 Jan 22 peter 252
4129 19 Jan 22 peter 253
4129 19 Jan 22 peter 254   void MatrixMutable::move_assign(MatrixMutable&& other)
4129 19 Jan 22 peter 255   {
4129 19 Jan 22 peter 256     copy_assign(other);
4129 19 Jan 22 peter 257   }
4129 19 Jan 22 peter 258
4129 19 Jan 22 peter 259
4129 19 Jan 22 peter 260   void MatrixMutable::visit(detail::Mover& mover)
4129 19 Jan 22 peter 261   {
4129 19 Jan 22 peter 262     mover.copy_assign(*this);
4129 19 Jan 22 peter 263   }
4129 19 Jan 22 peter 264
4129 19 Jan 22 peter 265
4119 07 Nov 21 peter 266   void swap(MatrixMutable& a, MatrixMutable& b)
3657 13 Jul 17 peter 267   {
2689 27 Feb 12 peter 268     assert(a.gsl_matrix_p());
2689 27 Feb 12 peter 269     assert(b.gsl_matrix_p());
774 01 Mar 07 jari 270     int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p());
774 01 Mar 07 jari 271     if (status)
3743 12 Jul 18 peter 272       throw utility::GSL_error("swap(Matrix&,Matrix&)",status);
774 01 Mar 07 jari 273   }
774 01 Mar 07 jari 274
774 01 Mar 07 jari 275
680 11 Oct 06 jari 276 }}} // of namespace utility, yat and thep

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