yat/utility/BasicMatrix.cc

Code
Comments
Other
Rev Date Author Line
3652 28 Jun 17 peter 1 // $Id$
3652 28 Jun 17 peter 2
3652 28 Jun 17 peter 3 /*
4359 23 Aug 23 peter 4   Copyright (C) 2017, 2022 Peter Johansson
3652 28 Jun 17 peter 5
3652 28 Jun 17 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
3652 28 Jun 17 peter 7
3652 28 Jun 17 peter 8   The yat library is free software; you can redistribute it and/or
3652 28 Jun 17 peter 9   modify it under the terms of the GNU General Public License as
3652 28 Jun 17 peter 10   published by the Free Software Foundation; either version 3 of the
3652 28 Jun 17 peter 11   License, or (at your option) any later version.
3652 28 Jun 17 peter 12
3652 28 Jun 17 peter 13   The yat library is distributed in the hope that it will be useful,
3652 28 Jun 17 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
3652 28 Jun 17 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3652 28 Jun 17 peter 16   General Public License for more details.
3652 28 Jun 17 peter 17
3652 28 Jun 17 peter 18   You should have received a copy of the GNU General Public License
3652 28 Jun 17 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
3652 28 Jun 17 peter 20 */
3652 28 Jun 17 peter 21
3652 28 Jun 17 peter 22 #include <config.h>
3652 28 Jun 17 peter 23
3652 28 Jun 17 peter 24 #include "BasicMatrix.h"
3652 28 Jun 17 peter 25 #include "Exception.h"
3652 28 Jun 17 peter 26
3652 28 Jun 17 peter 27 #include <gsl/gsl_blas.h>
3652 28 Jun 17 peter 28 #include <gsl/gsl_matrix.h>
3652 28 Jun 17 peter 29
3652 28 Jun 17 peter 30 #include <cassert>
3652 28 Jun 17 peter 31 #include <sstream>
3652 28 Jun 17 peter 32
3652 28 Jun 17 peter 33 namespace theplu {
3652 28 Jun 17 peter 34 namespace yat {
3652 28 Jun 17 peter 35 namespace utility {
3652 28 Jun 17 peter 36 namespace detail {
3652 28 Jun 17 peter 37
3652 28 Jun 17 peter 38   gsl_matrix* allocate_matrix(size_t row, size_t col)
3652 28 Jun 17 peter 39   {
3652 28 Jun 17 peter 40     gsl_matrix* result = NULL;
3652 28 Jun 17 peter 41     if (row && col) {
3652 28 Jun 17 peter 42       result = gsl_matrix_alloc(row, col);
3652 28 Jun 17 peter 43       if (!result)
3652 28 Jun 17 peter 44         throw GSL_error("gsl_matrix_alloc failed");
3652 28 Jun 17 peter 45     }
3652 28 Jun 17 peter 46     else if (row || col) {
3652 28 Jun 17 peter 47       std::stringstream msg;
3652 28 Jun 17 peter 48       msg << "allocate_matrix(" << row << ", " << col << "): "
3652 28 Jun 17 peter 49           << "incorrect dimensions\n";
3652 28 Jun 17 peter 50       throw GSL_error(msg.str());
3652 28 Jun 17 peter 51     }
3652 28 Jun 17 peter 52     return result;
3652 28 Jun 17 peter 53   }
3652 28 Jun 17 peter 54
3652 28 Jun 17 peter 55
3652 28 Jun 17 peter 56   void assign(gsl_matrix* dest, const gsl_matrix* src)
3652 28 Jun 17 peter 57   {
3652 28 Jun 17 peter 58     assert(dest);
3652 28 Jun 17 peter 59     assert(src);
3652 28 Jun 17 peter 60     assert(rows(dest) == rows(src));
3652 28 Jun 17 peter 61     assert(columns(dest) == columns(src));
3652 28 Jun 17 peter 62     if (gsl_matrix_memcpy(dest, src))
3652 28 Jun 17 peter 63       throw utility::GSL_error("copy Matrix dimension mis-match");
3652 28 Jun 17 peter 64   }
3652 28 Jun 17 peter 65
3652 28 Jun 17 peter 66
3652 28 Jun 17 peter 67   size_t columns(const gsl_matrix* p)
3652 28 Jun 17 peter 68   {
3652 28 Jun 17 peter 69     return (p ? p->size2 : 0);
3652 28 Jun 17 peter 70   }
3652 28 Jun 17 peter 71
3652 28 Jun 17 peter 72
3652 28 Jun 17 peter 73   void copy(gsl_matrix*& dest, const gsl_matrix* src)
3652 28 Jun 17 peter 74   {
3652 28 Jun 17 peter 75     if (!src) {
3652 28 Jun 17 peter 76       deallocate(dest);
3652 28 Jun 17 peter 77       return;
3652 28 Jun 17 peter 78     }
3652 28 Jun 17 peter 79     reallocate(dest, rows(src), columns(src));
3652 28 Jun 17 peter 80     assign(dest, src);
3652 28 Jun 17 peter 81   }
3652 28 Jun 17 peter 82
3652 28 Jun 17 peter 83
3652 28 Jun 17 peter 84   gsl_matrix* create_copy(const gsl_matrix* p)
3652 28 Jun 17 peter 85   {
3652 28 Jun 17 peter 86     if (!p)
3652 28 Jun 17 peter 87       return NULL;
3652 28 Jun 17 peter 88     gsl_matrix* result = NULL;
3652 28 Jun 17 peter 89     copy(result, p);
3652 28 Jun 17 peter 90     return result;
3652 28 Jun 17 peter 91   }
3652 28 Jun 17 peter 92
3652 28 Jun 17 peter 93
3652 28 Jun 17 peter 94   void deallocate(gsl_matrix*& p)
3652 28 Jun 17 peter 95   {
3907 09 May 20 peter 96     gsl_matrix_free(p);
3907 09 May 20 peter 97     p = NULL;
3652 28 Jun 17 peter 98   }
3652 28 Jun 17 peter 99
3652 28 Jun 17 peter 100
4124 13 Jan 22 peter 101   void move_if_owner(gsl_matrix*& dest, gsl_matrix*& src)
4124 13 Jan 22 peter 102   {
4124 13 Jan 22 peter 103     // Note that dest and src can be one of 1) nullptr, 2) a view, or
4124 13 Jan 22 peter 104     // 3) owner of its data, which means we have nine different cases.
4124 13 Jan 22 peter 105
4124 13 Jan 22 peter 106     // If src is nullptr copy is fast anyway
4124 13 Jan 22 peter 107     //
4124 13 Jan 22 peter 108     // We cannot move from a view bcs that create a very unexpected
4124 13 Jan 22 peter 109     // link between dest and the owner of the src data.
4124 13 Jan 22 peter 110     //
4124 13 Jan 22 peter 111     // If dest is a view, we must copy becasue we cannot break the
4124 13 Jan 22 peter 112     // link between the view and its owner.
4124 13 Jan 22 peter 113     if (!src || !src->owner || (dest && !dest->owner))
4124 13 Jan 22 peter 114       copy(dest, src);
4124 13 Jan 22 peter 115     // that leaves two cases:
4124 13 Jan 22 peter 116     // a) owner -> owner
4124 13 Jan 22 peter 117     // b) owner -> nullptr
4124 13 Jan 22 peter 118     // where a) is straightforward; for b) the dest nullptr from a
4124 13 Jan 22 peter 119     // Matrix (or similar class owning its data) - it cannot be from a
4124 13 Jan 22 peter 120     // MatrixView because a MatrixView is not allowed to alter its
4124 13 Jan 22 peter 121     // size.
4124 13 Jan 22 peter 122     else {
4124 13 Jan 22 peter 123       assert(src && src->owner);
4200 19 Aug 22 peter 124       assert(dest==nullptr || dest->owner);
4124 13 Jan 22 peter 125       std::swap(dest, src);
4124 13 Jan 22 peter 126     }
4124 13 Jan 22 peter 127   }
4124 13 Jan 22 peter 128
4124 13 Jan 22 peter 129
3652 28 Jun 17 peter 130   void reallocate(gsl_matrix*& p, size_t row, size_t col)
3652 28 Jun 17 peter 131   {
4129 19 Jan 22 peter 132     assert(p==nullptr || p->owner);
4124 13 Jan 22 peter 133     assert(static_cast<bool>(row) == static_cast<bool>(col));
3652 28 Jun 17 peter 134     if (row!=rows(p) || col!=columns(p)) {
3652 28 Jun 17 peter 135       gsl_matrix* tmp = allocate_matrix(row, col);
3652 28 Jun 17 peter 136       // allocation above might throw, so wait modifying \a p until
3652 28 Jun 17 peter 137       // we know no exception was thrown.
3652 28 Jun 17 peter 138       deallocate(p);
3652 28 Jun 17 peter 139       p = tmp;
3652 28 Jun 17 peter 140     }
3652 28 Jun 17 peter 141   }
3652 28 Jun 17 peter 142
3652 28 Jun 17 peter 143
3652 28 Jun 17 peter 144   size_t rows(const gsl_matrix* p)
3652 28 Jun 17 peter 145   {
3652 28 Jun 17 peter 146     return (p ? p->size1 : 0);
3652 28 Jun 17 peter 147   }
3652 28 Jun 17 peter 148
3652 28 Jun 17 peter 149
4129 19 Jan 22 peter 150   const double* end(const gsl_matrix* p)
4129 19 Jan 22 peter 151   {
4129 19 Jan 22 peter 152     assert(p);
4129 19 Jan 22 peter 153     assert(p->block);
4129 19 Jan 22 peter 154     return p->data + p->block->size;
4129 19 Jan 22 peter 155   }
4129 19 Jan 22 peter 156
4129 19 Jan 22 peter 157
4129 19 Jan 22 peter 158   bool overlap(const gsl_matrix* lhs, const gsl_matrix* rhs)
4129 19 Jan 22 peter 159   {
4129 19 Jan 22 peter 160     if (!lhs || !rhs)
4129 19 Jan 22 peter 161       return false;
4129 19 Jan 22 peter 162
4129 19 Jan 22 peter 163     return (rhs->data >= lhs->data && rhs->data < end(lhs)) ||
4129 19 Jan 22 peter 164       (lhs->data >= rhs->data && lhs->data < end(rhs));
4129 19 Jan 22 peter 165   }
4129 19 Jan 22 peter 166
4129 19 Jan 22 peter 167
3652 28 Jun 17 peter 168 }}}} // of namespace detail, utility, yat and thep

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