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