G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsAPALM.hpp
Go to the documentation of this file.
1
14#pragma once
15
16namespace gismo
17{
18
19template <class T>
20#ifdef GISMO_WITH_MPI
21gsAPALM<T>::gsAPALM( gsALMBase<T> * ALM,
22 const gsAPALMData<T,solution_t> & Data,
23 const gsMpiComm & comm )
24#else
25gsAPALM<T>::gsAPALM( gsALMBase<T> * ALM,
26 const gsAPALMData<T,solution_t> & Data,
27 const gsSerialComm & comm )
28#endif
29:
30m_ALM(ALM),
31m_dataEmpty(Data),
32m_comm(comm)
33{
34 GISMO_ASSERT(m_dataEmpty.empty(),"gsAPALMData must be empty; it will be used to define the options!");
35 this->_defaultOptions();
36 this->_getOptions();
37
38 //Get size and rank of the processor
39 m_proc_count = m_comm.size();
40 m_rank = m_comm.rank();
41}
42
43template <class T>
44gsAPALM<T>::gsAPALM( gsALMBase<T> * ALM,
45 const gsAPALMData<T,solution_t> & Data)
46:
47m_ALM(ALM),
48m_dataEmpty(Data),
49#ifdef GISMO_WITH_MPI
50 m_comm_dummy(new gsMpiComm()),
51#else
52 m_comm_dummy(new gsSerialComm()),
53#endif
54m_comm(*m_comm_dummy)
55{
56 GISMO_ASSERT(m_dataEmpty.empty(),"gsAPALMData must be empty; it will be used to define the options!");
57 this->_defaultOptions();
58 this->_getOptions();
59
60 m_proc_count = 1;
61 m_rank = 0;
62}
63
64template <class T>
65void gsAPALM<T>::_defaultOptions()
66{
67 m_options.addInt("MaxIt","Maximum number of steps",1000);
68 m_options.addSwitch("Verbose","Verbosity",false);
69 m_options.addInt("SubIntervals","Number of subintervals",2);
70 m_options.addSwitch("SingularPoint","Enable singular point detection",false);
71 m_options.addReal("BranchLengthMultiplier","Multiplier for the length of branches after bifurcation",1);
72 m_options.addReal("BifLengthMultiplier","Multiplier for the length of the first interval after detection of a bifurcation point",0.1);
73}
74
75template <class T>
76void gsAPALM<T>::_getOptions()
77{
78 m_maxIterations = m_options.getInt("MaxIt");
79 m_verbose = m_options.getSwitch("Verbose");
80 m_subIntervals = m_options.getInt("SubIntervals");
81 m_singularPoint = m_options.getSwitch("SingularPoint");
82 m_branchLengthMult = m_options.getReal("BranchLengthMultiplier");
83 m_bifLengthMult = m_options.getReal("BifLengthMultiplier");
84}
85
86template <class T>
87void gsAPALM<T>::initialize()
88{
89 this->_getOptions();
90
91 // Initialize solutions
92 if (m_starts.size()==0)
93 {
94 T L0;
95 gsVector<T> U0(m_ALM->numDofs());
96 U0.setZero();
97 L0 = 0.0;
98 this->_initStart(U0,L0,m_ALM->getLength());
99 }
100}
101
102template <class T>
103void gsAPALM<T>::_initStart(const gsVector<T> & Ustart, const T & Lstart, const T & dL)
104{
105 // Initialize solutions
106 std::vector<gsVector<T>> Ustarts; Ustarts.push_back(Ustart);
107 std::vector<T> Lstarts; Lstarts.push_back(Lstart);
108 std::vector<T> dLs; dLs .push_back(dL);
109 this->_initStart(Ustarts,Lstarts,dLs);
110}
111
112template <class T>
113void gsAPALM<T>::_initStart(const std::vector<gsVector<T>> & Ustarts, const std::vector<T> & Lstarts, const std::vector<T> & dLs)
114{
115 GISMO_ASSERT(Ustarts.size()==Lstarts.size(),"Size mismatch in start solutions");
116 GISMO_ASSERT(Ustarts.size()==dLs.size(),"Size mismatch in start solutions");
117
118 for (size_t k=0; k!=Ustarts.size(); k++)
119 m_starts.push(std::make_tuple(std::make_pair(Ustarts[k],Lstarts[k]),dLs[k],false));
120}
121
122template <class T>
123void gsAPALM<T>::serialSolve(index_t Nsteps)
124{
125#ifdef GISMO_WITH_MPI
126 if (m_rank!=0)
127 m_comm.barrier();
128 else
129 {
130#endif
131 GISMO_ASSERT(m_starts.size()>0,"No start point is created. Call initialize first?");
132
133 T Lold, Lprev, L0;
134 gsMatrix<T> Uold, Uprev, U0;
135 T dL, dL0;
136 bool bisected, finished;
137 std::vector<solution_t> solutions;
138 std::vector<T> times;
139 std::vector<index_t> levels;
140 while (!m_starts.empty())
141 {
142 finished = false;
143 solutions.clear();
144 times.clear();
145 levels.clear();
146 // Initialize solutions
147 U0 = Uold = Uprev = std::get<0>(m_starts.front()).first;
148 L0 = Lold = Lprev = std::get<0>(m_starts.front()).second;
149 dL = dL0 = std::get<1>(m_starts.front());
150 bisected = std::get<2>(m_starts.front());
151 m_starts.pop();
152 T s = 0;
153
154 // If start point comes after a bisection, then multiply it by the corresponding multiplier.
155 if (bisected)
156 dL *= m_bifLengthMult;
157 else
158 {
159 // Store initial solution
160 solutions.push_back({U0,L0});
161 times.push_back(s);
162 levels.push_back(0);
163 }
164
165 index_t k=1;
166 std::vector<solution_t> tmpSolutions;
167 while (k<Nsteps && !finished)
168 {
169 if (m_verbose) gsMPIInfo(m_rank)<<"Load step "<< k<<"\t"<<"dL = "<<m_ALM->getLength()<<"; curve time = "<<s<<"\n";
170 solution_t start = std::make_pair(Uold,Lold);
171 solution_t prev = std::make_pair(Uprev,Lprev);
172
173 std::tuple<index_t, T, solution_t, solution_t> dataEntry = std::make_tuple(k,dL,start,prev);
174
175 this->_initiation(dataEntry,s,dL,tmpSolutions,finished);
176
177 // Update curve length
178 s += dL;
179
180 // Update solutions
181 solutions.push_back(tmpSolutions[0]);
182 times.push_back(s);
183 levels.push_back(0);
184
185 // Update previous and old solution
186 Uprev = Uold;
187 Lprev = Lold;
188 Uold = tmpSolutions[0].first;
189 Lold = tmpSolutions[0].second;
190
191 // Store new branch as start point
192 if (finished)
193 m_starts.push(std::make_tuple(tmpSolutions[1],dL0*m_branchLengthMult,true));
194
195 dL = dL0; // set the length to the length that we want after the first iteration
196
197 // Update step counter
198 k++;
199 } // end of steps
200
201 // Store data of the branch
202 gsAPALMData<T,solution_t> data = m_dataEmpty; // create new data set
203 data.setData(times,solutions); // initialize the dataset with the newly computed data
204 data.init(); // initialize the dataset
205 m_data.add(data); // add the dataset to the tree
206
207 // Store the solutions, times and levels
208 m_solutions.push_back(solutions);
209 m_times.push_back(times);
210 m_levels.push_back(levels);
211
212 } // end of start points
213
214#ifdef GISMO_WITH_MPI
215 m_comm.barrier();
216 }
217#endif
218}
219
220
221// NOTE: This does not make new branches!
222template <class T>
223void gsAPALM<T>::parallelSolve()
224{
225#ifdef GISMO_WITH_MPI
226 if (m_comm.size()==1)
227 this->parallelSolve_impl<false>();
228 else
229 this->parallelSolve_impl<true>();
230#else
231 this->parallelSolve_impl<false>();
232#endif
233}
234
235template <class T>
236void gsAPALM<T>::solve(index_t Nsteps)
237{
238#ifdef GISMO_WITH_MPI
239 if (m_comm.size()==1)
240 this->_solve_impl<false>(Nsteps);
241 else
242 this->_solve_impl<true>(Nsteps);
243#else
244 this->_solve_impl<false>(Nsteps);
245#endif
246}
247
248template <class T>
249template <bool _hasWorkers>
250typename std::enable_if< _hasWorkers, void>::type
251gsAPALM<T>::parallelSolve_impl()
252{
253#ifdef GISMO_WITH_MPI
254 GISMO_ASSERT(m_comm.size()>1,"Something went wrong. This implementation only works when >1 processes are availbale, but nprocesses = "<<m_comm.size());
255
256 bool stop = false; // stop signal
257 solution_t reference;
258 index_t ID;
259 index_t it = 0;
260 index_t branch;
261
262 std::vector<T> distances;
263 std::vector<solution_t> stepSolutions;
264 T lowerDistance, upperDistance;
265 std::pair<T,T> dataInterval;
266
267 index_t dataLevel;
268
269 std::tuple<index_t, T , solution_t, solution_t> dataEntry;
270
271 //----------------------------------------------------------------------------------------
272 //MAIN------------------------------------------------------------------------------------
273 //----------------------------------------------------------------------------------------
274 if (m_rank==0)
275 {
276 index_t njobs = 0;
277
278 gsMPIInfo(m_rank)<<"Adding workers...\n";
279 for (index_t w = 1; w!=m_proc_count; w++)
280 m_workers.push(w);
281
282 while (!m_data.empty() && it < m_maxIterations && !m_workers.empty())
283 {
284 gsMPIInfo(m_rank)<<"Available workers:\t";
285 std::queue<index_t> workers_copy = m_workers;
286 while (!workers_copy.empty())
287 {
288 gsMPIInfo(m_rank)<<workers_copy.front()<<"\t";
289 workers_copy.pop();
290 }
291 gsMPIInfo(m_rank)<<"\n";
292
293
294 branch = m_data.getFirstNonEmptyBranch();
295 gsMPIInfo(m_rank)<<"There are "<<m_data.branch(branch).nActive()<<" active jobs and "<<m_data.branch(branch).nWaiting()<<" jobs in the queue of branch "<<branch<<"\n";
296
297 dataEntry = m_data.branch(branch).pop();
298 ID = std::get<0>(dataEntry);
299 dataLevel = m_data.branch(branch).jobLevel(ID);
300 bool success = m_data.branch(branch).getReferenceByID(ID,reference);
301 GISMO_ENSURE(success,"Reference not found");
302
303 this->_sendMainToWorker(m_workers.front(),stop);
304 this->_sendMainToWorker(m_workers.front(),
305 branch,
306 ID,
307 dataLevel);
308 this->_sendMainToWorker(m_workers.front(),
309 dataEntry,
310 m_data.branch(branch).jobTimes(ID),
311 reference
312 );
313
314 m_workers.pop();
315 it++;
316 njobs++;
317 }
318
319 while (njobs > 0)
320 {
321 gsMPIInfo(m_rank)<<njobs<<" job(s) running\n";
322
323 index_t source = MPI_ANY_SOURCE;
324 this->_recvWorkerToMain(source,
325 branch,
326 ID);
327 this->_recvWorkerToMain(source,
328 distances,
329 stepSolutions,
330 upperDistance,
331 lowerDistance);
332
333 m_data.branch(branch).submit(ID,distances,stepSolutions,upperDistance,lowerDistance);
334 m_data.branch(branch).finishJob(ID);
335
336 // Remove job
337 njobs--;
338 // Add worker to pool
339 m_workers.push(source);
340
341 gsMPIInfo(m_rank)<<"Available workers:\t";
342 std::queue<index_t> workers_copy = m_workers;
343 while (!workers_copy.empty())
344 {
345 gsInfo<<workers_copy.front()<<"\t";
346 workers_copy.pop();
347 }
348 gsInfo<<"\n";
349
350 // As long as the queue is not empty, the iterations are lower than the max number and the pool of workers is not empty
351 while (!m_data.empty() && it < m_maxIterations && !m_workers.empty())
352 {
353 // SAME WHILE LOOP AS ABOVE!!!!!!!!
354 branch = m_data.getFirstNonEmptyBranch();
355 gsMPIInfo(m_rank)<<"There are "<<m_data.branch(branch).nActive()<<" active jobs and "<<m_data.branch(branch).nWaiting()<<" jobs in the queue of branch "<<branch<<"\n";
356
357 dataEntry = m_data.branch(branch).pop();
358 ID = std::get<0>(dataEntry);
359 dataLevel = m_data.branch(branch).jobLevel(ID);
360 bool success = m_data.branch(branch).getReferenceByID(ID,reference);
361 GISMO_ENSURE(success,"Reference not found");
362 this->_sendMainToWorker(m_workers.front(),stop);
363 this->_sendMainToWorker(m_workers.front(),
364 branch,
365 ID,
366 dataLevel);
367 this->_sendMainToWorker(m_workers.front(),
368 dataEntry,
369 m_data.branch(branch).jobTimes(ID),
370 reference
371 );
372 m_workers.pop();
373
374 it++;
375 njobs++;
376 }
377 }
378 stop = true;
379 this->_sendMainToAll(stop);
380
381 this->_finalize();
382 }
383 else
384 {
385 while (!stop) // until loop breaks due to stop signal
386 {
387 this->_recvMainToWorker(0,stop);
388 if (stop)
389 break;
390
391 this->_recvMainToWorker(0,
392 branch,
393 ID,
394 dataLevel);
395 this->_recvMainToWorker(0,
396 dataEntry,
397 dataInterval,
398 reference
399 );
400
401 this->_correction(dataEntry,
402 dataInterval,
403 dataLevel,
404 reference,
405 distances,
406 stepSolutions,
407 upperDistance,
408 lowerDistance);
409
410 this->_sendWorkerToMain(0,
411 branch,
412 ID);
413 this->_sendWorkerToMain(0,
414 distances,
415 stepSolutions,
416 upperDistance,
417 lowerDistance);
418 }
419 }
420
421#else
422 GISMO_NO_IMPLEMENTATION;
423#endif
424
425 // #ifdef GISMO_WITH_MPI
426 // gsMPIInfo(m_rank)<<"[MPI Process "<<m_rank<<" of "<<m_comm.size()<<"] Hi!\n";
427 // #endif
428}
429
430template <class T>
431template <bool _hasWorkers>
432typename std::enable_if<!_hasWorkers, void>::type
433gsAPALM<T>::parallelSolve_impl()
434{
435 solution_t reference;
436 index_t ID;
437 index_t it = 0;
438 std::vector<T> distances;
439 std::vector<solution_t> stepSolutions;
440 T lowerDistance, upperDistance;
441 index_t branch;
442 index_t dataLevel;
443
444 m_data.print();
445 std::tuple<index_t, T , solution_t, solution_t> dataEntry;
446 while (!m_data.empty() && it < m_maxIterations)
447 {
448 branch = m_data.getFirstNonEmptyBranch();
449 gsMPIInfo(m_rank)<<"There are "<<m_data.branch(branch).nActive()<<" active jobs and "<<m_data.branch(branch).nWaiting()<<" jobs in the queue of branch "<<branch<<"\n";
450
451 dataEntry = m_data.branch(branch).pop();
452 ID = std::get<0>(dataEntry);
453 dataLevel = m_data.branch(branch).jobLevel(ID);
454 bool success = m_data.branch(branch).getReferenceByID(ID,reference);
455 GISMO_ENSURE(success,"Reference not found");
456 this->_correction(dataEntry,
457 m_data.branch(branch).jobTimes(ID),
458 dataLevel,
459 reference,
460 distances,
461 stepSolutions,
462 upperDistance,
463 lowerDistance
464 );
465
466 m_data.branch(branch).submit(ID,distances,stepSolutions,upperDistance,lowerDistance);
467 m_data.branch(branch).finishJob(ID);
468
469 it++;
470 }
471 this->_finalize();
472}
473
474template <class T>
475void gsAPALM<T>::_finalize()
476{
477 std::vector<solution_t> solutions;
478 std::vector<T> times;
479 std::vector<index_t> levels;
480 index_t nBranches = m_data.nBranches();
481 m_solutions.resize(nBranches);
482 m_times.resize(nBranches);
483 m_levels.resize(nBranches);
484 m_lvlSolutions.resize(nBranches);
485 m_lvlTimes.resize(nBranches);
486 for (index_t b=0; b!=nBranches; b++)
487 {
488 std::tie(times,solutions,levels) = m_data.branch(b).getFlatSolution();
489 m_times.at(b) = times;
490 m_solutions.at(b) = solutions;
491 m_levels.at(b) = levels;
492
493 index_t maxLevel = m_data.branch(b).maxLevel();
494 m_lvlSolutions[b].resize(maxLevel+1);
495 m_lvlTimes[b].resize(maxLevel+1);
496 for (size_t k=0; k!=m_solutions.size(); k++)
497 {
498 GISMO_ASSERT(m_lvlSolutions[b].size() > (size_t)m_levels[b][k],"level mismatch, maxLevel = " << maxLevel << "level = "<<m_levels[b][k]);
499 GISMO_ASSERT(m_lvlTimes[b].size() > (size_t)m_levels[b][k],"level mismatch, maxLevel = " << maxLevel << "level = "<<m_levels[b][k]);
500 m_lvlSolutions[b][m_levels[b][k]].push_back(&(m_solutions[b][k]));
501 m_lvlTimes[b][m_levels[b][k]].push_back(&(m_times[b][k]));
502 }
503 }
504}
505
506template <class T>
507template <bool _hasWorkers>
508typename std::enable_if< _hasWorkers, void>::type
509gsAPALM<T>::_solve_impl(index_t Nsteps)
510{
511#ifdef GISMO_WITH_MPI
512 GISMO_ASSERT(m_comm.size()>1,"Something went wrong. This implementation only works when >1 processes are availbale, but nprocesses = "<<m_comm.size());
513
514 bool stop = false; // stop signal
515 solution_t reference;
516 index_t ID;
517 index_t it = 0;
518 index_t branch;
519 index_t dataLevel;
520 std::tuple<index_t, T , solution_t, solution_t> dataEntry;
521
522 // Correction
523 std::pair<T,T> dataInterval;
524 std::vector<T> distances;
525 std::vector<solution_t> stepSolutions;
526 T lowerDistance, upperDistance;
527
528 // Initiation
529 T tstart;
530 T distance;
531 std::vector<solution_t> solutions;
532 bool bifurcation;
533
534
535 //----------------------------------------------------------------------------------------
536 //MAIN------------------------------------------------------------------------------------
537 //----------------------------------------------------------------------------------------
538 if (m_rank==0)
539 {
540 index_t njobs = 0;
541
542 // K keeps track of the iterations per branch
543 std::vector<index_t> K(1);
544 K[0] = 1;
545
546 // Fill the start interval list
547 T L0;
548 gsMatrix<T> U0;
549 T dL;
550 while (!m_starts.empty())
551 {
552 // Initialize solutions
553 U0 = std::get<0>(m_starts.front()).first;
554 L0 = std::get<0>(m_starts.front()).second;
555 dL = std::get<1>(m_starts.front());
556 m_starts.pop();
557
558 gsAPALMData<T,solution_t> data = m_dataEmpty;
559 branch = m_data.add(data);
560 m_data.branch(branch).addStartPoint(T(0),std::make_pair(U0,L0));
561 m_data.branch(branch).setLength(dL); // sets the default length that is used when started from a Point
562 }
563
564 gsMPIInfo(m_rank)<<"Adding workers...\n";
565 for (index_t w = 1; w!=m_proc_count; w++)
566 m_workers.push(w);
567
568 while (!m_data.empty() && it < m_maxIterations && !m_workers.empty())
569 {
570 gsMPIInfo(m_rank)<<"Available workers:\t";
571 std::queue<index_t> workers_copy = m_workers;
572 while (!workers_copy.empty())
573 {
574 gsMPIInfo(m_rank)<<workers_copy.front()<<"\t";
575 workers_copy.pop();
576 }
577 gsMPIInfo(m_rank)<<"\n";
578
579
580 branch = m_data.getFirstNonEmptyBranch();
581 gsMPIInfo(m_rank)<<"There are "<<m_data.branch(branch).nActive()<<" active jobs and "<<m_data.branch(branch).nWaiting()<<" jobs in the queue of branch "<<branch<<"\n";
582
583 dataEntry = m_data.branch(branch).pop();
584 ID = std::get<0>(dataEntry);
585 dataLevel = m_data.branch(branch).jobLevel(ID);
586 // Initialization intervals start at level 0
587 if (dataLevel==0)
588 {
589 this->_sendMainToWorker(m_workers.front(),stop);
590 this->_sendMainToWorker(m_workers.front(),
591 branch,
592 ID,
593 dataLevel);
594 this->_sendMainToWorker(m_workers.front(),
595 dataEntry,
596 m_data.branch(branch).jobStartTime(ID)
597 );
598 }
599 // Correction intervals have level > 0
600 else
601 {
602 bool success = m_data.branch(branch).getReferenceByID(ID,reference);
603 GISMO_ENSURE(success,"Reference not found");
604
605 this->_sendMainToWorker(m_workers.front(),stop);
606 this->_sendMainToWorker(m_workers.front(),
607 branch,
608 ID,
609 dataLevel);
610 this->_sendMainToWorker(m_workers.front(),
611 dataEntry,
612 m_data.branch(branch).jobTimes(ID),
613 reference
614 );
615 }
616 m_workers.pop();
617 it++;
618 njobs++;
619 }
620
621 while (njobs > 0)
622 {
623 gsMPIInfo(m_rank)<<njobs<<" job(s) running\n";
624 index_t source = MPI_ANY_SOURCE;
625 // Receive meta-data
626 this->_recvWorkerToMain(source,
627 branch,
628 ID);
629
630 // Initialization intervals start at level 0
631 dataLevel = m_data.branch(branch).jobLevel(ID);
632 if (dataLevel == 0)
633 {
634 this->_recvWorkerToMain(source,
635 distance,
636 solutions,
637 bifurcation);
638
639 T tstart = m_data.branch(branch).jobStartTime(ID);
640 m_data.branch(branch).appendData(tstart+distance,solutions[0],false);
641 m_data.branch(branch).finishJob(ID);
642
643 if (K[branch]++ < Nsteps-1) // add a new point
644 {
645 if (!bifurcation)
646 {
647 GISMO_ASSERT(solutions.size()==1,"There must be one solution, but solutions.size() = "<<solutions.size());
648 m_data.branch(branch).appendPoint(true);
649 // sets the default length that is used when started from a Point.
650 // After bifurcation, it's the original one times the branch length multiplier times the bifurcation length multiplier
651 if (branch!=0 && ID==0)
652 m_data.branch(branch).setLength(m_data.branch(branch).getLength()/m_bifLengthMult);
653 }
654 else
655 {
656 GISMO_ASSERT(solutions.size()==2,"There must be two solutions!");
657 gsAPALMData<T,solution_t> data = m_dataEmpty;
658 branch = m_data.add(data);
659 K.push_back(1);
660 m_data.branch(branch).addStartPoint(T(0),solutions[1],true);
661 // Sets the default length that is used when started from a Point.
662 // After bifurcation, it's the original one times the branch length multiplier times the bifurcation length multiplier
663 m_data.branch(branch).setLength(m_ALM->getLength()*m_branchLengthMult*m_bifLengthMult);
664 }
665 }
666 }
667 // Correction intervals have level > 0
668 else
669 {
670 this->_recvWorkerToMain(source,
671 distances,
672 stepSolutions,
673 upperDistance,
674 lowerDistance);
675
676 m_data.branch(branch).submit(ID,distances,stepSolutions,upperDistance,lowerDistance);
677 m_data.branch(branch).finishJob(ID);
678 }
679
680 // Remove job
681 njobs--;
682 // Add worker to pool
683 m_workers.push(source);
684
685 gsMPIInfo(m_rank)<<"Available workers:\t";
686 std::queue<index_t> workers_copy = m_workers;
687 while (!workers_copy.empty())
688 {
689 gsInfo<<workers_copy.front()<<"\t";
690 workers_copy.pop();
691 }
692 gsInfo<<"\n";
693
694 // As long as the queue is not empty, the iterations are lower than the max number and the pool of workers is not empty
695 while (!m_data.empty() && it < m_maxIterations && !m_workers.empty())
696 {
697 // // SAME WHILE LOOP AS ABOVE!!!!!!!!
698 branch = m_data.getFirstNonEmptyBranch();
699 gsMPIInfo(m_rank)<<"There are "<<m_data.branch(branch).nActive()<<" active jobs and "<<m_data.branch(branch).nWaiting()<<" jobs in the queue of branch "<<branch<<"\n";
700
701 dataEntry = m_data.branch(branch).pop();
702 ID = std::get<0>(dataEntry);
703
704 // Initialization intervals start at level 0
705 dataLevel = m_data.branch(branch).jobLevel(ID);
706 if (dataLevel==0)
707 {
708 this->_sendMainToWorker(m_workers.front(),stop);
709 this->_sendMainToWorker(m_workers.front(),
710 branch,
711 ID,
712 dataLevel);
713
714 this->_sendMainToWorker(m_workers.front(),
715 dataEntry,
716 m_data.branch(branch).jobStartTime(ID)
717 );
718 }
719 // Correction intervals have level > 0
720 else
721 {
722 bool success = m_data.branch(branch).getReferenceByID(ID,reference);
723 GISMO_ENSURE(success,"Reference not found");
724
725 this->_sendMainToWorker(m_workers.front(),stop);
726 this->_sendMainToWorker(m_workers.front(),
727 branch,
728 ID,
729 dataLevel);
730 this->_sendMainToWorker(m_workers.front(),
731 dataEntry,
732 m_data.branch(branch).jobTimes(ID),
733 reference
734 );
735 }
736 m_workers.pop();
737 it++;
738 njobs++;
739 }
740 }
741 stop = true;
742 this->_sendMainToAll(stop);
743
744 this->_finalize();
745 }
746 else
747 {
748 while (!stop) // until loop breaks due to stop signal
749 {
750 this->_recvMainToWorker(0,stop);
751 if (stop)
752 break;
753
754 // Receive meta-data
755 this->_recvMainToWorker(0,
756 branch,
757 ID,
758 dataLevel);
759
760 // Initialization intervals start at level 0
761 if (dataLevel==0)
762 {
763 this->_recvMainToWorker(0,
764 dataEntry,
765 tstart
766 );
767
768 this->_initiation(dataEntry,
769 tstart,
770 distance,
771 solutions,
772 bifurcation
773 );
774
775 this->_sendWorkerToMain(0,
776 branch,
777 ID);
778 this->_sendWorkerToMain(0,
779 distance,
780 solutions,
781 bifurcation);
782 }
783 // Correction intervals have level > 0
784 else
785 {
786 this->_recvMainToWorker(0,
787 dataEntry,
788 dataInterval,
789 reference
790 );
791
792 this->_correction(dataEntry,
793 dataInterval,
794 dataLevel,
795 reference,
796 distances,
797 stepSolutions,
798 upperDistance,
799 lowerDistance);
800
801 this->_sendWorkerToMain(0,
802 branch,
803 ID);
804 this->_sendWorkerToMain(0,
805 distances,
806 stepSolutions,
807 upperDistance,
808 lowerDistance);
809 }
810 }
811 }
812
813#else
814 GISMO_NO_IMPLEMENTATION;
815#endif
816
817 // #ifdef GISMO_WITH_MPI
818 // gsMPIInfo(m_rank)<<"[MPI Process "<<m_rank<<" of "<<m_comm.size()<<"] Hi!\n";
819 // #endif
820}
821
822template <class T>
823template <bool _hasWorkers>
824typename std::enable_if<!_hasWorkers, void>::type
825gsAPALM<T>::_solve_impl(index_t Nsteps)
826{
827 index_t ID;
828 index_t it = 0;
829 std::vector<index_t> K(1);
830 K[0] = 1;
831
832 index_t branch;
833 index_t dataLevel;
834
835 T L0;
836 gsMatrix<T> U0;
837 T dL;
838 while (!m_starts.empty())
839 {
840 // Initialize solutions
841 U0 = std::get<0>(m_starts.front()).first;
842 L0 = std::get<0>(m_starts.front()).second;
843 dL = std::get<1>(m_starts.front());
844 m_starts.pop();
845
846 gsAPALMData<T,solution_t> data = m_dataEmpty;
847 branch = m_data.add(data);
848 m_data.branch(branch).addStartPoint(T(0),std::make_pair(U0,L0));
849 m_data.branch(branch).setLength(dL); // sets the default length that is used when started from a Point
850 }
851
852 m_data.print();
853 std::tuple<index_t, T , solution_t, solution_t> dataEntry;
854 while (!m_data.empty() && it < m_maxIterations)
855 {
856 branch = m_data.getFirstNonEmptyBranch();
857 gsMPIInfo(m_rank)<<"There are "<<m_data.branch(branch).nActive()<<" active jobs and "<<m_data.branch(branch).nWaiting()<<" jobs in the queue of branch "<<branch<<"\n";
858
859 dataEntry = m_data.branch(branch).pop();
860 ID = std::get<0>(dataEntry);
861
862 dataLevel = m_data.branch(branch).jobLevel(ID);
863 if (dataLevel==0)
864 {
865 if (m_verbose) gsMPIInfo(m_rank)<<"Initialization\n";
866 T startTime = m_data.branch(branch).jobStartTime(ID);
867 T distance;
868 std::vector<solution_t> solutions;
869 bool bifurcation;
870 this->_initiation(dataEntry,
871 startTime,
872 distance,
873 solutions,
874 bifurcation
875 );
876
877 m_data.branch(branch).appendData(startTime+distance,solutions[0],false);
878 m_data.branch(branch).finishJob(ID);
879 // if no bifurcation is hit, we mark the end of the interval as the start point
880
881 if (K[branch]++ < Nsteps - 1) // add a new point
882 {
883 if (!bifurcation)
884 {
885 GISMO_ASSERT(solutions.size()==1,"There must be one solution!");
886 m_data.branch(branch).appendPoint(true);
887 // sets the default length that is used when started from a Point.
888 // After bifurcation, it's the original one times the branch length multiplier times the bifurcation length multiplier
889 if (branch!=0 && ID==0)
890 m_data.branch(branch).setLength(m_data.branch(branch).getLength()/m_bifLengthMult);
891 }
892 else
893 {
894 GISMO_ASSERT(solutions.size()==2,"There must be two solutions!");
895 gsAPALMData<T,solution_t> data = m_dataEmpty;
896 branch = m_data.add(data);
897 K.push_back(1);
898 m_data.branch(branch).addStartPoint(T(0),solutions[1],true); // use the second solution that is given in solutions
899 // Sets the default length that is used when started from a Point.
900 // After bifurcation, it's the original one times the branch length multiplier times the bifurcation length multiplier
901 m_data.branch(branch).setLength(m_ALM->getLength()*m_branchLengthMult*m_bifLengthMult);
902 }
903 }
904 }
905 else
906 {
907 if (m_verbose) gsMPIInfo(m_rank)<<"Correction\n";
908 solution_t reference;
909 std::vector<T> distances;
910 std::vector<solution_t> stepSolutions;
911 T lowerDistance, upperDistance;
912 bool success = m_data.branch(branch).getReferenceByID(ID,reference);
913 GISMO_ENSURE(success,"Reference not found");
914 this->_correction(dataEntry,
915 m_data.branch(branch).jobTimes(ID),
916 dataLevel,
917 reference,
918 distances,
919 stepSolutions,
920 upperDistance,
921 lowerDistance
922 );
923
924 m_data.branch(branch).submit(ID,distances,stepSolutions,upperDistance,lowerDistance);
925 m_data.branch(branch).finishJob(ID);
926
927 it++;
928 }
929 }
930 this->_finalize();
931}
932
933// NOTE: This does not make new branches!
934template <class T>
935void gsAPALM<T>::_initiation( const std::tuple<index_t, T , solution_t, solution_t> & dataEntry,
936 const T & startTime,
937 T & distance,
938 std::vector<solution_t>&solutions,
939 bool & bifurcation )
940{
941 solution_t start, prev;
942 index_t ID;
943 T tstart = startTime;
944 T dL, dL0;
945 gsVector<T> DeltaU;
946 T DeltaL;
947
948 T Lprev, Lold;
949 gsMatrix<T> Uprev, Uold;
950
952 solutions.clear();
953 std::tie(ID,dL,start,prev) = dataEntry;
954 std::tie(Uold,Lold) = start;
955 std::tie(Uprev,Lprev) = prev;
956 dL0 = dL;
957
958 m_ALM->setLength(dL);
959 m_ALM->setSolution(Uold,Lold);
960 // m_ALM->resetStep();
961 m_ALM->setPrevious(Uprev,Lprev);
962 gsMPIDebug(m_rank)<<"Start - ||u|| = "<<Uold.norm()<<", L = "<<Lold<<"\n";
963 gsMPIDebug(m_rank)<<"Prev - ||u|| = "<<Uprev.norm()<<", L = "<<Lprev<<"\n";
964
965 bool diverged = true;
966 bifurcation = false;
967 while (diverged)
968 {
969 gsMPIInfo(m_rank)<<"Starting with ID "<<ID<<" from (|U|,L) = ("<<Uold.norm()<<","<<Lold<<"), curve time = "<<tstart<<", arc-length = "<<dL<<"\n";
970 // Set a step
971 gsStatus status = m_ALM->step();
972 diverged = (status!=gsStatus::Success);
973 if (status==gsStatus::NotConverged || status==gsStatus::AssemblyError)
974 {
975 if (m_verbose) gsMPIInfo(m_rank)<<"Error: Loop terminated, arc length method did not converge.\n";
976 dL = dL / 2;
977 m_ALM->setLength(dL);
978 m_ALM->setSolution(Uold,Lold);
979 continue;
980 }
981 }
982 dL = dL0;
983 if (m_singularPoint)
984 {
985 m_ALM->computeStability(m_ALM->options().getSwitch("Quasi"));
986 if (m_ALM->stabilityChange())
987 {
988 gsMPIInfo(m_rank)<<"Bifurcation spotted!"<<"\n";
989 if (m_ALM->isBifurcation(false))
990 {
991 m_ALM->computeSingularPoint(Uold,Lold,false,false,false);
992 bifurcation = true;
993 }
994 }
995 }
996
997 DeltaU = m_ALM->solutionU() - Uold;
998 DeltaL = m_ALM->solutionL() - Lold;
999 distance = m_ALM->distance(DeltaU,DeltaL);
1000 solutions.push_back(std::make_pair(m_ALM->solutionU(),m_ALM->solutionL()));
1001
1002 this->serialStepOutput(solutions[0],tstart,ID);
1003
1004 if (bifurcation)
1005 {
1006 m_ALM->switchBranch();
1007 // add the extra point after bifurcation to the export
1008 solutions.push_back(std::make_pair(m_ALM->solutionU(),m_ALM->solutionL()));
1009 }
1010}
1011
1012// NOTE: This does not make new branches!
1013template <class T>
1014void gsAPALM<T>::_correction( const std::tuple<index_t, T , solution_t, solution_t> & dataEntry,
1015 const std::pair<T,T> & dataInterval,
1016 const index_t & dataLevel,
1017 const solution_t & dataReference,
1018 std::vector<T> & distances,
1019 std::vector<solution_t>&stepSolutions,
1020 T & upperDistance,
1021 T & lowerDistance )
1022{
1023 solution_t start, prev, reference;
1024 index_t ID;
1025 T tstart = 0;
1026 T tend = 0;
1027 T dL, dL0;
1028 gsVector<T> DeltaU;
1029 T DeltaL;
1030 index_t Nintervals = m_subIntervals;
1031 bool bisected = false;
1032 T dL_rem = 0;
1033
1034 T Lprev, Lold;
1035 gsMatrix<T> Uprev, Uold;
1036
1038 stepSolutions.resize(Nintervals);
1039 std::vector<T> stepTimes(Nintervals);
1040 distances.resize(Nintervals+1);
1041
1042 std::tie(ID,dL0,start,prev) = dataEntry;
1043 std::tie(Uold,Lold) = start;
1044 std::tie(Uprev,Lprev) = prev;
1045 std::tie(tstart,tend) = dataInterval;
1046
1047 gsMPIDebug(m_rank)<<"tstart = "<<tstart<<" , "<<"tend = "<<tend<<"\n";
1048
1049 gsMatrix<T> Uori = Uold;
1050 T Lori = Lold;
1051
1052 dL0 = dL0 / Nintervals;
1053 dL = dL0;
1054
1055 m_ALM->setLength(dL);
1056 m_ALM->setSolution(Uold,Lold);
1057 // m_ALM->resetStep();
1058 m_ALM->setPrevious(Uprev,Lprev);
1059 gsMPIDebug(m_rank)<<"Start - ||u|| = "<<Uold.norm()<<", L = "<<Lold<<"\n";
1060 gsMPIDebug(m_rank)<<"Prev - ||u|| = "<<Uprev.norm()<<", L = "<<Lprev<<"\n";
1061
1062 T s = 0;
1063 T time = tstart;
1064
1065 gsMPIInfo(m_rank)<<"Starting with ID "<<ID<<" from (|U|,L) = ("<<Uold.norm()<<","<<Lold<<"), curve time = "<<time<<"\n";
1066 for (index_t k = 0; k!=Nintervals; k++)
1067 {
1068 gsMPIDebug(m_rank)<<"Interval "<<k+1<<" of "<<Nintervals<<"\n";
1069 gsMPIDebug(m_rank)<<"Start - ||u|| = "<<Uold.norm()<<", L = "<<Lold<<"\n";
1070
1071 gsStatus status = m_ALM->step();
1072 if (status==gsStatus::NotConverged || status==gsStatus::AssemblyError)
1073 {
1074 gsMPIInfo(m_rank)<<"Error: Loop terminated, arc length method did not converge.\n";
1075 dL = dL / 2.;
1076 dL_rem += dL; // add the remainder of the interval to dL_rem
1077 m_ALM->setLength(dL);
1078 m_ALM->setSolution(Uold,Lold);
1079 bisected = true;
1080 k -= 1;
1081 continue;
1082 }
1083
1084 std::pair<gsVector<T>,T> pair = std::make_pair(m_ALM->solutionU(),m_ALM->solutionL());
1085 stepSolutions.at(k) = pair;
1086 DeltaU = m_ALM->solutionU() - Uold;
1087 DeltaL = m_ALM->solutionL() - Lold;
1088
1089 distances.at(k) = m_ALM->distance(DeltaU,DeltaL);
1090
1091 s += distances.at(k);
1092 time += distances.at(k);
1093 stepTimes.at(k) = time;
1094
1095 Uold = m_ALM->solutionU();
1096 Lold = m_ALM->solutionL();
1097
1098 if (!bisected) // if bisected = false
1099 dL = dL0;
1100 else
1101 {
1102 // if the current interval has finished, but was refined before.
1103 // The next interval should have the remaining length.
1104 // Also, Nintervals should increase
1105 //
1106 dL = dL_rem;
1107 Nintervals++;
1108 stepSolutions.resize(Nintervals);
1109 stepTimes.resize(Nintervals);
1110 distances.resize(Nintervals+1);
1111 }
1112
1113 this->parallelStepOutput(pair,time,k);
1114
1115 m_ALM->setLength(dL);
1116 dL_rem = 0;
1117 bisected = false;
1118 }
1119
1120 gsMPIDebug(m_rank)<<"Ref - ||u|| = "<<dataReference.first.norm()<<", L = "<<dataReference.second<<"\n";
1121
1122 GISMO_ASSERT(dataReference.first.normalized().size()==m_ALM->solutionU().normalized().size(),"Reference solution and current solution have different size! ref.size() = "<<dataReference.first.normalized().size()<<"; sol.size() = "<<m_ALM->solutionU().normalized().size());
1123 GISMO_ASSERT((dataReference.first.normalized()).dot(m_ALM->solutionU().normalized())>-0.8,
1124 "Reference is almost in the opposite direction of the solution. Are branches mirrored? result:" << (reference.first.normalized()).dot(m_ALM->solutionU().normalized()));
1125
1126 DeltaU = dataReference.first - m_ALM->solutionU();
1127 DeltaL = dataReference.second - m_ALM->solutionL();
1128 distances.back() = m_ALM->distance(DeltaU,DeltaL);
1129
1130 DeltaU = dataReference.first - Uori;
1131 DeltaL = dataReference.second - Lori;
1132 upperDistance = m_ALM->distance(DeltaU,DeltaL);
1133
1134 DeltaU = stepSolutions.back().first - Uori;
1135 DeltaL = stepSolutions.back().second - Lori;
1136 lowerDistance = m_ALM->distance(DeltaU,DeltaL);
1137
1138 std::vector<solution_t> stepSolutionsExport(Nintervals+2);
1139 std::vector<T> stepTimesExport(Nintervals+2);
1140 // Export parallel interval output
1141 // Solutions
1142 std::pair<gsVector<T>,T> front = std::make_pair(start.first,start.second);
1143 stepSolutionsExport.front() = front;
1144 std::pair<gsVector<T>,T> back = std::make_pair(dataReference.first,dataReference.second);
1145 stepSolutionsExport.back() = back;
1146 // Times
1147 stepTimesExport.front() = tstart;
1148 stepTimesExport.back() = tend;
1149
1150 // Temporarily swap the interval solutions and times
1151 for (index_t k=0; k!=Nintervals; k++)
1152 {
1153 std::swap(stepSolutions.at(k),stepSolutionsExport.at(k+1));
1154 std::swap(stepTimes.at(k),stepTimesExport.at(k+1));
1155 }
1156
1157 this->parallelIntervalOutput(stepSolutionsExport,stepTimesExport,dataLevel,ID);
1158
1159 // And swap them back
1160 for (index_t k=0; k!=Nintervals; k++)
1161 {
1162 std::swap(stepSolutionsExport.at(k+1),stepSolutions.at(k));
1163 std::swap(stepTimesExport.at(k+1),stepTimes.at(k));
1164 }
1165}
1166
1167// -----------------------------------------------------------------------------------------------------
1168// MPI functions
1169// -----------------------------------------------------------------------------------------------------
1170#ifdef GISMO_WITH_MPI
1171
1172template <class T>
1173void gsAPALM<T>::_sendMainToWorker( const index_t & workerID,
1174 const index_t & branch,
1175 const index_t & jobID,
1176 const index_t & dataLevel)
1177{
1178 gsMPIInfo(m_rank)<<"Sending data from "<<m_rank<<" to "<<workerID<<"\n";
1179 MPI_Request req[3];
1180 m_comm.isend(&branch ,1,workerID,&req[0] ,0 );
1181 m_comm.isend(&jobID ,1,workerID,&req[1] ,1 );
1182 m_comm.isend(&dataLevel ,1,workerID,&req[2] ,2 );
1183 MPI_Waitall( 3, req, MPI_STATUSES_IGNORE );
1184}
1185
1186template <class T>
1187void gsAPALM<T>::_sendMainToWorker( const index_t & workerID,
1188 const std::tuple<index_t, T , solution_t, solution_t> & dataEntry,
1189 const std::pair<T,T> & dataInterval,
1190 const solution_t & dataReference )
1191{
1192 gsMPIInfo(m_rank)<<"Sending data from "<<m_rank<<" to "<<workerID<<"\n";
1193
1194 index_t ID;
1195 T dL0;
1196 solution_t start, prev;
1197 gsVector<T> startU, prevU, refU;
1198 index_t startSize,prevSize, refSize;
1199 T startL, prevL, refL;
1200
1201 T tstart, tend;
1202
1203 std::tie(ID,dL0,start,prev) = dataEntry;
1204 std::tie(tstart,tend) = dataInterval;
1205
1206 std::tie(startU,startL) = start;
1207 startSize = startU.size();
1208 std::tie(prevU,prevL) = prev;
1209 prevSize = prevU.size();
1210 std::tie(refU,refL) = dataReference;
1211 refSize = refU.size();
1212
1213 // Put all data 0-14 in a struct and send a single struct
1214 MPI_Request req[10];
1215 m_comm.isend(&ID ,1,workerID,&req[0] ,0 );
1216 m_comm.isend(&dL0 ,1,workerID,&req[1] ,1 );
1217 m_comm.isend(&tstart ,1,workerID,&req[2] ,2 );
1218 m_comm.isend(&tend ,1,workerID,&req[3] ,3 );
1219
1220 m_comm.isend(&startSize ,1,workerID,&req[4] ,4 );
1221 m_comm.isend(&startL ,1,workerID,&req[5] ,5 );
1222
1223 m_comm.isend(&prevSize ,1,workerID,&req[6] ,6 );
1224 m_comm.isend(&prevL ,1,workerID,&req[7] ,7 );
1225
1226 m_comm.isend(&refSize ,1,workerID,&req[8] ,8 );
1227 m_comm.isend(&refL ,1,workerID,&req[9] ,9 );
1228
1229 MPI_Waitall( 10, req , MPI_STATUSES_IGNORE );
1230 MPI_Request req_data[3];
1231 m_comm.isend(startU.data(), startSize,workerID,&req_data[0],10);
1232 m_comm.isend(prevU.data(), prevSize, workerID,&req_data[1],11);
1233 m_comm.isend(refU.data(), refSize, workerID,&req_data[2],12);
1234
1235 MPI_Waitall( 3 , req_data, MPI_STATUSES_IGNORE );
1236}
1237
1238template <class T>
1239void gsAPALM<T>::_sendMainToWorker( const index_t & workerID,
1240 const std::tuple<index_t, T , solution_t, solution_t> & dataEntry,
1241 const T & startTime)
1242{
1243 gsMPIInfo(m_rank)<<"Sending data from "<<m_rank<<" to "<<workerID<<"\n";
1244
1245 index_t ID;
1246 T dL0;
1247 solution_t start, prev;
1248 gsVector<T> startU, prevU;
1249 index_t startSize,prevSize;
1250 T startL, prevL;
1251
1252 T tstart;
1253
1254 std::tie(ID,dL0,start,prev) = dataEntry;
1255 tstart = startTime;
1256
1257 std::tie(startU,startL) = start;
1258 startSize = startU.size();
1259 std::tie(prevU,prevL) = prev;
1260 prevSize = prevU.size();
1261
1262 // Put all data 0-10 in a struct and send a single struct
1263 MPI_Request req[7];
1264 m_comm.isend(&ID ,1,workerID,&req[0] ,0 );
1265 m_comm.isend(&dL0 ,1,workerID,&req[1] ,1 );
1266 m_comm.isend(&tstart ,1,workerID,&req[2] ,2 );
1267
1268 m_comm.isend(&startSize ,1,workerID,&req[3] ,3 );
1269 m_comm.isend(&startL ,1,workerID,&req[4] ,4 );
1270
1271 m_comm.isend(&prevSize ,1,workerID,&req[5] ,5 );
1272 m_comm.isend(&prevL ,1,workerID,&req[6] ,6 );
1273 MPI_Waitall( 7, req , MPI_STATUSES_IGNORE );
1274
1275 MPI_Request req_data[2];
1276 m_comm.isend(startU.data(), startSize,workerID,&req_data[0],7);
1277 m_comm.isend(prevU.data(), prevSize, workerID,&req_data[1],8);
1278
1279 MPI_Waitall( 2, req_data, MPI_STATUSES_IGNORE );
1280}
1281
1282template <class T>
1283void gsAPALM<T>::_sendMainToWorker( const index_t & workerID,
1284 const bool & stop )
1285{
1286 // NOTE: when we handle the stop signal with isend/irecv, we get a deadlock
1287 // TODO: handle stop with broadcast
1288 gsMPIInfo(m_rank)<<"Sending stop signal from "<<m_rank<<" to "<<workerID<<"\n";
1289 m_comm.send(&stop,1,workerID,99);
1290}
1291
1292template <class T>
1293void gsAPALM<T>::_sendMainToAll( const bool & stop )
1294{
1295 for (index_t w = 1; w!=m_proc_count; w++)
1296 this->_sendMainToWorker(w,stop);
1297}
1298
1299template <class T>
1300void gsAPALM<T>::_recvMainToWorker( const index_t & sourceID,
1301 index_t & branch,
1302 index_t & jobID,
1303 index_t & dataLevel)
1304{
1305 gsMPIInfo(m_rank)<<"Receiving data on "<<m_rank<<" from "<<sourceID<<"\n";
1306 MPI_Request req[3];
1307 m_comm.irecv(&branch ,1,sourceID,&req[0] ,0 );
1308 m_comm.irecv(&jobID ,1,sourceID,&req[1] ,1 );
1309 m_comm.irecv(&dataLevel ,1,sourceID,&req[2] ,2 );
1310 MPI_Waitall( 3, req, MPI_STATUSES_IGNORE );
1311}
1312
1313
1314template <class T>
1315void gsAPALM<T>::_recvMainToWorker( const index_t & sourceID,
1316 std::tuple<index_t, T , solution_t, solution_t> & dataEntry,
1317 std::pair<T,T> & dataInterval,
1318 solution_t & dataReference)
1319{
1320 gsMPIInfo(m_rank)<<"Receiving data on "<<m_rank<<" from "<<sourceID<<"\n";
1321
1322 index_t ID;
1323 T dL0;
1324 solution_t start, prev;
1325 gsVector<T> startU, prevU, refU;
1326 index_t startSize,prevSize, refSize;
1327 T startL, prevL, refL;
1328
1329 T tstart, tend;
1330
1331 // Put all data 0-14 in a struct and recv a single struct
1332 MPI_Request req[10];
1333 m_comm.irecv(&ID ,1,sourceID,&req[0] ,0 );
1334 m_comm.irecv(&dL0 ,1,sourceID,&req[1] ,1 );
1335 m_comm.irecv(&tstart ,1,sourceID,&req[2] ,2 );
1336 m_comm.irecv(&tend ,1,sourceID,&req[3] ,3 );
1337
1338 m_comm.irecv(&startSize ,1,sourceID,&req[4] ,4 );
1339 m_comm.irecv(&startL ,1,sourceID,&req[5] ,5 );
1340
1341 m_comm.irecv(&prevSize ,1,sourceID,&req[6] ,6 );
1342 m_comm.irecv(&prevL ,1,sourceID,&req[7] ,7 );
1343
1344 m_comm.irecv(&refSize ,1,sourceID,&req[8] ,8 );
1345 m_comm.irecv(&refL ,1,sourceID,&req[9] ,9 );
1346
1347 MPI_Waitall( 10, req, MPI_STATUSES_IGNORE );
1348 dataInterval = std::make_pair(tstart,tend);
1349
1350 startU.resize(startSize);
1351 prevU .resize(prevSize);
1352 refU .resize(refSize);
1353
1354 MPI_Request req_data[3];
1355 m_comm.irecv(startU.data(), startSize,sourceID,&req_data[0],10);
1356 m_comm.irecv(prevU.data(), prevSize, sourceID,&req_data[1],11);
1357 m_comm.irecv(refU.data(), refSize, sourceID,&req_data[2],12);
1358
1359 MPI_Waitall( 3, req_data, MPI_STATUSES_IGNORE );
1360
1361 start = std::make_pair(startU,startL);
1362 prev = std::make_pair(prevU ,prevL);
1363 dataReference = std::make_pair(refU ,refL);
1364
1365 dataEntry = std::make_tuple(ID,dL0,start,prev);
1366}
1367
1368template <class T>
1369void gsAPALM<T>::_recvMainToWorker( const index_t & sourceID,
1370 std::tuple<index_t, T , solution_t, solution_t> & dataEntry,
1371 T & startTime)
1372{
1373 gsMPIInfo(m_rank)<<"Receiving data on "<<m_rank<<" from "<<sourceID<<"\n";
1374
1375 index_t ID;
1376 T dL0;
1377 solution_t start, prev;
1378 gsVector<T> startU, prevU;
1379 index_t startSize,prevSize;
1380 T startL, prevL;
1381
1382 // Put all data 0-14 in a struct and recv a single struct
1383 MPI_Request req[7];
1384 m_comm.irecv(&ID ,1,sourceID,&req[0] ,0 );
1385 m_comm.irecv(&dL0 ,1,sourceID,&req[1] ,1 );
1386 m_comm.irecv(&startTime ,1,sourceID,&req[2] ,2 );
1387
1388 m_comm.irecv(&startSize ,1,sourceID,&req[3] ,3 );
1389 m_comm.irecv(&startL ,1,sourceID,&req[4] ,4 );
1390
1391 m_comm.irecv(&prevSize ,1,sourceID,&req[5] ,5 );
1392 m_comm.irecv(&prevL ,1,sourceID,&req[6] ,6 );
1393
1394 MPI_Waitall( 7, req, MPI_STATUSES_IGNORE );
1395
1396 startU.resize(startSize);
1397 prevU .resize(prevSize);
1398
1399 MPI_Request req_data[2];
1400 m_comm.irecv(startU.data(), startU.size(),sourceID,&req_data[0],7);
1401 m_comm.irecv(prevU.data(), prevU.size(), sourceID,&req_data[1],8);
1402
1403 MPI_Waitall( 2, req_data, MPI_STATUSES_IGNORE );
1404
1405 start = std::make_pair(startU,startL);
1406 prev = std::make_pair(prevU ,prevL);
1407
1408 dataEntry = std::make_tuple(ID,dL0,start,prev);
1409}
1410
1411template <class T>
1412void gsAPALM<T>::_recvMainToWorker( const index_t & sourceID,
1413 bool & stop)
1414{
1415 // when we handle the stop signal with isend/irecv, we get a deadlock
1416 m_comm.recv(&stop,1,sourceID,99);
1417}
1418
1419// Sends metadata
1420template <class T>
1421void gsAPALM<T>::_sendWorkerToMain( const index_t & mainID,
1422 const index_t & branch,
1423 const index_t & jobID)
1424{
1425 gsMPIInfo(m_rank)<<"Sending data from "<<m_rank<<" to "<<mainID<<"\n";
1426
1427 index_t tag = 0;
1428
1429 MPI_Request req_meta[2];
1430 m_comm.isend(&branch ,1, mainID, &req_meta[0], tag++);
1431 m_comm.isend(&jobID ,1, mainID, &req_meta[1], tag++);
1432 MPI_Waitall( 2, req_meta, MPI_STATUSES_IGNORE );
1433}
1434
1435// Receives metadata
1436template <class T>
1437void gsAPALM<T>::_recvWorkerToMain( index_t & sourceID,
1438 index_t & branch,
1439 index_t & jobID)
1440{
1441
1442 index_t tag = 0;
1443
1444 MPI_Status status;
1445 MPI_Request req_meta[2];
1446 m_comm.irecv(&branch ,1, sourceID, &req_meta[0], tag++);
1447 MPI_Wait ( &req_meta[0], &status );
1448 if (sourceID==MPI_ANY_SOURCE) sourceID = status.MPI_SOURCE; // Overwrite source ID to be unique
1449
1450 // Request metadata
1451 m_comm.irecv(&jobID ,1, sourceID, &req_meta[1], tag++);
1452 MPI_Waitall( 2, req_meta, MPI_STATUSES_IGNORE );
1453}
1454
1455template <class T>
1456void gsAPALM<T>::_sendWorkerToMain( const index_t & mainID,
1457 const std::vector<T> & distances,
1458 const std::vector<solution_t> & stepSolutions,
1459 const T & upperDistance,
1460 const T & lowerDistance )
1461{
1462 gsMPIInfo(m_rank)<<"Sending data from "<<m_rank<<" to "<<mainID<<"\n";
1463
1464 index_t tag = 0;
1465 index_t size = stepSolutions.size();
1466 index_t vectorSize;
1467 T load;
1468
1469 MPI_Request req_meta[2];
1470 m_comm.isend(&upperDistance ,1, mainID, &req_meta[0], tag++);
1471 m_comm.isend(&lowerDistance ,1, mainID, &req_meta[1], tag++);
1472
1473 MPI_Request req_solSize;
1474 m_comm.isend(&size ,1, mainID, &req_solSize, tag++);
1475
1476 MPI_Request req_vecSizes[size];
1477 MPI_Request req_loads[size];
1478 MPI_Request req_distances[size+1];
1479 for (index_t k=0; k!=size; k++)
1480 {
1481 m_comm.isend(&distances.at(k) ,1, mainID, &req_distances[k] , tag++);
1482 load = stepSolutions.at(k).second;
1483 m_comm.isend(&load ,1, mainID, &req_loads[k] , tag++);
1484 vectorSize= stepSolutions.at(k).first.size();
1485 m_comm.isend(&vectorSize ,1, mainID, &req_vecSizes[k] , tag++);
1486 }
1487 m_comm.isend(&(distances.at(size)),1, mainID, &req_distances[size], tag++);
1488
1489 MPI_Request req_data[size];
1490 for (index_t k=0; k!=size; k++)
1491 m_comm.isend( stepSolutions.at(k).first.data(),
1492 stepSolutions.at(k).first.size(), mainID, &req_data[k], tag++);
1493
1494 MPI_Wait(&req_solSize,MPI_STATUS_IGNORE);
1495 MPI_Waitall( size, req_vecSizes, MPI_STATUSES_IGNORE );
1496 MPI_Waitall( size, req_data, MPI_STATUSES_IGNORE );
1497 MPI_Waitall( size, req_loads, MPI_STATUSES_IGNORE );
1498 MPI_Waitall( size+1, req_distances, MPI_STATUSES_IGNORE );
1499 MPI_Waitall( 2, req_meta, MPI_STATUSES_IGNORE );
1500}
1501
1502template <class T>
1503void gsAPALM<T>::_sendWorkerToMain( const index_t & mainID,
1504 const T & distance,
1505 const std::vector<solution_t> & solutions,
1506 const bool & bifurcation)
1507{
1508 gsMPIInfo(m_rank)<<"Sending data from "<<m_rank<<" to "<<mainID<<"\n";
1509
1510 index_t tag = 0;
1511 index_t size = solutions.size();
1512 index_t vectorSize;
1513 T load;
1514
1515 MPI_Request req_meta[2];
1516 m_comm.isend(&distance ,1, mainID, &req_meta[0], tag++);
1517 m_comm.isend(&bifurcation ,1, mainID, &req_meta[1], tag++);
1518
1519 MPI_Request req_solSize;
1520 m_comm.isend(&size ,1, mainID, &req_solSize, tag++);
1521
1522 MPI_Request req_vecSizes[size];
1523 MPI_Request req_loads[size];
1524 for (index_t k=0; k!=size; k++)
1525 {
1526 load = solutions.at(k).second;
1527 m_comm.isend(&load ,1, mainID, &req_loads[k] , tag++);
1528 vectorSize= solutions.at(k).first.size();
1529 m_comm.isend(&vectorSize ,1, mainID, &req_vecSizes[k] , tag++);
1530 }
1531
1532 MPI_Request req_data[size];
1533 for (index_t k=0; k!=size; k++)
1534 m_comm.isend( solutions.at(k).first.data(),
1535 solutions.at(k).first.size(), mainID, &req_data[k], tag++);
1536
1537 MPI_Wait(&req_solSize,MPI_STATUS_IGNORE);
1538 MPI_Waitall( size, req_vecSizes, MPI_STATUSES_IGNORE );
1539 MPI_Waitall( size, req_data, MPI_STATUSES_IGNORE );
1540 MPI_Waitall( size, req_loads, MPI_STATUSES_IGNORE );
1541 MPI_Waitall( 2, req_meta, MPI_STATUSES_IGNORE );
1542}
1543
1544template <class T>
1545void gsAPALM<T>::_recvWorkerToMain( index_t & sourceID,
1546 std::vector<T> & distances,
1547 std::vector<solution_t> & stepSolutions,
1548 T & upperDistance,
1549 T & lowerDistance)
1550{
1551
1552 index_t tag = 0;
1553 index_t size;
1554
1555 MPI_Status status;
1556 MPI_Request req_meta[2];
1557 m_comm.irecv(&upperDistance ,1, sourceID, &req_meta[0], tag++);
1558 MPI_Wait ( &req_meta[0], &status );
1559 if (sourceID==MPI_ANY_SOURCE) sourceID = status.MPI_SOURCE; // Overwrite source ID to be unique
1560 gsMPIInfo(m_rank)<<"Receiving data on "<<m_rank<<" from "<<sourceID<<"\n";
1561
1562 // Request metadata
1563 m_comm.irecv(&lowerDistance ,1, sourceID, &req_meta[1], tag++);
1564
1565 // Request the number of solution intervals
1566 MPI_Request req_solSize;
1567 m_comm.irecv(&size ,1, sourceID, &req_solSize, tag++);
1568
1569 MPI_Wait(&req_solSize,MPI_STATUS_IGNORE);
1570
1571 // temporary objects to store solutions and the sizes of the solution vectors
1572 std::vector<gsVector<T>> Us(size);
1573 std::vector<T > Ls(size);
1574 std::vector<index_t > vectorSizes(size);
1575 // final containers for distances and step solutions
1576 distances.resize(size+1);
1577 stepSolutions.resize(size);
1578
1579 // Collect distances, solution vector sizes and loads
1580 MPI_Request req_vecSizes[size];
1581 MPI_Request req_loads[size];
1582 MPI_Request req_distances[size+1];
1583 for (index_t k=0; k!=size; k++)
1584 {
1585 m_comm.irecv(&distances.at(k) ,1, sourceID, &req_distances[k] , tag++);
1586 m_comm.irecv(&Ls.at(k) ,1, sourceID, &req_loads[k] , tag++);
1587 m_comm.irecv(&vectorSizes.at(k) ,1, sourceID, &req_vecSizes[k] , tag++);
1588 }
1589 m_comm.irecv(&(distances.at(size)),1, sourceID, &req_distances[size], tag++);
1590
1591 // When all solution vector sizes are collected, resize the solution vectors
1592 MPI_Waitall( size, req_vecSizes, MPI_STATUSES_IGNORE );
1593 for (index_t k=0; k!=size; k++)
1594 Us[k].resize(vectorSizes.at(k));
1595
1596 // Collect the solution data
1597 MPI_Request req_data[size];
1598 for (index_t k=0; k!=size; k++)
1599 m_comm.irecv(Us.at(k).data() ,vectorSizes.at(k) , sourceID, &req_data[k], tag++);
1600
1601 // When the solution data and the loads are collected, put them in the stepSolutions container.
1602 MPI_Waitall( size, req_data, MPI_STATUSES_IGNORE );
1603 MPI_Waitall( size, req_loads, MPI_STATUSES_IGNORE );
1604 for (index_t k=0; k!=size; k++)
1605 stepSolutions.at(k) = std::make_pair(Us.at(k),Ls.at(k));
1606
1607 // Wait for the distances to finish
1608 MPI_Waitall( size+1, req_distances, MPI_STATUSES_IGNORE );
1609 // Wait for the meta data to finish
1610 MPI_Waitall( 2, req_meta, MPI_STATUSES_IGNORE );
1611}
1612
1613template <class T>
1614void gsAPALM<T>::_recvWorkerToMain( index_t & sourceID,
1615 T & distance,
1616 std::vector<solution_t> & solutions,
1617 bool & bifurcation)
1618{
1619 index_t tag = 0;
1620 index_t size;
1621
1622 MPI_Status status;
1623 MPI_Request req_meta[2];
1624 m_comm.irecv(&distance ,1, sourceID, &req_meta[0], tag++);
1625 MPI_Wait ( &req_meta[0], &status );
1626 if (sourceID==MPI_ANY_SOURCE) sourceID = status.MPI_SOURCE; // Overwrite source ID to be unique
1627 gsMPIInfo(m_rank)<<"Receiving data on "<<m_rank<<" from "<<sourceID<<"\n";
1628
1629 // Request metadata
1630 m_comm.irecv(&bifurcation ,1, sourceID, &req_meta[1], tag++);
1631
1632 // Request the number of solution intervals
1633 MPI_Request req_solSize;
1634 m_comm.irecv(&size ,1, sourceID, &req_solSize, tag++);
1635
1636 MPI_Wait(&req_solSize,MPI_STATUS_IGNORE);
1637
1638 // temporary objects to store solutions and the sizes of the solution vectors
1639 std::vector<gsVector<T>> Us(size);
1640 std::vector<T > Ls(size);
1641 std::vector<index_t > vectorSizes(size);
1642 // final containers for distances and step solutions
1643 solutions.resize(size);
1644
1645 // Collect distances, solution vector sizes and loads
1646 MPI_Request req_vecSizes[size];
1647 MPI_Request req_loads[size];
1648 for (index_t k=0; k!=size; k++)
1649 {
1650 m_comm.irecv(&Ls.at(k) ,1, sourceID, &req_loads[k] , tag++);
1651 m_comm.irecv(&vectorSizes.at(k) ,1, sourceID, &req_vecSizes[k] , tag++);
1652 }
1653
1654 // When all solution vector sizes are collected, resize the solution vectors
1655 MPI_Waitall( size, req_vecSizes, MPI_STATUSES_IGNORE );
1656 for (index_t k=0; k!=size; k++)
1657 Us[k].resize(vectorSizes.at(k));
1658
1659 // Collect the solution data
1660 MPI_Request req_data[size];
1661 for (index_t k=0; k!=size; k++)
1662 m_comm.irecv(Us.at(k).data() ,vectorSizes.at(k) , sourceID, &req_data[k], tag++);
1663
1664 // When the solution data and the loads are collected, put them in the solutions container.
1665 MPI_Waitall( size, req_data, MPI_STATUSES_IGNORE );
1666 MPI_Waitall( size, req_loads, MPI_STATUSES_IGNORE );
1667 for (index_t k=0; k!=size; k++)
1668 solutions.at(k) = std::make_pair(Us.at(k),Ls.at(k));
1669
1670 // Wait for the meta data to finish
1671 MPI_Waitall( 2, req_meta, MPI_STATUSES_IGNORE );
1672}
1673
1674
1675#endif
1676
1677// -----------------------------------------------------------------------------------------------------
1678} // namespace gismo
index_t
#define index_t
Definition gsConfig.h:32
GISMO_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition gsDebug.h:129
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition gsDebug.h:102
gsInfo
#define gsInfo
Definition gsDebug.h:43
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::gsStatus
gsStatus
Definition gsStructuralAnalysisTypes.h:21
gismo::distance
T distance(gsMatrix< T > const &A, gsMatrix< T > const &B, index_t i=0, index_t j=0, bool cols=false)
compute a distance between the point number in the set and the point number <j> in the set ; by def...