[Debian-med-packaging] MP patch for PhyML

Diego Darriba López ddarriba at udc.es
Tue Feb 15 22:38:16 UTC 2011


Dear Andreas, 


My OpenMP implementation perform the MP just using pragma directives, so if the -fopenmp flag is not set, the compiler will just ignore these lines and it will perform a single processor build. I checked both builds with gcc. We would be in a different situation if I had used OpenMP API calls. In fact, the #include<omp.h> I included in my patch it is not necessary at all (sorry). This way, the MP or sequential build just depends on the use of the -fopenmp flag. 


Keep in mind that this flag differs between compilares (e.g, -openmp for icc, -mp for openUH, etc.) if you want to bring support for different compilers in the autoconf file. 


Kind regards, 


Diego. 


De: "Andreas Tille" <andreas at an3as.eu> 
Para: "Diego Darriba" <ddarriba at udc.es> 
CC: "s guindon" <s.guindon at auckland.ac.nz>, "Stephane Guindon" <guindon at stat.auckland.ac.nz>, "Olivier Gascuel" <gascuel at lirmm.fr>, "Jean-Francois Dufayard" <jeanfrancois.dufayard at gmail.com>, "David Posada" <dposada at uvigo.es>, "Debian Med Packaging Team" <debian-med-packaging at lists.alioth.debian.org> 
Enviados: Martes, 15 de Febrero 2011 21:25:51 
Asunto: Re: MP patch for PhyML 

Dear Diego, 

I admit I have no experience with MP builds but I wonder whether you 
should not include the changes into #ifdef USEMP / #else / #endif 
statements to enable both - single processor and MP builds. Or does the 
patch work with both and it is just the build option which controls 
everything? 

Kind regards 

Andreas. 

