MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
schsolman.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbmath/schsolman.cc,v 1.60 2017/01/12 14:43:54 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  * Paolo Mantegazza <mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  */
31 
32 /* SchurMatrixHandler
33  SchurVectorHandler
34  SchurSolutionManager */
35 
36 /*
37  * Copyright 1999-2017 Giuseppe Quaranta <quaranta@aero.polimi.it>
38  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
39  */
40 
41 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
42 
43 #include <iomanip>
44 
45 #ifdef USE_MPI
46 #include "solman.h"
47 #include "parser.h"
48 #include "linsol.h"
49 #include "schsolman.h"
50 #include "mysleep.h"
51 
52 /* FIXME: we must move the test to libmbmath */
53 #include "nonlin.h"
54 
55 /* NOTE: define to use Wait/Waitall instead of Test/Testall
56  * Apparently, this results in far better performances,
57  * so we might want to extend it to all other communications */
58 #define USE_MPI_WAIT
59 
60 /*
61  * Le occorrenze della routine 'mysleep' vanno commentate quando
62  * si opera su di una macchina dove le comunicazioni sono efficienti
63  * e non ci sono altri utenti operanti che esauriscono le risorse.
64  * In caso contrario la si utilizza per ottenere delle misure significative
65  * di tempi di CPU, anche se i tempi solari ottenibili non sono proporzionati
66  */
67 
68 #ifdef MPI_PROFILING
69 extern "C" {
70 #include <mpe.h>
71 #include <stdio.h>
72 };
73 #endif /* MPI_PROFILING */
74 
75 #undef min
76 #undef max
77 #include <algorithm>
78 
80  integer iBlocks,
81  integer* pLocalDofs,
82  int iDim1,
83  integer* pInterfDofs,
84  int iDim2,
85  SolutionManager* pLSM,
86  LinSol &ls)
87 :
88 SolvCommSize(0),
89 iPrbmSize(iSize),
90 iPrbmBlocks(iBlocks),
91 iBlkSize(iSize),
92 pLocDofs(pLocalDofs),
93 pIntDofs(pInterfDofs),
94 iLocVecDim(iDim1),
95 iIntVecDim(iDim2),
96 pRecvDim(NULL),
97 pDispl(NULL),
98 pDofsRecvdList(NULL),
99 pSchurDofs(NULL),
100 iSchurIntDim(0),
101 pGlbToLoc(NULL),
102 pSchGlbToLoc(NULL),
103 pBlockLenght(NULL),
104 pTypeDsp(NULL),
105 ppNewTypes(NULL),
106 pBuffer(NULL),
107 pMH(NULL),
108 pRVH(NULL),
109 pSolVH(NULL),
110 pSchMH(NULL),
111 pSchVH(NULL),
112 pSolSchVH(NULL),
113 pBMH(NULL),
114 pdCM(NULL),
115 prVH(NULL),
116 pgVH(NULL),
117 pSolrVH(NULL),
118 pGSReq(NULL),
119 pGRReq(NULL),
120 pLocalSM(pLSM),
121 pInterSM(0),
122 bNewMatrix(false)
123 {
124  DEBUGCOUT("Entering SchurSolutionManager::SchurSolutionManager()"
125  << std::endl);
126 
127  ASSERT(iPrbmSize > 0);
128  ASSERT(pLocDofs != NULL);
129 
130  /* inizializzo il communicator */
131  SolvComm = MBDynComm.Dup();
132  SolvCommSize = SolvComm.Get_size();
133  MyRank = SolvComm.Get_rank();
134 
135  if (SolvCommSize <= 1) {
136  silent_cerr("SchurSolutionManager: "
137  "invalid communicator size " << SolvCommSize
138  << std::endl);
140  }
141 
142 #ifdef DEBUG
143  ASSERT(pIntDofs != NULL);
144 #endif /* DEBUG */
145 
146  DEBUGCOUT("Solution Communicator Size: " << SolvCommSize << std::endl);
147 
148  /* determina le dimensioni locali e il vettore
149  * di trasferimento Globale -> Locale;
150  * sul master (MyRank == 0) la matrice C e' quella
151  * di tutte le interfacce */
152 
153  /* initializes iSchurIntDim among other stuff... */
154  InitializeComm();
155 
156  /* utilizza iWorkSpaceSize come un coefficiente moltiplicativo */
157 
158  /* NOTE: iWorkSpaceSize is likely to be 0, since now most of the
159  * linear solvers don't need any; this is not a big deal, resulting
160  * in GetSolutionManager() being called with a 0 iWorkSpaceSize arg */
161  iWorkSpaceSize = ls.iGetWorkSpaceSize();
162  integer iIntWorkSpaceSize = iWorkSpaceSize*(iSchurIntDim*iSchurIntDim)/(iPrbmSize*iPrbmSize);
163 
164  if (MyRank == 0) {
165  pInterSM = ls.GetSolutionManager(iSchurIntDim,
166  iIntWorkSpaceSize);
167  }
168 
169  /* estrae i puntatori alle matrici e ai vettori creati
170  * dai solutori per l'assemblaggio degli SchurHandlers */
171  pBMH = pLocalSM->pMatHdl();
172  prVH = pLocalSM->pResHdl();
173  pSolrVH = pLocalSM->pSolHdl();
174 
175  if (MyRank == 0) {
176  pSchMH = pInterSM->pMatHdl();
177  pSchVH = pInterSM->pResHdl();
178  pSolSchVH = pInterSM->pSolHdl();
179  }
180 
181  /* Definisce le matrici globali */
182  if (prVH->pdGetVec() == pSolrVH->pdGetVec()) {
184  SchurMatrixHandler(iLocVecDim,
185  iIntVecDim, pBMH, pGlbToLoc));
186 
187  } else {
189  SchurMatrixHandlerUm(iLocVecDim,
190  iIntVecDim, pBMH, pGlbToLoc));
191  }
192 
194  SchurVectorHandler(iLocVecDim,
195  iIntVecDim, prVH, pGlbToLoc));
196 
197  pdCM = pMH->GetCMat();
198 
199  /* NOTE: from my understanding, pgVH is the interface portion
200  * of the pRVH (the SchurVectorhandler of the residual), which
201  * is allocated internally by the constructor; it is passed
202  * to pSolVH (the SchurVectorhandler of the solution) as well,
203  * to act as interface portion of the solution, so the residual
204  * and the solution share the same memory area. */
205  pgVH = pRVH->GetIVec();
206 
208  SchurVectorHandler(iLocVecDim,
209  iIntVecDim, pSolrVH, pgVH, pGlbToLoc));
210 
211  /* Creazione di NewType per la trasmissione
212  * del vettore soluzione calcolato */
213  if (MyRank == 0) {
214  SAFENEWARR(pBlockLenght, int, pDispl[SolvCommSize]);
215  InitializeList(pBlockLenght, pDispl[SolvCommSize], 1);
216  SAFENEWARR(pTypeDsp, MPI::Aint, pDispl[SolvCommSize]);
217 
218  doublereal *pdm1 = pSolSchVH->pdGetVec() - 1;
219  MPI::Aint DispTmp = MPI::Get_address(pdm1 + 1);
220  for (int i = 0; i < pDispl[SolvCommSize]; i++) {
221  int j = i%iBlkSize;
222  int blk = i/iBlkSize;
223  pTypeDsp[i] = MPI::Get_address(&pdm1[pSchGlbToLoc[pDofsRecvdList[j] + blk*iBlkSize]]) - DispTmp;
224  }
225 
226  SAFENEWARR(ppNewTypes, MPI::Datatype*, SolvCommSize);
227  MPI::Aint* pCurrDispl = pTypeDsp;
228 
229  for (int i = 0; i < SolvCommSize; i++) {
230  ppNewTypes[i] = NULL;
231  SAFENEWWITHCONSTRUCTOR(ppNewTypes[i], MPI::Datatype,
232  MPI::Datatype(MPI::DOUBLE.Create_hindexed(pRecvDim[i],
233  pBlockLenght, pCurrDispl)));
234  ppNewTypes[i]->Commit();
235  pCurrDispl += pRecvDim[i];
236  }
237  }
238 
239  silent_cout("Local dimension on process " << MyRank
240  << ": " << iLocVecDim << std::endl);
241 
242  if (MyRank == 0) {
243  silent_cout("Interface dimension: " << iSchurIntDim
244  << std::endl);
245  }
246 
247 #ifdef DEBUG
248  if (MyRank == 0) {
249  silent_cout("Interface Dofs " << std::endl);
250  for (int i = 0; i < iSchurIntDim; i++) {
251  silent_cout(pSchurDofs[i] << " ");
252  }
253  silent_cout(std::endl);
254  }
255 
256  IsValid();
257 #endif /* DEBUG */
258 }
259 
261 {
262  DEBUGCOUT("Entering SchurSolutionManager::~SchurSolutionManager()"
263  << endl);
264 
265 #ifdef DEBUG
266  IsValid();
267 #endif /* DEBUG */
268 
269  if (pRecvDim != NULL) {
271  }
272 
273  if (pDispl != NULL) {
274  SAFEDELETEARR(pDispl);
275  }
276 
277  if (pDofsRecvdList != NULL) {
279  }
280 
281  if (pSchurDofs != NULL) {
283  }
284 
285  if (pGlbToLoc != NULL) {
287  }
288 
289  if (pSchGlbToLoc != NULL) {
291  }
292 
293  if (pBlockLenght != NULL) {
295  }
296 
297  if (pTypeDsp != NULL) {
298  SAFEDELETEARR(pTypeDsp);
299  }
300 
301  if (pBuffer != NULL) {
303  }
304 
305  if (ppNewTypes != NULL) {
306  for (int i = 0; i < SolvCommSize; i++) {
307  ppNewTypes[i]->Free();
308  SAFEDELETE(ppNewTypes[i]);
309  }
310 
311  SAFEDELETEARR(ppNewTypes);
312  }
313 
314  if (pMH != NULL) {
315  SAFEDELETE(pMH);
316  }
317 
318  if (pRVH != NULL) {
319  SAFEDELETE(pRVH);
320  }
321 
322  if (pSolVH != NULL) {
324  }
325 
326  if (pInterSM != NULL) {
328  }
329 
330  if (pGSReq != NULL) {
331  SAFEDELETEARR(pGSReq);
332  }
333 
334  if (pGRReq != NULL) {
335  SAFEDELETEARR(pGRReq);
336  }
337 }
338 
339 #ifdef DEBUG
340 void
341 SchurSolutionManager::IsValid(void) const
342 {
343  NO_OP;
344 }
345 #endif /* DEBUG */
346 
347 /* Inizializza il gestore delle matrici */
348 
349 void
351 {
352 #ifdef DEBUG
353  IsValid();
354 #endif /* DEBUG */
355 
356  pLocalSM->MatrReset();
357  pMH->MatEFCReset();
358  if (MyRank == 0) {
359  pInterSM->MatrReset();
360  }
361  bNewMatrix = true;
362 }
363 
364 /* Inizializzatore "speciale" */
365 void
367 {
369  pMH->MatEFCReset();
370  if (MyRank == 0) {
372  }
373  bNewMatrix = true;
374 }
375 
376 /* Risolve i blocchi */
377 
378 void
380 {
381  DEBUGCOUT("Entering SchurSolutionManager::Solve()" << endl);
382 #ifdef DEBUG
383  IsValid();
384 #endif /* DEBUG */
385 #ifdef MPI_PROFILING
386  MPE_Log_event(31, 0, "start");
387 #endif /* MPI_PROFILING */
388 
389  /* Fattorizzazione matrice B */
390  pLocalSM->Solve();
391 
392  /* NOTE: we need to set the matrix handler after each solution,
393  * because when the local SM is using an optimized sparse matrix,
394  * it may switch from a generic representation, e.g. SpMapMH,
395  * to a compressed representation, e.g. CCMH or DirMH
396  */
398 
399 #ifdef MPI_PROFILING
400  MPE_Log_event(32, 0, "end");
401 #endif /* MPI_PROFILING */
402 
403  /* prodotto g = g - F*r */
404  pMH->CompNewg(*pgVH, *pSolrVH);
405 
406 #ifdef MPI_PROFILING
407  MPE_Log_event(13, 0, "start");
408  MPE_Log_send(0, G_TAG, iIntVecDim);
409 #endif /* MPI_PROFILING */
410 
411  /* on master node... */
412  if (MyRank == 0) {
413 #ifdef MPI_PROFILING
414  MPE_Log_event(19, 0, "start");
415 #endif /* MPI_PROFILING */
416 
417  /* receive from remote only... */
418  for (int i = 1; i < SolvCommSize; i++){
419  pGRReq[i] = SolvComm.Irecv(&pBuffer[pDispl[i]],
420  pRecvDim[i], MPI::DOUBLE, i, G_TAG);
421  }
422 
423  /* emulate transmission...*/
424  (void)memmove(&pBuffer[pDispl[0]], pgVH->pdGetVec(),
425  sizeof(double)*pRecvDim[0]);
426 
427  /* on remote nodes... */
428  } else {
429  /* Inizializza Trasmissione di g */
430  pGSReq[0] = SolvComm.Isend(pgVH->pdGetVec(), iIntVecDim,
431  MPI::DOUBLE, 0, G_TAG);
432  }
433 
434  /* Calcolo di E'. Va fatto solo se la matrice e' stata rifattorizzata */
435  if (bNewMatrix) {
436  for (int i = 0; i < iIntVecDim; i++) {
439  /* fa solo la back substitution perche'
440  * la fattorizzazione e' gia' stata lanciata */
441  pLocalSM->Solve();
442  }
443  }
444 
447 
448  /* on master node... */
449  if (MyRank == 0) {
450 #ifdef USE_MPI_WAIT
451  MPI::Request::Waitall(SolvCommSize - 1, &pGRReq[1]);
452 #else /* ! USE_MPI_WAIT */
453  while (true) {
454  /* testing remote only... */
455  if (MPI::Request::Testall(SolvCommSize - 1, &pGRReq[1]))
456  {
457 #ifdef MPI_PROFILING
458  MPE_Log_event(20, 0, "end");
459  for (int i = 1; i < SolvCommSize; i++){
460  MPE_Log_receive(i, G_TAG, pRecvDim[i]);
461  }
462 #endif /* MPI_PROFILING */
463  break;
464  }
465  MYSLEEP(1500);
466  }
467 #endif /* ! USE_MPI_WAIT */
468 
469  /* Assemblaggio del vettore delle incognite g
470  * sulla macchina master */
471 
472 #ifdef DEBUG_SCHUR
473  std::cout << "# Solve 19:";
474  for (int i = 0; i < pDispl[SolvCommSize]; i++) {
475  std::cout << " " << pBuffer[i];
476  }
477  std::cout << std::endl;
478 #endif /* DEBUG_SCHUR */
479 
480  /* Assembla pSchVH */
481  pSchVH->Reset();
482  if (iPrbmBlocks == 1) {
483  for (int i = 0; i < pDispl[SolvCommSize]; i++) {
485  }
486 
487  } else {
488  for (int i = 0; i < pDispl[SolvCommSize]; i++) {
489  int j = i%iBlkSize;
490  int blk = i/iBlkSize;
491 
493  }
494  }
495 
496  /* on remote nodes... */
497  } else {
498 #ifdef USE_MPI_WAIT
499  pGSReq[0].Wait();
500 #else /* ! USE_MPI_WAIT */
501  /* verifica completamento ricezioni e trasmissione vettore g */
502  while (true) {
503  if (pGSReq[0].Test()) {
504 #ifdef MPI_PROFILING
505  MPE_Log_event(14, 0, "end");
506 #endif /* MPI_PROFILING */
507  break;
508  }
509  MYSLEEP(150);
510  }
511 #endif /* ! USE_MPI_WAIT */
512  }
513 
514  /* assembla le matrici di schur locali se e' stata
515  * appena assemblata una nuova matrice di partenza */
516  if (bNewMatrix) {
517  AssSchur();
518  bNewMatrix = false;
519  }
520 
521  /* on master node... */
522  if (MyRank == 0) {
523  /* risoluzione schur pb sul master */
524 
525 #ifdef MPI_PROFILING
526  MPE_Log_event(35, 0, "start");
527 #endif /* MPI_PROFILING */
528 
529  pInterSM->Solve();
530 
531 #ifdef MPI_PROFILING
532  MPE_Log_event(36, 0, "end");
533 #endif /* MPI_PROFILING */
534 
535  /* invia la soluzione calcolata */
536 #ifdef MPI_PROFILING
537  MPE_Log_event(13, 0, "start");
538 #endif /* MPI_PROFILING */
539  doublereal *pd = pSolSchVH->pdGetVec();
540  for (int i = 1; i < SolvCommSize; i++) {
541  /* sending to remote only... */
542  pGRReq[i] = SolvComm.Isend(pd, 1,
543  *ppNewTypes[i], i, G_TAG);
544 #ifdef MPI_PROFILING
545  MPE_Log_send(i, G_TAG, pRecvDim[i]);
546 #endif /* MPI_PROFILING */
547  }
548 
549  /* emulate local transmission... */
550  memmove(pgVH->pdGetVec(), pd, sizeof(double)*pRecvDim[0]);
551 
552 #ifdef USE_MPI_WAIT
553  MPI::Request::Waitall(SolvCommSize - 1, &pGRReq[1]);
554 #else /* ! USE_MPI_WAIT */
555  /* Verifica completamento trasmissioni */
556  while (true) {
557  /* testing remote only... */
558  if (MPI::Request::Testall(SolvCommSize - 1, &pGRReq[1]))
559  {
560 #ifdef MPI_PROFILING
561  MPE_Log_event(14, 0, "end");
562 #endif /* MPI_PROFILING */
563  break;
564  }
565  MYSLEEP(150);
566  }
567 #endif /* ! USE_MPI_WAIT */
568 
569  /* on other nodes... */
570  } else {
571 #ifdef MPI_PROFILING
572  MPE_Log_event(19, 0, "start");
573 #endif /* MPI_PROFILING */
574 
575  pGSReq[0] = SolvComm.Irecv(pgVH->pdGetVec(), iIntVecDim,
576  MPI::DOUBLE, 0, G_TAG);
577 
578 #ifdef USE_MPI_WAIT
579  pGSReq[0].Wait();
580 #else /* ! USE_MPI_WAIT */
581  while (true) {
582  if (pGSReq[0].Test()) {
583 #ifdef MPI_PROFILING
584  MPE_Log_event(20, 0, "end");
585  MPE_Log_receive(0, G_TAG, iIntVecDim);
586 #endif /* MPI_PROFILING */
587  break;
588  }
589  MYSLEEP(1500);
590  }
591 #endif /* ! USE_MPI_WAIT */
592  }
593 
594  /* calcolo la soluzione corretta per x = Solr - E'*g */
595  pMH->CompNewf(*pSolrVH, *pgVH);
596 
597 #ifdef DEBUG_SCHUR
598  std::cout << "# Solve 20a:";
599  for (int i = 1; i <= iIntVecDim; i++) {
600  std::cout << " " << pgVH->operator()(i);
601  }
602  std::cout << std::endl;
603  std::cout << "# Solve 20b:";
604  for (int i = 1; i <= iLocVecDim; i++) {
605  std::cout << " " << pSolrVH->operator()(i);
606  }
607  std::cout << std::endl;
608 #endif /* DEBUG_SCHUR */
609 }
610 
611 void
613 {
614 #ifdef MPI_PROFILING
615  MPE_Log_event(41, 0, "start");
616 #endif /* MPI_PROFILING */
617 
618  DEBUGCOUT("Entering SchurSolutionManager::AssSchur()" << endl);
619 
620 #ifdef DEBUG
621  IsValid();
622 #endif /* DEBUG */
623 
624  pMH->CompLocSchur();
625 
626  /* on local node... */
627  if (MyRank == 0) {
628 #ifdef MPI_PROFILING
629  MPE_Log_event(19, 0, "start");
630 #endif /* MPI_PROFILING */
631 
632  /* remote nodes contribution... */
633  int iOffset = pRecvDim[0]*pRecvDim[0];
634  for (int i = 1; i < SolvCommSize; i++) {
635  int iShift = pRecvDim[i]*pRecvDim[i];
636 
637  pGRReq[i] = SolvComm.Irecv(&pBuffer[iOffset],
638  iShift, MPI::DOUBLE, i, S_TAG);
639 
640  iOffset += iShift;
641  }
642 
643  /* local node contribution... */
644  (void)memmove(&pBuffer[pDispl[0]], pdCM,
645  sizeof(double)*pRecvDim[0]*pRecvDim[0]);
646 
647 #ifdef USE_MPI_WAIT
648  MPI::Request::Waitall(SolvCommSize - 1, &pGRReq[1]);
649 #else /* ! USE_MPI_WAIT */
650  while (true) {
651  if (MPI::Request::Testall(SolvCommSize - 1, &pGRReq[1]))
652  {
653 #ifdef MPI_PROFILING
654  MPE_Log_event(20, 0, "end");
655  for (int i = 1; i < SolvCommSize; i++){
656  MPE_Log_receive(i, S_TAG, pRecvDim[i]*pRecvDim[i]);
657  }
658 #endif /* MPI_PROFILING */
659  break;
660  }
661  MYSLEEP(1500);
662  }
663 #endif /* ! USE_MPI_WAIT */
664 
665  /* Assembla Schur matrix */
666  pSchMH->Reset();
667  doublereal *pbmd = &pBuffer[-pDispl[0]];
668  for (int i = 0; i < SolvCommSize; i++) {
669  if (iPrbmBlocks == 1) {
670  for (int j = pDispl[i]; j < pDispl[i + 1]; j++) {
671  int idxC = pSchGlbToLoc[pDofsRecvdList[j]];
672 
673  for (int k = pDispl[i]; k < pDispl[i + 1]; k++) {
674  int idxR = pSchGlbToLoc[pDofsRecvdList[k]];
675 
676  pSchMH->IncCoef(idxR, idxC, pbmd[k]);
677  }
678 
679  pbmd += pRecvDim[i];
680  }
681 
682  } else {
683  for (int j = pDispl[i]; j < pDispl[i + 1]; j++) {
684 
685  int ic = j%iBlkSize;
686  int bc = j/iBlkSize;
687  int idxC = pSchGlbToLoc[pDofsRecvdList[ic] + bc*iBlkSize];
688 
689  for (int k = pDispl[i]; k < pDispl[i + 1]; k++) {
690  int ir = k%iBlkSize;
691  int br = k/iBlkSize;
692  int idxR = pSchGlbToLoc[pDofsRecvdList[ir] + br*iBlkSize];
693 
694  pSchMH->IncCoef(idxR, idxC, pbmd[k]);
695  }
696 
697  pbmd += pRecvDim[i];
698  }
699  }
700 
701  /* NOTE: this is required to allow the innermost loop
702  * on k to be between pDispl[i] and pDispl[i + 1];
703  * as a consequence, pbmd must be decreased
704  * by pDispl[i] at each i.
705  * Note that, in general, pDispl[i + 1] - pDispl[i]
706  * is equal to pRecvDim[i] */
707  pbmd -= pDispl[i + 1] - pDispl[i];
708  }
709 
710  /* on other nodes... */
711  } else {
712  /* Trasmette le Schur locali */
713 #ifdef MPI_PROFILING
714  MPE_Log_event(13, 0, "start");
715  MPE_Log_send(0, S_TAG, iIntVecDim*iIntVecDim);
716 #endif /* MPI_PROFILING */
717 
718  pGSReq[0] = SolvComm.Isend(pdCM, iIntVecDim*iIntVecDim,
719  MPI::DOUBLE, 0, S_TAG);
720 
721 #ifdef USE_MPI_WAIT
722  pGSReq[0].Wait();
723 #else /* !USE_MPI_WAIT */
724  /* verifica completamento ricezioni e trasmissione */
725  while (true) {
726  if (pGSReq[0].Test()) {
727 #ifdef MPI_PROFILING
728  MPE_Log_event(14, 0, "end");
729 #endif /* MPI_PROFILING */
730  break;
731  }
732  MYSLEEP(150);
733  }
734 #endif /* !USE_MPI_WAIT */
735  }
736 
737 #ifdef MPI_PROFILING
738  MPE_Log_event(42, 0, "end");
739 #endif /* MPI_PROFILING */
740 }
741 
742 void
744 {
745  DEBUGCOUT("Entering SchurSolutionManager::InitializeComm()" << endl);
746 
747  /* il master riceve informazioni sulle dimensioni
748  * dei vettori scambiati */
749  if (MyRank == 0) {
750  SAFENEWARR(pRecvDim, int, SolvCommSize);
751  SAFENEWARR(pDispl, int, SolvCommSize + 1);
752  pDispl[0] = 0;
753  }
754 
755  /* trasmissione dimensioni interfacce locali */
756  SolvComm.Gather(&iIntVecDim, 1, MPI::INT, pRecvDim, 1, MPI::INT, 0);
757  if (MyRank == 0) {
758  for (int i = 0; i < SolvCommSize; i++){
759  pDispl[i + 1] = pDispl[i] + pRecvDim[i];
760  }
761  }
762 
763  if (MyRank == 0 && pDispl[SolvCommSize] == 0) {
764  silent_cerr("SchurSolutionManager::InitializeComm(): "
765  "empty problem" << std::endl);
767  }
768 
769  if (MyRank == 0) {
770  SAFENEWARR(pSchurDofs, integer, pDispl[SolvCommSize]);
771  SAFENEWARR(pDofsRecvdList, integer, pDispl[SolvCommSize]);
772  }
773 
774  SolvComm.Gatherv(pIntDofs, iIntVecDim, MPI::INT,
775  pDofsRecvdList, pRecvDim, pDispl, MPI::INT, 0);
776 
777  if (MyRank == 0) {
778  for (int i = 0; i < pDispl[SolvCommSize]; i++) {
779  pSchurDofs[i] = pDofsRecvdList[i];
780  }
781 
782  /* ordina gli indici dei residui locali */
783  std::sort(pSchurDofs, pSchurDofs + pDispl[SolvCommSize]);
784 
785  /* elimina le ripetizioni */
786  integer* p = std::unique(pSchurDofs,
787  pSchurDofs + pDispl[SolvCommSize]);
788 
789  /* dimensione effettiva del residuo locale */
790  iSchurIntDim = p - pSchurDofs;
791  }
792 
793  /* Vettore di trasformazione locale-globale */
796  for (int i = 0; i < iPrbmSize*iPrbmBlocks + 1; i++) {
797  pGlbToLoc[i] = 0;
798  pSchGlbToLoc[i] = 0;
799  }
800 
801  for (int i = 0; i < iLocVecDim; i++) {
802  for (int j = 0; j < iPrbmBlocks; j++) {
803  pGlbToLoc[pLocDofs[i] + j*iBlkSize] = i + j*iBlkSize + 1;
804  }
805  }
806 
807  /* per distinguerli gli indici dei nodi di interfaccia sono negativi */
808  for (int i = 0; i < iIntVecDim; i++) {
809  for (int j = 0; j < iPrbmBlocks; j++) {
810  pGlbToLoc[pIntDofs[i] + j*iBlkSize] = -(i + j*iBlkSize + 1);
811  }
812  }
813 
814  /* Global to local per la matrice di schur */
815  if (MyRank == 0) {
816  for (int i = 0; i < iSchurIntDim; i++) {
817  for (int j = 0; j < iPrbmBlocks; j++) {
818  pSchGlbToLoc[pSchurDofs[i]+ j*iBlkSize] = i + j*iBlkSize + 1;
819  }
820  }
821  }
822 
823  iLocVecDim *= iPrbmBlocks;
824  iIntVecDim *= iPrbmBlocks;
826  iSchurIntDim *= iPrbmBlocks;
827 
828  /* creo i buffer per la ricezione e trasmissione dei messaggi */
829  if (MyRank == 0) {
830  /* i messaggi + grandi sono legati alla ricezione
831  * delle matrici di schur */
832  integer iTmpTot = 0;
833  for (int i = 0; i < SolvCommSize; i++) {
834  pRecvDim[i] = pRecvDim[i]*iPrbmBlocks;
835  iTmpTot += pRecvDim[i]*pRecvDim[i];
836  }
837 
838  for(int i = 0; i < SolvCommSize; i++){
839  pDispl[i + 1] = pDispl[i] + pRecvDim[i];
840  }
841 
842  /* buffer di ricezione */
843  SAFENEWARR(pBuffer, doublereal, iTmpTot);
844 
845  } else {
846  /* il messaggi + grandi sono le ricezioni
847  * dei valori di interfaccia */
848  SAFENEWARR(pBuffer, doublereal, iIntVecDim);
849  }
850 
851  if (MyRank == 0){
852  SAFENEWARRNOFILL(pGSReq, MPI::Request, SolvCommSize);
853  SAFENEWARRNOFILL(pGRReq, MPI::Request, SolvCommSize);
854 
855  } else {
856  SAFENEWARRNOFILL(pGSReq, MPI::Request, 1);
857  SAFENEWARRNOFILL(pGRReq, MPI::Request, 1);
858  }
859 }
860 
861 /* sposta il puntatore al vettore del residuo */
862 doublereal *
864 {
865  silent_cerr("SchurSolutionManager::pdSetResVec(): "
866  "you should not be here!! "
867  "Aborting..." << std::endl);
869 }
870 
871 /* sposta il puntatore al vettore del residuo */
872 doublereal *
874 {
875  silent_cerr("SchurSolutionManager::pdSetSolVec(): "
876  "you should not be here!! "
877  "Aborting..." << std::endl);
879 }
880 
881 /* Rende disponibile l'handler per la matrice */
884 {
885  ASSERT(pMH != NULL);
886  return pMH;
887 }
888 
889 /* Rende disponibile l'handler per il termine noto */
892 {
893  ASSERT(pRVH != NULL);
894  return pRVH;
895 }
896 
897 /* Rende disponibile l'handler per la soluzione */
900 {
901  ASSERT(pSolVH != NULL);
902  return pSolVH;
903 }
904 
905 void
907 {
908  DEBUGCOUT("Entering SchurSolutionManager::StartExchIntRes()" << endl);
909 
910  /* Inizializza Trasmissione di g */
911 
912 #ifdef MPI_PROFILING
913  MPE_Log_event(13, 0, "start");
914  MPE_Log_send(0, G_TAG, iIntVecDim);
915 #endif /* MPI_PROFILING */
916 
917  if (MyRank == 0) {
918 #ifdef MPI_PROFILING
919  MPE_Log_event(19, 0, "start");
920 #endif /* MPI_PROFILING */
921 
922  /* collect remote contributions... */
923  for (int i = 1; i < SolvCommSize; i++) {
924  pGRReq[i] = SolvComm.Irecv(&pBuffer[pDispl[i]],
925  pRecvDim[i], MPI::DOUBLE, i, G_TAG);
926  }
927 
928  /* set up local contribution... */
929  (void)memmove(&pBuffer[pDispl[0]], pgVH->pdGetVec(),
930  sizeof(double)*pRecvDim[0]);
931 
932  } else {
933  pGSReq[0] = SolvComm.Isend(pgVH->pdGetVec(), iIntVecDim,
934  MPI::DOUBLE, 0, G_TAG);
935  }
936 }
937 
938 void
940  const NonlinearSolverTest* t)
941 {
942  DEBUGCOUT("Entering SchurSolutionManager::ComplExchIntRes()" << endl);
943 
944  /* NOTE: right now, all we transmit is the partial result
945  * of the test, as computed by the caller of this function;
946  * the master node is in charge of computing the contribution
947  * of the interface residual */
948  const size_t DBLMSGSIZE = 1;
949  doublereal d[DBLMSGSIZE];
950  d[0] = dRes;
951 
952  if (MyRank == 0) {
953 #ifdef USE_MPI_WAIT
954  MPI::Request::Waitall(SolvCommSize - 1, &pGRReq[1]);
955 #else /* ! USE_MPI_WAIT */
956  while (true) {
957  if (MPI::Request::Testall(SolvCommSize - 1, &pGRReq[1]))
958  {
959 #ifdef MPI_PROFILING
960  MPE_Log_event(20, 0, "end");
961  for (int i = 1; i < SolvCommSize; i++){
962  MPE_Log_receive(i, G_TAG, pRecvDim[i]);
963  }
964 #endif /* MPI_PROFILING */
965  break;
966  }
967  MYSLEEP(1500);
968  }
969 #endif /* ! USE_MPI_WAIT */
970 
971  /* Assembla pSVH */
972  pSchVH->Reset();
973  if (iPrbmBlocks == 1) {
974  for (int i = 0; i < pDispl[SolvCommSize]; i++) {
975  pSchVH->IncCoef(pSchGlbToLoc[pDofsRecvdList[i]], pBuffer[i]);
976  }
977 
978  } else {
979  for (int i = 0; i < pDispl[SolvCommSize]; i++) {
980  int j = i%iBlkSize;
981  int blk = i/iBlkSize;
982 
983  pSchVH->IncCoef(pSchGlbToLoc[pDofsRecvdList[j] + blk*iBlkSize],
984  pBuffer[i]);
985  }
986  }
987 
988  /* interface contribution to error */
989  for (int iCntp1 = 1; iCntp1 <= iSchurIntDim; iCntp1++) {
990  t->TestOne(d[0], *pSchVH, iCntp1);
991  }
992 
993  for (int i = 1; i < SolvCommSize; i++){
994  pGRReq[i] = SolvComm.Irecv(&pBuffer[DBLMSGSIZE*i],
995  DBLMSGSIZE, MPI::DOUBLE,
996  i, G_TAG + 100);
997  }
998 
999 #ifdef USE_MPI_WAIT
1000  MPI::Request::Waitall(SolvCommSize - 1, &pGRReq[1]);
1001 #else /* ! USE_MPI_WAIT */
1002  while (true) {
1003  if (MPI::Request::Testall(SolvCommSize - 1, &pGRReq[1]))
1004  {
1005  break;
1006  }
1007  MYSLEEP(1500);
1008  }
1009 #endif /* ! USE_MPI_WAIT */
1010 
1011  for (int i = 1; i < SolvCommSize; i++) {
1012  t->TestMerge(d[0], pBuffer[DBLMSGSIZE*i]);
1013  }
1014 
1015  for (int i = 1; i < SolvCommSize; i++){
1016  SolvComm.Send(d, DBLMSGSIZE, MPI::DOUBLE, i, G_TAG);
1017  }
1018 
1019  } else {
1020  /* send the residual test value */
1021  SolvComm.Send(d, DBLMSGSIZE, MPI::DOUBLE, 0, G_TAG + 100);
1022 
1023 #ifdef USE_MPI_WAIT
1024  pGSReq[0].Wait();
1025 #else /* ! USE_MPI_WAIT */
1026  /* verifica completamento ricezioni e trasmissione */
1027  while (true) {
1028  if (pGSReq[0].Test()) {
1029 #ifdef MPI_PROFILING
1030  MPE_Log_event(14, 0, "end");
1031 #endif /* MPI_PROFILING */
1032  break;
1033  }
1034  MYSLEEP(150);
1035  }
1036 #endif /* ! USE_MPI_WAIT */
1037 
1038  /* wait for the cumulative residual test */
1039  pGRReq[0] = SolvComm.Irecv(d, DBLMSGSIZE, MPI::DOUBLE,
1040  0, G_TAG);
1041 
1042 #ifdef USE_MPI_WAIT
1043  pGRReq[0].Wait();
1044 #else /* ! USE_MPI_WAIT */
1045  while (true) {
1046  if (pGRReq[0].Test()) {
1047  break;
1048  }
1049  MYSLEEP(150);
1050  }
1051 #endif /* ! USE_MPI_WAIT */
1052  }
1053 
1054  dRes = d[0];
1055 }
1056 
1057 void
1059 {
1060  DEBUGCOUT("Entering SchurSolutionManager::StartExchIntSol()" << endl);
1061 
1062  /* FIXME: the solution interface should have alread been exchanged
1063  * during the solution process, so we don't start anything... */
1064  return;
1065 }
1066 
1067 void
1069  const NonlinearSolverTest* t)
1070 {
1071  DEBUGCOUT("Entering SchurSolutionManager::ComplExchIntSol()" << endl);
1072 
1073  /* right now, all we transmit is the partial result of the test,
1074  * as computed by the caller of this function */
1075  const size_t DBLMSGSIZE = 1;
1076  doublereal d[DBLMSGSIZE];
1077  d[0] = dSol;
1078 
1079  if (MyRank == 0) {
1080  /* Note: the interface contribution should have already
1081  * been transmitted during Solve(), and stored in pgVH. */
1082 
1083  /* interface contribution to error */
1084  for (int iCntp1 = 1; iCntp1 <= iSchurIntDim; iCntp1++) {
1085  t->TestOne(d[0], *pgVH, iCntp1);
1086  }
1087 
1088  for (int i = 1; i < SolvCommSize; i++){
1089  pGRReq[i] = SolvComm.Irecv(&pBuffer[DBLMSGSIZE*i],
1090  DBLMSGSIZE, MPI::DOUBLE,
1091  i, G_TAG + 100);
1092  }
1093 
1094 #ifdef USE_MPI_WAIT
1095  MPI::Request::Waitall(SolvCommSize - 1, &pGRReq[1]);
1096 #else /* ! USE_MPI_WAIT */
1097  while (true) {
1098  if (MPI::Request::Testall(SolvCommSize - 1, &pGRReq[1]))
1099  {
1100  break;
1101  }
1102  MYSLEEP(1500);
1103  }
1104 #endif /* ! USE_MPI_WAIT */
1105 
1106  for (int i = 1; i < SolvCommSize; i++) {
1107  t->TestMerge(d[0], pBuffer[DBLMSGSIZE*i]);
1108  }
1109 
1110  for (int i = 1; i < SolvCommSize; i++){
1111  SolvComm.Send(d, DBLMSGSIZE, MPI::DOUBLE, i, G_TAG);
1112  }
1113 
1114  } else {
1115  /* send the residual test value */
1116  SolvComm.Send(d, DBLMSGSIZE, MPI::DOUBLE, 0, G_TAG + 100);
1117 
1118 #ifdef USE_MPI_WAIT
1119  pGSReq[0].Wait();
1120 #else /* ! USE_MPI_WAIT */
1121  /* verifica completamento ricezioni e trasmissione */
1122  while (true) {
1123  if (pGSReq[0].Test()) {
1124 #ifdef MPI_PROFILING
1125  MPE_Log_event(14, 0, "end");
1126 #endif /* MPI_PROFILING */
1127  break;
1128  }
1129  MYSLEEP(150);
1130  }
1131 #endif /* ! USE_MPI_WAIT */
1132 
1133  /* wait for the cumulative residual test */
1134  pGRReq[0] = SolvComm.Irecv(d, DBLMSGSIZE, MPI::DOUBLE,
1135  0, G_TAG);
1136 
1137 #ifdef USE_MPI_WAIT
1138  pGRReq[0].Wait();
1139 #else /* ! USE_MPI_WAIT */
1140  while (true) {
1141  if (pGRReq[0].Test()) {
1142  break;
1143  }
1144  MYSLEEP(150);
1145  }
1146 #endif /* ! USE_MPI_WAIT */
1147  }
1148 
1149  dSol = d[0];
1150 }
1151 
1152 /* SchurSolutionManager - End */
1153 
1154 #endif /* USE_MPI */
1155 
SchurVectorHandler * pRVH
Definition: schsolman.h:114
virtual void Reset(void)=0
void MatrInitialize(void)
virtual doublereal * GetECol(const integer iCol) const
Definition: schurmh.h:144
void ComplExchIntRes(doublereal &d, const NonlinearSolverTest *t)
virtual void MatEFCReset(void)
Definition: schurmh.h:156
virtual void IncCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: mh.cc:374
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
virtual doublereal * pdGetVec(void) const =0
doublereal * pdSetSolVec(doublereal *pd)
Definition: solman.cc:101
integer * pSchurDofs
Definition: schsolman.h:96
integer * pIntDofs
Definition: schsolman.h:89
doublereal * pdCM
Definition: schsolman.h:123
SolutionManager *const GetSolutionManager(integer iNLD, integer iLWS=0) const
Definition: linsol.cc:455
#define SAFEDELETEARR(pnt)
Definition: mynewmem.h:713
SchurSolutionManager(integer iSize, integer iBlocks, integer *pLocalDofs, int iDim1, integer *pInterfDofs, int iDim2, SolutionManager *pLSM, LinSol &ls)
virtual void IncCoef(integer iRow, const doublereal &dCoef)=0
virtual VectorHandler & CompNewf(VectorHandler &f, const VectorHandler &g) const
Definition: schurmh.h:452
#define MYSLEEP(t)
Definition: mysleep.h:49
virtual void CompLocSchur(void)
Definition: schurmh.h:436
virtual VectorHandler & CompNewg(VectorHandler &g, const VectorHandler &f) const
Definition: schurmh.h:422
#define NO_OP
Definition: myassert.h:74
VectorHandler * pgVH
Definition: schsolman.h:125
void StartExchIntSol(void)
SolutionManager * pLocalSM
Definition: schsolman.h:133
virtual doublereal * GetEColSol(const integer iCol) const
Definition: schurmh.h:150
virtual void TestMerge(doublereal &dResCurr, const doublereal &dResNew) const =0
Definition: linsol.h:39
virtual void Reset(void)=0
integer * pGlbToLoc
Definition: schsolman.h:99
integer * pDofsRecvdList
Definition: schsolman.h:93
const int S_TAG
Definition: schsolman.h:66
VectorHandler * pResHdl(void) const
MatrixHandler * pSchMH
Definition: schsolman.h:118
virtual MatrixHandler * pMatHdl(void) const =0
Definition: mbdyn.h:76
#define DEBUGCOUT(msg)
Definition: myassert.h:232
doublereal * pdSetResVec(doublereal *pd)
Definition: solman.cc:93
virtual void MatrReset(void)=0
VectorHandler * pSolVH
Definition: schsolman.h:116
virtual void Solve(void)=0
SchurMatrixHandler * pMH
Definition: schsolman.h:113
VectorHandler * pSolSchVH
Definition: schsolman.h:120
#define ASSERT(expression)
Definition: colamd.c:977
virtual void MatrInitialize(void)
Definition: solman.cc:72
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
const int G_TAG
Definition: schsolman.h:65
integer * pSchGlbToLoc
Definition: schsolman.h:100
doublereal * pBuffer
Definition: schsolman.h:108
MatrixHandler * pMatHdl(void) const
void StartExchIntRes(void)
VectorHandler * pSolHdl(void) const
#define SAFENEWARRNOFILL(pnt, item, sz)
Definition: mynewmem.h:704
doublereal * pdSetResVec(doublereal *pRes)
integer iGetWorkSpaceSize(void) const
Definition: linsol.cc:339
void ComplExchIntSol(doublereal &d, const NonlinearSolverTest *t)
void SetBMat(MatrixHandler *pBM)
Definition: schurmh.cc:69
#define SAFENEWARR(pnt, item, sz)
Definition: mynewmem.h:701
doublereal * pdSetSolVec(doublereal *pSol)
virtual void TestOne(doublereal &dRes, const VectorHandler &Vec, const integer &iIndex, doublereal dCoef) const =0
void InitializeList(T *list, integer dim, T value)
Definition: schsolman.h:56
virtual ~SchurSolutionManager(void)
VectorHandler * prVH
Definition: schsolman.h:124
double doublereal
Definition: colamd.c:52
void InitializeComm(void)
long int integer
Definition: colamd.c:51
VectorHandler * pSolrVH
Definition: schsolman.h:126
integer * pLocDofs
Definition: schsolman.h:88
SolutionManager * pInterSM
Definition: schsolman.h:134
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710
VectorHandler * pSchVH
Definition: schsolman.h:119