The p.adjust function in R allows you to put in p values and offers several adjustment methods. For many p values, see https://www.amazon.com/Large-Scale-Inference-Estimation-Prediction-Mathematical/dp/0521192498/187-5605237-0781513?ie=UTF8&qid=1282261412&ref_=sr_1_4&s=books&sr=1-4 among others.

Thanks Daniel,

Probably i should give some more context. In this particular case I have a 35x35 matrix, which implies 595 pairwise correlations. Buy doing 10000 simulations, the lowest possible value for p will be 0.0001. If i apply a bonferroni correction here I will obtain p=0.0595 which will be insufficient to reject the null-hypothesis in all cases.

With >11900 simulations i begin to have the possibility of rejecting the null-hypothesis, but then the problem becomes too computationally intensive (the simulations takes forever if i put more simulations than that).

I'm inclined to just reject the null hypothesis for these multiple comparisons if the observed value falls outside the simulated ones. In other words, for 10000 sims, i would only reject if p=0.0001. This seems awfully stringent criteria, so I was wondering if there is any better way of doing it.

I'll plug a couple of sites on adjusting sets of p-values:

Discussion and SAS code:

http://www.biostathandbook.com/multiplecomparisons.html

R code:

http://rcompanion.org/rcompanion/f_01.html

I don't believe that there is a fully satisfactory answer to your problem. The approach you describe in your last comment seems the best, with the addition of a note describing the challenges the situation presents. A priori comparisons based on theory or prior experiemental work is always preferable if possible, of course. You are to be commended for appreciating the problems raised by large numbers of a posteriori comparisons, too many do not.t.