HermesCommon  2.0
umfpack_solver.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).
4 // Email: hpfem-group@unr.edu, home page: http://hpfem.org/.
5 //
6 // Hermes 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 // Hermes 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 Hermes; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
22 #include "config.h"
23 #ifdef WITH_UMFPACK
24 #include "umfpack_solver.h"
25 
26 extern "C"
27 {
28 #include <umfpack.h>
29 }
30 
31 namespace Hermes
32 {
33  namespace Algebra
34  {
35  static int find_position(int *Ai, int Alen, int idx)
36  {
37  assert(Ai != NULL);
38  assert(Alen > 0);
39  assert(idx >= 0);
40 
41  register int lo = 0, hi = Alen - 1, mid;
42 
43  while (true)
44  {
45  mid = (lo + hi) >> 1;
46 
47  if(idx < Ai[mid]) hi = mid - 1;
48  else if(idx > Ai[mid]) lo = mid + 1;
49  else break;
50 
51  // Sparse matrix entry not found (raise an error when trying to add
52  // value to this position, return 0 when obtaining value there).
53  if(lo > hi)
54  {
55  mid = -1;
56  break;
57  }
58  }
59  return mid;
60  }
61 
62  template<typename Scalar>
63  CSCMatrix<Scalar>::CSCMatrix()
64  {
65  this->size = 0; nnz = 0;
66  Ap = NULL;
67  Ai = NULL;
68  Ax = NULL;
69  }
70 
71  template<typename Scalar>
72  CSCMatrix<Scalar>::CSCMatrix(unsigned int size)
73  {
74  this->size = size;
75  this->alloc();
76  }
77 
78  template<typename Scalar>
79  CSCMatrix<Scalar>::~CSCMatrix()
80  {
81  free();
82  }
83 
84  template<typename Scalar>
85  void CSCMatrix<Scalar>::multiply_with_vector(Scalar* vector_in, Scalar* vector_out)
86  {
87  int n = this->size;
88  for (int j = 0; j<n; j++) vector_out[j] = 0;
89  for (int j = 0; j<n; j++)
90  {
91  for (int i = Ap[j]; i < Ap[j + 1]; i++)
92  {
93  vector_out[Ai[i]] += vector_in[j]*Ax[i];
94  }
95  }
96  }
97 
98  template<typename Scalar>
99  void CSCMatrix<Scalar>::multiply_with_Scalar(Scalar value)
100  {
101  for (unsigned int i = 0; i < this->nnz; i++) Ax[i] *= value;
102  }
103 
104  template<typename Scalar>
105  void CSCMatrix<Scalar>::alloc()
106  {
107  assert(this->pages != NULL);
108 
109  // initialize the arrays Ap and Ai
110  Ap = new int[this->size + 1];
111  int aisize = this->get_num_indices();
112  Ai = new int[aisize];
113 
114  // sort the indices and remove duplicities, insert into Ai
115  unsigned int i;
116  int pos = 0;
117  for (i = 0; i < this->size; i++)
118  {
119  Ap[i] = pos;
120  pos += this->sort_and_store_indices(this->pages[i], Ai + pos, Ai + aisize);
121  }
122  Ap[i] = pos;
123 
124  delete [] this->pages;
125  this->pages = NULL;
126 
127  nnz = Ap[this->size];
128 
129  Ax = new Scalar[nnz];
130  memset(Ax, 0, sizeof(Scalar) * nnz);
131  }
132 
133  template<typename Scalar>
134  void CSCMatrix<Scalar>::free()
135  {
136  nnz = 0;
137  if(Ap != NULL)
138  {
139  delete [] Ap;
140  Ap = NULL;
141  }
142  if(Ai != NULL)
143  {
144  delete [] Ai;
145  Ai = NULL;
146  }
147  if(Ax != NULL)
148  {
149  delete [] Ax;
150  Ax = NULL;
151  }
152  }
153 
154  template<typename Scalar>
155  Scalar CSCMatrix<Scalar>::get(unsigned int m, unsigned int n)
156  {
157  // Find m-th row in the n-th column.
158  int mid = find_position(Ai + Ap[n], Ap[n + 1] - Ap[n], m);
159 
160  if(mid < 0) // if the entry has not been found
161  return 0.0;
162  else
163  return Ax[Ap[n] + mid];
164  }
165 
166  template<typename Scalar>
167  void CSCMatrix<Scalar>::zero()
168  {
169  memset(Ax, 0, sizeof(Scalar) * nnz);
170  }
171 
172  template<>
173  void CSCMatrix<double>::add(unsigned int m, unsigned int n, double v)
174  {
175  if(v != 0.0) // ignore zero values.
176  {
177  // Find m-th row in the n-th column.
178  int pos = find_position(Ai + Ap[n], Ap[n + 1] - Ap[n], m);
179  // Make sure we are adding to an existing non-zero entry.
180  if(pos < 0)
181  {
182  this->info("CSCMatrix<Scalar>::add(): i = %d, j = %d.", m, n);
183  throw Hermes::Exceptions::Exception("Sparse matrix entry not found: [%i, %i]", m, n);
184  }
185 
186 #pragma omp atomic
187  Ax[Ap[n] + pos] += v;
188  }
189  }
190 
191  template<>
192  void CSCMatrix<std::complex<double> >::add(unsigned int m, unsigned int n, std::complex<double> v)
193  {
194  if(v != 0.0) // ignore zero values.
195  {
196  // Find m-th row in the n-th column.
197  int pos = find_position(Ai + Ap[n], Ap[n + 1] - Ap[n], m);
198  // Make sure we are adding to an existing non-zero entry.
199  if(pos < 0)
200  {
201  this->info("CSCMatrix<Scalar>::add(): i = %d, j = %d.", m, n);
202  throw Hermes::Exceptions::Exception("Sparse matrix entry not found: [%i, %i]", m, n);
203  }
204 
205 #pragma omp critical
206  Ax[Ap[n] + pos] += v;
207  }
208  }
209 
210  template<typename Scalar>
211  void CSCMatrix<Scalar>::add_to_diagonal_blocks(int num_stages, CSCMatrix<Scalar>* mat_block)
212  {
213  int ndof = mat_block->get_size();
214  if(this->get_size() != (unsigned int) num_stages * ndof)
215  throw Hermes::Exceptions::Exception("Incompatible matrix sizes in CSCMatrix<Scalar>::add_to_diagonal_blocks()");
216 
217  for (int i = 0; i < num_stages; i++)
218  {
219  this->add_as_block(ndof*i, ndof*i, mat_block);
220  }
221  }
222 
223  template<typename Scalar>
224  void CSCMatrix<Scalar>::add_sparse_to_diagonal_blocks(int num_stages, SparseMatrix<Scalar>* mat)
225  {
226  add_to_diagonal_blocks(num_stages, static_cast<CSCMatrix<Scalar>*>(mat));
227  }
228 
229  template<typename Scalar>
230  unsigned int CSCMatrix<Scalar>::get_nnz() const
231  {
232  return this->nnz;
233  }
234 
235  template<typename Scalar>
236  void CSCMatrix<Scalar>::add_as_block(unsigned int offset_i, unsigned int offset_j, CSCMatrix<Scalar>* mat)
237  {
238  UMFPackIterator<Scalar> mat_it(mat);
239  UMFPackIterator<Scalar> this_it(this);
240 
241  // Sanity check.
242  bool this_not_empty = this_it.init();
243  if(!this_not_empty) throw Hermes::Exceptions::Exception("Empty matrix detected in CSCMatrix<Scalar>::add_as_block().");
244 
245  // Iterate through the small matrix column by column and add all nonzeros
246  // to the large one.
247  bool mat_not_finished = mat_it.init();
248  if(!mat_not_finished) throw Hermes::Exceptions::Exception("Empty matrix detected in CSCMatrix<Scalar>::add_as_block().");
249 
250  int mat_i, mat_j;
251  Scalar mat_val;
252  while(mat_not_finished)
253  {
254  mat_it.get_current_position(mat_i, mat_j, mat_val);
255  bool found = this_it.move_to_position(mat_i + offset_i, mat_j + offset_j);
256  if(!found)
257  throw Hermes::Exceptions::Exception("Nonzero matrix entry at %d, %d not found in CSCMatrix<Scalar>::add_as_block().",
258  mat_i + offset_i, mat_j + offset_j);
259  this_it.add_to_current_position(mat_val);
260  mat_not_finished = mat_it.move_ptr();
261  }
262  }
263 
264  template<typename Scalar>
265  void CSCMatrix<Scalar>::add_matrix(CSCMatrix<Scalar>* mat)
266  {
267  assert(this->get_size() == mat->get_size());
268  // Create iterators for both matrices.
269  UMFPackIterator<Scalar> mat_it(mat);
270  UMFPackIterator<Scalar> this_it(this);
271  int mat_i, mat_j;
272  Scalar mat_val;
273  int this_i, this_j;
274  Scalar this_val;
275 
276  bool mat_not_finished = mat_it.init();
277  bool this_not_finished = this_it.init();
278  while(mat_not_finished && this_not_finished)
279  {
280  mat_it.get_current_position(mat_i, mat_j, mat_val);
281  //printf("mat: current position %d %d %g\n", mat_i, mat_j, mat_val);
282  this_it.get_current_position(this_i, this_j, this_val);
283  //printf("this: current position %d %d %g\n", this_i, this_j, this_val);
284  while(mat_i != this_i || mat_j != this_j)
285  {
286  //printf("SHOULD NOT BE HERE\n");
287  this_not_finished = this_it.move_ptr();
288  if(!this_not_finished)
289  {
290  printf("Entry %d %d does not exist in the matrix to which it is contributed.\n", mat_i, mat_j);
291  throw Hermes::Exceptions::Exception("Incompatible matrices in add_umfpack_matrix().");
292  }
293  this_it.get_current_position(this_i, this_j, this_val);
294  }
295  this_it.add_to_current_position(mat_val);
296  mat_not_finished = mat_it.move_ptr();
297  this_not_finished = this_it.move_ptr();
298  if(mat_not_finished && !this_not_finished)
299  throw Hermes::Exceptions::Exception("Incompatible matrices in add_umfpack_matrix().");
300  }
301  }
302 
303  template<typename Scalar>
304  void CSCMatrix<Scalar>::add_to_diagonal(Scalar v)
305  {
306  for (unsigned int i = 0; i<this->size; i++)
307  {
308  add(i, i, v);
309  }
310  };
311 
312  template<typename Scalar>
313  void CSCMatrix<Scalar>::add(unsigned int m, unsigned int n, Scalar **mat, int *rows, int *cols)
314  {
315  for (unsigned int i = 0; i < m; i++) // rows
316  for (unsigned int j = 0; j < n; j++) // cols
317  if(rows[i] >= 0 && cols[j] >= 0) // not Dir. dofs.
318  add(rows[i], cols[j], mat[i][j]);
319  }
320 
321  double inline real(double x)
322  {
323  return x;
324  }
325 
326  double inline imag(double x)
327  {
328  return 0;
329  }
330 
331  double inline real(std::complex<double> x)
332  {
333  return x.real();
334  }
335 
336  double inline imag(std::complex<double> x)
337  {
338  return x.imag();;
339  }
340 
341  template<>
342  bool CSCMatrix<double>::dump(FILE *file, const char *var_name, EMatrixDumpFormat fmt, char* number_format)
343  {
344  switch (fmt)
345  {
346  case DF_MATLAB_SPARSE:
347  fprintf(file, "%% Size: %dx%d\n%% Nonzeros: %d\ntemp = zeros(%d, 3);\ntemp =[\n",
348  this->size, this->size, nnz, nnz);
349  for (unsigned int j = 0; j < this->size; j++)
350  for (int i = Ap[j]; i < Ap[j + 1]; i++)
351  {
352  fprintf(file, "%d %d ", Ai[i] + 1, j + 1);
353  Hermes::Helpers::fprint_num(file, Ax[i], number_format);
354  fprintf(file, "\n");
355  }
356  fprintf(file, "];\n%s = spconvert(temp);\n", var_name);
357 
358  return true;
359 
360  case DF_MATRIX_MARKET:
361  {
362  fprintf(file, "%%%%Matrix<Scalar>Market matrix coordinate real symmetric\n");
363  int nnz_sym = 0;
364  for (unsigned int j = 0; j < this->size; j++)
365  for (int i = Ap[j]; i < Ap[j + 1]; i++)
366  if((int)j <= Ai[i]) nnz_sym++;
367  fprintf(file, "%d %d %d\n", this->size, this->size, nnz_sym);
368  for (unsigned int j = 0; j < this->size; j++)
369  for (int i = Ap[j]; i < Ap[j + 1]; i++)
370  // The following line was replaced with the one below, because it gave a warning
371  // to cause code abort at runtime.
372  //if(j <= Ai[i]) fprintf(file, "%d %d %24.15e\n", Ai[i] + 1, j + 1, Ax[i]);
373  if((int)j <= Ai[i])
374  {
375  fprintf(file, "%d %d ", Ai[i] + 1, (int)j + 1);
376  Hermes::Helpers::fprint_num(file, Ax[i], number_format);
377  fprintf(file, "\n");
378  }
379 
380  return true;
381  }
382 
383  case DF_HERMES_BIN:
384  {
385  hermes_fwrite("HERMESX\001", 1, 8, file);
386  int ssize = sizeof(double);
387  hermes_fwrite(&ssize, sizeof(int), 1, file);
388  hermes_fwrite(&this->size, sizeof(int), 1, file);
389  hermes_fwrite(&nnz, sizeof(int), 1, file);
390  hermes_fwrite(Ap, sizeof(int), this->size + 1, file);
391  hermes_fwrite(Ai, sizeof(int), nnz, file);
392  hermes_fwrite(Ax, sizeof(double), nnz, file);
393  return true;
394  }
395 
396  case DF_PLAIN_ASCII:
397  exit(1);
398  {
399  const double zero_cutoff = 1e-10;
400  double *ascii_entry_buff = new double[nnz];
401  int *ascii_entry_i = new int[nnz];
402  int *ascii_entry_j = new int[nnz];
403  int k = 0;
404 
405  // If real or imaginary part of Scalar entry is below zero_cutoff
406  // it's not included in ascii file, and number of non-zeros is reduced by one.
407  for (unsigned int j = 0; j < size; j++)
408  {
409  for (int i = Ap[j]; i < Ap[j + 1]; i++)
410  {
411  if(real(Ax[i]) > zero_cutoff || imag(Ax[i]) > zero_cutoff)
412  {
413  ascii_entry_buff[k] = Ax[i];
414  ascii_entry_i[k] = Ai[i];
415  ascii_entry_j[k] = j;
416  k++;
417  }
418  else
419  nnz -= 1;
420  }
421  }
422 
423  fprintf(file, "%d\n", size);
424  fprintf(file, "%d\n", nnz);
425  for (unsigned int k = 0; k < nnz; k++)
426  fprintf(file, "%d %d %f\n", ascii_entry_i[k], ascii_entry_j[k], ascii_entry_buff[k]);
427 
428  //Free memory
429  delete [] ascii_entry_buff;
430  delete [] ascii_entry_i;
431  delete [] ascii_entry_j;
432 
433  //Clear pointer
434  ascii_entry_buff = NULL;
435  ascii_entry_i = NULL;
436  ascii_entry_j = NULL;
437 
438  return true;
439  }
440 
441  default:
442  return false;
443  }
444  }
445 
446  template<>
447  bool CSCMatrix<std::complex<double> >::dump(FILE *file, const char *var_name, EMatrixDumpFormat fmt, char* number_format)
448  {
449  switch (fmt)
450  {
451  case DF_MATLAB_SPARSE:
452  fprintf(file, "%% Size: %dx%d\n%% Nonzeros: %d\ntemp = zeros(%d, 3);\ntemp =[\n",
453  this->size, this->size, nnz, nnz);
454  for (unsigned int j = 0; j < this->size; j++)
455  for (int i = Ap[j]; i < Ap[j + 1]; i++)
456  {
457  fprintf(file, "%d %d ", Ai[i] + 1, j + 1);
458  Hermes::Helpers::fprint_num(file, Ax[i], number_format);
459  fprintf(file, "\n");
460  }
461  fprintf(file, "];\n%s = spconvert(temp);\n", var_name);
462 
463  return true;
464 
465  case DF_MATRIX_MARKET:
466  {
467  fprintf(file, "%%%%Matrix<Scalar>Market matrix coordinate real symmetric\n");
468  int nnz_sym = 0;
469  for (unsigned int j = 0; j < this->size; j++)
470  for (int i = Ap[j]; i < Ap[j + 1]; i++)
471  if((int)j <= Ai[i]) nnz_sym++;
472  fprintf(file, "%d %d %d\n", this->size, this->size, nnz_sym);
473  for (unsigned int j = 0; j < this->size; j++)
474  for (int i = Ap[j]; i < Ap[j + 1]; i++)
475  // The following line was replaced with the one below, because it gave a warning
476  // to cause code abort at runtime.
477  //if(j <= Ai[i]) fprintf(file, "%d %d %24.15e\n", Ai[i] + 1, j + 1, Ax[i]);
478  if((int)j <= Ai[i])
479  {
480  fprintf(file, "%d %d ", Ai[i] + 1, (int)j + 1);
481  Hermes::Helpers::fprint_num(file, Ax[i], number_format);
482  fprintf(file, "\n");
483  }
484 
485  return true;
486  }
487 
488  case DF_HERMES_BIN:
489  {
490  hermes_fwrite("HERMESX\001", 1, 8, file);
491  int ssize = sizeof(std::complex<double>);
492  hermes_fwrite(&ssize, sizeof(int), 1, file);
493  hermes_fwrite(&this->size, sizeof(int), 1, file);
494  hermes_fwrite(&nnz, sizeof(int), 1, file);
495  hermes_fwrite(Ap, sizeof(int), this->size + 1, file);
496  hermes_fwrite(Ai, sizeof(int), nnz, file);
497  hermes_fwrite(Ax, sizeof(std::complex<double>), nnz, file);
498  return true;
499  }
500 
501  case DF_PLAIN_ASCII:
502  exit(1);
503  {
504  const double zero_cutoff = 1e-10;
505  std::complex<double> *ascii_entry_buff = new std::complex<double>[nnz];
506  int *ascii_entry_i = new int[nnz];
507  int *ascii_entry_j = new int[nnz];
508  int k = 0;
509 
510  // If real or imaginary part of Scalar entry is below zero_cutoff
511  // it's not included in ascii file, and number of non-zeros is reduced by one.
512  for (unsigned int j = 0; j < size; j++)
513  {
514  for (int i = Ap[j]; i < Ap[j + 1]; i++)
515  {
516  if(real(Ax[i]) > zero_cutoff || imag(Ax[i]) > zero_cutoff)
517  {
518  ascii_entry_buff[k] = Ax[i];
519  ascii_entry_i[k] = Ai[i];
520  ascii_entry_j[k] = j;
521  k++;
522  }
523  else
524  nnz -= 1;
525  }
526  }
527 
528  fprintf(file, "%d\n", size);
529  fprintf(file, "%d\n", nnz);
530  for (unsigned int k = 0; k < nnz; k++)
531  fprintf(file, "%d %d %E %E\n", ascii_entry_i[k], ascii_entry_j[k], ascii_entry_buff[k].real(), ascii_entry_buff[k].imag());
532 
533  //Free memory
534  delete [] ascii_entry_buff;
535  delete [] ascii_entry_i;
536  delete [] ascii_entry_j;
537 
538  //Clear pointer
539  ascii_entry_buff = NULL;
540  ascii_entry_i = NULL;
541  ascii_entry_j = NULL;
542 
543  return true;
544  }
545 
546  default:
547  return false;
548  }
549  }
550 
551  template<typename Scalar>
552  unsigned int CSCMatrix<Scalar>::get_matrix_size() const
553  {
554  return this->size;
555  }
556 
557  template<typename Scalar>
558  double CSCMatrix<Scalar>::get_fill_in() const
559  {
560  return nnz / (double) (this->size * this->size);
561  }
562 
563  template<typename Scalar>
564  void CSCMatrix<Scalar>::create(unsigned int size, unsigned int nnz, int* ap, int* ai, Scalar* ax)
565  {
566  this->nnz = nnz;
567  this->size = size;
568  this->Ap = new int[this->size + 1]; assert(this->Ap != NULL);
569  this->Ai = new int[nnz]; assert(this->Ai != NULL);
570  this->Ax = new Scalar[nnz]; assert(this->Ax != NULL);
571  for (unsigned int i = 0; i < this->size + 1; i++) this->Ap[i] = ap[i];
572  for (unsigned int i = 0; i < nnz; i++)
573  {
574  this->Ax[i] = ax[i];
575  this->Ai[i] = ai[i];
576  }
577  }
578 
579  template<typename Scalar>
580  CSCMatrix<Scalar>* CSCMatrix<Scalar>::duplicate()
581  {
582  CSCMatrix<Scalar>* new_matrix = new CSCMatrix<Scalar>();
583  new_matrix->create(this->get_size(), this->get_nnz(), this->get_Ap(), this->get_Ai(), this->get_Ax());
584  return new_matrix;
585  }
586 
587  template<typename Scalar>
588  int *CSCMatrix<Scalar>::get_Ap()
589  {
590  return this->Ap;
591  }
592 
593  template<typename Scalar>
594  int *CSCMatrix<Scalar>::get_Ai()
595  {
596  return this->Ai;
597  }
598 
599  template<typename Scalar>
600  Scalar *CSCMatrix<Scalar>::get_Ax()
601  {
602  return this->Ax;
603  }
604 
605  template<typename Scalar>
606  UMFPackVector<Scalar>::UMFPackVector()
607  {
608  v = NULL;
609  this->size = 0;
610  }
611 
612  template<typename Scalar>
613  UMFPackVector<Scalar>::UMFPackVector(unsigned int size)
614  {
615  v = NULL;
616  this->size = size;
617  this->alloc(size);
618  }
619 
620  template<typename Scalar>
621  UMFPackVector<Scalar>::~UMFPackVector()
622  {
623  free();
624  }
625 
626  template<typename Scalar>
627  void UMFPackVector<Scalar>::alloc(unsigned int n)
628  {
629  free();
630  this->size = n;
631  v = new Scalar[n];
632  this->zero();
633  }
634 
635  template<typename Scalar>
636  void UMFPackVector<Scalar>::zero()
637  {
638  memset(v, 0, this->size * sizeof(Scalar));
639  }
640 
641  template<typename Scalar>
642  void UMFPackVector<Scalar>::change_sign()
643  {
644  for (unsigned int i = 0; i < this->size; i++) v[i] *= -1.;
645  }
646 
647  template<typename Scalar>
648  void UMFPackVector<Scalar>::free()
649  {
650  delete [] v;
651  v = NULL;
652  this->size = 0;
653  }
654 
655  template<typename Scalar>
656  void UMFPackVector<Scalar>::set(unsigned int idx, Scalar y)
657  {
658  v[idx] = y;
659  }
660 
661  template<>
662  void UMFPackVector<double>::add(unsigned int idx, double y)
663  {
664 #pragma omp atomic
665  v[idx] += y;
666  }
667 
668  template<>
669  void UMFPackVector<std::complex<double> >::add(unsigned int idx, std::complex<double> y)
670  {
671 #pragma omp critical(UMFPackVector_add)
672  v[idx] += y;
673  }
674 
675  template<typename Scalar>
676  void UMFPackVector<Scalar>::add(unsigned int n, unsigned int *idx, Scalar *y)
677  {
678  for (unsigned int i = 0; i < n; i++)
679  v[idx[i]] += y[i];
680  }
681 
682  template<typename Scalar>
683  Scalar UMFPackVector<Scalar>::get(unsigned int idx)
684  {
685  return v[idx];
686  }
687 
688  template<typename Scalar>
689  void UMFPackVector<Scalar>::extract(Scalar *v) const
690  {
691  memcpy(v, this->v, this->size * sizeof(Scalar));
692  }
693 
694  template<typename Scalar>
695  void UMFPackVector<Scalar>::add_vector(Vector<Scalar>* vec)
696  {
697  assert(this->length() == vec->length());
698  for (unsigned int i = 0; i < this->length(); i++) this->v[i] += vec->get(i);
699  }
700 
701  template<typename Scalar>
702  void UMFPackVector<Scalar>::add_vector(Scalar* vec)
703  {
704  for (unsigned int i = 0; i < this->length(); i++) this->v[i] += vec[i];
705  }
706 
707  template<typename Scalar>
708  Scalar *UMFPackVector<Scalar>::get_c_array()
709  {
710  return this->v;
711  }
712 
713  template<>
714  bool UMFPackVector<double>::dump(FILE *file, const char *var_name, EMatrixDumpFormat fmt, char* number_format)
715  {
716  switch (fmt)
717  {
718  case DF_MATLAB_SPARSE:
719  fprintf(file, "%% Size: %dx1\n%s =[\n", this->size, var_name);
720  for (unsigned int i = 0; i < this->size; i++)
721  {
722  Hermes::Helpers::fprint_num(file, v[i], number_format);
723  fprintf(file, "\n");
724  }
725  fprintf(file, " ];\n");
726  return true;
727 
728  case DF_HERMES_BIN:
729  {
730  hermes_fwrite("HERMESR\001", 1, 8, file);
731  int ssize = sizeof(double);
732  hermes_fwrite(&ssize, sizeof(int), 1, file);
733  hermes_fwrite(&this->size, sizeof(int), 1, file);
734  hermes_fwrite(v, sizeof(double), this->size, file);
735  return true;
736  }
737 
738  case DF_PLAIN_ASCII:
739  {
740  fprintf(file, "\n");
741  for (unsigned int i = 0; i < size; i++)
742  {
743  Hermes::Helpers::fprint_num(file, v[i], number_format);
744  fprintf(file, "\n");
745  }
746 
747  return true;
748  }
749 
750  default:
751  return false;
752  }
753  }
754 
755  template<>
756  bool UMFPackVector<std::complex<double> >::dump(FILE *file, const char *var_name, EMatrixDumpFormat fmt, char* number_format)
757  {
758  switch (fmt)
759  {
760  case DF_MATLAB_SPARSE:
761  fprintf(file, "%% Size: %dx1\n%s =[\n", this->size, var_name);
762  for (unsigned int i = 0; i < this->size; i++)
763  {
764  Hermes::Helpers::fprint_num(file, v[i], number_format);
765  fprintf(file, "\n");
766  }
767  fprintf(file, " ];\n");
768  return true;
769 
770  case DF_HERMES_BIN:
771  {
772  hermes_fwrite("HERMESR\001", 1, 8, file);
773  int ssize = sizeof(std::complex<double>);
774  hermes_fwrite(&ssize, sizeof(int), 1, file);
775  hermes_fwrite(&this->size, sizeof(int), 1, file);
776  hermes_fwrite(v, sizeof(std::complex<double>), this->size, file);
777  return true;
778  }
779 
780  case DF_PLAIN_ASCII:
781  {
782  fprintf(file, "\n");
783  for (unsigned int i = 0; i < size; i++)
784  {
785  fprintf(file, "%E %E\n", v[i].real(), v[i].imag());
786  }
787 
788  return true;
789  }
790 
791  default:
792  return false;
793  }
794  }
795 
796  template class HERMES_API CSCMatrix<double>;
797  template class HERMES_API CSCMatrix<std::complex<double> >;
798  template class HERMES_API UMFPackMatrix<double>;
799  template class HERMES_API UMFPackMatrix<std::complex<double> >;
800  template class HERMES_API UMFPackVector<double>;
801  template class HERMES_API UMFPackVector<std::complex<double> >;
802  }
803 
804  namespace Solvers
805  {
806  template<typename Scalar>
807  bool UMFPackIterator<Scalar>::init()
808  {
809  if(this->size == 0 || this->nnz == 0) return false;
810  this->Ap_pos = 0;
811  this->Ai_pos = 0;
812  return true;
813  }
814 
815  template<typename Scalar>
816  UMFPackIterator<Scalar>::UMFPackIterator(CSCMatrix<Scalar>* mat)
817  {
818  this->size = mat->get_size();
819  this->nnz = mat->get_nnz();
820  this->Ai = mat->get_Ai();
821  this->Ap = mat->get_Ap();
822  this->Ax = mat->get_Ax();
823  this->Ai_pos = 0;
824  this->Ap_pos = 0;
825  }
826 
827  template<typename Scalar>
828  void UMFPackIterator<Scalar>::get_current_position(int& i, int& j, Scalar& val)
829  {
830  i = Ai[Ai_pos];
831  j = Ap_pos;
832  val = Ax[Ai_pos];
833  }
834 
835  template<typename Scalar>
836  bool UMFPackIterator<Scalar>::move_to_position(int i, int j)
837  {
838  int ii, jj;
839  Scalar val;
840  get_current_position(ii, jj, val);
841  while (!(ii == i && jj == j))
842  {
843  if(!this->move_ptr()) return false;
844  get_current_position(ii, jj, val);
845  }
846  return true;
847  }
848 
849  template<typename Scalar>
850  bool UMFPackIterator<Scalar>::move_ptr()
851  {
852  if(Ai_pos >= nnz - 1) return false; // It is no longer possible to find next element.
853  if(Ai_pos + 1 >= Ap[Ap_pos + 1])
854  {
855  Ap_pos++;
856  }
857  Ai_pos++;
858  return true;
859  }
860 
861  template<typename Scalar>
862  void UMFPackIterator<Scalar>::add_to_current_position(Scalar val)
863  {
864  this->Ax[this->Ai_pos] += val;
865  }
866 
867  template<>
868  bool UMFPackLinearMatrixSolver<double>::setup_factorization()
869  {
870  // Perform both factorization phases for the first time.
871  int eff_fact_scheme;
872  if(factorization_scheme != HERMES_FACTORIZE_FROM_SCRATCH && symbolic == NULL && numeric == NULL)
873  eff_fact_scheme = HERMES_FACTORIZE_FROM_SCRATCH;
874  else
875  eff_fact_scheme = factorization_scheme;
876 
877  int status;
878  switch(eff_fact_scheme)
879  {
880  case HERMES_FACTORIZE_FROM_SCRATCH:
881  if(symbolic != NULL) umfpack_di_free_symbolic(&symbolic);
882 
883  //debug_log("Factorizing symbolically.");
884  status = umfpack_di_symbolic(m->get_size(), m->get_size(), m->get_Ap(), m->get_Ai(), m->get_Ax(), &symbolic, NULL, NULL);
885  if(status != UMFPACK_OK)
886  {
887  check_status("umfpack_di_symbolic", status);
888  return false;
889  }
890  if(symbolic == NULL)
891  throw Exceptions::Exception("umfpack_di_symbolic error: symbolic == NULL");
892 
893  case HERMES_REUSE_MATRIX_REORDERING:
894  case HERMES_REUSE_MATRIX_REORDERING_AND_SCALING:
895  if(numeric != NULL) umfpack_di_free_numeric(&numeric);
896 
897  //debug_log("Factorizing numerically.");
898  status = umfpack_di_numeric(m->get_Ap(), m->get_Ai(), m->get_Ax(), symbolic, &numeric, NULL, NULL);
899  if(status != UMFPACK_OK)
900  {
901  check_status("umfpack_di_numeric", status);
902  return false;
903  }
904  if(numeric == NULL)
905  throw Exceptions::Exception("umfpack_di_numeric error: numeric == NULL");
906  }
907 
908  return true;
909  }
910 
911  template<typename Scalar>
912  UMFPackLinearMatrixSolver<Scalar>::UMFPackLinearMatrixSolver(UMFPackMatrix<Scalar> *m, UMFPackVector<Scalar> *rhs)
913  : DirectSolver<Scalar>(HERMES_FACTORIZE_FROM_SCRATCH), m(m), rhs(rhs), symbolic(NULL), numeric(NULL)
914  {
915  }
916 
917  template<typename Scalar>
918  UMFPackLinearMatrixSolver<Scalar>::~UMFPackLinearMatrixSolver()
919  {
920  free_factorization_data();
921  }
922 
923  template<typename Scalar>
924  int UMFPackLinearMatrixSolver<Scalar>::get_matrix_size()
925  {
926  return m->get_size();
927  }
928 
929  template<>
930  bool UMFPackLinearMatrixSolver<std::complex<double> >::setup_factorization()
931  {
932  // Perform both factorization phases for the first time.
933  int eff_fact_scheme;
934  if(factorization_scheme != HERMES_FACTORIZE_FROM_SCRATCH && symbolic == NULL && numeric == NULL)
935  eff_fact_scheme = HERMES_FACTORIZE_FROM_SCRATCH;
936  else
937  eff_fact_scheme = factorization_scheme;
938 
939  int status;
940  switch(eff_fact_scheme)
941  {
942  case HERMES_FACTORIZE_FROM_SCRATCH:
943  if(symbolic != NULL)
944  umfpack_zi_free_symbolic(&symbolic);
945 
946  status = umfpack_zi_symbolic(m->get_size(), m->get_size(), m->get_Ap(), m->get_Ai(), (double *)m->get_Ax(), NULL, &symbolic, NULL, NULL);
947  if(status != UMFPACK_OK)
948  {
949  check_status("umfpack_di_symbolic", status);
950  return false;
951  }
952  if(symbolic == NULL)
953  throw Exceptions::Exception("umfpack_di_symbolic error: symbolic == NULL");
954 
955  case HERMES_REUSE_MATRIX_REORDERING:
956  case HERMES_REUSE_MATRIX_REORDERING_AND_SCALING:
957  if(numeric != NULL)
958  umfpack_zi_free_numeric(&numeric);
959 
960  status = umfpack_zi_numeric(m->get_Ap(), m->get_Ai(), (double *) m->get_Ax(), NULL, symbolic, &numeric, NULL, NULL);
961  if(status != UMFPACK_OK)
962  {
963  check_status("umfpack_di_numeric", status);
964  return false;
965  }
966  if(numeric == NULL)
967  throw Exceptions::Exception("umfpack_di_numeric error: numeric == NULL");
968  }
969 
970  return true;
971  }
972 
973  template<>
974  void UMFPackLinearMatrixSolver<double>::free_factorization_data()
975  {
976  if(symbolic != NULL) umfpack_di_free_symbolic(&symbolic);
977  symbolic = NULL;
978  if(numeric != NULL) umfpack_di_free_numeric(&numeric);
979  numeric = NULL;
980  }
981 
982  template<>
983  void UMFPackLinearMatrixSolver<std::complex<double> >::free_factorization_data()
984  {
985  if(symbolic != NULL) umfpack_zi_free_symbolic(&symbolic);
986  symbolic = NULL;
987  if(numeric != NULL) umfpack_zi_free_numeric(&numeric);
988  numeric = NULL;
989  }
990 
991  template<>
992  bool UMFPackLinearMatrixSolver<double>::solve()
993  {
994  assert(m != NULL);
995  assert(rhs != NULL);
996  assert(m->get_size() == rhs->length());
997 
998  this->tick();
999 
1000  if( !setup_factorization() )
1001  throw Exceptions::LinearMatrixSolverException("LU factorization could not be completed.");
1002 
1003  if(sln != NULL)
1004  delete [] sln;
1005  sln = new double[m->get_size()];
1006  memset(sln, 0, m->get_size() * sizeof(double));
1007  int status = umfpack_di_solve(UMFPACK_A, m->get_Ap(), m->get_Ai(), m->get_Ax(), sln, rhs->get_c_array(), numeric, NULL, NULL);
1008  if(status != UMFPACK_OK)
1009  {
1010  check_status("umfpack_di_solve", status);
1011  return false;
1012  }
1013 
1014  this->tick();
1015  time = this->accumulated();
1016 
1017  return true;
1018  }
1019 
1020  template<>
1021  bool UMFPackLinearMatrixSolver<std::complex<double> >::solve()
1022  {
1023  assert(m != NULL);
1024  assert(rhs != NULL);
1025  assert(m->get_size() == rhs->length());
1026 
1027  this->tick();
1028  if( !setup_factorization() )
1029  {
1030  this->warn("LU factorization could not be completed.");
1031  return false;
1032  }
1033 
1034  if(sln)
1035  delete [] sln;
1036  sln = new std::complex<double>[m->get_size()];
1037  memset(sln, 0, m->get_size() * sizeof(std::complex<double>));
1038  int status = umfpack_zi_solve(UMFPACK_A, m->get_Ap(), m->get_Ai(), (double *)m->get_Ax(), NULL, (double*) sln, NULL, (double *)rhs->get_c_array(), NULL, numeric, NULL, NULL);
1039  if(status != UMFPACK_OK)
1040  {
1041  check_status("umfpack_di_solve", status);
1042  return false;
1043  }
1044 
1045  this->tick();
1046  time = this->accumulated();
1047 
1048  return true;
1049  }
1050 
1051  template<typename Scalar>
1052  void UMFPackLinearMatrixSolver<Scalar>::check_status(const char *fn_name, int status)
1053  {
1054  switch (status)
1055  {
1056  case UMFPACK_OK: break;
1057  case UMFPACK_WARNING_singular_matrix: this->warn("%s: singular matrix!", fn_name); break;
1058  case UMFPACK_ERROR_out_of_memory: this->warn("%s: out of memory!", fn_name); break;
1059  case UMFPACK_ERROR_argument_missing: this->warn("%s: argument missing", fn_name); break;
1060  case UMFPACK_ERROR_invalid_Symbolic_object: this->warn("%s: invalid Symbolic object", fn_name); break;
1061  case UMFPACK_ERROR_invalid_Numeric_object: this->warn("%s: invalid Numeric object", fn_name); break;
1062  case UMFPACK_ERROR_different_pattern: this->warn("%s: different pattern", fn_name); break;
1063  case UMFPACK_ERROR_invalid_system: this->warn("%s: invalid system", fn_name); break;
1064  case UMFPACK_ERROR_n_nonpositive: this->warn("%s: n nonpositive", fn_name); break;
1065  case UMFPACK_ERROR_invalid_matrix: this->warn("%s: invalid matrix", fn_name); break;
1066  case UMFPACK_ERROR_internal_error: this->warn("%s: internal error", fn_name); break;
1067  default: this->warn("%s: unknown error (%d)", fn_name, status); break;
1068  }
1069  }
1070 
1071  template class HERMES_API UMFPackLinearMatrixSolver<double>;
1072  template class HERMES_API UMFPackLinearMatrixSolver<std::complex<double> >;
1073  }
1074 }
1075 #endif