22 #ifndef __HERMES_COMMON_MATRIX_H
23 #define __HERMES_COMMON_MATRIX_H
46 namespace DenseMatrixOperations
55 T **vec = (T **)
new char[
sizeof(T *) * m +
sizeof(T) * m * n];
56 memset(vec, 0,
sizeof(T *) * m +
sizeof(T) * m * n);
57 T *row = (T *) (vec + m);
58 for (
unsigned int i = 0; i < m; i++, row += n) vec[i] = row;
63 T **new_matrix_malloc(
unsigned int m,
unsigned int n = 0)
66 T **vec = (T **) malloc(
sizeof(T *) * m +
sizeof(T) * m * n);
67 memset(vec, 0,
sizeof(T *) * m +
sizeof(T) * m * n);
68 T *row = (T *) (vec + m);
69 for (
unsigned int i = 0; i < m; i++, row += n) vec[i] = row;
76 void copy_matrix(T** dest, T** src,
unsigned int m,
unsigned int n = 0)
79 for(
unsigned int i = 0; i < m; i++)
81 memcpy(dest[i], src[i], n*
sizeof(T));
92 void save_matrix_octave(
const std::string& matrix_name, T** matrix,
unsigned int m,
unsigned int n = 0,
const std::string& filename = std::string())
97 std::string fname = filename;
99 fname = matrix_name +
".mat";
102 std::ofstream fout(fname.c_str());
110 fout << std::string(
"# name: ") << matrix_name << std::endl;
111 fout << std::string(
"# type: matrix") << std::endl;
112 fout << std::string(
"# rows: ") << m << std::endl;
113 fout << std::string(
"# columns: ") << n << std::endl;
116 for(
unsigned int i = 0; i < m; i++)
118 for(
unsigned int k = 0; k < n; k++)
119 fout <<
' ' << matrix[i][k];
130 unsigned int m,
const std::string& filename = std::string())
133 std::string fname = filename;
135 fname = matrix_name +
".mat";
138 std::ofstream fout(fname.c_str());
146 fout << std::string(
"# name: ") << matrix_name << std::endl;
147 fout << std::string(
"# type: sparse matrix") << std::endl;
148 fout << std::string(
"# nnz: ") << Ap[m] << std::endl;
149 fout << std::string(
"# rows: ") << m << std::endl;
150 fout << std::string(
"# columns: ") << m << std::endl;
153 for (
int j = 0; j < m; j++)
154 for (
int i = Ap[j]; i < Ap[j + 1]; i++)
155 fout << j + 1 <<
" " << Ai[i] + 1 <<
" " << Ax[i] << std::endl;
164 void transpose(T **matrix,
unsigned int m,
unsigned int n)
166 unsigned int min = std::min(m, n);
167 for (
unsigned int i = 0; i < min; i++)
168 for (
unsigned int j = i + 1; j < min; j++)
169 std::swap(matrix[i][j], matrix[j][i]);
172 for (
unsigned int i = 0; i < m; i++)
173 for (
unsigned int j = m; j < n; j++)
174 matrix[j][i] = matrix[i][j];
176 for (
unsigned int i = n; i < m; i++)
177 for (
unsigned int j = 0; j < n; j++)
178 matrix[j][i] = matrix[i][j];
183 void chsgn(T **matrix,
unsigned int m,
unsigned int n)
185 for (
unsigned int i = 0; i < m; i++)
186 for (
unsigned int j = 0; j < n; j++)
187 matrix[i][j] = -matrix[i][j];
196 HERMES_API
void ludcmp(
double **a,
int n,
int *indx,
double *d);
206 void lubksb(
double **a,
int n,
int *indx, T *b)
211 for (i = 0; i < n; i++)
216 for (j = 0; j < i; j++) sum -= a[i][j]*b[j];
219 for (i = n-1; i >= 0; i--)
222 for (j = i + 1; j < n; j++) sum -= a[i][j]*b[j];
223 b[i] = sum / a[i][i];
231 HERMES_API
void choldc(
double **a,
int n,
double p[]);
241 void cholsl(
double **a,
int n,
double p[], T b[], T x[])
246 for (i = 0; i < n; i++)
250 while (--k >= 0) sum -= a[i][k] * x[k];
254 for (i = n-1; i >= 0; i--)
258 while (++k < n) sum -= a[k][i] * x[k];
280 template<
typename Scalar>
290 Matrix(
unsigned int size) { this->size = size;};
294 Matrix() { this->size = 0;};
297 virtual void alloc() = 0;
300 virtual void free() = 0;
306 virtual Scalar
get(
unsigned int m,
unsigned int n) = 0;
309 virtual void zero() = 0;
312 virtual void add_to_diagonal(Scalar v) = 0;
319 virtual void add(
unsigned int m,
unsigned int n, Scalar v) = 0;
328 virtual void add(
unsigned int m,
unsigned int n, Scalar **mat,
int *rows,
int *cols) = 0;
339 virtual unsigned int get_matrix_size()
const = 0;
347 template<
typename Scalar>
359 virtual void prealloc(
unsigned int n);
365 virtual void pre_add_ij(
unsigned int row,
unsigned int col);
370 virtual unsigned int get_size() {
return this->size; }
401 virtual void extract_row_copy(
unsigned int row,
unsigned int len,
402 unsigned int &n_entries,
double *vals,
403 unsigned int *idxs) { }
418 virtual void extract_col_copy(
unsigned int col,
unsigned int len,
419 unsigned int &n_entries,
double *vals,
420 unsigned int *idxs) { }
436 virtual double get_fill_in()
const = 0;
438 unsigned row_storage:1;
439 unsigned col_storage:1;
450 static const int PAGE_SIZE = 62;
471 int sort_and_store_indices(
Page *page,
int *buffer,
int *max);
474 int get_num_indices();
481 template<
typename Scalar>
490 virtual void alloc(
unsigned int ndofs) = 0;
492 virtual void free() = 0;
499 virtual Scalar
get(
unsigned int idx) = 0;
503 virtual void extract(Scalar *v)
const = 0;
506 virtual void zero() = 0;
509 virtual void change_sign() = 0;
515 virtual void set(
unsigned int idx, Scalar y) = 0;
521 virtual void add(
unsigned int idx, Scalar y) = 0;
526 virtual void add_vector(Scalar* vec) = 0;
533 virtual void add(
unsigned int n,
unsigned int *idx, Scalar *y) = 0;
536 unsigned int length()
const {
return this->size;}
543 virtual bool dump(FILE *file,
const char *var_name,
553 template<
typename Scalar> HERMES_API
558 template<
typename Scalar> HERMES_API