On Tue, Feb 15, 2011 at 06:58:19PM +0100, Diego Darriba wrote: 
> Dear Stephan et al, 
> 
> I send you the patch for the current phyml trunk. I checked the results 
> with several input data, and I think I solved correctly some 
> dependencies in the parallelized loops from my computer science point of 
> view. In the configure.ac it is just necessary to include the *-fopenmp* 
> flag for compiling (if using gcc compiler). I would modify it, but I 
> prefer you to do so. That way you can create a conditional Makefile, or 
> include the flag as a commentary, or whatever. 
> 
> If the compilation flag is not set, the application will compile with no 
> error, but it will run as before (i.e, sequentially). 
> 
> Please, let me know if you have any doubt or any problem. 
> 
> Kind Regards, 
> 
> Diego. 
> 
> El 15/02/11 02:51, Stephane Guindon escribió: 
>> Dear Andreas, 
>> 
>> Thank you very much for contributing to PhyML. 
>> Yes, it makes perfect sense to try and include this functionality 
>> into the latest release of the program. Please feel free to take 
>> a look at the most recent version of the sources on our svn server 
>> (http://code.google.com/p/phyml/source/browse/#svn%2Ftrunk%2Fsrc) 
>> and let me know how to amend the configure.ac, *.c and *.h file 
>> in an appropriate manner. 
>> 
>> Regards, 
>> 
>> -Stephane- 
>> 
>> On 15/02/11 03:26, Andreas Tille wrote: 
>>> Hi PhyML authors, 
>>> 
>>> when discssing with Diego about a program which uses PhyML he told me 
>>> about his MP patch which would enable making use of hybrid memory when 
>>> running PhyML on clusters. I was considering to port this patch to the 
>>> Debian packaged version where I intended to support two versions of 
>>> the package: one conventional phylip and one phylip-mp. 
>>> 
>>> When inspecting the patch I noticed that it was done against an older 
>>> version of PhyML (20100123). I would consider it a reasonable idea if 
>>> this patch would be included into the upstream version by just wrapping 
>>> it into "#ifdef"s properly to enable switching between both versions by 
>>> simply providing different options to make. If the patch would be taken 
>>> over into the official PhyML code this would simplify our work as 
>>> distributor (as well as Diego's work) and it might be of extra value for 
>>> other users. 
>>> 
>>> What do you think about this idea? 
>>> 
>>> Kind regards 
>>> 
>>> Andreas. 
>>> 
> 

> Index: src/lk.c 
> =================================================================== 
> --- src/lk.c (revisi??n: 582) 
> +++ src/lk.c (copia de trabajo) 
> @@ -11,8 +11,8 @@ 
> */ 
> 
> #include "lk.h" 
> +#include <omp.h> 
> 
> - 
> /* int LIM_SCALE; */ 
> /* phydbl LIM_SCALE_VAL; */ 
> /* phydbl BIG; */ 
> @@ -320,8 +320,9 @@ 
> 
> phydbl Lk(t_tree *tree) 
> { 
> - int br; 
> + int br, site; 
> int n_patterns; 
> + phydbl c_lnL; 
> 
> tree->old_lnL = tree->c_lnL; 
> 
> @@ -349,10 +350,17 @@ 
> 
> tree->c_lnL = .0; 
> tree->sum_min_sum_scale = .0; 
> - For(tree->curr_site,n_patterns) 
> + 
> + c_lnL = .0; 
> + #pragma omp parallel for reduction(+:c_lnL) 
> + For(site,n_patterns) 
> { 
> - if(tree->data->wght[tree->curr_site] > SMALL) Lk_Core(tree->noeud[0]->b[0],tree); 
> + if(tree->data->wght[site] > SMALL) { 
> + Lk_Core(tree->noeud[0]->b[0],tree,site); 
> + c_lnL += tree->c_lnL_sorted[site]; 
> + } 
> } 
> + tree->c_lnL = c_lnL; 
> 
> /* Qksort(tree->c_lnL_sorted,NULL,0,n_patterns-1); */ 
> 
> @@ -372,7 +380,8 @@ 
> 
> phydbl Lk_At_Given_Edge(t_edge *b_fcus, t_tree *tree) 
> { 
> - int n_patterns; 
> + int n_patterns,site; 
> + phydbl c_lnL; 
> 
> n_patterns = tree->n_pattern; 
> 
> @@ -396,10 +405,17 @@ 
> 
> tree->c_lnL = .0; 
> tree->sum_min_sum_scale = .0; 
> - For(tree->curr_site,n_patterns) 
> + 
> + c_lnL = .0; 
> + #pragma omp parallel for reduction(+:c_lnL) 
> + For(site,n_patterns) 
> { 
> - if(tree->data->wght[tree->curr_site] > SMALL) Lk_Core(b_fcus,tree); 
> + if(tree->data->wght[site] > SMALL) { 
> + Lk_Core(b_fcus,tree,site); 
> + c_lnL += tree->c_lnL_sorted[site]; 
> + } 
> } 
> + tree->c_lnL = c_lnL; 
> 
> /* Qksort(tree->c_lnL_sorted,NULL,0,n_patterns-1); */ 
> 
> @@ -422,7 +438,7 @@ 
> 
> #ifndef USE_OLD_LK 
> 
> -phydbl Lk_Core(t_edge *b, t_tree *tree) 
> +phydbl Lk_Core(t_edge *b, t_tree *tree, int site) 
> { 
> phydbl log_site_lk; 
> phydbl site_lk_cat, site_lk; 
> @@ -431,7 +447,7 @@ 
> phydbl max_sum_scale; 
> phydbl sum; 
> int ambiguity_check,state; 
> - int catg,ns,k,l,site; 
> + int catg,ns,k,l; 
> int dim1,dim2,dim3; 
> int *sum_scale_left_cat,*sum_scale_rght_cat; 
> phydbl multiplier; 
> @@ -439,6 +455,7 @@ 
> phydbl tmp; 
> phydbl logbig; 
> phydbl inv_site_lk; 
> + phydbl tree_site_lk_cat[tree->mod->n_catg]; 
> 
> logbig = LOG((phydbl)BIG); 
> 
> @@ -451,7 +468,6 @@ 
> site_lk_cat = .0; 
> ambiguity_check = -1; 
> state = -1; 
> - site = tree->curr_site; 
> ns = tree->mod->ns; 
> 
> if((b->rght->tax) && (!tree->mod->s_opt->greedy)) 
> @@ -537,7 +553,7 @@ 
> } 
> } 
> } 
> - tree->site_lk_cat[catg] = site_lk_cat; 
> + tree_site_lk_cat[catg] = site_lk_cat; 
> } 
> 
> max_sum_scale = (phydbl)BIG; 
> @@ -562,7 +578,7 @@ 
> Warn_And_Exit("\n"); 
> } 
> 
> - tmp = sum + (logbig - LOG(tree->site_lk_cat[catg]))/(phydbl)LOG2; 
> + tmp = sum + (logbig - LOG(tree_site_lk_cat[catg]))/(phydbl)LOG2; 
> if(tmp < max_sum_scale) max_sum_scale = tmp; /* min of the maxs */ 
> } 
> 
> @@ -575,7 +591,7 @@ 
> { 
> exponent = -(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale; 
> 
> - site_lk_cat = tree->site_lk_cat[catg]; 
> + site_lk_cat = tree_site_lk_cat[catg]; 
> 
> if(exponent > 0) 
> { 
> @@ -605,12 +621,12 @@ 
> PhyML_Printf("\n+ site=%4d cat=%4d site_lk_cat=%G sum_scale=%d max=%G fact=%d expo=%d dbl=%G", 
> tree->curr_site, 
> catg, 
> - tree->site_lk_cat[catg], 
> + tree_site_lk_cat[catg], 
> sum_scale_left_cat[catg]+sum_scale_rght_cat[catg], 
> max_sum_scale, 
> fact_sum_scale, 
> -(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale, 
> - (double)tree->site_lk_cat[catg] * pow(2.,-(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale)); 
> + (double)tree_site_lk_cat[catg] * pow(2.,-(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale)); 
> 
> site_lk_cat = BIG / 10; 
> } 
> @@ -623,12 +639,12 @@ 
> PhyML_Printf("\n- site=%4d cat=%4d site_lk_cat=%G sum_scale=%d max=%G fact=%d expo=%d dbl=%G", 
> tree->curr_site, 
> catg, 
> - tree->site_lk_cat[catg], 
> + tree_site_lk_cat[catg], 
> sum_scale_left_cat[catg]+sum_scale_rght_cat[catg], 
> max_sum_scale, 
> fact_sum_scale, 
> -(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale, 
> - (double)tree->site_lk_cat[catg] * pow(2.,-(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale)); 
> + (double)tree_site_lk_cat[catg] * pow(2.,-(sum_scale_left_cat[catg]+sum_scale_rght_cat[catg])+fact_sum_scale)); 
> 
> Exit("\n"); 
> } 
> @@ -636,13 +652,13 @@ 
> site_lk_cat = .0; 
> } 
> 
> - tree->site_lk_cat[catg] = site_lk_cat; 
> + tree_site_lk_cat[catg] = site_lk_cat; 
> } 
> 
> site_lk = .0; 
> For(catg,tree->mod->n_catg) 
> { 
> - site_lk += tree->site_lk_cat[catg] * tree->mod->gamma_r_proba[catg]; 
> + site_lk += tree_site_lk_cat[catg] * tree->mod->gamma_r_proba[catg]; 
> } 
> 
> /* The substitution model does include invariable sites */ 
> @@ -677,7 +693,7 @@ 
> 
> log_site_lk = LOG(site_lk) - (phydbl)LOG2 * fact_sum_scale; 
> 
> - For(catg,tree->mod->n_catg) tree->log_site_lk_cat[catg][site] = LOG(tree->site_lk_cat[catg]) - (phydbl)LOG2 * fact_sum_scale; 
> + For(catg,tree->mod->n_catg) tree->log_site_lk_cat[catg][site] = LOG(tree_site_lk_cat[catg]) - (phydbl)LOG2 * fact_sum_scale; 
> 
> if(isinf(log_site_lk) || isnan(log_site_lk)) 
> { 
> @@ -694,7 +710,6 @@ 
> /* Multiply log likelihood by the number of times this site pattern is found in the data */ 
> tree->c_lnL_sorted[site] = tree->data->wght[site]*log_site_lk; 
> 
> - tree->c_lnL += tree->data->wght[site]*log_site_lk; 
> /* tree->sum_min_sum_scale += (int)tree->data->wght[site]*min_sum_scale; */ 
> 
> return log_site_lk; 
> @@ -1118,6 +1133,7 @@ 
> Pij2 = d->b[dir2]->Pij_rr; 
> 
> /* For every site in the alignment */ 
> + #pragma omp parallel for private(state_v1,state_v2,ambiguity_check_v1,ambiguity_check_v2,catg,i,p1_lk1,p2_lk2) 
> For(site,n_patterns) 
> { 
> state_v1 = state_v2 = -1; 
> @@ -1399,6 +1415,7 @@ 
> Pij2 = d->b[dir2]->Pij_rr; 
> 
> /* For every site in the alignment */ 
> + #pragma omp parallel for private(state_v1,state_v2,ambiguity_check_v1,ambiguity_check_v2,catg,i,p1_lk1,p2_lk2) 
> For(site,n_patterns) 
> { 
> state_v1 = state_v2 = -1; 
> @@ -2234,13 +2251,13 @@ 
> sides of t_edge *b are up-to-date. Calculate the log-likelihood at one site. 
> */ 
> #ifdef USE_OLD_LK 
> -phydbl Lk_Core(t_edge *b, t_tree *tree) 
> +phydbl Lk_Core(t_edge *b, t_tree *tree, int site) 
> { 
> phydbl log_site_lk, site_lk, site_lk_cat; 
> phydbl scale_left, scale_rght; 
> phydbl sum; 
> int ambiguity_check,state; 
> - int catg,ns,k,l,site; 
> + int catg,ns,k,l; 
> int dim1,dim2,dim3; 
> 
> 
> @@ -2253,7 +2270,6 @@ 
> site_lk_cat = .0; 
> ambiguity_check = -1; 
> state = -1; 
> - site = tree->curr_site; 
> ns = tree->mod->ns; 
> 
> tree->old_site_lk[site] = tree->cur_site_lk[site]; 
> @@ -2383,7 +2399,7 @@ 
> tree->cur_site_lk[site] = log_site_lk; 
> /* Multiply log likelihood by the number of times this site pattern is found oin the data */ 
> tree->c_lnL_sorted[site] = tree->data->wght[site]*log_site_lk; 
> - tree->c_lnL += tree->c_lnL_sorted[site]; 
> + 
> return log_site_lk; 
> } 
> 
> @@ -2492,6 +2508,7 @@ 
> Pij2 = d->b[dir2]->Pij_rr; 
> 
> /* For every site in the alignment */ 
> + #pragma omp parallel for private(scale_v1,scale_v2,state_v1,state_v2,ambiguity_check_v1,ambiguity_check_v2,catg,i,j,p1_lk1,p2_lk2,max_p_lk) 
> For(site,n_patterns) 
> { 
> /* Current scaling values at that site */ 
> Index: src/lk.h 
> =================================================================== 
> --- src/lk.h (revisi??n: 582) 
> +++ src/lk.h (copia de trabajo) 
> @@ -51,7 +51,7 @@ 
> phydbl Update_Lk_At_Given_Edge(t_edge *b_fcus, t_tree *tree); 
> void Update_P_Lk_Greedy(t_tree *tree, t_edge *b_fcus, t_node *n); 
> void Get_All_Partial_Lk_Scale_Greedy(t_tree *tree, t_edge *b_fcus, t_node *a, t_node *d); 
> -phydbl Lk_Core(t_edge *b, t_tree *tree); 
> +phydbl Lk_Core(t_edge *b, t_tree *tree, int site); 
> phydbl Lk_Triplet(t_node *a, t_node *d, t_tree *tree); 
> void Print_Lk_Given_Edge_Recurr(t_node *a, t_node *d, t_edge *b, t_tree *tree); 
> phydbl *Post_Prob_Rates_At_Given_Edge(t_edge *b, phydbl *post_prob, t_tree *tree); 


-- 
http://fam-tille.de 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.alioth.debian.org/pipermail/debian-med-packaging/attachments/20110215/681a99e2/attachment-0001.htm>


More information about the Debian-med-packaging mailing list