1 /*
2 * MBDyn (C) is a multibody analysis code.
3 * http://www.mbdyn.org
4 *
5 * Copyright (C) 1996-2017
6 *
7 * Pierangelo Masarati <masarati@aero.polimi.it>
8 * Paolo Mantegazza <mantegazza@aero.polimi.it>
9 *
10 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
11 * via La Masa, 34 - 20156 Milano, Italy
12 * http://www.aero.polimi.it
13 *
14 * Changing this copyright notice is forbidden.
15 *
16 * This program is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation (version 2 of the License).
19 *
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program; if not, write to the Free Software
28 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29 */
30 /*
31 * Copyright (C) 2009
32 *
33 * Marco Morandini
34 *
35 */
36
37 /*****************************************************************************
38 * *
39 * ParNaive C++ WRAPPER *
40 * *
41 *****************************************************************************/
42
43 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
44
45 #ifdef USE_NAIVE_MULTITHREAD
46
47 /* FIXME: incompatible with RTAI at present */
48 #ifndef USE_RTAI
49
50 #include <sys/types.h>
51 #include <sys/stat.h>
52 #include <fcntl.h>
53 #include <sys/ioctl.h>
54 #include <cstdio>
55
56 #include <unistd.h>
57 #include <signal.h>
58
59 #include "parnaivewrap.h"
60 #include "mthrdslv.h"
61 #include "task2cpu.h"
62
63 #include "pmthrdslv.h"
64
65 /* ParNaiveSolver - begin */
66
67 /* Costruttore: si limita ad allocare la memoria */
ParNaiveSolver(unsigned nt,const integer & size,const doublereal & dMP,NaiveMatrixHandler * const a)68 ParNaiveSolver::ParNaiveSolver(unsigned nt, const integer &size,
69 const doublereal& dMP,
70 NaiveMatrixHandler *const a)
71 : LinearSolver(0),
72 iSize(size),
73 dMinPiv(dMP),
74 ppril(0),
75 pnril(0),
76 A(a),
77 nThreads(nt),
78 thread_data(0)
79 {
80 ASSERT(iN > 0);
81
82 piv.resize(iSize);
83 fwd.resize(iSize);
84 todo.resize(iSize);
85 row_locks.resize(iSize + 2);
86 col_locks.resize(iSize, AO_TS_INITIALIZER);
87
88 pthread_mutex_init(&thread_mutex, NULL);
89 pthread_cond_init(&thread_cond, NULL);
90
91 SAFENEWARR(ppril, integer *, iSize);
92 ppril[0] = 0;
93 SAFENEWARR(ppril[0], integer, iSize*iSize);
94 for (integer i = 1; i < iSize; i++) {
95 ppril[i] = ppril[i - 1] + iSize;
96 }
97 SAFENEWARR(pnril, integer, iSize);
98 #ifdef HAVE_MEMSET_H
99 memset(pnril, 0, sizeof(integer)*iSize);
100 #else /* ! HAVE_MEMSET_H */
101 for (integer row = 0; row < iSize; row++) {
102 pnril[row] = 0;
103 }
104 #endif /* ! HAVE_MEMSET_H */
105
106
107 SAFENEWARRNOFILL(thread_data, thread_data_t, nThreads);
108
109 (void)mbdyn_task2cpu(nThreads - 1);
110
111 for (unsigned t = 0; t < nThreads; t++) {
112 thread_data[t].pSLUS = this;
113 thread_data[t].threadNumber = t;
114 thread_data[t].retval = 0;
115
116 sem_init(&thread_data[t].sem, 0, 0);
117 if (pthread_create(&thread_data[t].thread, NULL, thread_op, &thread_data[t]) != 0) {
118 silent_cerr("ParNaiveSolver: pthread_create() failed "
119 "for thread " << t
120 << " of " << nThreads << std::endl);
121 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
122 }
123 }
124 }
125
126 /* Distruttore */
~ParNaiveSolver(void)127 ParNaiveSolver::~ParNaiveSolver(void)
128 {
129
130 thread_operation = ParNaiveSolver::EXIT;
131 thread_count = nThreads;
132
133 for (unsigned i = 0; i < nThreads; i++) {
134 sem_post(&thread_data[i].sem);
135 }
136
137 /* thread cleanup func? */
138
139 pthread_mutex_lock(&thread_mutex);
140 if (thread_count > 0) {
141 pthread_cond_wait(&thread_cond, &thread_mutex);
142 }
143 pthread_mutex_unlock(&thread_mutex);
144
145 for (unsigned i = 1; i < nThreads; i++) {
146 sem_destroy(&thread_data[i].sem);
147 }
148
149 pthread_mutex_destroy(&thread_mutex);
150 pthread_cond_destroy(&thread_cond);
151
152 if (ppril) {
153 if (ppril[0]) {
154 SAFEDELETEARR(ppril[0]);
155 }
156 SAFEDELETEARR(ppril);
157 }
158 if (pnril) {
159 SAFEDELETEARR(pnril);
160 }
161
162 /* other cleanup... */
163 }
164
165 void *
thread_op(void * arg)166 ParNaiveSolver::thread_op(void *arg)
167 {
168 thread_data_t *td = (thread_data_t *)arg;
169
170 silent_cout("ParNaiveSolver: thread " << td->threadNumber
171 << " [self=" << pthread_self()
172 << ",pid=" << getpid() << "]"
173 << " starting..." << std::endl);
174
175 /* deal with signals ... */
176 sigset_t newset /* , oldset */ ;
177 sigemptyset(&newset);
178 sigaddset(&newset, SIGTERM);
179 sigaddset(&newset, SIGINT);
180 sigaddset(&newset, SIGHUP);
181 pthread_sigmask(SIG_BLOCK, &newset, /* &oldset */ NULL);
182
183 (void)mbdyn_task2cpu(td->threadNumber);
184
185 bool bKeepGoing(true);
186
187 while (bKeepGoing) {
188 sem_wait(&td->sem);
189
190 switch (td->pSLUS->thread_operation) {
191 case ParNaiveSolver::FACTOR:
192 td->retval = pnaivfct(td->pSLUS->A->ppdRows,
193 td->pSLUS->A->iSize,
194 td->pSLUS->A->piNzr,
195 td->pSLUS->A->ppiRows,
196 td->pSLUS->A->piNzc,
197 td->pSLUS->A->ppiCols,
198 td->pSLUS->pnril,
199 td->pSLUS->ppril,
200 td->pSLUS->A->ppnonzero,
201 &td->pSLUS->piv[0],
202 &td->pSLUS->todo[0],
203 td->pSLUS->dMinPiv,
204 &td->pSLUS->row_locks[0],
205 &td->pSLUS->col_locks[0],
206 td->threadNumber,
207 td->pSLUS->nThreads
208 );
209
210 if (td->retval) {
211 if (td->retval & NAIVE_ENULCOL) {
212 silent_cerr("NaiveSolver: NAIVE_ENULCOL("
213 << (td->retval & ~NAIVE_ENULCOL) << ")" << std::endl);
214 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
215 }
216
217 if (td->retval & NAIVE_ENOPIV) {
218 silent_cerr("NaiveSolver: NAIVE_ENOPIV("
219 << (td->retval & ~NAIVE_ENOPIV) << ")" << std::endl);
220 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
221 }
222
223 /* default */
224 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
225 }
226 break;
227
228 case ParNaiveSolver::SOLVE:
229 pnaivslv(td->pSLUS->A->ppdRows,
230 td->pSLUS->A->iSize,
231 td->pSLUS->A->piNzc,
232 td->pSLUS->A->ppiCols,
233 td->pSLUS->LinearSolver::pdRhs,
234 &td->pSLUS->piv[0],
235 &td->pSLUS->fwd[0],
236 td->pSLUS->LinearSolver::pdSol,
237 &td->pSLUS->row_locks[0],
238 td->threadNumber,
239 td->pSLUS->nThreads
240 );
241 break;
242
243 case ParNaiveSolver::EXIT:
244 bKeepGoing = false;
245 break;
246
247 default:
248 silent_cerr("ParNaiveSolver: unhandled op"
249 << std::endl);
250 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
251 }
252
253 td->pSLUS->EndOfOp();
254 }
255
256 pthread_exit(NULL);
257 }
258
259 void
EndOfOp(void)260 ParNaiveSolver::EndOfOp(void)
261 {
262 bool last;
263
264 /* decrement the thread counter */
265 pthread_mutex_lock(&thread_mutex);
266 thread_count--;
267 last = (thread_count == 0);
268
269 /* if last thread, signal to restart */
270 if (last) {
271 #if 0
272 pthread_cond_broadcast(&thread_cond);
273 #else
274 pthread_cond_signal(&thread_cond);
275 #endif
276 }
277 pthread_mutex_unlock(&thread_mutex);
278 }
279
280 #ifdef DEBUG
281 void
IsValid(void) const282 ParNaiveSolver::IsValid(void) const
283 {
284 ASSERT(Aip != NULL);
285 ASSERT(App != NULL);
286 ASSERT(Axp != NULL);
287 ASSERT(iN > 0);
288 }
289 #endif /* DEBUG */
290
291 /* Fattorizza la matrice */
292 void
Factor(void)293 ParNaiveSolver::Factor(void)
294 {
295 #ifdef DEBUG
296 IsValid();
297 #endif /* DEBUG */
298
299 ASSERT(iNonZeroes > 0);
300
301 thread_operation = ParNaiveSolver::FACTOR;
302 thread_count = nThreads;
303
304 for (int i = 0; i < iSize; i++) {
305 piv[i] = -1;
306 todo[i] = -1;
307 row_locks[i] = 0;
308 pnril[0] = 0;
309 }
310
311 /* NOTE: these are for pivot_lock and sync_lock */
312 /* FIXME: bottleneck? */
313 row_locks[iSize] = 0;
314 row_locks[iSize + 1] = 0;
315
316 /* NOTE: no need to reset col_locks because they're
317 * always left equal to zero after use */
318
319 for (unsigned t = 0; t < nThreads; t++) {
320 thread_data[t].retval = 0;
321 sem_post(&thread_data[t].sem);
322 }
323
324 pthread_mutex_lock(&thread_mutex);
325 if (thread_count > 0) {
326 pthread_cond_wait(&thread_cond, &thread_mutex);
327 }
328 pthread_mutex_unlock(&thread_mutex);
329
330 }
331
332 /* Risolve */
333 void
Solve(void) const334 ParNaiveSolver::Solve(void) const
335 {
336 #ifdef DEBUG
337 IsValid();
338 #endif /* DEBUG */
339
340 if (bHasBeenReset) {
341 const_cast<ParNaiveSolver *>(this)->Factor();
342 bHasBeenReset = false;
343 }
344
345 for (int i = 0; i < iSize + 2; i++) {
346 row_locks[i] = 0;
347 }
348
349 thread_operation = ParNaiveSolver::SOLVE;
350 thread_count = nThreads;
351
352 for (unsigned t = 0; t < nThreads; t++) {
353 thread_data[t].retval = 0;
354 sem_post(&thread_data[t].sem);
355 }
356
357
358 pthread_mutex_lock(&thread_mutex);
359 if (thread_count > 0) {
360 pthread_cond_wait(&thread_cond, &thread_mutex);
361 }
362 pthread_mutex_unlock(&thread_mutex);
363 }
364
365 void
SetMat(NaiveMatrixHandler * const a)366 ParNaiveSolver::SetMat(NaiveMatrixHandler *const a)
367 {
368 A = a;
369 }
370
371 /* ParNaiveSolver - end */
372
373
374 /* ParNaiveSparseSolutionManager - begin: code */
375
376 /* Costruttore */
ParNaiveSparseSolutionManager(unsigned nt,const integer Dim,const doublereal dMP)377 ParNaiveSparseSolutionManager::ParNaiveSparseSolutionManager(unsigned nt,
378 const integer Dim, const doublereal dMP)
379 : A(0),
380 VH(Dim),
381 XH(Dim)
382 {
383 ASSERT(Dim > 0);
384 ASSERT((dMP >= 0.0) && (dMP <= 1.0));
385
386 SAFENEWWITHCONSTRUCTOR(A, NaiveMatrixHandler, NaiveMatrixHandler(Dim));
387 SAFENEWWITHCONSTRUCTOR(SolutionManager::pLS,
388 ParNaiveSolver,
389 ParNaiveSolver(nt, Dim, dMP, A));
390
391 pLS->pdSetResVec(VH.pdGetVec());
392 pLS->pdSetSolVec(XH.pdGetVec());
393 pLS->SetSolutionManager(this);
394
395 #ifdef DEBUG
396 IsValid();
397 #endif /* DEBUG */
398 }
399
400
401 /* Distruttore; verificare la distruzione degli oggetti piu' complessi */
~ParNaiveSparseSolutionManager(void)402 ParNaiveSparseSolutionManager::~ParNaiveSparseSolutionManager(void)
403 {
404 #ifdef DEBUG
405 IsValid();
406 #endif /* DEBUG */
407
408 /* Dealloca roba, ... */
409 if (A != 0) {
410 SAFEDELETE(A);
411 A = 0;
412 }
413
414 /* ... tra cui i thread */
415 }
416
417 #ifdef DEBUG
418 /* Test di validita' del manager */
419 void
IsValid(void) const420 ParNaiveSparseSolutionManager::IsValid(void) const
421 {
422 ASSERT(iMatSize > 0);
423
424 #ifdef DEBUG_MEMMANAGER
425 ASSERT(defaultMemoryManager.fIsPointerToBlock(VH.pdGetVec()));
426 ASSERT(defaultMemoryManager.fIsPointerToBlock(XH.pdGetVec()));
427 ASSERT(defaultMemoryManager.fIsPointerToBlock(pLS));
428 #endif /* DEBUG_MEMMANAGER */
429
430 ASSERT((VH.IsValid(), 1));
431 ASSERT((XH.IsValid(), 1));
432 ASSERT((pLS->IsValid(), 1));
433 }
434 #endif /* DEBUG */
435
436 /* Inizializza il gestore delle matrici */
437 void
MatrReset(void)438 ParNaiveSparseSolutionManager::MatrReset(void)
439 {
440 #ifdef DEBUG
441 IsValid();
442 #endif /* DEBUG */
443
444 pLS->Reset();
445 }
446
447
448 /* Risolve il problema */
449 void
Solve(void)450 ParNaiveSparseSolutionManager::Solve(void)
451 {
452 #ifdef DEBUG
453 IsValid();
454 #endif /* DEBUG */
455
456
457 pLS->Solve();
458
459 }
460
461 /* Rende disponibile l'handler per la matrice */
462 MatrixHandler*
pMatHdl(void) const463 ParNaiveSparseSolutionManager::pMatHdl(void) const
464 {
465 return A;
466 }
467
468 VectorHandler*
pResHdl(void) const469 ParNaiveSparseSolutionManager::pResHdl(void) const
470 {
471 #ifdef DEBUG
472 VH.IsValid();
473 #endif /* DEBUG */
474 return &VH;
475 }
476
477 /* Rende disponibile l'handler per la soluzione */
478 VectorHandler*
pSolHdl(void) const479 ParNaiveSparseSolutionManager::pSolHdl(void) const
480 {
481 #ifdef DEBUG
482 XH.IsValid();
483 #endif /* DEBUG */
484 return &XH;
485 }
486
487 /* ParNaiveSparseSolutionManager - end */
488
489
490 /* ParNaiveSparsePermSolutionManager - begin */
491
492 extern "C" {
493 #include "colamd.h"
494 }
495
ParNaiveSparsePermSolutionManager(unsigned nt,const integer Dim,const doublereal dMP)496 ParNaiveSparsePermSolutionManager::ParNaiveSparsePermSolutionManager(
497 unsigned nt,
498 const integer Dim,
499 const doublereal dMP)
500 : ParNaiveSparseSolutionManager(nt, Dim, dMP),
501 dMinPiv(dMP),
502 ePermState(PERM_NO)
503 {
504 perm.resize(Dim, 0);
505 invperm.resize(Dim, 0);
506
507 SAFEDELETE(A);
508 A = 0;
509 SAFENEWWITHCONSTRUCTOR(A, NaivePermMatrixHandler, NaivePermMatrixHandler(Dim, perm, invperm));
510
511 dynamic_cast<ParNaiveSolver *>(pLS)->SetMat(A);
512
513 MatrInitialize();
514 }
515
~ParNaiveSparsePermSolutionManager(void)516 ParNaiveSparsePermSolutionManager::~ParNaiveSparsePermSolutionManager(void)
517 {
518 NO_OP;
519 }
520
521 void
MatrReset(void)522 ParNaiveSparsePermSolutionManager::MatrReset(void)
523 {
524 if (ePermState == PERM_INTERMEDIATE) {
525 ePermState = PERM_READY;
526
527 pLS->pdSetResVec(VH.pdGetVec());
528 pLS->pdSetSolVec(XH.pdGetVec());
529
530 pLS->SetSolutionManager(this);
531 }
532
533 pLS->Reset();
534 }
535
536 void
ComputePermutation(void)537 ParNaiveSparsePermSolutionManager::ComputePermutation(void)
538 {
539 std::vector<integer> Ai;
540 A->MakeCCStructure(Ai, invperm);
541 doublereal knobs [COLAMD_KNOBS];
542 integer stats [COLAMD_STATS];
543 integer Alen = mbdyn_colamd_recommended (Ai.size(), A->iGetNumRows(), A->iGetNumCols());
544 Ai.resize(Alen);
545 mbdyn_colamd_set_defaults(knobs);
546 if (!mbdyn_colamd(A->iGetNumRows(), A->iGetNumCols(), Alen,
547 &(Ai[0]), &(invperm[0]), knobs, stats)) {
548 silent_cerr("colamd permutation failed" << std::endl);
549 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
550 }
551 for (integer i = 0; i < A->iGetNumRows(); i++) {
552 perm[invperm[i]] = i;
553 }
554 ePermState = PERM_INTERMEDIATE;
555 }
556
557 void
BackPerm(void)558 ParNaiveSparsePermSolutionManager::BackPerm(void)
559 {
560 for (integer i = 0; i < A->iGetNumCols(); i++) {
561 XH(invperm[i] + 1) = VH(i + 1);
562 }
563 }
564
565
566 /* Risolve il sistema Fattorizzazione + Backward Substitution */
567 void
Solve(void)568 ParNaiveSparsePermSolutionManager::Solve(void)
569 {
570 if (ePermState == PERM_NO) {
571 ComputePermutation();
572
573 } else {
574 pLS->pdSetSolVec(VH.pdGetVec());
575 }
576
577 pLS->Solve();
578
579 if (ePermState == PERM_READY) {
580 BackPerm();
581 pLS->pdSetSolVec(XH.pdGetVec());
582 }
583 }
584
585 /* Inizializzatore "speciale" */
586 void
MatrInitialize()587 ParNaiveSparsePermSolutionManager::MatrInitialize()
588 {
589 ePermState = PERM_NO;
590
591 for (integer i = 0; i < A->iGetNumRows(); i++) {
592 perm[i] = i;
593 invperm[i] = i;
594 }
595
596 MatrReset();
597 }
598
599 /* ParParNaiveSparsePermSolutionManager - end */
600
601
602 #endif /* ! USE_RTAI */
603
604 #endif /* USE_NAIVE_MULTITHREAD */
605