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-2022 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
36namespace 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 {
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 {
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 * row_start()
Access to C-style row_start array.
Definition: matrices.h:1060
int * column_index()
Access to C-style column index array.
Definition: matrices.h:1072
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
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
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
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.
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
double * values_pt()
access function to the underlying values
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 const OomphCommunicator * comm_pt() const
Get function for comm 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 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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...