HermesCommon  3.0
dense_matrix_operations.h
Go to the documentation of this file.
1 // This file is part of HermesCommon
2 //
3 // Copyright (c) 2009 hp-FEM group at the University of Nevada, Reno (UNR).
5 //
6 // Hermes2D is free software; you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published
8 // by the Free Software Foundation; either version 2 of the License,
9 // or (at your option) any later version.
10 //
11 // Hermes2D is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with Hermes2D; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
22 #ifndef __HERMES_COMMON_DENSE_MATRIX_OPERATIONS_H
23 #define __HERMES_COMMON_DENSE_MATRIX_OPERATIONS_H
24
25 #include "common.h"
26 #include "util/compat.h"
27 #include "util/memory_handling.h"
28 #include "exceptions.h"
29
30 namespace Hermes
31 {
33  namespace Algebra
34  {
36  namespace DenseMatrixOperations
37  {
41  template<typename T>
42  T **new_matrix(unsigned int m, unsigned int n = 0)
43  {
44  if (!n) n = m;
45  T **vec = (T**)calloc_with_check_direct_size<char>(sizeof(T *)* m + sizeof(T)* m * n);
46  T *row = (T *)(vec + m);
47  for (unsigned int i = 0; i < m; i++, row += n)
48  vec[i] = row;
49  return vec;
50  }
51
52  template<typename T>
53  T **new_matrix_malloc(unsigned int m, unsigned int n = 0)
54  {
55  if (!n) n = m;
56  T **vec = (T **)malloc(sizeof(T *)* m + sizeof(T)* m * n);
57  memset(vec, 0, sizeof(T *)* m + sizeof(T)* m * n);
58  T *row = (T *)(vec + m);
59  for (unsigned int i = 0; i < m; i++, row += n) vec[i] = row;
60  return vec;
61  }
62
65  template<typename T>
66  void copy_matrix(T** dest, T** src, unsigned int m, unsigned int n = 0)
67  {
68  if (n == 0) n = m;
69  for (unsigned int i = 0; i < m; i++)
70  {
71  memcpy(dest[i], src[i], n*sizeof(T));
72  }
73  }
74
81  template<typename T>
82  void save_matrix_octave(const std::string& matrix_name, T** matrix, unsigned int m, unsigned int n = 0, const std::string& filename = std::string())
83  {
84  if (n == 0) n = m;
85
86  //create filename
87  std::string fname = filename;
88  if (fname.empty())
89  fname = matrix_name + ".mat";
90
91  //open file
92  std::ofstream fout(fname.c_str());
93  if (!fout.is_open())
94  {
95  throw Hermes::Exceptions::Exception("Unable to save a matrix to a file \"%s\"", fname.c_str());
96  return;
97  }
98
100  fout << std::string("# name: ") << matrix_name << std::endl;
101  fout << std::string("# type: matrix") << std::endl;
102  fout << std::string("# rows: ") << m << std::endl;
103  fout << std::string("# columns: ") << n << std::endl;
104
105  //write contents
106  for (unsigned int i = 0; i < m; i++)
107  {
108  for (unsigned int k = 0; k < n; k++)
109  fout << ' ' << matrix[i][k];
110  fout << std::endl;
111  }
112
113  //finish
114  fout.close();
115  }
116
118  template<typename T>
119  void save_sparse_matrix_octave(const std::string& matrix_name, const T* Ax, const int* Ap, const int* Ai,
120  unsigned int m, const std::string& filename = std::string())
121  {
122  // create filename
123  std::string fname = filename;
124  if (fname.empty())
125  fname = matrix_name + ".mat";
126
127  // open file
128  std::ofstream fout(fname.c_str());
129  if (!fout.is_open())
130  {
131  throw Hermes::Exceptions::Exception("Unable to save a matrix to a file \"%s\"", fname.c_str());
132  return;
133  }
134
136  fout << std::string("# name: ") << matrix_name << std::endl;
137  fout << std::string("# type: sparse matrix") << std::endl;
138  fout << std::string("# nnz: ") << Ap[m] << std::endl;
139  fout << std::string("# rows: ") << m << std::endl;
140  fout << std::string("# columns: ") << m << std::endl;
141
142  // write contents
143  for (int j = 0; j < m; j++)
144  for (int i = Ap[j]; i < Ap[j + 1]; i++)
145  fout << j + 1 << " " << Ai[i] + 1 << " " << Ax[i] << std::endl;
146
147  // finish
148  fout.close();
149  }
150
153  template<typename T>
154  void transpose(T **matrix, unsigned int m, unsigned int n)
155  {
156  unsigned int min = std::min(m, n);
157  for (unsigned int i = 0; i < min; i++)
158  for (unsigned int j = i + 1; j < min; j++)
159  std::swap(matrix[i][j], matrix[j][i]);
160
161  if (m < n)
162  for (unsigned int i = 0; i < m; i++)
163  for (unsigned int j = m; j < n; j++)
164  matrix[j][i] = matrix[i][j];
165  else if (n < m)
166  for (unsigned int i = n; i < m; i++)
167  for (unsigned int j = 0; j < n; j++)
168  matrix[j][i] = matrix[i][j];
169  }
170
173  template<typename T>
174  void transpose(T *matrix, unsigned int m, unsigned int n, unsigned int size)
175  {
176  unsigned int min = std::min(m, n);
177  for (unsigned int i = 0; i < min; i++)
178  {
179  for (unsigned int j = i + 1; j < min; j++)
180  {
181  std::swap(matrix[i * size + j], matrix[j * size + i]);
182  }
183  }
184
185  if (m < n)
186  for (unsigned int i = 0; i < m; i++)
187  for (unsigned int j = m; j < n; j++)
188  matrix[j * size + i] = matrix[i * size + j];
189
190  else if (n < m)
191  for (unsigned int i = n; i < m; i++)
192  for (unsigned int j = 0; j < n; j++)
193  matrix[j * size + i] = matrix[i * size + j];
194  }
195
197  template<typename T>
198  void change_sign(T **matrix, unsigned int m, unsigned int n)
199  {
200  for (unsigned int i = 0; i < m; i++)
201  for (unsigned int j = 0; j < n; j++)
202  matrix[i][j] = -matrix[i][j];
203  }
205  template<typename T>
206  void change_sign(T *matrix, unsigned int m, unsigned int n, unsigned int size)
207  {
208  for (unsigned int i = 0; i < m; i++)
209  for (unsigned int j = 0; j < n; j++)
210  {
211  int local_matrix_index_array = i * size + j;
212  matrix[local_matrix_index_array] = -matrix[local_matrix_index_array];
213  }
214  }
215
222  template<typename T, typename Int>
223  HERMES_API void ludcmp(T **a, Int n, Int *indx, double *d);
224
232  template<typename T, typename S, typename Int>
233  HERMES_API void lubksb(T **a, Int n, Int *indx, S *b);
234
239  template<typename T, typename Int>
240  HERMES_API void choldc(T **a, Int n, T p[]);
241
249  template<typename T, typename Int>
250  void cholsl(double **a, Int n, double p[], T b[], T x[])
251  {
252  int i, k;
253  T sum;
254
255  for (i = 0; i < n; i++)
256  {
257  sum = b[i];
258  k = i;
259  while (--k >= 0) sum -= a[i][k] * x[k];
260  x[i] = sum / p[i];
261  }
262
263  for (i = n - 1; i >= 0; i--)
264  {
265  sum = x[i];
266  k = i;
267  while (++k < n) sum -= a[k][i] * x[k];
268  x[i] = sum / p[i];
269  }
270  }
271  }
272  }
273 }
274 #endif
HERMES_API void ludcmp(T **a, Int n, Int *indx, double *d)
void save_sparse_matrix_octave(const std::string &matrix_name, const T *Ax, const int *Ap, const int *Ai, unsigned int m, const std::string &filename=std::string())
Saves MxM sparse matrix to a octave file format.
General namespace for the Hermes library.
Exception interface Basically a std::exception, but with a constructor with string and with print_msg...
Definition: exceptions.h:49
void transpose(T **matrix, unsigned int m, unsigned int n)
T ** new_matrix(unsigned int m, unsigned int n=0)
File containing common definitions, and basic global enums etc. for HermesCommon. ...
void cholsl(double **a, Int n, double p[], T b[], T x[])
void copy_matrix(T **dest, T **src, unsigned int m, unsigned int n=0)
File containing definition of exceptions classes.
File containing platform compatibility layer, especially for Win / MSVC.
void change_sign(T **matrix, unsigned int m, unsigned int n)
Changes the sign of a matrix.
HERMES_API void choldc(T **a, Int n, T p[])
File containing common definitions, and basic global enums etc. for HermesCommon. ...
HERMES_API void lubksb(T **a, Int n, Int *indx, S *b)
void save_matrix_octave(const std::string &matrix_name, T **matrix, unsigned int m, unsigned int n=0, const std::string &filename=std::string())
Saves a dense matrix to a octave file format.