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