HermesCommon  3.0
cs_matrix.cpp
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 #include "cs_matrix.h"
23 #include "util/memory_handling.h"
24
25 namespace Hermes
26 {
27  namespace Algebra
28  {
29  double inline real(double x)
30  {
31  return x;
32  }
33
34  double inline imag(double x)
35  {
36  return 0;
37  }
38
39  double inline real(std::complex<double> x)
40  {
41  return x.real();
42  }
43
44  double inline imag(std::complex<double> x)
45  {
46  return x.imag();
47  }
48
49  template<typename Scalar>
50  int CSMatrix<Scalar>::find_position(int *Ai, int Alen, unsigned int idx)
51  {
52  register int lo = 0, hi = Alen - 1, mid;
53
54  while (true)
55  {
56  mid = (lo + hi) >> 1;
57
58  if (idx < Ai[mid])
59  hi = mid - 1;
60  else if (idx > Ai[mid])
61  lo = mid + 1;
62  else break;
63
65  // value to this position, return 0 when obtaining value there).
66  if (lo > hi)
67  {
68  mid = -1;
69  break;
70  }
71  }
72  return mid;
73  }
74
75  template<typename Scalar>
76  CSMatrix<Scalar>::CSMatrix() : SparseMatrix<Scalar>(), nnz(0), Ap(nullptr), Ai(nullptr), Ax(nullptr)
77  {
78  }
79
80  template<typename Scalar>
81  CSMatrix<Scalar>::CSMatrix(unsigned int size)
82  {
83  this->size = size;
84  this->alloc();
85  }
86
87  template<typename Scalar>
89  {
90  free();
91  }
92
93  template<typename Scalar>
95  {
96  for (unsigned int i = 0; i < this->nnz; i++)
97  Ax[i] *= value;
98  }
99
100  template<typename Scalar>
102  {
103  // initialize the arrays Ap and Ai
104  Ap = malloc_with_check<CSMatrix<Scalar>, int>(this->size + 1, this);
105  int aisize = this->get_num_indices();
106  Ai = malloc_with_check<CSMatrix<Scalar>, int>(aisize, this);
107
108  // sort the indices and remove duplicities, insert into Ai
109  unsigned int i;
110  int pos = 0;
111  for (i = 0; i < this->size; i++)
112  {
113  Ap[i] = pos;
114  pos += this->sort_and_store_indices(&this->pages[i], Ai + pos, Ai + aisize);
115  }
116  Ap[this->size] = pos;
117
118  free_with_check(this->pages);
119  free_with_check(this->next_pages);
120
121  nnz = Ap[this->size];
122
123  this->alloc_data();
124  }
125
126  template<typename Scalar>
128  {
129  Ax = calloc_with_check<CSMatrix<Scalar>, Scalar>(nnz, this);
130  }
131
132  template<typename Scalar>
134  {
135  nnz = 0;
136  free_with_check(Ap);
137  free_with_check(Ai);
138  free_with_check(Ax);
139  }
140
141  template<typename Scalar>
142  void CSMatrix<Scalar>::set_row_zero(unsigned int n)
143  {
144  for (int i = 0; i < Ap[n + 1] - Ap[n]; i++)
145  Ax[Ap[n] + i] = Scalar(0);
146  }
147
148  template<typename Scalar>
150  {
151  memset(Ax, 0, sizeof(Scalar)* nnz);
152  }
153
154  template<typename Scalar>
155  unsigned int CSMatrix<Scalar>::get_nnz() const
156  {
157  return this->nnz;
158  }
159
160  template<typename Scalar>
162  {
163  return nnz / (double)(this->size * this->size);
164  }
165
166  template<typename Scalar>
167  void CSMatrix<Scalar>::create(unsigned int size, unsigned int nnz, int* ap, int* ai, Scalar* ax)
168  {
169  this->nnz = nnz;
170  this->size = size;
171  this->Ap = malloc_with_check<CSMatrix<Scalar>, int>(this->size + 1, this);
172  this->Ai = malloc_with_check<CSMatrix<Scalar>, int>(nnz, this);
173  this->Ax = malloc_with_check<CSMatrix<Scalar>, Scalar>(nnz, this);
174  memcpy(this->Ap, ap, (this->size + 1) * sizeof(int));
175  memcpy(this->Ai, ai, this->nnz * sizeof(int));
176  memcpy(this->Ax, ax, this->nnz * sizeof(Scalar));
177  }
178
179  template<typename Scalar>
181  {
182  // The variable names are so to reflect CSC -> CSR direction.
183  // From the "Ap indexed by columns" to "Ap indexed by rows".
184  int* tempAp = malloc_with_check<CSMatrix<Scalar>, int>(this->size + 1, this);
185  int* tempAi = malloc_with_check<CSMatrix<Scalar>, int>(nnz, this);
186  Scalar* tempAx = malloc_with_check<CSMatrix<Scalar>, Scalar>(nnz, this);
187
188  int run_i = 0;
189  for (int target_row = 0; target_row < this->size; target_row++)
190  {
191  tempAp[target_row] = run_i;
192  for (int src_column = 0; src_column < this->size; src_column++)
193  {
194  for (int src_row = this->Ap[src_column]; src_row < this->Ap[src_column + 1]; src_row++)
195  {
196  if (this->Ai[src_row] == target_row)
197  {
198  tempAi[run_i] = src_column;
199  tempAx[run_i++] = this->Ax[src_row];
200  }
201  }
202  }
203  }
204
205  tempAp[this->size] = this->nnz;
206  memcpy(this->Ai, tempAi, sizeof(int)* nnz);
207  memcpy(this->Ap, tempAp, sizeof(int)* (this->size + 1));
208  memcpy(this->Ax, tempAx, sizeof(Scalar)* nnz);
209  free_with_check(tempAi);
210  free_with_check(tempAx);
211  free_with_check(tempAp);
212  }
213
214  template<typename Scalar>
216  {
217  return this->Ap;
218  }
219
220  template<typename Scalar>
222  {
223  return this->Ai;
224  }
225
226  template<typename Scalar>
227  Scalar *CSMatrix<Scalar>::get_Ax() const
228  {
229  return this->Ax;
230  }
231
232  template<typename Scalar>
233  void CSMatrix<Scalar>::add_as_block(unsigned int offset_i, unsigned int offset_j, SparseMatrix<Scalar>* mat)
234  {
235  if ((this->get_size() < offset_i + mat->get_size()) || (this->get_size() < offset_j + mat->get_size()))
236  throw Hermes::Exceptions::Exception("Incompatible matrix sizes in SparseMatrix<Scalar>::add_as_block()");
237
238  CSMatrix<Scalar>* csMatrix = dynamic_cast<CSMatrix<Scalar>*>(mat);
239  if (!mat)
240  {
242  }
243  else
244  {
245  for (unsigned short i = 0; i < csMatrix->get_size(); i++)
246  {
247  int index = csMatrix->Ap[i];
248  for (int j = 0; j < csMatrix->Ap[i + 1] - index; j++)
249  {
250  this->add(offset_i + csMatrix->Ai[index + j], offset_j + i, csMatrix->Ax[index + j]);
251  }
252  }
253  }
254  }
255
256  template<>
257  void CSMatrix<double>::add(unsigned int m, unsigned int n, double v)
258  {
259  if (v != 0.0) // ignore zero values.
260  {
261  // Find m-th row in the n-th column.
262  int pos = find_position(Ai + Ap[n], Ap[n + 1] - Ap[n], m);
263  // Make sure we are adding to an existing non-zero entry.
264  if (pos < 0)
265  {
266  this->info("CSMatrix<Scalar>::add(): i = %d, j = %d.", m, n);
268  }
269
270 #pragma omp atomic
271  Ax[Ap[n] + pos] += v;
272  }
273  }
274
275  template<>
276  void CSMatrix<std::complex<double> >::add(unsigned int m, unsigned int n, std::complex<double> v)
277  {
278  if (v != 0.0) // ignore zero values.
279  {
280  // Find m-th row in the n-th column.
281  int pos = find_position(Ai + Ap[n], Ap[n + 1] - Ap[n], m);
282  // Make sure we are adding to an existing non-zero entry.
283  if (pos < 0)
284  {
285  this->info("CSMatrix<Scalar>::add(): i = %d, j = %d.", m, n);
287  }
288
290  Ax[Ap[n] + pos] += v;
291  }
292  }
293
294  template<typename Scalar>
295  Scalar CSMatrix<Scalar>::get(unsigned int m, unsigned int n) const
296  {
297  // Find m-th row in the n-th column.
298  int mid = find_position(Ai + Ap[n], Ap[n + 1] - Ap[n], m);
299
300  if (mid < 0) // if the entry has not been found
301  return 0.0;
302  else
303  return Ax[Ap[n] + mid];
304  }
305
306  template<typename Scalar>
308  {
309  }
310
311  template<typename Scalar>
312  CSCMatrix<Scalar>::CSCMatrix(unsigned int size) : CSMatrix<Scalar>(size)
313  {
314  }
315
316  template<typename Scalar>
318  {
319  }
320
321  template<>
322  void CSCMatrix<double>::add(unsigned int m, unsigned int n, double v)
323  {
325  }
326
327  template<>
328  void CSCMatrix<std::complex<double> >::add(unsigned int m, unsigned int n, std::complex<double> v)
329  {
331  }
332
333  template<typename Scalar>
334  Scalar CSCMatrix<Scalar>::get(unsigned int m, unsigned int n) const
335  {
336  return CSMatrix<Scalar>::get(m, n);
337  }
338
339  template<typename Scalar>
340  void CSCMatrix<Scalar>::multiply_with_vector(Scalar* vector_in, Scalar*& vector_out, bool vector_out_initialized) const
341  {
342  if (!vector_out_initialized)
343  vector_out = malloc_with_check<Scalar>(this->size);
344  memset(vector_out, 0, sizeof(Scalar)* this->size);
345  {
346  for (int i = 0; i < this->size; i++)
347  {
348  for (int j = 0; j < this->Ap[i + 1] - this->Ap[i]; j++)
349  vector_out[this->Ai[this->Ap[i] + j]] += this->Ax[this->Ap[i] + j] * vector_in[i];
350  }
351  }
352  }
353
354  static int i_coordinate(int i, int j, bool invert)
355  {
356  if (invert)
357  return i;
358  else
359  return j;
360  }
361
362  static int j_coordinate(int i, int j, bool invert)
363  {
364  if (invert)
365  return j;
366  else
367  return i;
368  }
369
370  template<typename Scalar>
371  void CSMatrix<Scalar>::export_to_file(const char *filename, const char *var_name, MatrixExportFormat fmt, char* number_format, bool invert_storage)
372  {
373  switch (fmt)
374  {
376  {
377  FILE* file = fopen(filename, "w");
378  if (!file)
379  throw Exceptions::IOException(Exceptions::IOException::Write, filename);
381  fprintf(file, "%%%%MatrixMarket matrix coordinate real general\n");
382  else
383  fprintf(file, "%%%%MatrixMarket matrix coordinate complex general\n");
384
385  fprintf(file, "%d %d %d\n", this->size, this->size, this->nnz);
386
387  if (invert_storage)
388  this->switch_orientation();
389  for (unsigned int j = 0; j < this->size; j++)
390  {
391  for (int i = Ap[j]; i < Ap[j + 1]; i++)
392  {
393  Hermes::Helpers::fprint_coordinate_num(file, i_coordinate(Ai[i] + 1, j + 1, invert_storage), j_coordinate(Ai[i] + 1, j + 1, invert_storage), Ax[i], number_format);
394  fprintf(file, "\n");
395  }
396  }
397  if (invert_storage)
398  this->switch_orientation();
399
400  fclose(file);
401  }
402  break;
403
405  {
406 #ifdef WITH_MATIO
407  mat_sparse_t sparse;
408  sparse.nzmax = this->nnz;
409  if (invert_storage)
410  this->switch_orientation();
411
412  sparse.nir = this->nnz;
413  sparse.ir = this->Ai;
414  sparse.njc = this->size + 1;
415  sparse.jc = (int *)this->Ap;
416  sparse.ndata = this->nnz;
417
418  size_t dims[2];
419  dims[0] = this->size;
420  dims[1] = this->size;
421
422  mat_t *mat = Mat_CreateVer(filename, "", MAT_FT_MAT5);
423
424  matvar_t *matvar;
425
426  // For complex. No allocation here.
427  double* Ax_re = nullptr;
428  double* Ax_im = nullptr;
429
430  // For real.
432  {
433  sparse.data = Ax;
434  matvar = Mat_VarCreate(var_name, MAT_C_SPARSE, MAT_T_DOUBLE, 2, dims, &sparse, MAT_F_DONT_COPY_DATA);
435  }
436  else
437  {
438  // For complex.
439  Ax_re = malloc_with_check<CSMatrix<Scalar>, double>(this->nnz, this);
440  Ax_im = malloc_with_check<CSMatrix<Scalar>, double>(this->nnz, this);
441  struct mat_complex_split_t z = { Ax_re, Ax_im };
442
443  for (int i = 0; i < this->nnz; i++)
444  {
445  Ax_re[i] = ((std::complex<double>)(this->Ax[i])).real();
446  Ax_im[i] = ((std::complex<double>)(this->Ax[i])).imag();
447  sparse.data = &z;
448  }
449  matvar = Mat_VarCreate(var_name, MAT_C_SPARSE, MAT_T_DOUBLE, 2, dims, &sparse, MAT_F_DONT_COPY_DATA | MAT_F_COMPLEX);
450  }
451
452  if (matvar)
453  {
454  Mat_VarWrite(mat, matvar, MAT_COMPRESSION_ZLIB);
455  Mat_VarFree(matvar);
456  }
457  if (invert_storage)
458  this->switch_orientation();
459  free_with_check(Ax_re);
460  free_with_check(Ax_im);
461  Mat_Close(mat);
462
463  if (!matvar)
464  throw Exceptions::IOException(Exceptions::IOException::Write, filename);
465 #endif
466  }
467  break;
468
470  {
471  FILE* file = fopen(filename, "w");
472  if (!file)
473  throw Exceptions::IOException(Exceptions::IOException::Write, filename);
474
475  if (invert_storage)
476  this->switch_orientation();
477  for (unsigned int j = 0; j < this->size; j++)
478  {
479  for (int i = Ap[j]; i < Ap[j + 1]; i++)
480  {
481  Helpers::fprint_coordinate_num(file, i_coordinate(Ai[i], j, invert_storage), j_coordinate(Ai[i], j, invert_storage), Ax[i], number_format);
482  fprintf(file, "\n");
483  }
484  }
485  if (invert_storage)
486  this->switch_orientation();
487
488  fclose(file);
489  }
490  break;
491
492  case EXPORT_FORMAT_MATLAB_SIMPLE:
493  {
494  FILE* file = fopen(filename, "w");
495  if (invert_storage)
496  this->switch_orientation();
497  fprintf(file, "%% Size: %dx%d\n%% Nonzeros: %d\ntemp = zeros(%d, 3);\ntemp =[\n",
498  this->size, this->size, nnz, nnz);
499  for (unsigned int j = 0; j < this->size; j++)
500  for (int i = Ap[j]; i < Ap[j + 1]; i++)
501  {
502  fprintf(file, "%d %d ", Ai[i] + 1, j + 1);
503  Hermes::Helpers::fprint_num(file, Ax[i], number_format);
504  fprintf(file, "\n");
505  }
506  fprintf(file, "];\n");
507  if (invert_storage)
508  this->switch_orientation();
509  fclose(file);
510  }
511  break;
512 #ifdef WITH_BSON
513  case EXPORT_FORMAT_BSON:
514  {
515  // Init bson
516  bson bw;
517  bson_init(&bw);
518
519  // Matrix size.
520  bson_append_int(&bw, "size", this->size);
521  // Nonzeros.
522  bson_append_int(&bw, "nnz", this->nnz);
523
524  if (invert_storage)
525  this->switch_orientation();
526
527  bson_append_start_array(&bw, "Ap");
528  for (unsigned int i = 0; i < this->size; i++)
529  bson_append_int(&bw, "p", this->Ap[i]);
530  bson_append_finish_array(&bw);
531
532  bson_append_start_array(&bw, "Ai");
533  for (unsigned int i = 0; i < this->nnz; i++)
534  bson_append_int(&bw, "i", this->Ai[i]);
535  bson_append_finish_array(&bw);
536
537  bson_append_start_array(&bw, "Ax");
538  for (unsigned int i = 0; i < this->nnz; i++)
539  bson_append_double(&bw, "x", real(this->Ax[i]));
540  bson_append_finish_array(&bw);
541
543  {
544  bson_append_start_array(&bw, "Ax-imag");
545  for (unsigned int i = 0; i < this->nnz; i++)
546  bson_append_double(&bw, "x-i", imag(this->Ax[i]));
547  bson_append_finish_array(&bw);
548  }
549  bson_append_finish_array(&bw);
550
551  if (invert_storage)
552  this->switch_orientation();
553
554  // Done.
555  bson_finish(&bw);
556
557  // Write to disk.
558  FILE *fpw;
559  fpw = fopen(filename, "wb");
560  const char *dataw = (const char *)bson_data(&bw);
561  fwrite(dataw, bson_size(&bw), 1, fpw);
562  fclose(fpw);
563
564  bson_destroy(&bw);
565  }
566  break;
567 #endif
568  }
569  }
570
571  template<typename Scalar>
572  void CSCMatrix<Scalar>::export_to_file(const char *filename, const char *var_name, MatrixExportFormat fmt, char* number_format)
573  {
574  CSMatrix<Scalar>::export_to_file(filename, var_name, fmt, number_format, false);
575  }
576
577  template<typename Scalar>
578  void CSRMatrix<Scalar>::export_to_file(const char *filename, const char *var_name, MatrixExportFormat fmt, char* number_format)
579  {
580  CSMatrix<Scalar>::export_to_file(filename, var_name, fmt, number_format, true);
581  }
582
583  template<typename Scalar>
584  void CSMatrix<Scalar>::import_from_file(const char *filename, const char *var_name, MatrixExportFormat fmt, bool invert_storage)
585  {
586  this->free();
587
588  switch (fmt)
589  {
591  throw Exceptions::MethodNotOverridenException("CSMatrix<Scalar>::import_from_file - Matrix Market");
592  break;
594  {
595 #ifdef WITH_MATIO
596  mat_t *matfp;
597  matvar_t *matvar;
598
599  matfp = Mat_Open(filename, MAT_ACC_RDONLY);
600
601  if (!matfp)
603
605
606  if (matvar)
607  {
608  mat_sparse_t *sparse = (mat_sparse_t *)matvar->data;
609
610  this->nnz = sparse->nir;
611  this->Ax = malloc_with_check<CSMatrix<Scalar>, Scalar>(this->nnz, this);
612  this->Ai = malloc_with_check<CSMatrix<Scalar>, int>(this->nnz, this);
613  this->size = sparse->njc;
614  this->Ap = malloc_with_check<CSMatrix<Scalar>, int>(this->size + 1, this);
615
616  void* data = nullptr;
618  data = sparse->data;
619  else
620  {
621  std::complex<double>* complex_data = malloc_with_check<CSMatrix<Scalar>, std::complex<double> >(this->nnz, this);
622  double* real_array = (double*)((mat_complex_split_t*)sparse->data)->Re;
623  double* imag_array = (double*)((mat_complex_split_t*)sparse->data)->Im;
624  for (int i = 0; i < this->nnz; i++)
625  complex_data[i] = std::complex<double>(real_array[i], imag_array[i]);
626  data = (void*)complex_data;
627  }
628  memcpy(this->Ax, data, this->nnz * sizeof(Scalar));
630  free_with_check(data);
631  memcpy(this->Ap, sparse->jc, this->size * sizeof(int));
632  this->Ap[this->size] = this->nnz;
633  memcpy(this->Ai, sparse->ir, this->nnz * sizeof(int));
634
635  if (invert_storage)
636  this->switch_orientation();
637  }
638
639  Mat_Close(matfp);
640
641  if (!matvar)
643 #else
644  throw Exceptions::Exception("MATIO not included.");
645 #endif
646  }
647  break;
648
650  throw Exceptions::MethodNotOverridenException("CSMatrix<Scalar>::import_from_file - Simple format");
651
652 #ifdef WITH_BSON
653  case EXPORT_FORMAT_BSON:
654  {
655  FILE *fpr;
656  fpr = fopen(filename, "rb");
657
658  // file size:
659  fseek(fpr, 0, SEEK_END);
660  int size = ftell(fpr);
661  rewind(fpr);
662
663  // allocate memory to contain the whole file:
664  char *datar = malloc_with_check<char>(size);
666  fclose(fpr);
667
668  bson br;
669  bson_init_finished_data(&br, datar, 0);
670
671  bson_iterator it;
672  bson sub;
673  bson_find(&it, &br, "size");
674  this->size = bson_iterator_int(&it);
675  bson_find(&it, &br, "nnz");
676  this->nnz = bson_iterator_int(&it);
677
678  this->Ap = malloc_with_check<CSMatrix<Scalar>, int>(this->size + 1, this);
679  this->Ai = malloc_with_check<CSMatrix<Scalar>, int>(nnz, this);
680  this->Ax = malloc_with_check<CSMatrix<Scalar>, Scalar>(nnz, this);
681
682  // coeffs
683  bson_iterator it_coeffs;
684  bson_find(&it_coeffs, &br, "Ap");
685  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
686  bson_iterator_init(&it, &sub);
687  int index_coeff = 0;
688  while (bson_iterator_next(&it))
689  this->Ap[index_coeff++] = bson_iterator_int(&it);
690  this->Ap[this->size] = this->nnz;
691
692  bson_find(&it_coeffs, &br, "Ai");
693  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
694  bson_iterator_init(&it, &sub);
695  index_coeff = 0;
696  while (bson_iterator_next(&it))
697  this->Ai[index_coeff++] = bson_iterator_int(&it);
698
699  bson_find(&it_coeffs, &br, "Ax");
700  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
701  bson_iterator_init(&it, &sub);
702  index_coeff = 0;
703  while (bson_iterator_next(&it))
704  this->Ax[index_coeff++] = bson_iterator_double(&it);
705
707  {
708  bson_find(&it_coeffs, &br, "Ax-imag");
709  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
710  bson_iterator_init(&it, &sub);
711  index_coeff = 0;
712  while (bson_iterator_next(&it))
713  ((std::complex<double>)this->Ax[index_coeff++]).imag(bson_iterator_double(&it));
714  }
715
716  if (invert_storage)
717  this->switch_orientation();
718
719  bson_destroy(&br);
720  free_with_check(datar);
721  }
722 #endif
723  break;
724  }
725  }
726
727  template<typename Scalar>
728  void CSCMatrix<Scalar>::import_from_file(const char *filename, const char *var_name, MatrixExportFormat fmt)
729  {
730  CSMatrix<Scalar>::import_from_file(filename, var_name, fmt, false);
731  }
732
733  template<typename Scalar>
734  void CSRMatrix<Scalar>::import_from_file(const char *filename, const char *var_name, MatrixExportFormat fmt)
735  {
736  CSMatrix<Scalar>::import_from_file(filename, var_name, fmt, true);
737  }
738
739  template<typename Scalar>
741  {
743  new_matrix->create(this->get_size(), this->get_nnz(), this->get_Ap(), this->get_Ai(), this->get_Ax());
744  return new_matrix;
745  }
746
747  template<typename Scalar>
749  {
750  }
751
752  template<typename Scalar>
753  CSRMatrix<Scalar>::CSRMatrix(unsigned int size) : CSMatrix<Scalar>(size)
754  {
755  }
756
757  template<typename Scalar>
759  {
760  }
761
762  template<>
763  void CSRMatrix<double>::add(unsigned int m, unsigned int n, double v)
764  {
766  }
767
768  template<>
769  void CSRMatrix<std::complex<double> >::add(unsigned int m, unsigned int n, std::complex<double> v)
770  {
772  }
773
774  template<typename Scalar>
775  Scalar CSRMatrix<Scalar>::get(unsigned int m, unsigned int n) const
776  {
777  return CSMatrix<Scalar>::get(n, m);
778  }
779
780  template<typename Scalar>
781  void CSRMatrix<Scalar>::pre_add_ij(unsigned int row, unsigned int col)
782  {
783  if (this->pages[row].count >= SparseMatrix<Scalar>::PAGE_SIZE)
784  {
785  typename SparseMatrix<Scalar>::Page* final_page = &(this->pages[row]);
786  while (final_page->next != nullptr && final_page->count >= SparseMatrix<Scalar>::PAGE_SIZE)
787  final_page = final_page->next;
788
789  if (final_page->next == nullptr && final_page->count >= SparseMatrix<Scalar>::PAGE_SIZE)
790  {
791  final_page->next = new typename SparseMatrix<Scalar>::Page(true);
792  final_page = final_page->next;
793  }
794  final_page->idx[final_page->count++] = col;
795  }
796  else
797  this->pages[row].idx[this->pages[row].count++] = col;
798  }
799
800  template<typename Scalar>
802  {
804  new_matrix->create(this->get_size(), this->get_nnz(), this->get_Ap(), this->get_Ai(), this->get_Ax());
805  return new_matrix;
806  }
807  }
808 }
809
810 template class HERMES_API Hermes::Algebra::CSMatrix < double > ;
811 template class HERMES_API Hermes::Algebra::CSMatrix < std::complex<double> > ;
812
813 template class HERMES_API Hermes::Algebra::CSCMatrix < double > ;
814 template class HERMES_API Hermes::Algebra::CSCMatrix < std::complex<double> > ;
815
816 template class HERMES_API Hermes::Algebra::CSRMatrix < double > ;
817 template class HERMES_API Hermes::Algebra::CSRMatrix < std::complex<double> > ;
virtual unsigned int get_size() const
Definition: matrix.cpp:82
void switch_orientation()
Switches CSR / CSC arrays.
Definition: cs_matrix.cpp:180
General namespace for the Hermes library.
void import_from_file(const char *filename, const char *var_name, MatrixExportFormat fmt, bool invert_storage=false)
Definition: cs_matrix.cpp:584
Exception interface Basically a std::exception, but with a constructor with string and with print_msg...
Definition: exceptions.h:49
General CSR Matrix class. (can be used in umfpack, in that case use the CSCMatrix subclass...
Definition: cs_matrix.h:161
virtual void add(unsigned int m, unsigned int n, Scalar v)
int * Ap
Index to Ax/Ai, where each column / row starts.
Definition: cs_matrix.h:121
T ** new_matrix(unsigned int m, unsigned int n=0)
File containing common definitions, and basic global enums etc. for HermesCommon. ...
static int find_position(int *Ai, int Alen, unsigned int idx)
Finds the correct position to insert / retrieve elements.
Definition: cs_matrix.cpp:50
void multiply_with_vector(Scalar *vector_in, Scalar *&vector_out, bool vector_out_initialized) const
Multiply with a vector.
Definition: cs_matrix.cpp:340
General (abstract) sparse matrix representation in Hermes.
SparseMatrix< Scalar > * duplicate() const
Duplicates a matrix (including allocation).
Definition: cs_matrix.cpp:801
MatrixExportFormat
Format of file matrix and vector output.
virtual void add_as_block(unsigned int i, unsigned int j, SparseMatrix< Scalar > *mat)
Definition: cs_matrix.cpp:233
Basic cs (Compressed sparse) matrix classes and operations.
Plain ascii file lines contains row column and value.
virtual void zero()
Utility method.
Definition: cs_matrix.cpp:149
virtual void add(unsigned int m, unsigned int n, Scalar v)
virtual void export_to_file(const char *filename, const char *var_name, MatrixExportFormat fmt, char *number_format="%lf")
Definition: cs_matrix.cpp:572
General CSC Matrix class. (can be used in umfpack, in that case use the CSCMatrix subclass...
Definition: cs_matrix.h:131
IO exception. Internal. Exception occurs when something fails to be written to / read from the disk...
Definition: exceptions.h:76
void create(unsigned int size, unsigned int nnz, int *ap, int *ai, Scalar *ax)
Definition: cs_matrix.cpp:167
int idx[PAGE_SIZE]
buffer for storing indices
Definition: matrix.h:197
virtual void set_row_zero(unsigned int n)
Utility method.
Definition: cs_matrix.cpp:142
virtual Scalar get(unsigned int Ai_data_index, unsigned int Ai_index) const
Definition: cs_matrix.cpp:295
void export_to_file(const char *filename, const char *var_name, MatrixExportFormat fmt, char *number_format="%lf")
Definition: cs_matrix.cpp:578
virtual Scalar get(unsigned int m, unsigned int n) const
Definition: cs_matrix.cpp:334
void import_from_file(const char *filename, const char *var_name, MatrixExportFormat fmt)
Definition: cs_matrix.cpp:734
int * Ai
Row / Column indices of values in Ax.
Definition: cs_matrix.h:119
Type deductors.
Definition: common.h:276
SparseMatrix< Scalar > * duplicate() const
Duplicates a matrix (including allocation).
Definition: cs_matrix.cpp:740
unsigned char count
number of indices stored
Definition: matrix.h:193
Scalar * get_Ax() const
Definition: cs_matrix.cpp:227
void multiply_with_Scalar(Scalar value)
Multiplies matrix with a Scalar.
Definition: cs_matrix.cpp:94
Structure for storing indices in sparse matrix.
Definition: matrix.h:191
virtual double get_fill_in() const
Utility method.
Definition: cs_matrix.cpp:161
virtual Scalar get(unsigned int m, unsigned int n) const
Definition: cs_matrix.cpp:775
Method is not overriden and should be.
Definition: exceptions.h:201
virtual void import_from_file(const char *filename, const char *var_name, MatrixExportFormat fmt)
Definition: cs_matrix.cpp:728
General CS Matrix class. Either row- or column- specific (see subclassses).
Definition: cs_matrix.h:35
CSCMatrix()
Default constructor.
Definition: cs_matrix.cpp:307
void export_to_file(const char *filename, const char *var_name, MatrixExportFormat fmt, char *number_format="%lf", bool invert_storage=false)
Definition: cs_matrix.cpp:371
virtual void add_as_block(unsigned int i, unsigned int j, SparseMatrix< Scalar > *mat)
Definition: matrix.cpp:207
virtual void pre_add_ij(unsigned int row, unsigned int col)
Definition: cs_matrix.cpp:781
Page * next
pointer to next page
Definition: matrix.h:199
Matrix Market which can be read by pysparse library.
CSMatrix()
Default constructor.
Definition: cs_matrix.cpp:76
virtual void free()
Utility method.
Definition: cs_matrix.cpp:133
virtual unsigned int get_nnz() const
Utility method.
Definition: cs_matrix.cpp:155
virtual void alloc()
Allocate utility storage (row, column indices, etc.).
Definition: cs_matrix.cpp:101
CSRMatrix()
Default constructor.
Definition: cs_matrix.cpp:748
virtual void add(unsigned int Ai_data_index, unsigned int Ai_index, Scalar v)