frontal_solver.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 //Interface to HSL frontal solver (fortran)
31 
32 
33 #ifdef OOMPH_HAS_MPI
34 #include "mpi.h"
35 #endif
36 
37 //oomph-lib headers
38 #include "cfortran.h"
39 #include "frontal.h"
40 #include "frontal_solver.h"
41 #include "problem.h"
42 #include "mesh.h"
43 #include "matrices.h"
44 #include "assembly_handler.h"
45 
46 namespace oomph
47 {
48 
49 
50 //====================================================================
51 /// Special solver for problems with one DOF -- HSL_MA42 can't handle
52 /// that!
53 //====================================================================
54 void HSL_MA42::solve_for_one_dof(Problem* const &problem_pt,
55  DoubleVector &result)
56 {
57 
58  //Find the number of elements
59  unsigned long n_el = problem_pt->mesh_pt()->nelement();
60 
61 #ifdef PARANOID
62  //Find the number of dofs (variables)
63  int n_dof = problem_pt->ndof();
64  // Check
65  if (n_dof!=1)
66  {
67  throw OomphLibError(
68  "You can only call this if the problem has just one dof!",
69  OOMPH_CURRENT_FUNCTION,
70  OOMPH_EXCEPTION_LOCATION);
71  }
72 #endif
73 
74  // Global jac and residuals are scalars!
75  double global_jac=0.0;
76  double global_res=0.0;
77 
78  //Locally cache pointer to assembly handler
79  AssemblyHandler* const assembly_handler_pt =
80  problem_pt->assembly_handler_pt();
81 
82 
83  // Do the assembly
84  for(unsigned long e=0;e<n_el;e++)
85  {
86  //Get pointer to the element
87  GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
88 
89  //Get number of variables in element
90  int nvar = assembly_handler_pt->ndof(elem_pt);
91 
92  // Only assemble if there's something to be assembled
93  if (nvar>0)
94  {
95  //Set up matrices and Vector for call to get_residuals
96  Vector<double> residuals(nvar);
97  DenseMatrix<double> jacobian(nvar);
98 
99  //Get the residuals and jacobian
100  assembly_handler_pt->get_jacobian(elem_pt,residuals,jacobian);
101 
102  // Add contribution
103  global_jac+=jacobian(0,0);
104  global_res+=residuals[0];
105  }
106  }
107 
108  // "Solve"
109  result[0]=global_res/global_jac;
110 
111  //If we are storing the matrix, allocate the storage
112  if(Enable_resolve)
113  {
114  //If storage has been allocated delete it
115  if(W!=0) {delete[] W;}
116  //Now allocate new storage
117  W = new double[1];
118  //Store the global jacobian
119  W[0] = global_jac;
120  //Set the number of degrees of freedom
121  N_dof=1;
122  }
123  //Set the sign of the jacobian
124  problem_pt->sign_of_jacobian() =
125  static_cast<int>(std::fabs(global_jac)/global_jac);
126 }
127 
128 //====================================================================
129 /// Wrapper for HSL MA42 frontal solver
130 //====================================================================
131 void HSL_MA42::solve(Problem* const &problem_pt, DoubleVector &result)
132 {
133 #ifdef OOMPH_HAS_MPI
134  if(problem_pt->communicator_pt()->nproc() > 1)
135  {
136  std::ostringstream error_message;
137  error_message
138  << "HSL_MA42 solver cannot be used in parallel.\n"
139  << "Please change to another linear solver.\n"
140  << "If you want to use a frontal solver then try MumpsSolver\n";
141 
142  throw OomphLibError(error_message.str(),
143  OOMPH_CURRENT_FUNCTION,
144  OOMPH_EXCEPTION_LOCATION);
145  }
146 #endif
147 
148 
149  //Find the number of elements
150  unsigned long n_el = problem_pt->mesh_pt()->nelement();
151 
152  //Find the number of dofs (variables) and store for resolves
153  N_dof = problem_pt->ndof();
154 
155  //Build the distribution .. this is a serial solver
156  LinearAlgebraDistribution dist(problem_pt->communicator_pt(),
157  N_dof,false);
158  this->build_distribution(dist);
159 
160  //Cast the number of dofs into an integer for the HSL solver
161  int n_dof = N_dof;
162 
163  // Special case: Just one variable! MA42 can't handle that
164  if (n_dof==1)
165  {
166  solve_for_one_dof(problem_pt,result);
167  return;
168  }
169 
170  // Reorder?
171  if (this->Reorder_flag)
172  {
173  reorder_elements(problem_pt);
174  }
175 
176  //Set the control flags
177  int ifsize[5];
178  double cntl[2], rinfo[2];
179 
180  //Set up memory for last and for ndf
181  int ndf, nmaxe=2, nrhs=1, lrhs=1;
182 
183  //Memory for last and dx should be allocated dynamically from the heap
184  //rather than the stack, otherwise one gets a nasty segmentation fault
185  int *last = new int[n_dof];
186  //Set up memory for dx
187  double **dx = new double *;
188  *dx = new double[n_dof];
189 
190  // Provide storage for exact values required for lenbuf array
191  int exact_lenbuf[3]={0,0,0};
192  bool exact_lenbuf_available=false;
193 
194  // Storage for recommendations for lenbuf so we can check how
195  // close our factors are
196  int lenbuf0_recommendation=0;
197  int lenbuf1_recommendation=0;
198  int lenbuf2_recommendation=0;
199 
200  // Provide storage
201  int lenbuf[3];
202 
203 
204  // Provide storage for exact values required for lenfle array
205  int exact_lenfle[3]={0,0,0};
206  bool exact_lenfle_available=false;
207 
208  // Storage for recommendation for lenfle so we can check how
209  // close our factor is
210  int lenfle0_recommendation=0;
211  int lenfle1_recommendation=0;
212  int lenfle2_recommendation=0;
213 
214  // Provide storage
215  int lenfle[3];
216 
217 
218  // Storage for recommendations for front size so we can check how
219  // close our factors are
220  double front0_recommendation=0;
221  double front1_recommendation=0;
222 
223 
224  // Flag to indicate if exact, required values for nfront are available
225  // from previous unsucessful solve
226  bool exact_nfront_available=false;
227 
228  // Storage for front sizes
229  int nfront[2];
230 
231  //Locally cache pointer to assembly handler
232  AssemblyHandler* const assembly_handler_pt =
233  problem_pt->assembly_handler_pt();
234 
235 
236  // Loop over size allocation/solver until we've made all the arrays
237  //=================================================================
238  // big enough...
239  //==============
240  bool not_done=true;
241  while (not_done)
242  {
243 
244  //Initialise frontal solver
245  //=========================
246  MA42ID(Icntl,cntl,Isave);
247 
248 
249  //Loop over the elements to setup variables associated with elements
250  //==================================================================
251  for(unsigned long e=0;e<n_el;e++)
252  {
253  //Pointer to the element
254  GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
255 
256  // MH: changed array to allocateable
257  int* ivar;
258 
259  //Get number of variables in element
260  int nvar = assembly_handler_pt->ndof(elem_pt);
261 
262 
263  // Multiple equations
264  //-------------------
265  if (nvar!=1)
266  {
267  //Copy global equation numbers into local array
268  ivar=new int[nvar];
269 
270  for(int i=0;i<nvar;i++)
271  {
272  //Need to add one to equation numbers
273  ivar[i] = 1 + assembly_handler_pt->eqn_number(elem_pt,i);
274  }
275  }
276 
277  // Just one equation
278  //------------------
279  else
280  {
281 
282  //Copy global equation numbers into local array - minimum length: 2
283  ivar=new int[2];
284 
285  // Here's the number of the only equation
286  int only_eqn = assembly_handler_pt->eqn_number(elem_pt,0);
287 
288  //Need to add one to equation numbers
289  ivar[0] = 1 + only_eqn;
290 
291  // Now add a dummy eqn either before or after the current one
292  int dummy_eqn;
293  if (only_eqn!=(n_dof-1))
294  {
295  dummy_eqn=only_eqn+1;
296  }
297  else
298  {
299  dummy_eqn=only_eqn-1;
300  }
301  ivar[1]=1+dummy_eqn;
302 
303  nvar=2;
304  }
305 
306 
307  //Now call the frontal routine
308  //GB: added check to exclude case with no dofs
309  if (nvar)
310  {
311  MA42AD(nvar,ivar,ndf,last,n_dof,Icntl,Isave,Info);
312  //ALH : throw an exception if the frontal_solver fails
313  if(Info[0] < 0) {throw NewtonSolverError(true);}
314  }
315 
316  // Cleanup
317  delete[] ivar;
318  ivar=0;
319 
320  }
321 
322 
323  //Loop over the elements again for symbolic factorisation
324  //=======================================================
325  for(unsigned long e=0;e<n_el;e++)
326  {
327  //Pointer to the element
328  GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
329 
330  // MH: changed array to allocateable
331  int* ivar;
332 
333  //Get number of variables in element
334  int nvar = assembly_handler_pt->ndof(elem_pt);
335 
336  // Multiple equations
337  //-------------------
338  if (nvar!=1)
339  {
340  //Copy global equation numbers into local array
341  ivar=new int[nvar];
342 
343  for(int i=0;i<nvar;i++)
344  {
345  //Need to add one to equation numbers
346  ivar[i] = 1 + assembly_handler_pt->eqn_number(elem_pt,i);
347  }
348  }
349 
350  // Just one equation
351  //------------------
352  else
353  {
354 
355  //Copy global equation numbers into local array - minimum length: 2
356  ivar=new int[2];
357 
358  // Here's the number of the only equation
359  int only_eqn = assembly_handler_pt->eqn_number(elem_pt,0);
360 
361  //Need to add one to equation numbers
362  ivar[0] = 1 + only_eqn;
363 
364  // Now add a dummy eqn either before or after the current one
365  int dummy_eqn;
366  if (only_eqn!=(n_dof-1))
367  {
368  dummy_eqn=only_eqn+1;
369  }
370  else
371  {
372  dummy_eqn=only_eqn-1;
373  }
374  ivar[1]=1+dummy_eqn;
375 
376  nvar=2;
377 
378  }
379 
380  //GB: added check to exclude case with no dofs
381  if (nvar)
382  {
383  MA42JD(nvar,ivar,ndf,last,nmaxe,ifsize,Icntl,Isave,Info);
384 
385  //ALH : throw an exception if the frontal_solver fails
386  if(Info[0] < 0) {throw NewtonSolverError(true);}
387  }
388 
389  // Cleanup
390  delete[] ivar;
391  ivar=0;
392 
393  } // end of symbolic factorisation
394 
395 
396 
397  // Front sizes: "Somewhat larger than..."
398  //---------------------------------------
399  front0_recommendation=ifsize[0];
400  front1_recommendation=ifsize[1];
401  if (!exact_nfront_available)
402  {
403  nfront[0]=int(ceil(Front_factor*double(front0_recommendation)));
404  nfront[1]=int(ceil(Front_factor*double(front1_recommendation)));
405 
406  if(n_dof < nfront[0]) nfront[0] = n_dof;
407  if(n_dof < nfront[1]) nfront[1] = n_dof;
408  }
409  //We have a problem if the front size is cocked
410 
411 
412  // Sizes for lenbuf[]: "Somewhat larger than..."
413  //----------------------------------------------
414  lenbuf0_recommendation=ifsize[2] + ndf;
415  //If we are storing the matrix,
416  //we need to allocate storage for lenbuf_1
417  if(Enable_resolve) {lenbuf1_recommendation=ifsize[3];}
418  //Otherwise set it to zero
419  else {lenbuf1_recommendation=0;}
420  lenbuf2_recommendation=ifsize[4];
421 
422 
423  //Enable use of direct access files?
424  //----------------------------------
426  {
427  if (Doc_stats)
428  {
429  oomph_info << "Using direct access files" << std::endl;
430  }
431 
432  // Unit numbers for direct access files (middle one only used for
433  // re-solve; set to zero as this is not used here).
434  int istrm[3];
435  istrm[0]=17;
436  //If we are storing the matrix, set the stream
437  if(Enable_resolve) {istrm[1]=19;}
438  //otherwise, set it to zero
439  else {istrm[1] = 0;}
440  istrm[2]=18;
441 
442  //Set up the memory: Need memory "somewhat larger" than ...
443  double factor=1.1;
444  lenbuf[0] = int(ceil(factor*double(10*(nfront[0] + 1))));
445  //If we are storing the matrix, set the value
446  if(Enable_resolve)
447  {lenbuf[1] = int(ceil(factor*double(10*(nfront[0]))));}
448  //Otherwise, set to zero
449  else {lenbuf[1] = 0;}
450  lenbuf[2] = int(ceil(factor*double(10*(2*nfront[0]+5))));
451 
452 
453  //Initial size allocation: Need memory "somewhat larger" than ...
454  if (exact_lenfle_available)
455  {
456  lenfle[0] = exact_lenfle[0];
457  lenfle[1] = exact_lenfle[1];
458  lenfle[2] = exact_lenfle[2];
459  }
460  else
461  {
462  lenfle0_recommendation=ifsize[2]+ndf;
463  lenfle1_recommendation=ifsize[3];
464  lenfle2_recommendation=ifsize[4];
465  lenfle[0]=int(ceil(Lenfle_factor*double(lenfle0_recommendation)));
466  lenfle[1]=int(ceil(Lenfle_factor*double(lenfle1_recommendation)));
467  lenfle[2]=int(ceil(Lenfle_factor*double(lenfle2_recommendation)));
468  }
469 
470  //Setup direct access files
471  MA42PD(istrm,lenbuf,lenfle,Icntl,Isave,Info);
472  }
473  else
474  {
475  if (Doc_stats)
476  {
477  oomph_info << "Not using direct access files" << std::endl;
478  }
479 
480  //Initial size allocation: Need memory "somewhat larger" than ...
481  if (exact_lenbuf_available)
482  {
483  lenbuf[0] = exact_lenbuf[0];
484  lenbuf[1] = exact_lenbuf[1];
485  lenbuf[2] = exact_lenbuf[2];
486  }
487  else
488  {
489  lenbuf[0] = int(ceil(Lenbuf_factor0*double(lenbuf0_recommendation)));
490  //If we are storing the matrix, set the buffer size
491  if(Enable_resolve)
492  {lenbuf[1] = int(ceil(Lenbuf_factor1*double(lenbuf1_recommendation)));}
493  //Otherwise its zero
494  else {lenbuf[1] = 0;}
495  lenbuf[2] = int(ceil(Lenbuf_factor2*double(lenbuf2_recommendation)));
496  }
497  }
498 
499 
500 
501  if (Doc_stats)
502  {
503  oomph_info << "\n FRONT SIZE (min and actual): "
504  << ifsize[0] << " " << nfront[0] << std::endl;
505  }
506 
507 
508  // Initialise success
509  bool success=true;
510 
511  // Workspace arrays: calculate amount in local variables initiall
512  int lw = 1 + lenbuf[0] + lenbuf[1] + nfront[0]*nfront[1];
513  if(lrhs*nfront[0] > nrhs*nfront[1])
514  {
515  lw += lrhs*nfront[0];
516  }
517  else
518  {
519  lw += nrhs*nfront[1];
520  }
521 
522  int liw = lenbuf[2] + 2*nfront[0] + 4*nfront[1];
523 
524  //Set up memory for w and iw
525  //Again allocate dynamically from the heap, rather than the stack
526  //If the required amount of storage has changed, reallocate
527  //and store the values in the object member data
528  if(lw != Lw)
529  {
530  //Set the new value of Lw (member data)
531  Lw = lw;
532  //Delete the previously allocated storage
533  if(W) delete[] W;
534  //Now allocate new storage
535  W = new (std::nothrow) double [Lw];
536  if (!W)
537  {
538  throw OomphLibError("Out of memory in allocation of w\n",
539  OOMPH_CURRENT_FUNCTION,
540  OOMPH_EXCEPTION_LOCATION);
541  }
542  }
543 
544  //If the required amount of storage has changed, reallocate
545  if(liw != Liw)
546  {
547  //Set the new value of Liw (member data)
548  Liw = liw;
549  //Delete the previously allocated storgae
550  if(IW!=0) {delete[] IW;}
551  //Now allocate new storage
552  IW = new (std::nothrow) int [Liw];
553  if (!IW)
554  {
555  throw OomphLibError("Out of memory in allocation of iw\n",
556  OOMPH_CURRENT_FUNCTION,
557  OOMPH_EXCEPTION_LOCATION);
558  }
559  }
560 
561  //Doc
562  if (Doc_stats)
563  {
564  double temp = (lenbuf[0]*sizeof(double) + lenbuf[2]*sizeof(int))/
565  (1024.0*1024.0);
566  oomph_info << "\n ESTIMATED MEMORY USAGE " << temp << "Mb" << std::endl;
567  }
568 
569 
570  //Now do the actual frontal assembly and solve loop
571  //=================================================
572  for(unsigned long e=0;e<n_el;e++)
573  {
574  //Get pointer to the element
575  GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
576 
577  //Get number of variables in element
578  int nvar = assembly_handler_pt->ndof(elem_pt);
579 
580  // MH: changed array to allocatable
581  int* ivar;
582 
583  // Multiple equations
584  //-------------------
585  if (nvar!=1)
586  {
587  //Copy global equation numbers into local array
588  ivar=new int[nvar];
589 
590  for(int i=0;i<nvar;i++)
591  {
592  //Need to add one to equation numbers
593  ivar[i] = 1 + assembly_handler_pt->eqn_number(elem_pt,i);
594  }
595  nmaxe = nvar;
596  }
597 
598  // Just one equation
599  //------------------
600  else
601  {
602 
603  //Copy global equation numbers into local array - minimum length: 2
604  ivar=new int[2];
605 
606  // Here's the number of the only equation
607  int only_eqn = assembly_handler_pt->eqn_number(elem_pt,0);
608 
609  //Need to add one to equation numbers
610  ivar[0] = 1 + only_eqn;
611 
612  // Now add a dummy eqn either before or after the current one
613  int dummy_eqn;
614  if (only_eqn!=(n_dof-1))
615  {
616  dummy_eqn=only_eqn+1;
617  }
618  else
619  {
620  dummy_eqn=only_eqn-1;
621  }
622  ivar[1]=1+dummy_eqn;
623 
624  nmaxe = 2;
625  }
626 
627  //Set up matrices and Vector for call to get_residuals
628  Vector<double> residuals(nvar);
629  DenseMatrix<double> jacobian(nvar);
630 
631  //Get the residuals and jacobian
632  assembly_handler_pt->get_jacobian(elem_pt,residuals,jacobian);
633 
634  // Add dummy rows and columns if required
635  if (nvar==1)
636  {
637  double onlyjac=jacobian(0,0);
638  jacobian.resize(2,2);
639  jacobian(0,0)=onlyjac;
640  jacobian(1,0)=0.0;
641  jacobian(0,1)=0.0;
642  jacobian(1,1)=0.0;
643 
644  double onlyres=residuals[0];
645  residuals.resize(2);
646  residuals[0]=onlyres;
647  residuals[1]=0.0;
648 
649  nvar=2;
650  }
651 
652 
653  //GB: added check to exclude case with no dofs
654  if (nvar)
655  {
656  //Now translate the residuals and jacobian into form to be
657  //taken by stupid FORTRAN frontal solver
658  //double avar[nvar][nmaxe], rhs[1][nmaxe];
659  //Assign memory on the heap rather than the stack because the ndofs
660  //could get too large
661  double **avar = new double*[nvar];
662  double *tmp = new double[nvar*nmaxe];
663  for(int i=0;i<nvar;i++) {avar[i] = &tmp[i*nmaxe];}
664  double **rhs = new double*[1];
665  rhs[0] = new double[nmaxe];
666 
667 
668  for(int i=0;i<nmaxe;i++)
669  {
670  rhs[0][i] = residuals[i];
671  for(int j=0;j<nvar;j++)
672  {
673  //Note need to transpose here
674  avar[j][i] = jacobian(i,j);
675  }
676  }
677 
678  //Call the frontal solver
679  MA42BD(nvar,ivar,ndf,last,nmaxe,avar,nrhs,rhs,lrhs,n_dof,dx,nfront,
680  lenbuf,Lw,W,Liw,IW,Icntl,cntl,Isave,Info,rinfo);
681 
682 
683  // Error check:
684  if (Info[0] < 0)
685  {
686 
687  if (Doc_stats) oomph_info << "Error: Info[0]=" << Info[0] << std::endl;
688 
689  // Solve isn't going to be successful
690  success=false;
691 
692  // Error: Entries in lenbuf too small -- this is covered
693  if (Info[0]==-16)
694  {
695  if (Doc_stats) oomph_info << "lenbuf[] too small -- can recover..."
696  << std::endl;
697  }
698  else if (Info[0]==-12)
699  {
700  if (Doc_stats) oomph_info << "nfront[] too small -- can recover..."
701  << std::endl;
702  }
703  else if (Info[0]==-17)
704  {
705  if (Doc_stats) oomph_info << "lenfle[] too small -- can recover..."
706  << std::endl;
707  }
708  else
709  {
710  if (Doc_stats)
711  {
712  oomph_info<<"Can't recover from this error"<<std::endl;
713  }
714  throw NewtonSolverError(true);
715  }
716  }
717 
718  //Clean up the memory
719  delete[] rhs[0]; rhs[0] = 0;
720  delete[] rhs; rhs =0;
721  delete[] avar[0]; avar[0] = 0;
722  delete[] avar;
723  }
724 
725  // Cleanup
726  delete[] ivar;
727  ivar=0;
728 
729  } // end of loop over elements for assemble/solve
730 
731  //Set the sign of the jacobian matrix
732  problem_pt->sign_of_jacobian() = Info[1];
733 
734  //If we are not storing the matrix, then delete the assigned workspace
735  if(!Enable_resolve)
736  {
737  //Free the workspace memory assigned and set stored values to zero
738  delete[] IW; IW=0; Liw=0;
739  delete[] W; W=0; Lw=0;
740  //Reset the number of dofs
741  N_dof=0;
742  }
743 
744  // We've done the elements -- did we encounter an error
745  // that forces us to re-allocate workspace?
746  if (success)
747  {
748  // We're done!
749  not_done=false;
750  }
751  else
752  {
753  //Reset sizes for lenbuf
755  {
756  exact_lenbuf[0] = Info[4];
757  exact_lenbuf[1] = Info[5];
758  exact_lenbuf[2] = Info[6];
759  exact_lenbuf_available=true;
760  }
761 
762  // nfront[] has already been overwritten with required values -- don't
763  // need to update anything. Just indicate that new values are
764  // available so they don't have to be re-assigned.
765  exact_nfront_available=true;
766 
767  //Reset sizes for lenfle
768  exact_lenfle[0] = lenbuf[0]*Info[9];
769  exact_lenfle[1] = lenbuf[1]*Info[10];
770  exact_lenfle[2] = lenbuf[2]*Info[11];
771  exact_lenfle_available=true;
772 
773  }
774 
775  } // end of loop over shuffling of workspace array sizes until it all
776  // fits...
777 
778  result.build(&dist);
779  //Now load the values back into result
780  for(int i=0;i<n_dof;i++) result[i] = dx[0][i];
781 
782  //Print the actual memory used
783  if (Doc_stats)
784  {
786  {
787  double temp = (Info[4]*sizeof(double) + Info[6]*sizeof(int))/
788  (1024.0*1024.0);
789  oomph_info << " TOTAL MEMORY USED " << temp << "Mb" << std::endl;
790  }
791  oomph_info << std::endl;
792  oomph_info << "lenbuf[]: " << lenbuf[0] << " "
793  << lenbuf[1] << " " << lenbuf[2] << " " << std::endl;
794  oomph_info << "lenbuf[] factors required and initially specified:"
795  << std::endl;
796  oomph_info << "lenbuf[0]: " << Info[4]/double(lenbuf0_recommendation) << " "
797  << Lenbuf_factor0 << std::endl;
798  oomph_info << "lenbuf[1]: " << Info[5]/double(lenbuf1_recommendation)
799  << " "
800  << Lenbuf_factor1 << std::endl;
801 
802  oomph_info << "lenbuf[2]: " << Info[6]/double(lenbuf2_recommendation) << " "
803  << Lenbuf_factor2 << std::endl;
804  oomph_info << "nfront[] factors required and initially specified:"
805  << std::endl;
806  oomph_info << "nfront[0]: " << nfront[0]/double(front0_recommendation)
807  << " " << Front_factor << std::endl;
808  oomph_info << "nfront[1]: " << nfront[1]/double(front1_recommendation)
809  << " " << Front_factor << std::endl;
811  {
812  oomph_info << "lenfle[]: " << lenfle[0] << " "
813  << lenfle[1] << " " << lenfle[2] << " " << std::endl;
814  oomph_info << "lenfle[] factors required and initially specified:"
815  << std::endl;
816  oomph_info << "lenfle[0]: "
817  << Info[9]*lenbuf[0]/double(lenfle0_recommendation)
818  << " " << Lenfle_factor << std::endl;
819  oomph_info << "lenfle[1]: "
820  << Info[10]*lenbuf[1]/double(lenfle1_recommendation)
821  << " " << Lenfle_factor << std::endl;
822  oomph_info << "lenfle[2]: "
823  << Info[11]*lenbuf[2]/double(lenfle2_recommendation)
824  << " " << Lenfle_factor << std::endl;
825  }
826  }
827 
828  //Free the memory assigned
829  delete[] *dx;
830  delete dx;
831  delete[] last;
832 }
833 
834 //====================================================================
835 /// Wrapper for HSL MA42 frontal solver
836 //====================================================================
837 void HSL_MA42::resolve(const DoubleVector &rhs, DoubleVector &result)
838 {
839  //If we haven't stored the factors complain
840  if(W==0)
841  {
842  throw OomphLibError("Resolve called without first solving",
843  OOMPH_CURRENT_FUNCTION,
844  OOMPH_EXCEPTION_LOCATION);
845  }
846 
847  //Find the number of dofs (variables)
848  int n_dof = N_dof;
849 
850  //Check that the number of DOFs is equal to the size of the RHS vector
851 #ifdef PARANOID
852  if(n_dof != static_cast<int>(rhs.nrow()))
853  {
854  throw OomphLibError(
855  "RHS does not have the same dimension as the linear system",
856  OOMPH_CURRENT_FUNCTION,
857  OOMPH_EXCEPTION_LOCATION);
858  }
859 #endif
860 
861  //Special case: just one variable! MA42 can't handle it
862  //Solve using the stored jacobian (single double)
863  if(n_dof==1) {result[0]= rhs[0]/W[0]; return;}
864 
865  //Otherwise actually call the MA42 routine
866  //We are solving the original matrix (not the transpose)
867  int trans=0;
868  //There is only one rhs
869  int nrhs=1;
870  //The size of the vectors is the number of degrees of freedom in the problem
871  int lx = n_dof;
872 
873  //Allocate storage for the RHS
874  double **b = new double*;
875  *b = new double[n_dof];
876  //Load in the RHS vector
877  for(int i=0;i<n_dof;i++) {b[0][i] = rhs[i];}
878 
879  //Storage for the results
880  double **x = new double*;
881  *x = new double[n_dof];
882 
883  //Now call the resolver
884  MA42CD(trans,nrhs,lx,b,x,Lw,W,Liw,IW,Icntl,Isave,Info);
885 
886  //If there has been an error throw it
887  if(Info[0] < 0) {throw NewtonSolverError(true);}
888 
889  result.build(this->distribution_pt());
890  //Put the answer into the result array
891  for(int i=0;i<n_dof;++i) {result[i] = x[0][i];}
892 
893  //Delete the allocated storage
894  delete[] *x;
895  delete x;
896  delete[] *b;
897  delete b;
898 }
899 
900 
901 //=========================================================================
902 ///Function to reorder the elements according to Sloan's algorithm
903 //=========================================================================
904 void HSL_MA42::reorder_elements(Problem* const &problem_pt)
905 {
906  //Find the number of elements
907  unsigned long n_el = problem_pt->mesh_pt()->nelement();
908 
909  //Find the number of dofs
910  int n_dof = problem_pt->ndof();
911 
912  Vector<int> order(n_el);
913 
914  //Control parameters
915  int icntl[10], info[15], direct, nsup, vars[n_dof], svar[n_dof];
916  int liw, lw, *perm=0;
917  double wt[3], rinfo[6];
918 
919 
920  // Direct ordering: 1
921  direct=1;
922 
923  // Initial guess for workspaces (deliberately too small so the
924  // routine fails and returns the required sizes
925  liw=1;
926  lw=1;
927 
928  //Initialise things here
929  MC63ID(icntl);
930 
931 
932  // Suppress printing of error and warning messages?
933  if (!Doc_stats)
934  {
935  icntl[0]=-1;
936  icntl[1]=-1;
937  }
938 
939  //Set up iw and w
940  int *iw = new int [liw];
941  double *w = new double [lw];
942 
943  //Set up the vectors that hold the element connectivity information
944  //Can initialise the first entry of eltptr to 1
945  Vector<int> eltvar, eltptr(1,1);
946 
947  //Locally cache pointer to assembly handler
948  AssemblyHandler* const assembly_handler_pt =
949  problem_pt->assembly_handler_pt();
950 
951  //Now loop over all the elements and assemble eltvar and eltptr
952  for(unsigned long e=0;e<n_el;e++)
953  {
954  //Set up the default order
955  order[e] = e;
956 
957  //Get pointer to the element
958  GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
959 
960  //Get the number of variables in the element via the assembly handler
961  int nvar = assembly_handler_pt->ndof(elem_pt);
962 
963  // Multiple equations
964  //-------------------
965  if (nvar!=1)
966  {
967  //Copy the equation numbers into the global array
968  for(int i=0;i<nvar;i++)
969  {
970  //Need to add one to equation numbers
971  eltvar.push_back(1 + assembly_handler_pt->eqn_number(elem_pt,i));
972  }
973  eltptr.push_back(eltptr.back()+nvar);
974  }
975  // Just one equation
976  //------------------
977  else
978  {
979  // Here's the number of the only equation
980  int only_eqn=assembly_handler_pt->eqn_number(elem_pt,0);
981 
982  //Need to add one to equation numbers
983  eltvar.push_back(1 + only_eqn);
984 
985  // Now add a dummy eqn either before or after the current one
986  int dummy_eqn;
987  if (only_eqn!=(n_dof-1)) {dummy_eqn=only_eqn+1;}
988  else {dummy_eqn=only_eqn-1;}
989 
990  eltvar.push_back(1+dummy_eqn);
991 
992  eltptr.push_back(eltptr.back()+2);
993  }
994 
995  } //End of loop over the elements
996 
997  //Set the value of n_e (number of entries in the element variable list
998  int n_e = eltvar.size();
999 
1000  do
1001  {
1002  //Call the reordering routine
1003  MC63AD(direct,n_dof,n_el,n_e,&eltvar[0],&eltptr[0],&order[0],
1004  perm,nsup,vars,svar,wt,liw,iw,lw,w,icntl,info,rinfo);
1005 
1006  //If not enough space in iw or w, allocate enough and try again
1007  if(info[0] == -4)
1008  {
1009  if (Doc_stats) oomph_info << " Reallocating liw to "
1010  << info[4] << std::endl;
1011  delete[] iw;
1012  liw = info[4];
1013  iw = new int[liw];
1014  }
1015 
1016  //If not enough space in w, allocate enough and try again
1017  if(info[0] == -2)
1018  {
1019  if (Doc_stats) oomph_info << " Reallocating lw to "
1020  << info[5] << std::endl;
1021  delete[] w;
1022  lw = info[5];
1023  w = new double [lw];
1024  }
1025 
1026  } while(info[0] < 0);
1027 
1028 
1029  if (Doc_stats)
1030  {
1031  oomph_info <<"\n Initial wavefront details :\n max "
1032  << rinfo[0] << "\tmean "
1033  << rinfo[1] << "\tprofile " << rinfo[2];
1034  oomph_info <<"\n Reordered wavefront details:\n max "
1035  << rinfo[3] << "\tmean "
1036  << rinfo[4] << "\tprofile " << rinfo[5];
1037  oomph_info << std::endl;
1038  }
1039 
1040  //Free the memory allocated
1041  delete[] w;
1042  delete[] iw;
1043 
1044 
1045 
1046  // Store re-ordered elements in temporary vector
1047  Vector<GeneralisedElement*> tmp_element_pt(n_el);
1048  for (unsigned e=0;e<n_el;e++)
1049  {
1050  tmp_element_pt[e]=problem_pt->mesh_pt()->element_pt(order[e]-1);
1051  }
1052  for (unsigned e=0;e<n_el;e++)
1053  {
1054  problem_pt->mesh_pt()->element_pt(e)=tmp_element_pt[e];
1055  }
1056 
1057 }
1058 
1059 }
A Generalised Element class.
Definition: elements.h:76
A class to handle errors in the Newton solver.
Definition: problem.h:2820
double Lenbuf_factor2
Factor to increase storage for lenbuf[2]; see MA42 documentation for details.
int Icntl[8]
Control flag for MA42; see MA42 documentation for details.
bool Reorder_flag
Reorder elements with Sloan's algorithm?
bool Doc_stats
Doc the solver stats or stay quiet?
bool Use_direct_access_files
Use direct access files?
cstr elem_len * i
Definition: cfortran.h:607
unsigned long N_dof
Size of the linear system.
int Liw
Size of the integer workspace array.
unsigned nrow() const
access function to the number of global rows.
bool Enable_resolve
Boolean that indicates whether the matrix (or its factors, in the case of direct solver) should be st...
Definition: linear_solver.h:80
The Problem class.
Definition: problem.h:152
double Lenbuf_factor0
Factor to increase storage for lenbuf[0]; see MA42 documentation for details.
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1264
OomphInfo oomph_info
e
Definition: cfortran.h:575
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1221
double Lenbuf_factor1
Factor to increase storage for lenbuf[1]; see MA42 documentation for details.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
AssemblyHandler *& assembly_handler_pt()
Return a pointer to the assembly handler object.
Definition: problem.h:1502
virtual unsigned ndof(GeneralisedElement *const &elem_pt)
Return the number of degrees of freedom in the element elem_pt.
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.
int Info[23]
Control flag for MA42; see MA42 documentation for details.
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
int & sign_of_jacobian()
Access function for the sign of the global jacobian matrix. This will be set by the linear solver...
Definition: problem.h:2301
int * IW
Integer workspace storage for MA42.
void reorder_elements(Problem *const &problem_pt)
Function to reorder the elements based on Sloan's algorithm.
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1574
double * W
Workspace storage for MA42.
double Front_factor
Factor to increase storage for front size; see MA42 documentation for details.
A class that is used to define the functions used to assemble the elemental contributions to the resi...
void resolve(const DoubleVector &rhs, DoubleVector &result)
Return the solution to the linear system Ax = result, where A is the most recently factorised jacobia...
virtual void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable" for elem_pt.
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
virtual unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Return the global equation number of the local unknown ieqn_local in elem_pt.
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results Vector which contains the solution of the li...
double Lenfle_factor
Factor to increase size of direct access files; see MA42 documentation for details.
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 Isave[45]
Control flag for MA42; see MA42 documentation for details.
void solve_for_one_dof(Problem *const &problem_pt, DoubleVector &result)
Special solver for problems with 1 dof (MA42 can't handle this case so solve() forwards the "solve" t...
int Lw
Size of the workspace array, W.
void resize(const unsigned long &n)
Definition: matrices.h:511