general_purpose_preconditioners.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Config header generated by autoconfig
27 #ifdef HAVE_CONFIG_H
28 #include <oomph-lib-config.h>
29 #endif
30 
31 
32 // oomph-lib includes
34 
35 
36 namespace oomph
37 {
38  //=============================================================
39  /// Setup diagonal preconditioner: Store the inverse of the
40  /// diagonal entries from the fully
41  /// assembled matrix.
42  //=============================================================
44  {
45  // first attempt to cast to DistributableLinearAlgebraObject
46  DistributableLinearAlgebraObject* dist_matrix_pt =
48 
49  // if it is a distributable matrix
50  if (dist_matrix_pt != 0)
51  {
52  // cache the number of first_rows and nrow_local
53  unsigned nrow_local = dist_matrix_pt->nrow_local();
54  unsigned first_row = dist_matrix_pt->first_row();
55 
56  // resize the inverse diagonal storage
57  Inv_diag.resize(nrow_local);
58 
59  // Extract the diagonal entries
60  for (unsigned i = 0; i < nrow_local; i++)
61  {
62  unsigned index = i + first_row;
63 #ifdef PARANOID
64  if ((*matrix_pt())(i, index) == 0.0)
65  {
66  throw OomphLibError(
67  "Zero diagonal in matrix --> Cannot use diagonal preconditioner.",
68  OOMPH_CURRENT_FUNCTION,
69  OOMPH_EXCEPTION_LOCATION);
70  }
71 #endif
72  Inv_diag[i] = 1.0 / (*matrix_pt())(i, index);
73  }
74 
75  // store the distribution
76  this->build_distribution(dist_matrix_pt->distribution_pt());
77  }
78 
79  // else it is not a distributable matrix
80  else
81  {
82  // # of rows in the matrix
83  unsigned n_row = matrix_pt()->nrow();
84 
85  // Resize the Inv_diag vector to accommodate the # of
86  // diagonal entries
87  Inv_diag.resize(n_row);
88 
89  // Extract the diagonal entries
90  for (unsigned i = 0; i < n_row; i++)
91  {
92 #ifdef PARANOID
93  if ((*matrix_pt())(i, i) == 0.0)
94  {
95  throw OomphLibError(
96  "Zero diagonal in matrix --> Cannot use diagonal preconditioner.",
97  OOMPH_CURRENT_FUNCTION,
98  OOMPH_EXCEPTION_LOCATION);
99  }
100  else
101 #endif
102  {
103  Inv_diag[i] = 1.0 / (*matrix_pt())(i, i);
104  }
105  }
106 
107  // create the distribution
108  LinearAlgebraDistribution dist(comm_pt(), n_row, false);
109  this->build_distribution(dist);
110  }
111  }
112 
113 
114  //=============================================================================
115  /// Apply preconditioner: Multiply r by the inverse of the diagonal.
116  //=============================================================================
118  const DoubleVector& r, DoubleVector& z)
119  {
120 #ifdef PARANOID
121  if (*r.distribution_pt() != *this->distribution_pt())
122  {
123  std::ostringstream error_message_stream;
124  error_message_stream
125  << "The r vector must have the same distribution as the "
126  "preconditioner. "
127  << "(this is the same as the matrix passed to setup())";
128  throw OomphLibError(error_message_stream.str(),
129  OOMPH_CURRENT_FUNCTION,
130  OOMPH_EXCEPTION_LOCATION);
131  }
132  if (z.built())
133  {
134  if (*z.distribution_pt() != *this->distribution_pt())
135  {
136  std::ostringstream error_message_stream;
137  error_message_stream
138  << "The z vector distribution has been setup; it must have the "
139  << "same distribution as the r vector (and preconditioner).";
140  throw OomphLibError(error_message_stream.str(),
141  OOMPH_CURRENT_FUNCTION,
142  OOMPH_EXCEPTION_LOCATION);
143  }
144  }
145 #endif
146 
147  // if z has not been setup then rebuild it
148  if (!z.built())
149  {
150  z.build(this->distribution_pt(), 0.0);
151  }
152 
153  // apply the preconditioner
154  const double* r_values = r.values_pt();
155  double* z_values = z.values_pt();
156  unsigned nrow_local = this->nrow_local();
157  for (unsigned i = 0; i < nrow_local; i++)
158  {
159  z_values[i] = Inv_diag[i] * r_values[i];
160  }
161  }
162 
163  //=============================================================================
164  /// Setup the lumped preconditioner. Specialisation for CCDoubleMatrix.
165  //=============================================================================
166  template<>
168  {
169  // # of rows in the matrix
170  Nrow = matrix_pt()->nrow();
171 
172  // Create the vector for the inverse lumped
173  if (Inv_lumped_diag_pt != 0)
174  {
175  delete[] this->Inv_lumped_diag_pt;
176  }
177  Inv_lumped_diag_pt = new double[this->Nrow];
178 
179  // zero the vector
180  for (unsigned i = 0; i < Nrow; i++)
181  {
182  Inv_lumped_diag_pt[i] = 0.0;
183  }
184 
185  // cast the Double Base Matrix to Compressed Column Double Matrix
186  CCDoubleMatrix* cc_matrix_pt = dynamic_cast<CCDoubleMatrix*>(matrix_pt());
187 
188  // get the matrix
189  int* m_row_index = cc_matrix_pt->row_index();
190  double* m_value = cc_matrix_pt->value();
191  unsigned m_nnz = cc_matrix_pt->nnz();
192 
193  // intially set positive matrix to true
194  Positive_matrix = true;
195 
196  // lump the matrix
197  for (unsigned i = 0; i < m_nnz; i++)
198  {
199  // if the matrix contains negative coefficient the matrix not positive
200  if (m_value[i] < 0.0)
201  {
202  Positive_matrix = false;
203  }
204 
205  // computed lumped matrix - temporarily stored in Inv_lumped_diag_pt
206  Inv_lumped_diag_pt[m_row_index[i]] += m_value[i];
207  }
208 
209  // invert the lumped matrix
210  for (unsigned i = 0; i < Nrow; i++)
211  {
212  Inv_lumped_diag_pt[i] = 1.0 / Inv_lumped_diag_pt[i];
213  }
214 
215  // create the distribution
216  LinearAlgebraDistribution dist(comm_pt(), Nrow, false);
217  this->build_distribution(dist);
218  }
219 
220  //=============================================================================
221  /// Setup the lumped preconditioner. Specialisation for CRDoubleMatrix.
222  //=============================================================================
223  template<>
225  {
226  // first attempt to cast to CRDoubleMatrix
227  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
228 
229  // # of rows in the matrix
230  Nrow = cr_matrix_pt->nrow_local();
231 
232  // Create the vector for the inverse lumped
233  if (Inv_lumped_diag_pt != 0)
234  {
235  delete[] this->Inv_lumped_diag_pt;
236  }
237  Inv_lumped_diag_pt = new double[this->Nrow];
238 
239  // zero the vector
240  for (unsigned i = 0; i < Nrow; i++)
241  {
242  Inv_lumped_diag_pt[i] = 0.0;
243  }
244 
245  // get the matrix
246  int* m_row_start = cr_matrix_pt->row_start();
247  double* m_value = cr_matrix_pt->value();
248 
249  // intially set positive matrix to true
250  Positive_matrix = true;
251 
252  // lump and invert matrix
253  for (unsigned i = 0; i < Nrow; i++)
254  {
255  Inv_lumped_diag_pt[i] = 0.0;
256  for (int j = m_row_start[i]; j < m_row_start[i + 1]; j++)
257  {
258  // if the matrix contains negative coefficient the matrix not positive
259  if (m_value[j] < 0.0)
260  {
261  Positive_matrix = false;
262  }
263 
264  // computed lumped coef (i,i)- temporarily stored in Inv_lumped_diag_pt
265  Inv_lumped_diag_pt[i] += m_value[j];
266  }
267  // invert coef (i,i)
268  Inv_lumped_diag_pt[i] = 1.0 / Inv_lumped_diag_pt[i];
269  }
270 
271  // store the distribution
272  this->build_distribution(cr_matrix_pt->distribution_pt());
273  }
274 
275 
276  //=============================================================================
277  /// Apply preconditioner: Multiply r by the inverse of the lumped matrix
278  //=============================================================================
279  template<typename MATRIX>
281  const DoubleVector& r, DoubleVector& z)
282  {
283 #ifdef PARANOID
284  if (*r.distribution_pt() != *this->distribution_pt())
285  {
286  std::ostringstream error_message_stream;
287  error_message_stream
288  << "The r vector must have teh same distribution as the "
289  "preconditioner. "
290  << "(this is the same as the matrix passed to setup())";
291  throw OomphLibError(error_message_stream.str(),
292  OOMPH_CURRENT_FUNCTION,
293  OOMPH_EXCEPTION_LOCATION);
294  }
295  if (z.built())
296  {
297  if (*z.distribution_pt() != *this->distribution_pt())
298  {
299  std::ostringstream error_message_stream;
300  error_message_stream
301  << "The z vector distribution has been setup; it must have the "
302  << "same distribution as the r vector (and preconditioner).";
303  throw OomphLibError(error_message_stream.str(),
304  OOMPH_CURRENT_FUNCTION,
305  OOMPH_EXCEPTION_LOCATION);
306  }
307  }
308 #endif
309 
310  z.build(r.distribution_pt(), 0.0);
311  for (unsigned i = 0; i < Nrow; i++)
312  {
313  z[i] = Inv_lumped_diag_pt[i] * r[i];
314  }
315  }
316 
317 
318  // ensure the lumped preconditioner get built
321 
322 
323  //=============================================================================
324  /// setup ILU(0) preconditioner for Matrices of CCDoubleMatrix type
325  //=============================================================================
327  {
328  // cast the Double Base Matrix to Compressed Column Double Matrix
329  CCDoubleMatrix* cc_matrix_pt = dynamic_cast<CCDoubleMatrix*>(matrix_pt());
330 
331 #ifdef PARANOID
332  if (cc_matrix_pt == 0)
333  {
334  std::ostringstream error_msg;
335  error_msg << "Failed to conver matrix_pt to CCDoubleMatrix*.";
336  throw OomphLibError(
337  error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
338  }
339 #endif
340 
341  // number of rows in matrix
342  int n_row = cc_matrix_pt->nrow();
343 
344  // set the distribution
345  LinearAlgebraDistribution dist(comm_pt(), n_row, false);
346  this->build_distribution(dist);
347 
348  // declares variables to store number of non zero entires in L and U
349  int l_nz = 0;
350  int u_nz = 0;
351 
352  // create space for m matrix
353  int* m_column_start;
354  int* m_row_index;
355  double* m_value;
356 
357  // get the m matrix
358  m_column_start = cc_matrix_pt->column_start();
359  m_row_index = cc_matrix_pt->row_index();
360  m_value = cc_matrix_pt->value();
361 
362  // find number non zero entries in L and U
363  for (int i = 0; i < n_row; i++)
364  {
365  for (int j = m_column_start[i]; j < m_column_start[i + 1]; j++)
366  {
367  if (m_row_index[j] > i)
368  {
369  l_nz++;
370  }
371  else
372  {
373  u_nz++;
374  }
375  }
376  }
377 
378  // resize vectors to store the data for the lower prior to building the
379  // matrices
380  L_column_start.resize(n_row + 1);
381  L_row_entry.resize(l_nz);
382 
383  // and the upper matrix
384  U_column_start.resize(n_row + 1);
385  U_row_entry.resize(u_nz);
386 
387  // set first column pointers to zero
388  L_column_start[0] = 0;
389  U_column_start[0] = 0;
390 
391  // split the matrix into L and U
392  for (int i = 0; i < n_row; i++)
393  {
394  L_column_start[i + 1] = L_column_start[i];
395  U_column_start[i + 1] = U_column_start[i];
396  for (int j = m_column_start[i]; j < m_column_start[i + 1]; j++)
397  {
398  if (m_row_index[j] > i)
399  {
400  int k = L_column_start[i + 1]++;
401  L_row_entry[k].index() = m_row_index[j];
402  L_row_entry[k].value() = m_value[j];
403  }
404  else
405  {
406  int k = U_column_start[i + 1]++;
407  U_row_entry[k].index() = m_row_index[j];
408  U_row_entry[k].value() = m_value[j];
409  }
410  }
411  }
412 
413  // sort each row entry vector into row index order for each column
414 
415  // loop over the columns
416  for (unsigned i = 0; i < unsigned(n_row); i++)
417  {
418  // sort the columns of the L matrix
419  std::sort(L_row_entry.begin() + L_column_start[i],
420  L_row_entry.begin() + L_column_start[i + 1]);
421 
422  // sort the columns of the U matrix
423  std::sort(U_row_entry.begin() + U_column_start[i],
424  U_row_entry.begin() + U_column_start[i + 1]);
425  }
426 
427 
428  // factorise matrix
429  int i;
430  unsigned j, pn, qn, rn;
431  pn = 0;
432  qn = 0;
433  rn = 0;
434  double multiplier;
435  for (i = 0; i < n_row - 1; i++)
436  {
437  multiplier = U_row_entry[U_column_start[i + 1] - 1].value();
438  for (j = L_column_start[i]; j < L_column_start[i + 1]; j++)
439  L_row_entry[j].value() /= multiplier;
440  for (j = U_column_start[i + 1]; j < U_column_start[i + 2] - 1; j++)
441  {
442  multiplier = U_row_entry[j].value();
443  qn = j + 1;
444  rn = L_column_start[i + 1];
445  for (pn = L_column_start[U_row_entry[j].index()];
446  pn < L_column_start[U_row_entry[j].index() + 1] &&
447  static_cast<int>(L_row_entry[pn].index()) <= i + 1;
448  pn++)
449  {
450  while (qn < U_column_start[i + 2] &&
451  U_row_entry[qn].index() < L_row_entry[pn].index())
452  qn++;
453  if (qn < U_column_start[i + 2] &&
454  L_row_entry[pn].index() == U_row_entry[qn].index())
455  U_row_entry[qn].value() -= multiplier * L_row_entry[pn].value();
456  }
457  for (; pn < L_column_start[U_row_entry[j].index() + 1]; pn++)
458  {
459  while (rn < L_column_start[i + 2] &&
460  L_row_entry[rn].index() < L_row_entry[pn].index())
461  rn++;
462  if (rn < L_column_start[i + 2] &&
463  L_row_entry[pn].index() == L_row_entry[rn].index())
464  L_row_entry[rn].value() -= multiplier * L_row_entry[pn].value();
465  }
466  }
467  }
468  }
469 
470 
471  //=============================================================================
472  /// setup ILU(0) preconditioner for Matrices of CRDoubleMatrix Type
473  //=============================================================================
475  {
476  // cast the Double Base Matrix to Compressed Column Double Matrix
477  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
478 
479 #ifdef PARANOID
480  if (cr_matrix_pt == 0)
481  {
482  std::ostringstream error_msg;
483  error_msg << "Failed to conver matrix_pt to CRDoubleMatrix*.";
484  throw OomphLibError(
485  error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
486  }
487 #endif
488 
489  // if the matrix is distributed then build global version
490  bool built_global = false;
491  if (cr_matrix_pt->distributed())
492  {
493  // get the global matrix
494  CRDoubleMatrix* global_matrix_pt = cr_matrix_pt->global_matrix();
495 
496  // store at cr_matrix pointer
497  cr_matrix_pt = global_matrix_pt;
498 
499  // set the flag so we can delete later
500  built_global = true;
501  }
502 
503  // store the Distribution
504  this->build_distribution(cr_matrix_pt->distribution_pt());
505 
506  // number of rows in matrix
507  int n_row = cr_matrix_pt->nrow();
508 
509  // declares variables to store number of non zero entires in L and U
510  int l_nz = 0;
511  int u_nz = 0;
512 
513  // create space for m matrix
514  int* m_row_start;
515  int* m_column_index;
516  double* m_value;
517 
518  // get the m matrix
519  m_row_start = cr_matrix_pt->row_start();
520  m_column_index = cr_matrix_pt->column_index();
521  m_value = cr_matrix_pt->value();
522 
523  // find number non zero entries in L and U
524  for (int i = 0; i < n_row; i++)
525  {
526  for (int j = m_row_start[i]; j < m_row_start[i + 1]; j++)
527  {
528  if (m_column_index[j] < i)
529  {
530  l_nz++;
531  }
532  else
533  {
534  u_nz++;
535  }
536  }
537  }
538 
539  // resize vectors to store the data for the lower prior to building the
540  // matrices
541  L_row_start.resize(n_row + 1);
542  L_row_entry.resize(l_nz);
543 
544  // and the upper matrix
545  U_row_start.resize(n_row + 1);
546  U_row_entry.resize(u_nz);
547 
548  // set first column pointers to zero
549  L_row_start[0] = 0;
550  U_row_start[0] = 0;
551 
552  // split the matrix into L and U
553  for (int i = 0; i < n_row; i++)
554  {
555  L_row_start[i + 1] = L_row_start[i];
556  U_row_start[i + 1] = U_row_start[i];
557  for (int j = m_row_start[i]; j < m_row_start[i + 1]; j++)
558  {
559  if (m_column_index[j] < i)
560  {
561  int k = L_row_start[i + 1]++;
562  L_row_entry[k].value() = m_value[j];
563  L_row_entry[k].index() = m_column_index[j];
564  }
565  else
566  {
567  int k = U_row_start[i + 1]++;
568  U_row_entry[k].value() = m_value[j];
569  U_row_entry[k].index() = m_column_index[j];
570  }
571  }
572  }
573 
574 
575  // factorise matrix
576  unsigned i, j, pn, qn, rn;
577  pn = 0;
578  qn = 0;
579  rn = 0;
580  double multiplier;
581  for (i = 1; i < static_cast<unsigned>(n_row); i++)
582  {
583  for (j = L_row_start[i]; j < L_row_start[i + 1]; j++)
584  {
585  pn = U_row_start[L_row_entry[j].index()];
586  multiplier = (L_row_entry[j].value() /= U_row_entry[pn].value());
587  qn = j + 1;
588  rn = U_row_start[i];
589  for (pn++; pn < U_row_start[L_row_entry[j].index() + 1] &&
590  U_row_entry[pn].index() < i;
591  pn++)
592  {
593  while (qn < L_row_start[i + 1] &&
594  L_row_entry[qn].index() < U_row_entry[pn].index())
595  qn++;
596  if (qn < L_row_start[i + 1] &&
597  U_row_entry[pn].index() == L_row_entry[qn].index())
598  L_row_entry[qn].value() -= multiplier * U_row_entry[pn].value();
599  }
600  for (; pn < U_row_start[L_row_entry[j].index() + 1]; pn++)
601  {
602  while (rn < U_row_start[i + 1] &&
603  U_row_entry[rn].index() < U_row_entry[pn].index())
604  rn++;
605  if (rn < U_row_start[i + 1] &&
606  U_row_entry[pn].index() == U_row_entry[rn].index())
607  U_row_entry[rn].value() -= multiplier * U_row_entry[pn].value();
608  }
609  }
610  }
611 
612  // if we built the global matrix then delete it
613  if (built_global)
614  {
615  delete cr_matrix_pt;
616  }
617  }
618 
619 
620  //=============================================================================
621  /// Apply ILU(0) preconditioner for CCDoubleMatrix: Solve Ly=r then
622  /// Uz=y and return z
623  //=============================================================================
625  const DoubleVector& r, DoubleVector& z)
626  {
627  // # of rows in the matrix
628  int n_row = r.nrow();
629 
630  // store the distribution of z
631  LinearAlgebraDistribution* z_dist = 0;
632  if (z.built())
633  {
634  z_dist = new LinearAlgebraDistribution(z.distribution_pt());
635  }
636 
637  // copy r to z
638  z = r;
639 
640  // if z is distributed then change to global
641  if (z.distributed())
642  {
643  z.redistribute(this->distribution_pt());
644  }
645 
646  // solve Ly=r (note L matrix is unit and diagonal is not stored)
647  for (unsigned i = 0; i < static_cast<unsigned>(n_row); i++)
648  {
649  for (unsigned j = L_column_start[i]; j < L_column_start[i + 1]; j++)
650  {
651  z[L_row_entry[j].index()] =
652  z[L_row_entry[j].index()] - z[i] * L_row_entry[j].value();
653  }
654  }
655 
656  // solve Uz=y
657  double x;
658  for (int i = n_row - 1; i >= 0; i--)
659  {
660  x = z[i] / U_row_entry[U_column_start[i + 1] - 1].value();
661  z[i] = x;
662  for (unsigned j = U_column_start[i]; j < U_column_start[i + 1] - 1; j++)
663  {
664  z[U_row_entry[j].index()] =
665  z[U_row_entry[j].index()] - x * U_row_entry[j].value();
666  }
667  }
668 
669  // if the distribution of z was preset the redistribute to original
670  if (z_dist != 0)
671  {
672  z.redistribute(z_dist);
673  delete z_dist;
674  }
675  }
676 
677  //=============================================================================
678  /// Apply ILU(0) preconditioner for CRDoubleMatrix: Solve Ly=r then
679  /// Uz=y
680  /// and return z
681  //=============================================================================
683  const DoubleVector& r, DoubleVector& z)
684  {
685  // # of rows in the matrix
686  int n_row = r.nrow();
687 
688  // store the distribution of z
689  LinearAlgebraDistribution* z_dist = 0;
690  if (z.built())
691  {
692  z_dist = new LinearAlgebraDistribution(z.distribution_pt());
693  }
694 
695  // copy r to z
696  z = r;
697 
698  // if z is distributed then change to global
699  if (z.distributed())
700  {
701  z.redistribute(this->distribution_pt());
702  }
703 
704  // solve Ly=r (note L matrix is unit and diagonal is not stored)
705  double t;
706  for (int i = 0; i < n_row; i++)
707  {
708  t = 0;
709  for (unsigned j = L_row_start[i]; j < L_row_start[i + 1]; j++)
710  {
711  t = t + L_row_entry[j].value() * z[L_row_entry[j].index()];
712  }
713  z[i] = z[i] - t;
714  }
715 
716  // solve Uz=y
717  for (int i = n_row - 1; i >= 0; i--)
718  {
719  t = 0;
720  for (unsigned j = U_row_start[i] + 1; j < U_row_start[i + 1]; j++)
721  {
722  t = t + U_row_entry[j].value() * z[U_row_entry[j].index()];
723  }
724  z[i] = z[i] - t;
725  z[i] = z[i] / U_row_entry[U_row_start[i]].value();
726  }
727 
728  // if the distribution of z was preset the redistribute to original
729  if (z_dist != 0)
730  {
731  z.redistribute(z_dist);
732  delete z_dist;
733  }
734  }
735 } // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
//////////////////////////////////////////////////////////////// ////////////////////////////////////...
Definition: matrices.h:2791
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:2817
int * column_start()
Access to C-style column_start array.
Definition: matrices.h:2692
int * row_index()
Access to C-style row index array.
Definition: matrices.h:2704
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
int * column_index()
Access to C-style column index array.
Definition: matrices.h:1072
int * row_start()
Access to C-style row_start array.
Definition: matrices.h:1060
CRDoubleMatrix * global_matrix() const
if this matrix is distributed then a the equivalent global matrix is built using new and returned....
Definition: matrices.cc:2431
double * value()
Access to C-style value array.
Definition: matrices.h:1084
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:1002
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
bool distributed() const
distribution is serial or distributed
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
unsigned nrow() const
access function to the number of global rows.
unsigned nrow_local() const
access function for the num of local rows on this processor.
unsigned first_row() const
access function for the first row on this processor
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
double * values_pt()
access function to the underlying values
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
The contents of the vector are redistributed to match the new distribution. In a non-MPI rebuild this...
bool built() const
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to z, i.e. z=D^-1.
Vector< double > Inv_diag
Vector of inverse diagonal entries.
void setup()
Setup the preconditioner (store diagonal) from the fully assembled matrix.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to z, i.e. z=D^-1.
void setup()
Setup the preconditioner (store diagonal) from the fully assembled matrix. Problem pointer is ignored...
An OomphLibError object which should be thrown when an run-time error is encountered....
virtual DoubleMatrixBase * matrix_pt() const
Get function for matrix pointer.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
virtual const OomphCommunicator * comm_pt() const
Get function for comm pointer.
virtual void setup()=0
Setup the preconditioner. Pure virtual generic interface function.
T * value()
Access to C-style value array.
Definition: matrices.h:616
unsigned long nnz() const
Return the number of nonzero entries.
Definition: matrices.h:640
//////////////////////////////////////////////////////////////////// ////////////////////////////////...