src/mesh/cfMesh/utilities/containers/VRWGraph/VRWGraphSMPModifierTemplates.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     cfMesh is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "VRWGraphSMPModifier.H"
  27 #include "labelPair.H"
  28
  29 # ifdef USE_OMP
  30 #include <omp.h>
  31 # endif
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 template<class ListType>
  41 void VRWGraphSMPModifier::setSizeAndRowSize(const ListType& s)
  42 {
  43     graph_.rows_.setSize(s.size());
  44
  45     # ifdef USE_OMP
  46     label nThreads = 3 * omp_get_num_procs();
  47     if( s.size() < 1000 )
  48         nThreads = 1;
  49     # else
  50     const label nThreads(1);
  51     # endif
  52
  53     label nEntries(0);
  54     DynList<label> procEntries;
  55     procEntries.setSize(nThreads);
  56
  57     # ifdef USE_OMP
  58     # pragma omp parallel num_threads(nThreads)
  59     # endif
  60     {
  61         # ifdef USE_OMP
  62         label& nLocalEntries = procEntries[omp_get_thread_num()];
  63         # else
  64         label& nLocalEntries = procEntries[0];
  65         # endif
  66         nLocalEntries = 0;
  67
  68         # ifdef USE_OMP
  69         # pragma omp for schedule(static)
  70         # endif
  71         forAll(s, i)
  72             nLocalEntries += s[i];
  73
  74         # ifdef USE_OMP
  75         # pragma omp critical
  76         # endif
  77         nEntries += nLocalEntries;
  78
  79         # ifdef USE_OMP
  80         # pragma omp barrier
  81
  82         # pragma omp master
  83         # endif
  84         {
  85             graph_.data_.setSize(nEntries);
  86         }
  87
  88         # ifdef USE_OMP
  89         # pragma omp barrier
  90         # endif
  91
  92         label start(0);
  93         # ifdef USE_OMP
  94         for(label i=0;i<omp_get_thread_num();++i)
  95             start += procEntries[i];
  96         # endif
  97
  98         # ifdef USE_OMP
  99         # pragma omp for schedule(static)
 100         # endif
 101         forAll(s, i)
 102         {
 103             graph_.rows_[i].start() = start;
 104             graph_.rows_[i].size() = s[i];
 105             start += s[i];
 106         }
 107     }
 108 }
 109
 110 template<class GraphType>
 111 void VRWGraphSMPModifier::reverseAddressing(const GraphType& origGraph)
 112 {
 113     graph_.setSize(0);
 114     labelLongList nAppearances;
 115
 116     # ifdef USE_OMP
 117     label nThreads = 3 * omp_get_num_procs();
 118     if( origGraph.size() < 1000 )
 119         nThreads = 1;
 120     # else
 121     const label nThreads(1);
 122     # endif
 123
 124     label minRow(INT_MAX), maxRow(-1);
 125     List<List<LongList<labelPair> > > dataForOtherThreads(nThreads);
 126
 127     # ifdef USE_OMP
 128     # pragma omp parallel num_threads(nThreads)
 129     # endif
 130     {
 131         # ifdef USE_OMP
 132         const label threadI = omp_get_thread_num();
 133         # else
 134         const label threadI(0);
 135         # endif
 136
 137         List<LongList<labelPair> >& dot = dataForOtherThreads[threadI];
 138         dot.setSize(nThreads);
 139
 140         //- find min and max entry in the graph
 141         //- they are used for assigning ranges of values local for each process
 142         label localMinRow(INT_MAX), localMaxRow(-1);
 143         # ifdef USE_OMP
 144         # pragma omp for schedule(static)
 145         # endif
 146         forAll(origGraph, rowI)
 147         {
 148             forAll(origGraph[rowI], i)
 149             {
 150                 const label entryI = origGraph[rowI][i];
 151                 localMaxRow = Foam::max(localMaxRow, entryI);
 152                 localMinRow = Foam::min(localMinRow, entryI);
 153             }
 154         }
 155
 156         ++localMaxRow;
 157
 158         # ifdef USE_OMP
 159         # pragma omp critical
 160         # endif
 161         {
 162             minRow = Foam::min(minRow, localMinRow);
 163             maxRow = Foam::max(maxRow, localMaxRow);
 164
 165             nAppearances.setSize(maxRow);
 166         }
 167
 168         # ifdef USE_OMP
 169         # pragma omp barrier
 170
 171         //- initialise appearances
 172         # pragma omp for schedule(static)
 173         # endif
 174         for(label i=0;i<maxRow;++i)
 175             nAppearances[i] = 0;
 176
 177         # ifdef USE_OMP
 178         # pragma omp barrier
 179         # endif
 180
 181         const label range = (maxRow - minRow) / nThreads + 1;
 182         const label localMin = minRow + threadI * range;
 183         const label localMax = Foam::min(localMin + range, maxRow);
 184
 185         //- find the number of appearances of each element in the original graph
 186         # ifdef USE_OMP
 187         # pragma omp for schedule(static)
 188         # endif
 189         forAll(origGraph, rowI)
 190         {
 191             forAll(origGraph[rowI], j)
 192             {
 193                 const label entryI = origGraph[rowI][j];
 194
 195                 const label threadNo = (entryI - minRow) / range;
 196
 197                 if( threadNo == threadI )
 198                 {
 199                     ++nAppearances[entryI];
 200                 }
 201                 else
 202                 {
 203                     dot[threadNo].append(labelPair(entryI, rowI));
 204                 }
 205             }
 206         }
 207
 208         # ifdef USE_OMP
 209         # pragma omp barrier
 210         # endif
 211
 212         //- count the appearances which are not local to the processor
 213         for(label i=0;i<nThreads;++i)
 214         {
 215             const LongList<labelPair>& data =
 216                 dataForOtherThreads[i][threadI];
 217
 218             forAll(data, j)
 219                 ++nAppearances[data[j].first()];
 220         }
 221
 222         # ifdef USE_OMP
 223         # pragma omp barrier
 224
 225         //- allocate graph
 226         # pragma omp master
 227         # endif
 228         setSizeAndRowSize(nAppearances);
 229
 230         # ifdef USE_OMP
 231         # pragma omp barrier
 232         # endif
 233
 234         for(label i=localMin;i<localMax;++i)
 235         {
 236             nAppearances[i] = 0;
 237         }
 238
 239         //- start filling reverse addressing graph
 240         //- update data from processors with smaller labels
 241         for(label i=0;i<threadI;++i)
 242         {
 243             const LongList<labelPair>& data =
 244                 dataForOtherThreads[i][threadI];
 245
 246             forAll(data, j)
 247             {
 248                 const label entryI = data[j].first();
 249                 graph_(entryI, nAppearances[entryI]++) = data[j].second();
 250             }
 251         }
 252
 253         //- update data local to the processor
 254         # ifdef USE_OMP
 255         # pragma omp for schedule(static)
 256         # endif
 257         forAll(origGraph, rowI)
 258         {
 259             forAll(origGraph[rowI], j)
 260             {
 261                 const label entryI = origGraph[rowI][j];
 262
 263                 if( (entryI >= localMin) && (entryI < localMax) )
 264                     graph_(entryI, nAppearances[entryI]++) = rowI;
 265             }
 266         }
 267
 268         //- update data from the processors with higher labels
 269         for(label i=threadI+1;i<nThreads;++i)
 270         {
 271             const LongList<labelPair>& data =
 272                 dataForOtherThreads[i][threadI];
 273
 274             forAll(data, j)
 275             {
 276                 const label entryI = data[j].first();
 277                 graph_(entryI, nAppearances[entryI]++) = data[j].second();
 278             }
 279         }
 280     }
 281 }
 282
 283 template<class ListType, class GraphType>
 284 void VRWGraphSMPModifier::reverseAddressing
 285 (
 286     const ListType& mapper,
 287     const GraphType& origGraph
 288 )
 289 {
 290     ListType nAppearances;
 291
 292     # ifdef USE_OMP
 293     label nThreads = 3 * omp_get_num_procs();
 294     if( origGraph.size() < 1000 )
 295         nThreads = 1;
 296     # else
 297     const label nThreads(1);
 298     # endif
 299
 300     label minRow(INT_MAX), maxRow(-1);
 301     List<List<LongList<labelPair> > > dataForOtherThreads(nThreads);
 302
 303     # ifdef USE_OMP
 304     # pragma omp parallel num_threads(nThreads)
 305     # endif
 306     {
 307         # ifdef USE_OMP
 308         const label threadI = omp_get_thread_num();
 309         # else
 310         const label threadI(0);
 311         # endif
 312
 313         List<LongList<labelPair> >& dot = dataForOtherThreads[threadI];
 314         dot.setSize(nThreads);
 315
 316         //- find min and max entry in the graph
 317         //- they are used for assigning ranges of values local for each process
 318         label localMinRow(INT_MAX), localMaxRow(-1);
 319         # ifdef USE_OMP
 320         # pragma omp for schedule(static)
 321         # endif
 322         forAll(origGraph, rowI)
 323         {
 324             forAll(origGraph[rowI], i)
 325             {
 326                 const label entryI = mapper[origGraph[rowI][i]];
 327                 localMaxRow = Foam::max(localMaxRow, entryI);
 328                 localMinRow = Foam::min(localMinRow, entryI);
 329             }
 330         }
 331
 332         ++localMaxRow;
 333
 334         # ifdef USE_OMP
 335         # pragma omp critical
 336         # endif
 337         {
 338             minRow = Foam::min(minRow, localMinRow);
 339             maxRow = Foam::max(maxRow, localMaxRow);
 340             nAppearances.setSize(maxRow);
 341         }
 342
 343         # ifdef USE_OMP
 344         # pragma omp barrier
 345
 346         //- initialise appearances
 347         # pragma omp for schedule(static)
 348         # endif
 349         for(label i=0;i<maxRow;++i)
 350             nAppearances[i] = 0;
 351
 352         # ifdef USE_OMP
 353         # pragma omp barrier
 354         # endif
 355
 356         const label range = (maxRow - minRow) / nThreads + 1;
 357         const label localMin = minRow + threadI * range;
 358         const label localMax = Foam::min(localMin + range, maxRow);
 359
 360         //- find the number of appearances of each element in the original graph
 361         # ifdef USE_OMP
 362         # pragma omp for schedule(static)
 363         # endif
 364         forAll(origGraph, rowI)
 365         {
 366             forAll(origGraph[rowI], i)
 367             {
 368                 const label entryI = mapper[origGraph[rowI][i]];
 369
 370                 const label threadNo = (entryI - minRow) / range;
 371
 372                 if( threadNo == threadI )
 373                 {
 374                     ++nAppearances[entryI];
 375                 }
 376                 else
 377                 {
 378                     dot[threadNo].append(labelPair(entryI, rowI));
 379                 }
 380             }
 381         }
 382
 383         # ifdef USE_OMP
 384         # pragma omp barrier
 385         # endif
 386
 387         //- count the appearances which are not local to the processor
 388         for(label i=0;i<nThreads;++i)
 389         {
 390             const LongList<labelPair>& data =
 391                 dataForOtherThreads[i][threadI];
 392
 393             forAll(data, j)
 394                 ++nAppearances[data[j].first()];
 395         }
 396
 397         # ifdef USE_OMP
 398         # pragma omp barrier
 399
 400         //- allocate graph
 401         # pragma omp master
 402         # endif
 403         {
 404             setSizeAndRowSize(nAppearances);
 405         }
 406
 407         # ifdef USE_OMP
 408         # pragma omp barrier
 409         # endif
 410
 411         for(label i=localMin;i<localMax;++i)
 412         {
 413             nAppearances[i] = 0;
 414         }
 415
 416         //- start filling reverse addressing graph
 417         //- update data from processors with smaller labels
 418         for(label i=0;i<threadI;++i)
 419         {
 420             const LongList<labelPair>& data =
 421                 dataForOtherThreads[i][threadI];
 422
 423             forAll(data, j)
 424             {
 425                 const label entryI = data[j].first();
 426                 graph_(entryI, nAppearances[entryI]++) = data[j].second();
 427             }
 428         }
 429
 430         //- update data local to the processor
 431         # ifdef USE_OMP
 432         # pragma omp for schedule(static)
 433         # endif
 434         forAll(origGraph, rowI)
 435         {
 436             forAll(origGraph[rowI], j)
 437             {
 438                 const label entryI = mapper[origGraph[rowI][j]];
 439
 440                 if( (entryI >= localMin) && (entryI < localMax) )
 441                     graph_(entryI, nAppearances[entryI]++) = rowI;
 442             }
 443         }
 444
 445         //- update data from the processors with higher labels
 446         for(label i=threadI+1;i<nThreads;++i)
 447         {
 448             const LongList<labelPair>& data =
 449                 dataForOtherThreads[i][threadI];
 450
 451             forAll(data, j)
 452             {
 453                 const label entryI = data[j].first();
 454                 graph_(entryI, nAppearances[entryI]++) = data[j].second();
 455             }
 456         }
 457     }
 458 }
 459
 460 template<class ListType>
 461 void VRWGraphSMPModifier::reverseAddressing
 462 (
 463     const ListType& mapper,
 464     const VRWGraph& origGraph
 465 )
 466 {
 467     ListType nAppearances;
 468
 469     # ifdef USE_OMP
 470     label nThreads = 3 * omp_get_num_procs();
 471     if( origGraph.size() < 1000 )
 472         nThreads = 1;
 473     # else
 474     const label nThreads(1);
 475     # endif
 476
 477     label minRow(INT_MAX), maxRow(-1);
 478     List<List<LongList<labelPair> > > dataForOtherThreads(nThreads);
 479
 480     # ifdef USE_OMP
 481     # pragma omp parallel num_threads(nThreads)
 482     # endif
 483     {
 484         # ifdef USE_OMP
 485         const label threadI = omp_get_thread_num();
 486         # else
 487         const label threadI(0);
 488         # endif
 489
 490         List<LongList<labelPair> >& dot = dataForOtherThreads[threadI];
 491         dot.setSize(nThreads);
 492
 493         //- find min and max entry in the graph
 494         //- they are used for assigning ranges of values local for each process
 495         label localMinRow(INT_MAX), localMaxRow(-1);
 496         # ifdef USE_OMP
 497         # pragma omp for schedule(static)
 498         # endif
 499         forAll(origGraph, rowI)
 500         {
 501             forAllRow(origGraph, rowI, i)
 502             {
 503                 const label entryI = mapper[origGraph(rowI, i)];
 504                 localMaxRow = Foam::max(localMaxRow, entryI);
 505                 localMinRow = Foam::min(localMinRow, entryI);
 506             }
 507         }
 508
 509         ++localMaxRow;
 510
 511         # ifdef USE_OMP
 512         # pragma omp critical
 513         # endif
 514         {
 515             minRow = Foam::min(minRow, localMinRow);
 516             maxRow = Foam::max(maxRow, localMaxRow);
 517             nAppearances.setSize(maxRow);
 518         }
 519
 520         # ifdef USE_OMP
 521         # pragma omp barrier
 522
 523         //- initialise appearances
 524         # pragma omp for schedule(static)
 525         # endif
 526         for(label i=0;i<maxRow;++i)
 527             nAppearances[i] = 0;
 528
 529         # ifdef USE_OMP
 530         # pragma omp barrier
 531         # endif
 532
 533         const label range = (maxRow - minRow) / nThreads + 1;
 534         const label localMin = minRow + threadI * range;
 535         const label localMax = Foam::min(localMin + range, maxRow);
 536
 537         //- find the number of appearances of each element in the original graph
 538         # ifdef USE_OMP
 539         # pragma omp for schedule(static)
 540         # endif
 541         forAll(origGraph, rowI)
 542         {
 543             forAllRow(origGraph, rowI, i)
 544             {
 545                 const label entryI = mapper[origGraph(rowI, i)];
 546
 547                 const label threadNo = (entryI - minRow) / range;
 548
 549                 if( threadNo == threadI )
 550                 {
 551                     ++nAppearances[entryI];
 552                 }
 553                 else
 554                 {
 555                     dot[threadNo].append(labelPair(entryI, rowI));
 556                 }
 557             }
 558         }
 559
 560         # ifdef USE_OMP
 561         # pragma omp barrier
 562         # endif
 563
 564         //- count the appearances which are not local to the processor
 565         for(label i=0;i<nThreads;++i)
 566         {
 567             const LongList<labelPair>& data =
 568                 dataForOtherThreads[i][threadI];
 569
 570             forAll(data, j)
 571                 ++nAppearances[data[j].first()];
 572         }
 573
 574         # ifdef USE_OMP
 575         # pragma omp barrier
 576
 577         //- allocate graph
 578         # pragma omp master
 579         # endif
 580         {
 581             setSizeAndRowSize(nAppearances);
 582         }
 583
 584         # ifdef USE_OMP
 585         # pragma omp barrier
 586         # endif
 587         for(label i=localMin;i<localMax;++i)
 588             nAppearances[i] = 0;
 589
 590         //- start filling reverse addressing graph
 591         //- update data from processors with smaller labels
 592         for(label i=0;i<threadI;++i)
 593         {
 594             const LongList<labelPair>& data =
 595                 dataForOtherThreads[i][threadI];
 596
 597             forAll(data, j)
 598             {
 599                 const label entryI = data[j].first();
 600                 graph_(entryI, nAppearances[entryI]++) = data[j].second();
 601             }
 602         }
 603
 604         //- update data local to the processor
 605         # ifdef USE_OMP
 606         # pragma omp for schedule(static)
 607         # endif
 608         forAll(origGraph, rowI)
 609         {
 610             forAllRow(origGraph, rowI, j)
 611             {
 612                 const label entryI = mapper[origGraph(rowI, j)];
 613
 614                 if( (entryI >= localMin) && (entryI < localMax) )
 615                     graph_(entryI, nAppearances[entryI]++) = rowI;
 616             }
 617         }
 618
 619         //- update data from the processors with higher labels
 620         for(label i=threadI+1;i<nThreads;++i)
 621         {
 622             const LongList<labelPair>& data =
 623                 dataForOtherThreads[i][threadI];
 624
 625             forAll(data, j)
 626             {
 627                 const label entryI = data[j].first();
 628                 graph_(entryI, nAppearances[entryI]++) = data[j].second();
 629             }
 630         }
 631     }
 632 }
 633
 634 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 635
 636 } // End namespace Foam
 637
 638 // ************************************************************************* //