yat/utility/MatrixBase.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 "MatrixBase.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"
1017 01 Feb 08 peter 32 #include "VectorBase.h"
1027 02 Feb 08 peter 33 #include "VectorConstView.h"
680 11 Oct 06 jari 34 #include "utility.h"
273 14 Apr 05 peter 35
1792 11 Feb 09 peter 36 #include <gsl/gsl_blas.h>
1792 11 Feb 09 peter 37
2325 22 Sep 10 peter 38 #include <algorithm>
792 11 Mar 07 jari 39 #include <cassert>
1783 08 Feb 09 peter 40 #include <climits>
293 26 Apr 05 peter 41 #include <cmath>
4119 07 Nov 21 peter 42 #include <ostream>
1920 24 Apr 09 peter 43 #include <limits>
17 08 Jul 03 peter 44
42 26 Feb 04 jari 45 namespace theplu {
680 11 Oct 06 jari 46 namespace yat {
616 31 Aug 06 jari 47 namespace utility {
12 19 Jun 03 daniel 48
4119 07 Nov 21 peter 49   MatrixBase::~MatrixBase(void)
703 18 Dec 06 jari 50   {
42 26 Feb 04 jari 51   }
12 19 Jun 03 daniel 52
12 19 Jun 03 daniel 53
4119 07 Nov 21 peter 54   MatrixBase::const_iterator MatrixBase::begin(void) const
1098 18 Feb 08 jari 55   {
4136 21 Jan 22 peter 56     return const_iterator(*this, 0, 0);
1064 10 Feb 08 peter 57   }
1064 10 Feb 08 peter 58
1064 10 Feb 08 peter 59
4119 07 Nov 21 peter 60   MatrixBase::const_column_iterator MatrixBase::begin_column(size_t i) const
1064 10 Feb 08 peter 61   {
4133 21 Jan 22 peter 62     return const_column_iterator(&(*this)(0,i), gsl_matrix_p()->tda);
1064 10 Feb 08 peter 63   }
1064 10 Feb 08 peter 64
1064 10 Feb 08 peter 65
4119 07 Nov 21 peter 66   MatrixBase::const_row_iterator MatrixBase::begin_row(size_t i) const
1064 10 Feb 08 peter 67   {
1690 31 Dec 08 peter 68     return const_row_iterator(&(*this)(i,0), 1);
1064 10 Feb 08 peter 69   }
1064 10 Feb 08 peter 70
1064 10 Feb 08 peter 71
4119 07 Nov 21 peter 72   const VectorConstView MatrixBase::column_const_view(size_t col) const
792 11 Mar 07 jari 73   {
1098 18 Feb 08 jari 74     return VectorConstView(*this, col, false);
1098 18 Feb 08 jari 75   }
792 11 Mar 07 jari 76
792 11 Mar 07 jari 77
4119 07 Nov 21 peter 78   size_t MatrixBase::columns(void) const
716 25 Dec 06 jari 79   {
4124 13 Jan 22 peter 80     return detail::columns(gsl_matrix_p());
716 25 Dec 06 jari 81   }
716 25 Dec 06 jari 82
716 25 Dec 06 jari 83
4119 07 Nov 21 peter 84   MatrixBase::const_iterator MatrixBase::end(void) const
716 25 Dec 06 jari 85   {
4136 21 Jan 22 peter 86     return const_iterator(*this, rows(), 0);
1064 10 Feb 08 peter 87   }
1064 10 Feb 08 peter 88
1064 10 Feb 08 peter 89
4119 07 Nov 21 peter 90   MatrixBase::const_column_iterator MatrixBase::end_column(size_t i) const
1064 10 Feb 08 peter 91   {
4133 21 Jan 22 peter 92     return begin_column(i) + rows();
1064 10 Feb 08 peter 93   }
1064 10 Feb 08 peter 94
1064 10 Feb 08 peter 95
4119 07 Nov 21 peter 96   MatrixBase::const_row_iterator MatrixBase::end_row(size_t i) const
1064 10 Feb 08 peter 97   {
1103 18 Feb 08 peter 98     return const_row_iterator(&(*this)(i,0)+columns(), 1);
1064 10 Feb 08 peter 99   }
1064 10 Feb 08 peter 100
1064 10 Feb 08 peter 101
4119 07 Nov 21 peter 102   bool MatrixBase::equal(const MatrixBase& other, const double d) const
272 14 Apr 05 peter 103   {
788 10 Mar 07 jari 104     if (this==&other)
788 10 Mar 07 jari 105       return true;
272 14 Apr 05 peter 106     if (columns()!=other.columns() || rows()!=other.rows())
272 14 Apr 05 peter 107       return false;
420 02 Dec 05 jari 108     for (size_t i=0; i<rows(); i++)
1682 29 Dec 08 peter 109       if (!row_const_view(i).equal(other.row_const_view(i), d))
1682 29 Dec 08 peter 110         return false;
272 14 Apr 05 peter 111     return true;
272 14 Apr 05 peter 112   }
12 19 Jun 03 daniel 113
12 19 Jun 03 daniel 114
4119 07 Nov 21 peter 115   size_t MatrixBase::rows(void) const
42 26 Feb 04 jari 116   {
4124 13 Jan 22 peter 117     return detail::rows(gsl_matrix_p());
716 25 Dec 06 jari 118   }
716 25 Dec 06 jari 119
716 25 Dec 06 jari 120
4119 07 Nov 21 peter 121   const VectorConstView MatrixBase::row_const_view(size_t col) const
1009 01 Feb 08 peter 122   {
1027 02 Feb 08 peter 123     return VectorConstView(*this, col, true);
1009 01 Feb 08 peter 124   }
1009 01 Feb 08 peter 125
1009 01 Feb 08 peter 126
4119 07 Nov 21 peter 127   const double& MatrixBase::operator()(size_t row, size_t column) const
1009 01 Feb 08 peter 128   {
813 16 Mar 07 peter 129     assert(row<rows());
813 16 Mar 07 peter 130     assert(column<columns());
4124 13 Jan 22 peter 131     const double* d=gsl_matrix_const_ptr(gsl_matrix_p(), row, column);
755 18 Feb 07 jari 132     if (!d)
4119 07 Nov 21 peter 133       throw utility::GSL_error("MatrixBase::operator()",GSL_EINVAL);
755 18 Feb 07 jari 134     return *d;
716 25 Dec 06 jari 135   }
716 25 Dec 06 jari 136
716 25 Dec 06 jari 137
4119 07 Nov 21 peter 138   bool MatrixBase::operator==(const MatrixBase& other) const
716 25 Dec 06 jari 139   {
716 25 Dec 06 jari 140     return equal(other);
716 25 Dec 06 jari 141   }
421 02 Dec 05 jari 142
716 25 Dec 06 jari 143
4119 07 Nov 21 peter 144   bool MatrixBase::operator!=(const MatrixBase& other) const
716 25 Dec 06 jari 145   {
716 25 Dec 06 jari 146     return !equal(other);
716 25 Dec 06 jari 147   }
716 25 Dec 06 jari 148
716 25 Dec 06 jari 149
4119 07 Nov 21 peter 150   bool isnull(const MatrixBase& other)
12 19 Jun 03 daniel 151   {
792 11 Mar 07 jari 152     return gsl_matrix_isnull(other.gsl_matrix_p());
774 01 Mar 07 jari 153   }
774 01 Mar 07 jari 154
774 01 Mar 07 jari 155
4119 07 Nov 21 peter 156   double max(const MatrixBase& other)
774 01 Mar 07 jari 157   {
792 11 Mar 07 jari 158     return gsl_matrix_max(other.gsl_matrix_p());
774 01 Mar 07 jari 159   }
774 01 Mar 07 jari 160
774 01 Mar 07 jari 161
4119 07 Nov 21 peter 162   double min(const MatrixBase& other)
774 01 Mar 07 jari 163   {
792 11 Mar 07 jari 164     return gsl_matrix_min(other.gsl_matrix_p());
774 01 Mar 07 jari 165   }
774 01 Mar 07 jari 166
774 01 Mar 07 jari 167
4119 07 Nov 21 peter 168   void minmax_index(const MatrixBase& other,
4119 07 Nov 21 peter 169                     std::pair<size_t,size_t>& min,
4119 07 Nov 21 peter 170                     std::pair<size_t,size_t>& max)
774 01 Mar 07 jari 171   {
792 11 Mar 07 jari 172     gsl_matrix_minmax_index(other.gsl_matrix_p(), &min.first, &min.second,
774 01 Mar 07 jari 173                             &max.first, &max.second);
774 01 Mar 07 jari 174   }
774 01 Mar 07 jari 175
774 01 Mar 07 jari 176
4119 07 Nov 21 peter 177   std::ostream& operator<<(std::ostream& s, const MatrixBase& m)
774 01 Mar 07 jari 178   {
420 02 Dec 05 jari 179     s.setf(std::ios::dec);
1919 24 Apr 09 peter 180     std::streamsize precision = s.precision();
1920 24 Apr 09 peter 181     s.precision(std::numeric_limits<double>().digits10);
2748 18 Jun 12 peter 182     for(size_t i=0, j=0; i<m.rows(); i++) {
2748 18 Jun 12 peter 183       s << m(i,0);
2748 18 Jun 12 peter 184       for (j=1; j<m.columns(); ++j)
2748 18 Jun 12 peter 185         s << s.fill() << m(i,j);
2748 18 Jun 12 peter 186       if (i<m.rows()-1)
2748 18 Jun 12 peter 187         s << "\n";
2748 18 Jun 12 peter 188     }
1919 24 Apr 09 peter 189     s.precision(precision);
420 02 Dec 05 jari 190     return s;
42 26 Feb 04 jari 191   }
12 19 Jun 03 daniel 192
737 07 Jan 07 peter 193
4119 07 Nov 21 peter 194   double trace(const MatrixBase& m)
2735 29 May 12 peter 195   {
2735 29 May 12 peter 196     size_t n = std::min(m.rows(), m.columns());
2735 29 May 12 peter 197     double sum = 0.0;
2735 29 May 12 peter 198     for (size_t i=0; i<n; ++i)
2735 29 May 12 peter 199       sum += m(i,i);
2735 29 May 12 peter 200     return sum;
2735 29 May 12 peter 201   }
2735 29 May 12 peter 202
3651 28 Jun 17 peter 203
680 11 Oct 06 jari 204 }}} // of namespace utility, yat and thep

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