<html><head><meta http-equiv="Content-Type" content="text/html charset=us-ascii"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; ">Hi Asis,<div><br></div><div>I have cloned the git project and successfully installed it to R. The 'testCptsRealLive.R' does not work though. The 'load' command has problems reading the binary file. </div><div><br></div><div>Nevertheless, I took a look at your OpenMP code: </div><div><br></div><div><div>SEXP conditionalProbabilityTables( SEXP uniqueEdgeLengths, SEXP annos, SEXP</div><div> stringifiedAnnotations, SEXP annotsMutationProbTableList, SEXP</div><div> mutTblLengthColIndx, SEXP nThreads ) {</div><div> BEGIN_RCPP</div><div><br></div><div> //std::cout << "1" << "\n";</div><div> NumericVector numberThreads = NumericVector( nThreads );</div><div> //std::cout << "2" << "\n";</div><div> omp_set_num_threads( numberThreads(0) );</div><div> //std::cout << "3" << "\n";</div><div><br></div><div> NumericVector edgeLengths = NumericVector( uniqueEdgeLengths );</div><div> //std::cout << "4" << "\n";</div><div> CharacterVector edgeLengthsAsStrs = as<CharacterVector>( edgeLengths );</div><div> //std::cout << "5" << "\n";</div><div> List cpts = List( 0 );</div><div> /* std::map< std::string, std::vector<double> > cpts; */</div><div> //std::cout << "6" << "\n";</div><div><br></div><div> #pragma omp parallel for private( uniqueEdgeLengths, annos, stringifiedAnnotations, annotsMutationProbTableList, mutTblLengthColIndx, nThreads, edgeLengths ) shared( cpts )</div><div> for ( int i = 0; i < edgeLengths.size(); i++ ) </div><div> {</div><div> NumericVector currBranchLength( 1, edgeLengths( i ) );</div><div> NumericMatrix cpt = conditionalProbabilityTable( currBranchLength, annos,</div><div> stringifiedAnnotations, annotsMutationProbTableList,</div><div> mutTblLengthColIndx );</div><div> cpts.push_back( cpt, std::string( edgeLengthsAsStrs( i ) ) );</div><div> /* cpts[ std::string( edgeLengthsAsStrs( i ) ) ] = as<std::vector<double> >( cpt ); */</div><div> }</div><div> return( wrap( cpts ) );</div><div><br></div><div> END_RCPP</div><div>}</div><div><br></div><div>My comments: </div><div>1) Set always '++i' not 'i++'. The second directive involves a copy, whereas the first one doesn't. This makes the first one faster. </div><div><br></div><div style="font-size: 13px; ">2) You are using a function 'conditionalProbabilityTable' in your OpenMP-Loop. Have you tested this function?</div><div style="font-size: 13px; "><br></div><div style="font-size: 13px; ">3) The 'cpts' List-Object is shared in the loop, i.e. all threads can access it. Then you write-access the list not by index or name, but by push_back(). If this is thread-safe I cannot tell you, but I wouldn't use it here. </div><div style="font-size: 13px; ">This can cause undefined behavior, i.e. your program runs sometimes without errors and another time you get a segfault. Try to access the list at least with an index or name. Another possibility: Use the push_back</div><div style="font-size: 13px; ">command in a '#pragma omp critical' directive, which allows only one thread at a time access to the List (this can cause a longer runtime). </div><div style="font-size: 13px; "><br></div><div style="font-size: 13px; ">4) You are using in the called function 'conditionalProbabilityTable' another for-loop: </div><div style="font-size: 13px; "><br></div><div><span style="font-size: 13px;">SEXP conditionalProbabilityTable( SEXP branchLength, SEXP annos,</span></div><div><span style="font-size: 13px;"> SEXP stringifiedAnnotations, SEXP annotsMutationProbTables, SEXP</span></div><div><span style="font-size: 13px;"> mutTblLengthColIndx ) {</span></div><div><span style="font-size: 13px;"> BEGIN_RCPP</span></div><div><span style="font-size: 13px;"><br></span></div><div><span style="font-size: 13px;"> List annosLst( annos );</span></div><div><span style="font-size: 13px;"> CharacterVector annosAsStrs( stringifiedAnnotations );</span></div><div><span style="font-size: 13px;"> NumericMatrix cpt = NumericMatrix( annosLst.size(), annosLst.size() );</span></div><div><span style="font-size: 13px;"> cpt.attr( "dimnames" ) = List::create( annosAsStrs, annosAsStrs );</span></div><div><span style="font-size: 13px;"><br></span></div><div><span style="font-size: 13px;"> for ( int i = 0; i < annosLst.size(); i++ ) {</span></div><div><span style="font-size: 13px;"> CharacterVector compositeAnnotation = annosLst( i );</span></div><div><span style="font-size: 13px;"> double compAnnoMutProb = 1.0;</span></div><div><span style="font-size: 13px;"> std::string ua = "unknown";</span></div><div><span style="font-size: 13px;"> std::string caFirst = as<std::string>( compositeAnnotation( 0 ) );</span></div><div><span style="font-size: 13px;"> if ( ua != caFirst ) {</span></div><div><span style="font-size: 13px;"> compAnnoMutProb = mutationProbability( compositeAnnotation,</span></div><div><span style="font-size: 13px;"> branchLength, annotsMutationProbTables, mutTblLengthColIndx )( 0 );</span></div><div><span style="font-size: 13px;"> }</span></div><div><span style="font-size: 13px;"> double mutToOtherAnnosProb = compAnnoMutProb / ( annosLst.size() - 1 );</span></div><div><span style="font-size: 13px;"> NumericVector colmn( annosLst.size(), mutToOtherAnnosProb );</span></div><div><span style="font-size: 13px;"> colmn( i ) = 1.0 - compAnnoMutProb;</span></div><div><span style="font-size: 13px;"> cpt( _, i ) = colmn;</span></div><div><span style="font-size: 13px;"> }</span></div><div><span style="font-size: 13px;"><br></span></div><div><span style="font-size: 13px;"> return( wrap( cpt ) );</span></div><div><span style="font-size: 13px;"><br></span></div><div><span style="font-size: 13px;"> END_RCPP</span></div><div><span style="font-size: 13px;">}</span></div><div><span style="font-size: 13px;">In such cases, you should always think about which loop to parallelize. Usually you parallelize the loop with the bigger workload. Only in case this loop has only</span></div><div style="font-size: 13px; ">a few iterations it would make sense to parallelize the other loop. Another possibility: Create nested parallelized loops (advanced-level). For this you just use another '#pragma omp for' directive before the second for-loop. Set the variables OMP_NESTED to true. </div><div style="font-size: 13px; ">Set OMP_MAX_ACTIVE_LEVELS=2. Further make sure, that you have enough stack size for the threads: OMP_STACKSIZE=1000000. </div><div style="font-size: 13px; ">Do not even dream about an improvement, when working with two cores! You need at least 8 to see the work power of nesting!</div><div style="font-size: 13px; "><br></div><div style="font-size: 13px; ">5) Set the environment variable OMP_DYNAMIC=false. The dynamic chunking does in most cases not work and static chunking has a better performance. </div><div style="font-size: 13px; "><br></div><div style="font-size: 13px; ">6) The OpenMP library was intended to be used on very simple objects, i.e. arrays and vectors. One reason for this is also the socalled page placing, where you place the memory chunks, which are used by a certain core</div><div style="font-size: 13px; ">on the local memory of this core (NUMA-technology). If we have very complicated objects, it is not always clear, how the memory can be placed for an optimal access. </div><div style="font-size: 13px; ">So, to cut a long story short, I would always use very simple objects in an OpenMP directive. As NumericVectors and Armadillo-Vectors rely on an inner array for memory allocations </div><div style="font-size: 13px; ">(e.g. you can initialize both objects with an already existing stl::vector for memory reutilization), this objects can still be considered as very simple -> fast. Second, page placement</div><div style="font-size: 13px; ">is more easily achieved, if wanted/needed. Third, write-access can be done easily by index and is thread-safe if indices do not intersect or we have read-access only. </div><div style="font-size: 13px; ">Run over vectors and matrices and stack them later on into a list. </div><div style="font-size: 13px; "><br></div><div style="font-size: 13px; ">7) It seems sometimes to me, that you create a lot of objects. Check if these objects are all needed. For example, consider in ''conditionalProbabilityTable':</div><div><div style="font-size: medium; "><span style="font-size: 13px; ">colmn( i ) = 1.0 - compAnnoMutProb;</span></div><div style="font-size: medium; "><span style="font-size: 13px; ">cpt( _, i ) = colmn;</span></div><div style="font-size: medium; "><span style="font-size: 13px; ">I think, you could write this immediately into 'cpt'?</span></div><div style="font-size: medium; "><span style="font-size: 13px; "><br></span></div><div style="font-size: medium; "><span style="font-size: 13px; ">8) Use instruments for testing: A very good introduction and an overview can be found on the website of the HPC-Cluster in Aachen (to which by the way you could get access from the University Bonn), </span><span style="font-size: 13px; "><a href="http://www.rz.rwth-aachen.de/global/show_document.asp?id=aaaaaaaaaaclexs">http://www.rz.rwth-aachen.de/global/show_document.asp?id=aaaaaaaaaaclexs</a>. Identify with an instrument the </span></div><div style="font-size: medium; "><span style="font-size: 13px; ">big workload in your sequential program and try to parallelize this section. Then check another time where another workload is, and so on. I would further suggest to simplify your program, when you </span></div><div><span style="font-size: 13px;">encounter undefined behavior. This helps in most cases a lot. </span></div></div><div><span style="font-size: 13px;"><br></span></div><div><span style="font-size: 13px;">I hope I could help here. I already saw your new thread on the list and it seems, everything is running. So congrats and keep on!</span></div><div><span style="font-size: 13px;"><br></span></div><div><span style="font-size: 13px;"><br></span></div><div><span style="font-size: 13px;">Best</span></div><div><span style="font-size: 13px;"><br></span></div><div><span style="font-size: 13px;">Simon</span></div><div><span style="font-size: 13px;"><br></span></div><div><span style="font-size: 13px;">P.S. Of course we could meet in Bonn or Cologne (where I am living and working now). I am also very interested, what is needed in Bioinformatics, maybe we could interchange ideas. Let us keep in contact!</span></div><div style="font-size: 13px; "><br></div><div style="font-size: 13px; "><br></div><div><div>On Jun 1, 2013, at 8:53 PM, Asis Hallab <<a href="mailto:asis.hallab@gmail.com">asis.hallab@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><blockquote type="cite">Dear Simon,<br><br>thank you very much for your insights and offer to help.<br><br>I prepared a very small project on github:<br><a href="https://github.com/asishallab/PhyloFun_Rccp">https://github.com/asishallab/PhyloFun_Rccp</a><br><br>You'll find installation and test instructions as well as references<br>to the file containing the OpenMP using C++ code.<br><br>Your help is much appreciated.<br><br>Kind regards!<br></blockquote></div><br></div></body></html>