Flow123d  release_2.1.2-337-g6b7a56b
matrix.cpp
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file matrix.cpp
15  * @brief
16  */
17 
18 #include <cmath>
19 #include <cstring>
20 #include "system/exc_common.hh"
23 
24 using namespace mathfce;
25 
26 TMatrix::TMatrix(int size) {
27  nc = size;
28  nr = size;
29  elm = new double[ nc * nr ];
30 }
31 
32 TMatrix::TMatrix(int num_rows, int num_cols) {
33  nc = num_cols;
34  nr = num_rows;
35  elm = new double[ nc * nr ];
36 }
37 
39 {
40  nc = x.nc;
41  nr = x.nr;
42  elm = new double[nc*nr];
43  memcpy(elm,x.elm,nc*nr*sizeof(double));
44 }
45 
47  delete[] elm;
48 }
49 
50 TMVector::TMVector(int size) {
51  this->size = size;
52  elm = new double[size];
53 }
54 
56 {
57  size=x.size;
58  elm = new double[size];
59  memcpy(elm,x.elm,size*sizeof(double));
60 }
61 
63  delete[] elm;
64 }
65 
66 std::ostream & operator <<(std::ostream& stream, const TMatrix& M) {
67  for (int i = 1; i <= M.nr; i++) {
68  for (int j = 1; j <= M.nc; j++) {
69  stream << M.Get(i, j) << " ";
70  }
71  stream << "\n";
72  }
73  return stream;
74 }
75 
76 std::ostream & operator <<(std::ostream& stream, const TMVector& V) {
77  for (int i = 1; i <= V.size; i++) {
78  stream << V.elm[ i - 1 ] << "\n";
79  }
80  return stream;
81 }
82 
83 void TMatrix::Set(int row, int col, double value) {
84  if (row > nr)
85  THROW( ExcAssertMsg() << EI_Message( "Number of the row is greater than number of rows in the matrix.") );
86  if (col > nc)
87  THROW( ExcAssertMsg() << EI_Message( "Number of the column is greater than number of columns in the matrix.") );
88  elm[ (row - 1) * nc + col - 1 ] = value;
89  return;
90 }
91 
92 double TMatrix::Get(int row, int col) const {
93  if (row > nr)
94  THROW( ExcAssertMsg() << EI_Message( "Number of the row is greater than number of rows in the matrix.") );
95  if (col > nc)
96  THROW( ExcAssertMsg() << EI_Message( "Number of the column is greater than number of columns in the matrix.") );
97  return elm[ (row - 1) * nc + col - 1 ];
98 }
99 
100 void TMatrix::SwapRows(int r1, int r2) {
101  if (r1 > nr || r2 > nr) {
102  THROW( ExcAssertMsg() << EI_Message( "Number of the row is greater than number of rows in the matrix.") );
103  }
104 
105  for (int i = 1; i <= nc; i++) {
106  double tmp = Get(r1, i);
107  Set(r1, i, Get(r2, i));
108  Set(r2, i, tmp);
109  }
110 
111  return;
112 }
113 
114 void TMVector::SwapElements(int i1, int i2) {
115  if (i1 > size || i2 > size) {
116  THROW( ExcAssertMsg() << EI_Message( "Number of the element is greater than size of the vector.") );
117  }
118 
119  double tmp = elm[ i1 - 1 ];
120  elm[ i1 - 1 ] = elm[ i2 - 1 ];
121  elm[ i2 - 1 ] = tmp;
122 
123  return;
124 }
125 
126 double TMVector::Get(int i) {
127  if (i > size) {
128  THROW( ExcAssertMsg() << EI_Message("Number of the element is greater than size of the vector.") );
129  }
130 
131 
132  return elm[ i - 1 ];
133 }
134 
135 void TMVector::Set(int i, double value) {
136  if (i > size)
137  THROW( ExcAssertMsg() << EI_Message("Number of the element is greater than size of the vector.") );
138  elm[ i - 1 ] = value;
139  return;
140 }
141 
142 int TMatrix::NRows() const {
143  return nr;
144 }
145 
146 int TMatrix::NCols() const {
147  return nc;
148 }
149 
150 TNSolutions Gauss(const TMatrix& A, TMVector* X, const TMVector& B) {
151  TMatrix M(A);
152  TMVector b(B);
153 
154  for (int i = 1; i < M.NRows(); i++) {
155  double tmp = fabs(M.Get(i, i));
156  int row = -1;
157  for (int j = i + 1; j <= M.NRows(); j++)
158  if (fabs(M.Get(j, i)) > tmp) {
159  tmp = fabs(M.Get(j, i));
160  row = j;
161  }
162  if (tmp < epsilon) {
163  continue;
164  }
165 
166  if (row != -1) {
167  M.SwapRows(i, row);
168  b.SwapElements(i, row);
169  }
170 
171  for (int j = i + 1; j <= M.NRows(); j++) {
172  tmp = M.Get(j, i) / M.Get(i, i);
173  for (int k = i; k <= M.NCols(); k++) {
174  M.Set(j, k, M.Get(j, k) - tmp * M.Get(i, k));
175  }
176  b.Set(j, b.Get(j) - tmp * b.Get(i));
177  }
178  }
179 
180  for (int i = M.NRows(); i >= 1; i--) {
181  double tmp = b.Get(i);
182 
183  for (int k = i + 1; k <= M.NCols(); k++) {
184  tmp -= M.Get(i, k) * X->Get(k);
185  }
186 
187  if (IsZero(M.Get(i, i))) {
188  if (IsZero(tmp)) {
189  return inf_solutions;
190  } else {
191  return no_solution;
192  }
193  }
194 
195  X->Set(i, tmp / M.Get(i, i));
196  }
197 
198  return one_solution;
199 }
~TMVector()
Definition: matrix.cpp:62
double * elm
Definition: matrix.h:36
TNSolutions
Definition: matrix.h:23
int NCols() const
Definition: matrix.cpp:146
TMVector(int)
Definition: matrix.cpp:50
void Set(int, double)
Definition: matrix.cpp:135
int NRows() const
Definition: matrix.cpp:142
TNSolutions Gauss(const TMatrix &A, TMVector *X, const TMVector &B)
Definition: matrix.cpp:150
double Get(int)
Definition: matrix.cpp:126
bool IsZero(double)
Definition: mathfce.cpp:22
double Get(int, int) const
Definition: matrix.cpp:92
static constexpr bool value
Definition: json.hpp:87
~TMatrix()
Definition: matrix.cpp:46
void SwapElements(int, int)
Definition: matrix.cpp:114
void SwapRows(int, int)
Definition: matrix.cpp:100
Definition: matrix.h:31
const double epsilon
Definition: mathfce.h:23
TMatrix(int)
Definition: matrix.cpp:26
void Set(int, int, double)
Definition: matrix.cpp:83
int nc
Definition: matrix.h:33
double * elm
Definition: matrix.h:59
std::ostream & operator<<(std::ostream &stream, const TMatrix &M)
Definition: matrix.cpp:66
int nr
Definition: matrix.h:34
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
int size
Definition: matrix.h:57