Merge lp:~maddevelopers/mg5amcnlo/COLLIER into lp:~maddevelopers/mg5amcnlo/2.4.3

[Valentin Hirschi]

This branch introduces the longed-for interface to COLLIER!

All tests pass (incl. the comparison with the multitude of benchmark processes) and MG5aMC user can now effortlessly use COLLIER for all the applications that can involve MadLoop.

So please give it a try simply by running a couple of loop processes with different setup, and after having installed COLLIER with

MG5_aMC> install collier

Regarding performances, COLLIER holds up to its promises (when properly using its global caching system). I compiled a tex document that summarizes the profiling I did so far:

https://drive.google.com/open?id=0By4qXx23HZGESkU2c0RMQ0J1RGs

In summary, we observe that:

a) COLLIER is faster than CutTools but slower than Ninja, although one must keep in mind that Ninja's timings must be doubled in practice for the stability test.

b) COLLIER scales worse than integrand reduction with the rank of the loop *and its multiplicity* (and significantly so, something I had not appreciated before). This is striking when looking at the timings of the 2 > 6 process. But it is nice that the COLI branch of COLLIER has no hard limitation in rank or multiplicity.

c) COLLIER's stability is indeed astonishing! The numerical stability plots speak for themselves... And this was not even for the most difficult processes we can generate with MadLoop (see http://arxiv.org/pdf/1604.01363.pdf). However, COLLIER's internal stability assessment is quite a bit over-conservative.

Finally, there are a number of new COLLIER-related parameters that you can find at the end of the MadLoopParams.dat card (please have a look). Speaking of which, notice that I decided to change the default value of the following three existing MadLoop parameters starting from this version:

    a) NRotationsDP 1-> 0
       We have been using that for quite some time and it doubles
       MadLoop runtime. I think we now have enough confidence in the
       loop direction stability test to drop it.

    b) COLLIERUseInternalStabilityTest False -> True
       Even though for complicated processes this test seems to
       overshoot the actual stability by up to 5 digits, it seems to never
       undershoots it, so it is safe to rely on it and it comes in at no
       additional computational time.
       So we can use it and if at some point too many unstable PS point show
       (but that should almost never occur) with COLLIER, then we can turn this off
       and fall back on the loop-direction stability test instead.

    c) UseQPIntegrandForNinja False -> True
       For complicated processes where a lot of points need stability rescuing
       with Ninja, it turns out that quad. prec. in the integrand can
       significantly increases the accuracy (up to 4-5 digits), so it
       is actually better to leave this on, especially since we
       now have COLLIER first.

Let me know your comments/questions on the above, the interface and the code if you want to look at it.

Cheers

[Olivier Mattelaer (olivier-mattelaer)]

Hi Valentin,

This is so cool.
I did not have the time to check yet but this is some first comment:

1) I would refrain to include this in 2.4.1, mainly because I want to have that version as bug fix only and be released very very soon (in max two weeks).

2) I do not understand the reason "UseQPIntegrandForNinja False -> True". You said "especially since we now have COLLIER first." This is actually not really True 99% of the user will not install collier and your option should reflect that scenario (or you update your option according to the tool available which can be really interesting.

Cheers,

Olivier

[Valentin Hirschi (valentin-hirschi)]

Hi Olivier,

On Thu, May 26, 2016 at 2:03 AM, Olivier Mattelaer <
<email address hidden>> wrote:

> Hi Valentin,
>
> This is so cool.
> I did not have the time to check yet but this is some first comment:
>
> 1) I would refrain to include this in 2.4.1, mainly because I want to have
> that version as bug fix only and be released very very soon (in max two
> weeks).
>

I understand. But I think that this implementation is fairly robust as it
doesn't modify anything else except MadLoop's output (to the exception of
the 'unique_id' additional option to generate_subrpocess_v4).
But still, I understand if you would rather delay it, just for safety; I'm
ok with this.

> 2) I do not understand the reason "UseQPIntegrandForNinja False -> True".
> You said "especially since we now have COLLIER first." This is actually not
> really True 99% of the user will not install collier and your option should
> reflect that scenario (or you update your option according to the tool
> available which can be really interesting.
>

I will have COLLIER installed by default from a local copy that can then be
updated at any point by the user from the online tarball using the command
'install collier'. (I'll ask D&D for their consent on this first, event
though the tool is public, just to be politically correct). Basically, the
same strategy as I used for Ninja.
And COLLIER has *no dependence*, is in plain fortran90 and only needs a
minute or so to install. So yes, everybody will have it :).

Regarding the order, I think I want to still keep Ninja first, then
followed by COLLIER and finally CutTools (though it should almost never get
there anymore).

Cheers,

>
> Cheers,
>
> Olivier
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/COLLIER/+merge/295787
> Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/2.4.1.
>

--
Valentin

[Olivier Mattelaer (olivier-mattelaer)]

Hi Valentin,

I wanted to start this review today. But I was stop quickly:

Enter COLLIER directory and installing COLLIER
/Users/omatt/Documents/eclipse/COLLIER/HEPTools/HEPToolsInstallers/installCOLLIER.sh: line 14: cmake: command not found
make: *** No targets specified and no makefile found. Stop.
Copying the static library and modules generated to target directory
cp: libcollier.a: No such file or directory
cp: modules/*.mod: No such file or directory
Finished installing COLLIER

Cheers,

Olivier

review: Needs Fixing

[Valentin Hirschi (valentin-hirschi)]

Dear reviewers,

The points raised by Olivier and Rikkert (offline) have been resolved.
This branch is mature and I plan on merging it early next week, so please have a quick look/trial by then if you can.

Cheers,

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

No the tests are not working...
So you can not merge.

Cheers,

Olivier

review: Disapprove

[Valentin Hirschi (valentin-hirschi)]

Hi Olivier,

On Mon, Jun 27, 2016 at 7:36 AM, Olivier Mattelaer <
<email address hidden>> wrote:

> Review: Disapprove
>
> Hi,
>
> No the tests are not working...
> So you can not merge.
>

It goes without saying that I'll run and fix the tests before actually
performing the merge...
I like to do it at the very end when I'm really about to merge, and not
preemptively before each merging proposal (at least for these "simple"
branches; it avoids having to run the tests each time when you implement
something the reviewer asked for).

So can I take this as your green light from a structural/functional
standpoint?
I can also send you the test-fixed version for you to try just before I
merge if you want.

Cheers,

>
> Cheers,
>
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/COLLIER/+merge/298349
> Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/2.5.0.
>

--
Valentin

[Olivier Mattelaer (olivier-mattelaer)]

Hi Valentin,

It is always important that at least two people runs the test to be sure that all file are include correctly
and that you do not have test which are machine dependent (which means useless test).
So clearly you do not have a green light.

Cheers,

Olivier

> On Jun 27, 2016, at 17:30, Valentin Hirschi <email address hidden> wrote:
>
> Hi Olivier,
>
> On Mon, Jun 27, 2016 at 7:36 AM, Olivier Mattelaer <
> <email address hidden>> wrote:
>
>> Review: Disapprove
>>
>> Hi,
>>
>> No the tests are not working...
>> So you can not merge.
>>
>
> It goes without saying that I'll run and fix the tests before actually
> performing the merge...
> I like to do it at the very end when I'm really about to merge, and not
> preemptively before each merging proposal (at least for these "simple"
> branches; it avoids having to run the tests each time when you implement
> something the reviewer asked for).
>
> So can I take this as your green light from a structural/functional
> standpoint?
> I can also send you the test-fixed version for you to try just before I
> merge if you want.
>
> Cheers,
>
>
>>
>> Cheers,
>>
>>
>> Olivier
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/COLLIER/+merge/298349
>> Your team MadDevelopers is subscribed to branch
>> lp:~maddevelopers/mg5amcnlo/2.5.0.
>>
>
>
>
> --
> Valentin
>
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/COLLIER/+merge/298349
> You are reviewing the proposed merge of lp:~maddevelopers/mg5amcnlo/COLLIER into lp:~maddevelopers/mg5amcnlo/2.5.0.

[Valentin Hirschi (valentin-hirschi)]

I stressed "from a structural/functional standpoint" and also that I would
tell you to run the tests once I will have run/fix them myself.
So I'll let you know.

On Mon, Jun 27, 2016 at 1:00 PM, Olivier Mattelaer <
<email address hidden>> wrote:

> Hi Valentin,
>
> It is always important that at least two people runs the test to be sure
> that all file are include correctly
> and that you do not have test which are machine dependent (which means
> useless test).
> So clearly you do not have a green light.
>
> Cheers,
>
> Olivier
>
> > On Jun 27, 2016, at 17:30, Valentin Hirschi <email address hidden>
> wrote:
> >
> > Hi Olivier,
> >
> > On Mon, Jun 27, 2016 at 7:36 AM, Olivier Mattelaer <
> > <email address hidden>> wrote:
> >
> >> Review: Disapprove
> >>
> >> Hi,
> >>
> >> No the tests are not working...
> >> So you can not merge.
> >>
> >
> > It goes without saying that I'll run and fix the tests before actually
> > performing the merge...
> > I like to do it at the very end when I'm really about to merge, and not
> > preemptively before each merging proposal (at least for these "simple"
> > branches; it avoids having to run the tests each time when you implement
> > something the reviewer asked for).
> >
> > So can I take this as your green light from a structural/functional
> > standpoint?
> > I can also send you the test-fixed version for you to try just before I
> > merge if you want.
> >
> > Cheers,
> >
> >
> >>
> >> Cheers,
> >>
> >>
> >> Olivier
> >> --
> >>
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/COLLIER/+merge/298349
> >> Your team MadDevelopers is subscribed to branch
> >> lp:~maddevelopers/mg5amcnlo/2.5.0.
> >>
> >
> >
> >
> > --
> > Valentin
> >
> >
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/COLLIER/+merge/298349
> > You are reviewing the proposed merge of
> lp:~maddevelopers/mg5amcnlo/COLLIER into lp:~maddevelopers/mg5amcnlo/
> 2.5.0.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/COLLIER/+merge/298349
> Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/2.5.0.
>

--
Valentin

[Valentin Hirschi (valentin-hirschi)]

Hi Olivier,

So I fixed the usual default_configuration test and the IOTests. All the
others are passing except three acceptance tests, which are also failing in
2.4.3, and which are related to the fact that the command 'compute_width'
is broken when resorting to a numerical madevent integration.

This is made clear by the fact that 'test_width_computation' fails in
2.4.3. I suspected that it was related to the fact that your removed the
files 'results.dat' which is used by 'compute_width' and indeed the test
passes in rev.274 but not in rev.275 of 2.4.3, so could you fix this?
Once this is fixed, I expect all tests to pass in COLLIER branch. Can you
run them and confirm that you are ok to merge it in 2.5.0?

Cheers

On Mon, Jun 27, 2016 at 1:24 PM, Valentin Hirschi <
<email address hidden>> wrote:

> I stressed "from a structural/functional standpoint" and also that I would
> tell you to run the tests once I will have run/fix them myself.
> So I'll let you know.
>
> On Mon, Jun 27, 2016 at 1:00 PM, Olivier Mattelaer <
> <email address hidden>> wrote:
>
> > Hi Valentin,
> >
> > It is always important that at least two people runs the test to be sure
> > that all file are include correctly
> > and that you do not have test which are machine dependent (which means
> > useless test).
> > So clearly you do not have a green light.
> >
> > Cheers,
> >
> > Olivier
> >
> > > On Jun 27, 2016, at 17:30, Valentin Hirschi <
> <email address hidden>>
> > wrote:
> > >
> > > Hi Olivier,
> > >
> > > On Mon, Jun 27, 2016 at 7:36 AM, Olivier Mattelaer <
> > > <email address hidden>> wrote:
> > >
> > >> Review: Disapprove
> > >>
> > >> Hi,
> > >>
> > >> No the tests are not working...
> > >> So you can not merge.
> > >>
> > >
> > > It goes without saying that I'll run and fix the tests before actually
> > > performing the merge...
> > > I like to do it at the very end when I'm really about to merge, and not
> > > preemptively before each merging proposal (at least for these "simple"
> > > branches; it avoids having to run the tests each time when you
> implement
> > > something the reviewer asked for).
> > >
> > > So can I take this as your green light from a structural/functional
> > > standpoint?
> > > I can also send you the test-fixed version for you to try just before I
> > > merge if you want.
> > >
> > > Cheers,
> > >
> > >
> > >>
> > >> Cheers,
> > >>
> > >>
> > >> Olivier
> > >> --
> > >>
> >
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/COLLIER/+merge/298349
> > >> Your team MadDevelopers is subscribed to branch
> > >> lp:~maddevelopers/mg5amcnlo/2.5.0.
> > >>
> > >
> > >
> > >
> > > --
> > > Valentin
> > >
> > >
> >
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/COLLIER/+merge/298349
> > > You are reviewing the proposed merge of
> > lp:~maddevelopers/mg5amcnlo/COLLIER into lp:~maddevelopers/mg5amcnlo/
> > 2.5.0.
> >
> >
> > --
> >
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/COLLIER/+merge/298349
> > Your team MadDevelopers is subscribed to branch
> > lp:~maddevelopers/mg5amcnlo/2.5.0.
> >
>
>
>
> --
> Valentin
>
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/COLLIER/+merge/298349
> Your team MadDeve...

Unmerged revisions

437. By Valentin Hirschi on 2016-07-02

1. Reinstated the 'wrong' loop_MSSM for parallel tests, but only in the input_files.
2. Merged with latest version of 2.4.3.

436. By Valentin Hirschi on 2016-06-28

1. Fixed an issue in the compilation of SMWidth
2. Updated all IOTest to correspond to a setup without Ninja Quad prec. We will do so until clang releases a version supporting __float128 (and quadmath in general).

435. By Valentin Hirschi on 2016-06-28

1. Merged with latest 2.4.3

434. By Valentin Hirschi on 2016-06-28

1. Fixed IOTests and default config tests.

433. By Valentin Hirschi on 2016-06-24

1. Changed the default value of the ML reduction order to 6|7|1.

432. By Valentin Hirschi on 2016-06-24

1. Added comments and updated the Offline HEPTools installer.

431. By Valentin Hirschi on 2016-06-24

1. Merged with latest version 2.4.3

430. By Valentin Hirschi on 2016-06-24

1.

429. By Valentin Hirschi on 2016-06-24

1. Fixed BinothLHA.f

428. By Valentin Hirschi on 2016-06-24

1. Added LARS's fix to COLLIER for rank 7+
2. Completely isolated the OneLOop of IREGI and CutTools so as to prevent interferences with other reduction tools using OneLOop.

=== modified file 'Template/LO/Source/make_opts'
--- Template/LO/Source/make_opts	2016-07-01 22:45:47 +0000
+++ Template/LO/Source/make_opts	2016-07-02 20:37:49 +0000
@@ -1,9 +1,9 @@
 DEFAULT_F2PY_COMPILER=f2py
 DEFAULT_F_COMPILER=gfortran
-MACFLAG=-mmacosx-version-min=10.7
-DEFAULT_CPP_COMPILER=clang
-STDLIB=-lc++
-STDLIB_FLAG=-stdlib=libc++
+MACFLAG=
+DEFAULT_CPP_COMPILER=g++
+STDLIB=-lstdc++
+STDLIB_FLAG=
 #end_of_make_opts_variables
 # Rest of the makefile
 

=== modified file 'Template/NLO/SubProcesses/BinothLHA.f'
--- Template/NLO/SubProcesses/BinothLHA.f	2016-03-02 09:31:23 +0000
+++ Template/NLO/SubProcesses/BinothLHA.f	2016-07-02 20:37:49 +0000
@@ -85,6 +85,8 @@
 c itself again to perform stability check to make sure no unstable EPS
 c splips unnoticed.
          CALL FORCE_STABILITY_CHECK(.TRUE.)
+         CALL COLLIER_COMPUTE_UV_POLES(.FALSE.)
+         CALL COLLIER_COMPUTE_IR_POLES(.FALSE.)
          firsttime_run = .false.
       endif
       if (firsttime) then
@@ -102,6 +104,7 @@
 c Just set the accuracy found to a positive value as it is not specified
 c once the initial pole check is performed.
          if (mc_hel.eq.0) then
+
             call sloopmatrix_thres(p,virt_wgts,tolerance,accuracies
      $           ,ret_code)
             prec_found = accuracies(0)            
@@ -164,18 +167,19 @@
 c MadLoop initialization PS points.
       cpol=.false.
       if ((firsttime .or. mc_hel.eq.0) .and. mod(ret_code,100)/10.ne.3
-     $     .and. mod(ret_code,100)/10.ne.4) then 
+     $     .and. mod(ret_code,100)/10.ne.4) then
          call getpoles(p,QES2,madfks_double,madfks_single,fksprefact)
          avgPoleRes(1)=(single+madfks_single)/2.0d0
          avgPoleRes(2)=(double+madfks_double)/2.0d0
          PoleDiff(1)=dabs(single - madfks_single)
          PoleDiff(2)=dabs(double - madfks_double)
          if ((dabs(avgPoleRes(1))+dabs(avgPoleRes(2))).ne.0d0) then
-            cpol = .not. (((PoleDiff(1)+PoleDiff(2))/
+            cpol = .not. ((((PoleDiff(1)+PoleDiff(2))/
      $           (dabs(avgPoleRes(1))+dabs(avgPoleRes(2)))) .lt.
-     $           tolerance*10d0)
+     $           tolerance*10d0).or.(mod(ret_code,10).eq.7))
          else
-            cpol = .not.(PoleDiff(1)+PoleDiff(2).lt.tolerance*10d0)
+            cpol = .not.((PoleDiff(1)+PoleDiff(2).lt.tolerance*10d0).or.
+     $                   (mod(ret_code,10).eq.7))
          endif
          if (tolerance.lt.0.0d0) then
             cpol = .false.

=== modified file 'Template/NLO/SubProcesses/check_poles.f'
--- Template/NLO/SubProcesses/check_poles.f	2016-03-02 09:31:23 +0000
+++ Template/NLO/SubProcesses/check_poles.f	2016-07-02 20:37:49 +0000
@@ -116,6 +116,8 @@
 c Make sure that stability checks are always used by MadLoop, even for
 c initialization
       CALL FORCE_STABILITY_CHECK(.TRUE.)
+      CALL COLLIER_COMPUTE_UV_POLES(.TRUE.)
+      CALL COLLIER_COMPUTE_IR_POLES(.TRUE.)
 
 200   continue
           calculatedborn = .false.

=== modified file 'Template/NLO/SubProcesses/makefile_loop.inc'
--- Template/NLO/SubProcesses/makefile_loop.inc	2016-03-02 03:03:20 +0000
+++ Template/NLO/SubProcesses/makefile_loop.inc	2016-07-02 20:37:49 +0000
@@ -24,6 +24,7 @@
 		 $(patsubst %(dotf)s,%(doto)s,$(wildcard mp_coef_construction_*.f)) \
 		 $(patsubst %(dotf)s,%(doto)s,$(wildcard TIR_interface.f)) \
 		 $(patsubst %(dotf)s,%(doto)s,$(wildcard GOLEM_interface.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard COLLIER_interface.f)) \
 		 $(patsubst %(dotf)s,%(doto)s,$(wildcard compute_color_flows.f))
 
 # This is the core of madloop computationally wise, so make sure to turn optimizations on and bound checks off.

=== modified file 'Template/loop_material/StandAlone/Cards/MadLoopParams.dat'
--- Template/loop_material/StandAlone/Cards/MadLoopParams.dat	2016-04-06 10:19:13 +0000
+++ Template/loop_material/StandAlone/Cards/MadLoopParams.dat	2016-07-02 20:37:49 +0000
@@ -4,8 +4,8 @@
 
 !     
 #MLReductionLib
-!6|1|2
-! Default :: 6|1|2
+!6|7|1
+! Default :: 6|7|1
 !     The tensor integral reduction library.The current choices are:
 !     1 | CutTools
 !     2 | PJFry++
@@ -13,6 +13,7 @@
 !     4 | Golem95
 !     5 | Samurai
 !     6 | Ninja
+!     7 | COLLIER
 !     One can use the combinations to reduce integral,e.g.
 !     1|2|3 means first use CutTools, if it is not stable, use PJFry++, 
 !     if it is still unstable, use IREGI. If it failed, use QP of CutTools.
@@ -22,15 +23,15 @@
 !     When using quadruple precision with Ninja or CutTools, the reduction will 
 !     always be done in quadruple precision, but the parameters below allow you to
 !     chose if you want to also recompute the *integrand* in quadruple precision.
-!     Doing so is slow and typically does not improve the accuracy for Ninja, so it 
-!     is disabled by default in this case.
+!     Doing so is slow but might improve the accuracy in some situation.
 #UseQPIntegrandForCutTools
 !.TRUE.
 ! Default :: .TRUE.
 #UseQPIntegrandForNinja
-!.FALSE.
-! Default :: .FALSE.
+!.TRUE.
+! Default :: .TRUE.
 !
+
 !     =================================================================================
 !     The parameters below set the parameters for IREGI
 !     =================================================================================
@@ -38,9 +39,9 @@
 #IREGIMODE
 !2
 ! Default :: 2
-!     IMODE=0, IBP reduction
-!     IMODE=1, PaVe reduction
-!     IMODE=2, PaVe reduction with stablility improved by IBP reduction
+!     IREGIMODE=0, IBP reduction
+!     IREGIMODE=1, PaVe reduction
+!     IREGIMODE=2, PaVe reduction with stablility improved by IBP reduction
 
 #IREGIRECY
 !.TRUE.
@@ -59,7 +60,7 @@
 !       2   |  Double precision, loops reduced with propagator with reversed order                               
 !       4   |  Quadruple precision, loops reduced with propagator in original order                               
 !       5   |  Quadruple precision, loops reduced with propagator with reversed order                          
-!      -1   |  Exhaustive automated numerical stability checks. See below for details.
+!      -1   | Exhaustive automated numerical stability checks. See below for details.
 !     
 !     Due to the architecture of the program, you are better off
 !     rerunning the full PS point in quadruple precision than just a single loop 
@@ -84,8 +85,8 @@
 !     Below you can chose if you want to use zero, one or two rotations of the PS point 
 !     in QP.
 #NRotations_DP
-!1
-! Default :: 1
+!0
+! Default :: 0
 #NRotations_QP
 !0
 ! Default :: 0
@@ -114,6 +115,7 @@
 #ImprovePSPoint
 !2
 ! Default :: 2
+
 !     =================================================================================
 !     The parameters below set two CutTools internal parameters accessible to the user. 
 !     =================================================================================
@@ -226,4 +228,71 @@
 !.FALSE.
 ! Default :: .FALSE.
 
+!     =================================================================================
+!     The parameters below set the main parameters for COLLIER
+!     To edit more specific technical COLLIER parameters, modify directly the content
+!     of the subroutine 'INITCOLLIER' in the file 'MadLoopCommons.f'
+!     =================================================================================
+
+!     Decide if COLLIER must be computed multiple times to evaluate the UV pole residues
+!     (Withing a Monte-Carlo performed in MG5aMC, this is automatically disabled internally)
+#COLLIERComputeUVpoles
+!.TRUE.
+! Default :: .TRUE.
+
+!     Decide if COLLIER must be computed multiple times to evaluate the IR pole residues
+!     (Withing a Monte-Carlo performed in MG5aMC, this is automatically disabled internally)
+#COLLIERComputeIRpoles
+!.TRUE.
+! Default :: .TRUE.
+
+!     Decide if COLLIER must be computed multiple times to evaluate the IR pole residues
+#COLLIERRequiredAccuracy
+!1.0d-8
+! Default :: 1.0d-8
+!     A value of -1.0d0 means that it will be automatically set from MLStabThres.
+!     The default value of 1.0d-8 corresponds to the value for which COLLIER's authors
+!     have optimized the library.
+
+!     Decide whether to use COLLIER's internal stability test or the loop-direction
+!     switch test instead.
+#COLLIERUseInternalStabilityTest
+!.TRUE.
+! Default :: .TRUE.
+!     COLLIER's internal stability test is at no extra cost but not as reliable
+!     as the loop-direction switch test, which however doubles the reduction time.
+!     This parameter is only relevant when running MadLoop with CTModeRun=-1.
+!     If you find a large number of unstable points with COLLIER for complicated
+!     processes, set this parameter to .FALSE. to make sure the PS points flagged
+!     as unstable with COLLIER really are so.
+
+!     Set up to which N-loop to use the COLLIER global caching system.
+#COLLIERGlobalCache
+!-1
+! Default :: -1
+!     -1 : Enable the global cache for all loops 
+!      0 : Disable the global cache alltogether
+!      N : Enable the global cache but only for up to N-loops
+
+!     Use the global cache when evaluating the poles as well (more memory consuming)
+!     During a Monte-Carlo it is typically not useful anyway, because the pole  
+!     computation is automatically disabled for COLLIER, irrespectively of the value
+!     of the parameters COLLIERCompute<UV/IR>poles specified above.
+#COLLIERUseCacheForPoles
+!.FALSE.
+! Default :: .FALSE.
+
+!     Choose which branch(es) of COLLIER have to be used
+#COLLIERMode
+!1
+! Default :: 1
+!     COLLIERMode=1 : COLI branch 
+!     COLLIERMode=2 : DD branch
+!     COLLIERMode=3 : Both DD and COLI branch compared
+
+!     Decide if COLLIER can output its information in a log directory.
+#COLLIERCanOutput
+!.FALSE.
+! Default :: .FALSE.
+
 /* End of param file */

=== modified file 'Template/loop_material/StandAlone/SubProcesses/MadLoopCommons.inc'
--- Template/loop_material/StandAlone/SubProcesses/MadLoopCommons.inc	2016-02-08 23:58:00 +0000
+++ Template/loop_material/StandAlone/SubProcesses/MadLoopCommons.inc	2016-07-02 20:37:49 +0000
@@ -30,6 +30,86 @@
 
       end
 
+
+## if(collier_available){
+	  SUBROUTINE INITCOLLIER()
+C 
+C INITIALISATION OF COLLIER
+C
+C  
+C MODULE
+C
+	  USE COLLIER 
+C
+C CONSTANTS
+C
+	  character(len=*) no_folder 
+      parameter (no_folder='')
+	  character(len=*) folderOutput
+      parameter (folderOutput='COLLIER_output')
+
+C     Force COLLIER to completely reset from scratch
+      logical noreset
+      parameter (noreset=.FALSE.)
+C
+C LOCAL VARIABLES 
+C    
+      INTEGER N_CACHES
+C  
+C GLOBAL VARIABLES 
+C
+      include 'MadLoopParams.inc'
+C     Now obtain the overall maximum rank, maximum external lines
+C     and maximal number of processes.
+      include 'global_specs.inc'
+
+C ----------
+C BEGIN CODE
+C ----------
+
+C     Initialize Collier
+      IF (.NOT.COLLIERCanOutput) THEN
+        CALL Init_cll(MAXNEXTERNAL,OVERALLMAXRANK,no_folder,noreset)
+      ELSE
+        CALL Init_cll(MAXNEXTERNAL,OVERALLMAXRANK,folderOutput,noreset)
+	  ENDIF
+
+C     Set target accuracy, be conservative w.r.t user request
+C     Keep in mind that COLLIER has been optimized for 1d-8.
+      IF (COLLIERRequiredAccuracy.eq.-1.0d0) THEN
+        CALL SetReqAcc_cll(MLStabThres*1.0d-3)
+      ELSE
+        CALL SetReqAcc_cll(COLLIERRequiredAccuracy)
+	  ENDIF
+
+C     Set COLLIER mode
+      CALL SetMode_cll(COLLIERMode)
+
+C     Set the global cache strategy
+C     If we use the caches for the poles as well, then it is a total of
+C     4 caches per process.
+      IF (COLLIERUseCacheForPoles) THEN
+	    N_CACHES =4*NPROCS
+	  ELSE
+	    N_CACHES =NPROCS
+	  ENDIF
+      IF (COLLIERGlobalCache.eq.-1) THEN
+	    CALL InitCacheSystem_cll(N_CACHES,MAXNEXTERNAL)
+	  ELSEIF(COLLIERGlobalCache.gt.0) THEN
+	    CALL InitCacheSystem_cll(N_CACHES,COLLIERGlobalCache)
+	  ENDIF
+
+C     Make sure to start by first switching off all cache
+      IF (COLLIERGlobalCache.ne.0) THEN
+        call SwitchOffCacheSystem_cll()
+      ENDIF
+
+C     Specify below your other custom COLLIER parameter settings 
+C     [user_specific_COLLIER_settings]
+	  
+	  END
+## }
+
       SUBROUTINE SET_FORBID_HEL_DOUBLECHECK(ONOFF)
 C
 C Give the possibility to overwrite the value of MadLoopParams.dat
@@ -67,9 +147,10 @@
       DATA ML_INIT/.TRUE./
       common/ML_INIT/ML_INIT
 
-      LOGICAL CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT
-      DATA CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT/.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE./
-      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
+      LOGICAL CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT,COLLIERINIT
+      DATA CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT,COLLIERINIT/.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE./
+      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT, COLLIERINIT
+
 
       character(512) MLPath
       data MLPath/'[[NA]]'/      
@@ -163,6 +244,8 @@
 C         Stable with Samurai.
 C     U == 6
 C         Stable with Ninja in double precision.
+C     U == 7
+C         Stable with COLLIER.
 C     U == 8
 C         Stable with Ninja in quadruple precision.
 C     U == 9
@@ -201,10 +284,10 @@
             STOP 'Only CutTools and Ninja can use quardruple precision'
          ENDIF
       ENDIF
-      IF(MLRed.GE.1.AND.MLRed.LE.6)THEN
+      IF(MLRed.GE.1.AND.MLRed.LE.7)THEN
          SET_RET_CODE_U=MLRed
       ELSE
-         STOP 'Only CutTools, PJFry++, IREGI, Golem95, Samurai and Ninja are available'
+         STOP 'Only CutTools, PJFry++, IREGI, Golem95, Samurai, Ninja and COLLIER are available'
       ENDIF
       END
 
@@ -249,8 +332,11 @@
       ELSEIF(LIBNUM.EQ.6)THEN
 C Ninja 
          CALL DETECT_NINJA(NLOOPLINE,RANK,complex_mass,LPASS)
+      ELSEIF(LIBNUM.EQ.7)THEN
+C Collier 
+         CALL DETECT_COLLIER(NLOOPLINE,RANK,complex_mass,LPASS)
       ELSE
-         STOP 'ONLY CUTTOOLS,PJFry++,IREGI,Golem95 and Samurai are available'
+         STOP 'Only CutTools, PJFry++, IREGI, Golem95, Samurai, Ninja and COLLIER are available'
       ENDIF
       RETURN
       END
@@ -358,6 +444,65 @@
       RETURN
       END
 
+## if(collier_available){
+      SUBROUTINE DETECT_COLLIER(NLOOPLINE,RANK,complex_mass,LPASS)
+C
+C Detect whether Collier can be used or not
+C
+      USE COLLIER
+      IMPLICIT NONE
+C
+C CONSTANTS
+C
+C
+C ARGUMENTS
+C
+      INTEGER NLOOPLINE,RANK
+      LOGICAL complex_mass,LPASS
+C
+C LOCAL VARIABLES
+C
+      INTEGER CURRENT_COLLIERMODE
+C
+C GLOBAL VARIABLES
+C
+      LOGICAL CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT,COLLIERINIT
+      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT, COLLIERINIT
+      include 'MadLoopParams.inc'
+C ----------
+C BEGIN CODE
+C ----------
+      if (.NOT.COLLIERINIT) THEN
+        CALL GETMODE_CLL(CURRENT_COLLIERMODE)
+      ELSE
+        CURRENT_COLLIERMODE = COLLIERMODE
+      ENDIF
+      LPASS=.TRUE.
+	  IF (CURRENT_COLLIERMODE.NE.1) THEN
+C       The DD branch is used and it has limitations
+        IF((NLOOPLINE.GT.6).OR.(RANK.GT.NLOOPLINE)) then
+	      LPASS=.FALSE.
+        endif
+      ELSE
+C       Limitations of the COLI branch are academic.
+        LPASS=.TRUE. 
+	  ENDIF
+      RETURN
+      END
+## }else{
+      SUBROUTINE DETECT_COLLIER(NLOOPLINE,RANK,complex_mass,LPASS)
+C
+C ARGUMENTS
+C
+      INTEGER NLOOPLINE,RANK
+      LOGICAL complex_mass,LPASS
+C
+C     COLLIER is not available in this output. This subroutine is dummy.
+C
+      LPASS=.True.
+      END
+## }
+
       SUBROUTINE DETECT_PJFRY(NLOOPLINE,RANK,complex_mass,LPASS)
 C
 C DETECT whether PJFRY++ CAN BE USED OR NOT

=== modified file 'Template/loop_material/StandAlone/SubProcesses/MadLoopParamReader.f'
--- Template/loop_material/StandAlone/SubProcesses/MadLoopParamReader.f	2016-02-18 11:42:56 +0000
+++ Template/loop_material/StandAlone/SubProcesses/MadLoopParamReader.f	2016-07-02 20:37:49 +0000
@@ -33,6 +33,12 @@
                stop 'CTModeRun must be >= -1 and <=6.'
              endif 
 
+           else if (buff .eq. '#COLLIERGlobalCache') then
+             read(666,*,end=999) COLLIERGlobalCache
+             if (COLLIERGlobalCache .lt. -1) then
+               stop 'COLLIERGlobalCache must be >= -1'
+             endif
+
            else if (buff .eq. '#NRotations_DP') then
              read(666,*,end=999) NRotations_DP
              if (NRotations_DP .lt. 0 .or.
@@ -52,7 +58,14 @@
              if (MLStabThres.lt.0.0d0) then
                  stop 'MLStabThres must be >= 0'
              endif
-          
+
+           else if (buff .eq. '#COLLIERRequiredAccuracy') then
+             read(666,*,end=999) COLLIERRequiredAccuracy
+             if (COLLIERRequiredAccuracy.le.0.0d0.and.
+     &           COLLIERRequiredAccuracy.ne.-1.0d0) then
+                 stop 'COLLIERRequiredAccuracy must be > 0 or = -1.0'
+             endif
+
            else if (buff .eq. '#CTLoopLibrary') then
              read(666,*,end=999) CTLoopLibrary
              if (CTLoopLibrary.lt.2 .or.
@@ -90,6 +103,21 @@
                  stop 'MaxAttempts must be >= 1'
              endif
 
+          else if (buff .eq. '#COLLIERComputeUVpoles') then
+             read(666,*,end=999) COLLIERComputeUVpoles
+
+          else if (buff .eq. '#COLLIERComputeIRpoles') then
+             read(666,*,end=999) COLLIERComputeIRpoles
+
+		  else if (buff .eq. '#COLLIERUseInternalStabilityTest') then
+             read(666,*,end=999) COLLIERUseInternalStabilityTest
+
+          else if (buff .eq. '#COLLIERUseCacheForPoles') then
+             read(666,*,end=999) COLLIERUseCacheForPoles
+
+          else if (buff .eq. '#COLLIERCanOutput') then
+             read(666,*,end=999) COLLIERCanOutput
+
           else if (buff .eq. '#UseLoopFilter') then
              read(666,*,end=999) UseLoopFilter
 
@@ -127,7 +155,7 @@
 
           else if (buff .eq. '#MLReductionLib') then
              read(666,*,end=999) MLReductionLib_str
-             MLReductionLib(1:6)=0
+             MLReductionLib(1:7)=0
              MLReductionLib_str_save=MLReductionLib_str
              j=0
              DO
@@ -149,6 +177,8 @@
                    MLRed=5
                 ELSEIF(MLReductionLib_char.EQ.'6 ')THEN
                    MLRed=6
+                ELSEIF(MLReductionLib_char.EQ.'7 ')THEN
+                   MLRed=7
                 ELSE
                    PRINT *, 'MLReductionLib is wrong: '//
      $                  TRIM(MLReductionLib_str_save)
@@ -171,6 +201,12 @@
                    MLReductionLib_str=MLReductionLib_str(i+1:)
                 ENDIF
              ENDDO
+          else if (buff .eq. '#COLLIERMode') then
+             read(666,*,end=999) COLLIERMode
+             if (COLLIERMode .lt. 1 .or.
+     &            COLLIERMode .gt.3) then
+                stop 'COLLIERMode must be >=1 and <=3.'
+             endif
           else if (buff .eq. '#IREGIRECY') then
              read(666,*,end=999) IREGIRECY
           else if (buff .eq. '#IREGIMODE') then
@@ -218,8 +254,6 @@
      & '==============================================================='
       write(*,*) ' > MLReductionLib            = '
      $     //TRIM(MLReductionLib_str_save)
-      write(*,*) ' > IREGIMODE                 = ',IREGIMODE
-      write(*,*) ' > IREGIRECY                 = ',IREGIRECY
       write(*,*) ' > CTModeRun                 = ',CTModeRun
       write(*,*) ' > MLStabThres               = ',MLStabThres
       write(*,*) ' > NRotations_DP             = ',NRotations_DP
@@ -243,6 +277,19 @@
      &UseQPIntegrandForNinja
       write(*,*) ' > UseQPIntegrandForCutTools = ',
      &UseQPIntegrandForCutTools
+      write(*,*) ' > IREGIMODE                 = ',IREGIMODE
+      write(*,*) ' > IREGIRECY                 = ',IREGIRECY
+      write(*,*) ' > COLLIERMode               = ',COLLIERMode
+      write(*,*) ' > COLLIERRequiredAccuracy   = ',
+     $COLLIERRequiredAccuracy
+      write(*,*) ' > COLLIERCanOutput          = ',COLLIERCanOutput
+      write(*,*) ' > COLLIERComputeUVpoles     = ',COLLIERComputeUVpoles
+      write(*,*) ' > COLLIERComputeIRpoles     = ',COLLIERComputeIRpoles
+      write(*,*) ' > COLLIERGlobalCache        = ',COLLIERGlobalCache
+      write(*,*) ' > COLLIERUseCacheForPoles   = ',
+     &COLLIERUseCacheForPoles
+      write(*,*) ' > COLLIERUseInternalStabilityTest = ',
+     &COLLIERUseInternalStabilityTest
       write(*,*)
      & '==============================================================='
       paramPrinted=.TRUE.
@@ -259,12 +306,22 @@
       include "MadLoopParams.inc"
 
       MLReductionLib(1)=6
-      MLReductionLib(2:6)=0
+      MLReductionLib(2)=7
+      MLReductionLib(3)=1
+      MLReductionLib(4:7)=0
       IREGIMODE=2
       IREGIRECY=.TRUE.
+      COLLIERComputeIRpoles = .TRUE.
+      COLLIERComputeUVpoles = .TRUE.
+      COLLIERUseCacheForPoles = .FALSE.
+      COLLIERCanOutput = .FALSE.
+      COLLIERGlobalCache = -1
+      COLLIERMode=1
+      COLLIERRequiredAccuracy=1.0d-8
+      COLLIERUseInternalStabilityTest = .TRUE.
       CTModeInit=0
       CTModeRun=-1
-      NRotations_DP=1
+      NRotations_DP=0
       NRotations_QP=0
       MLStabThres=1.0d-3
       CTStabThres=1.0d-2
@@ -281,6 +338,6 @@
       OSThres=1.0d-13
       ImprovePSPoint=2
       UseQPIntegrandForCutTools=.True.
-      UseQPIntegrandForNinja=.False.
+      UseQPIntegrandForNinja=.True.
 
       end

=== modified file 'Template/loop_material/StandAlone/SubProcesses/MadLoopParams.inc'
--- Template/loop_material/StandAlone/SubProcesses/MadLoopParams.inc	2016-02-18 11:42:56 +0000
+++ Template/loop_material/StandAlone/SubProcesses/MadLoopParams.inc	2016-07-02 20:37:49 +0000
@@ -7,17 +7,24 @@
 !
 	  integer CTModeInit,CTModeRun,CheckCycle,MaxAttempts,
      &CTLoopLibrary,NRotations_DP,NRotations_QP,ImprovePSPoint,
-     &MLReductionLib(7),IREGIMODE,HelicityFilterLevel
-
-      real*8 MLStabThres,CTStabThres,ZeroThres,OSThres
-
-      logical UseLoopFilter,LoopInitStartOver,DoubleCheckHelicityFilter
+     &MLReductionLib(8),IREGIMODE,HelicityFilterLevel,COLLIERMode,
+     &COLLIERGlobalCache
+
+      real*8 MLStabThres,CTStabThres,ZeroThres,OSThres,COLLIERRequiredAccuracy
+
+      logical UseLoopFilter,LoopInitStartOver,DoubleCheckHelicityFilter,
+     &COLLIERComputeIRpoles,COLLIERComputeUVpoles,COLLIERCanOutput
 	  logical HelInitStartOver,IREGIRECY,WriteOutFilters
       logical UseQPIntegrandForNinja, UseQPIntegrandForCutTools
+	  logical COLLIERUseCacheForPoles,COLLIERUseInternalStabilityTest
 
       common /MADLOOP/CTModeInit,CTModeRun,NRotations_DP,NRotations_QP,
+     &COLLIERMode,COLLIERGlobalCache,
      &ImprovePSPoint,CheckCycle, MaxAttempts,UseLoopFilter,MLStabThres,
+     &COLLIERRequiredAccuracy,
      &CTStabThres,CTLoopLibrary,LoopInitStartOver,
+     &COLLIERComputeIRpoles,COLLIERComputeUVpoles,COLLIERCanOutput,
+     &COLLIERUseCacheForPoles,COLLIERUseInternalStabilityTest,
      &DoubleCheckHelicityFilter,ZeroThres,OSThres,HelInitStartOver,
      &MLReductionLib,IREGIMODE,HelicityFilterLevel,IREGIRECY,
      &WriteOutFilters,UseQPIntegrandForNinja,UseQPIntegrandForCutTools

=== modified file 'Template/loop_material/StandAlone/SubProcesses/makefile'
--- Template/loop_material/StandAlone/SubProcesses/makefile	2016-05-26 23:36:59 +0000
+++ Template/loop_material/StandAlone/SubProcesses/makefile	2016-07-02 20:37:49 +0000
@@ -45,6 +45,7 @@
  $(patsubst $(DOTF),$(DOTO),$(wildcard mp_coef_construction_*.f)) \
  $(patsubst $(DOTF),$(DOTO),$(wildcard TIR_interface.f)) \
  $(patsubst $(DOTF),$(DOTO),$(wildcard GOLEM_interface.f)) \
+ $(patsubst $(DOTF),$(DOTO),$(wildcard COLLIER_interface.f)) \
  $(patsubst $(DOTF),$(DOTO),$(wildcard compute_color_flows.f))
 
 OLP_PROCESS= MadLoopParamReader.o MadLoopCommons.o \
@@ -64,6 +65,7 @@
  $(patsubst $(DOTF),$(DOTO),$(wildcard $(LOOP_PREFIX)*/mp_coef_construction_*.f)) \
  $(patsubst $(DOTF),$(DOTO),$(wildcard $(LOOP_PREFIX)*/TIR_interface.f)) \
  $(patsubst $(DOTF),$(DOTO),$(wildcard $(LOOP_PREFIX)*/GOLEM_interface.f)) \
+ $(patsubst $(DOTF),$(DOTO),$(wildcard $(LOOP_PREFIX)*/COLLIER_interface.f)) \
  $(patsubst $(DOTF),$(DOTO),$(wildcard $(LOOP_PREFIX)*/compute_color_flows.f))
 
 $(PROG): check_sa.o $(PROCESS) makefile $(LIBS)

=== modified file 'input/.mg5_configuration_default.txt'
--- input/.mg5_configuration_default.txt	2016-04-28 19:07:40 +0000
+++ input/.mg5_configuration_default.txt	2016-07-02 20:37:49 +0000
@@ -180,12 +180,17 @@
 #! if samurai=/PATH/TO/samurai/lib: use that specific installation path for samurai
 # samurai = None
 
-#! Set the samurai directory containing ninja's library
-#! if auto: try to find it automatically on the system (default)
+#! Set the Ninja directory containing ninja's library
 #! if '' or None: disabling ninja 
 #! if ninja=/PATH/TO/ninja/lib: use that specific installation path for ninja 
 # ninja = ./HEPTools/lib
 
+#! Set the COLLIER directory containing COLLIER's library
+#! if '' or None: disabling COLLIER 
+#! if ninja=/PATH/TO/ninja/lib: use that specific installation path for COLLIER
+# Note that it is necessary that you have generated a static library for COLLIER
+# collier = ./HEPTools/lib 
+
 #! Set how MadLoop dependencies (such as CutTools) should be handled
 #!  > external : ML5 places a link to the MG5_aMC-wide libraries
 #!  > internal : ML5 copies all dependencies in the output so that it is independent

=== modified file 'madgraph/interface/amcatnlo_run_interface.py'
--- madgraph/interface/amcatnlo_run_interface.py	2016-07-02 10:13:48 +0000
+++ madgraph/interface/amcatnlo_run_interface.py	2016-07-02 20:37:49 +0000
@@ -1964,7 +1964,6 @@
         else:
             logger.info('The results of this run' + \
                         ' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
-        devnull.close()
 
     def combine_plots_HwU(self,jobs,out,normalisation=None):
         """Sums all the plots in the HwU format."""
@@ -2402,6 +2401,7 @@
                               4 : 'Golem95',
                               5 : 'Samurai',
                               6 : 'Ninja (double precision)',
+                              7 : 'COLLIER',
                               8 : 'Ninja (quadruple precision)',
                               9 : 'CutTools (quadruple precision)'}
         RetUnit_finder =re.compile(

=== modified file 'madgraph/interface/loop_interface.py'
--- madgraph/interface/loop_interface.py	2016-05-24 09:52:16 +0000
+++ madgraph/interface/loop_interface.py	2016-07-02 20:37:49 +0000
@@ -538,10 +538,11 @@
 
         # Fortran MadGraph5_aMC@NLO Standalone
         if self._export_format in self.supported_ML_format:
-            for me in matrix_elements:
+            for unique_id, me in enumerate(matrix_elements):
                 calls = calls + \
                         self._curr_exporter.generate_subprocess_directory_v4(\
-                            me, self._curr_fortran_model)
+                            me, self._curr_fortran_model, (unique_id+1))
+            self._curr_exporter.write_global_specs(matrix_elements)
             # If all ME's do not share the same maximum loop vertex rank and the
             # same loop maximum wavefunction size, we need to set the maximum
             # in coef_specs.inc of the HELAS Source. The SubProcesses/P* directory
@@ -554,6 +555,7 @@
                 if len(set(max_lwfspins))>1 or len(set(max_loop_vert_ranks))>1:
                     self._curr_exporter.fix_coef_specs(max(max_lwfspins),\
                                                        max(max_loop_vert_ranks))
+                self._curr_exporter.write_global_specs(matrix_elements)
 
         # Just the matrix.f files
         if self._export_format == 'matrix':

=== modified file 'madgraph/interface/madevent_interface.py'
--- madgraph/interface/madevent_interface.py	2016-06-28 11:12:01 +0000
+++ madgraph/interface/madevent_interface.py	2016-07-02 20:37:49 +0000
@@ -2540,7 +2540,7 @@
 #           helicity filters.
             self.MadLoopparam.set('CheckCycle',4, ifnotdefault=False)
             
-            # For now it is tricky to have eahc channel performing the helicity
+            # For now it is tricky to have each channel performing the helicity
             # double check. What we will end up doing is probably some kind
             # of new initialization round at the beginning of each launch
             # command, to reset the filters.    
@@ -2548,7 +2548,7 @@
                                                              ifnotdefault=False)
           
             # Thanks to TIR recycling, TIR is typically much faster for Loop-induced
-            # processes, so that we place OPP last.
+            # processes when not doing MC over helicities, so that we place OPP last.
             if not hasattr(self, 'run_card'):
                 run_card = banner_mod.RunCard(opt['run_card'])
             else:
@@ -2561,7 +2561,7 @@
     """You chose to set the preferred reduction technique in MadLoop to be OPP (see parameter MLReductionLib).
     Beware that this can bring significant slowdown; the optimal choice --when not MC over helicity-- being to first start with TIR reduction.""")
                 # We do not include GOLEM for now since it cannot recycle TIR coefs yet.
-                self.MadLoopparam.set('MLReductionLib','2|6|1', ifnotdefault=False)
+                self.MadLoopparam.set('MLReductionLib','7|6|1', ifnotdefault=False)
             else:
                 if 'MLReductionLib' in self.MadLoopparam.user_set and \
                     not (self.MadLoopparam.get('MLReductionLib').startswith('1') or
@@ -2569,7 +2569,7 @@
                     logger.warning(
     """You chose to set the preferred reduction technique in MadLoop to be different than OPP (see parameter MLReductionLib).
     Beware that this can bring significant slowdown; the optimal choice --when MC over helicity-- being to first start with OPP reduction.""")
-                self.MadLoopparam.set('MLReductionLib','6|1|2', ifnotdefault=False)
+                self.MadLoopparam.set('MLReductionLib','6|7|1', ifnotdefault=False)
 
             # Also TIR cache will only work when NRotations_DP=0 (but only matters
             # when not MC-ing over helicities) so it will be hard-reset by MadLoop

=== modified file 'madgraph/interface/madgraph_interface.py'
--- madgraph/interface/madgraph_interface.py	2016-06-28 11:12:01 +0000
+++ madgraph/interface/madgraph_interface.py	2016-07-02 20:37:49 +0000
@@ -910,8 +910,11 @@
         
         user_options = {'--energy':'1000','--split_orders':'-1',
                    '--reduction':'1|2|3|4|5|6','--CTModeRun':'-1',
-                   '--helicity':'-1','--seed':'-1'}
-        
+                   '--helicity':'-1','--seed':'-1','--collier_cache':'-1',
+                   '--collier_req_acc':'auto',
+                   '--collier_internal_stability_test':'False',
+                   '--collier_mode':'1'}  
+
         if args[0] in ['cms'] or args[0].lower()=='cmsoptions':
             # increase the default energy to 5000
             user_options['--energy']='5000'
@@ -2638,7 +2641,7 @@
                      'update', 'Delphes2', 'SysCalc', 'Golem95', 'PJFry',
                                                                       'QCDLoop']
     # The targets below are installed using the HEPToolsInstaller.py script
-    _advanced_install_opts = ['ninja']
+    _advanced_install_opts = ['ninja','collier']
     
     # The options below are commented for now but already available
 #    _advanced_install_opts += ['pythia8','zlib','boost','lhapdf6','lhapdf5','hepmc','mg5amc_py8_interface','oneloop']
@@ -2693,6 +2696,7 @@
                        'golem':'auto',
                        'samurai':None,
                        'ninja':'./HEPTools/lib',
+                       'collier':'./HEPTools/lib',
                        'lhapdf':'lhapdf-config',
                        'applgrid':'applgrid-config',
                        'amcfast':'amcfast-config',
@@ -3535,6 +3539,15 @@
                                        " must be an integer, not %s."%option[1])
             elif option[0]=='--reduction':
                 MLoptions['MLReductionLib']=[int(ir) for ir in option[1].split('|')]
+            elif option[0]=='--collier_mode':
+                MLoptions['COLLIERMode']=int(option[1])
+            elif option[0]=='--collier_cache':
+                MLoptions['COLLIERGlobalCache']=int(option[1])
+            elif option[0]=='--collier_req_acc':
+                if option[1]!='auto':
+                    MLoptions['COLLIERRequiredAccuracy']=float(option[1])
+            elif option[0]=='--collier_internal_stability_test':
+                MLoptions['COLLIERUseInternalStabilityTest']=eval(option[1])                
             elif option[0]=='--CTModeRun':
                 try:
                     MLoptions['CTModeRun']=int(option[1])  
@@ -3899,6 +3912,14 @@
                     logger_check.warning('Samurai not available on your system; it will be skipped.')
                     MLoptions["MLReductionLib"].remove(5)
         
+        if 'collier' in self.options and isinstance(self.options['collier'],str):
+            TIR_dir['collier_dir']=self.options['collier']
+        else:
+            if "MLReductionLib" in MLoptions:
+                if 7 in MLoptions["MLReductionLib"]:
+                    logger_check.warning('Collier not available on your system; it will be skipped.')
+                    MLoptions["MLReductionLib"].remove(7)
+        
         if 'ninja' in self.options and isinstance(self.options['ninja'],str):
             TIR_dir['ninja_dir']=self.options['ninja']
         else:
@@ -5398,7 +5419,9 @@
             self.options['mg5amc_py8_interface_path'] = \
                                  pjoin(MG5DIR,'HEPTools','MG5aMC_PY8_interface')
             self.exec_cmd('save options')      
-
+        elif tool == 'collier':
+            self.options['collier'] = pjoin(os.curdir,'HEPTools','lib')
+            self.exec_cmd('save options')
         elif tool == 'ninja':
             if not misc.get_ninja_quad_prec_support(pjoin(
                                               MG5DIR,'HEPTools','ninja','lib')):
@@ -5497,6 +5520,7 @@
                           'hepmc':['CPC 134 (2001) 41-46'],
                           'mg5amc_py8_interface':['arXiv:1410.3012','arXiv:XXXX.YYYYY'],
                           'ninja':['arXiv:1203.0291','arXiv:1403.1229','arXiv:1604.01363'],
+                          'collier':['arXiv:1604.06792'],
                           'oneloop':['arXiv:1007.4716']}
 
         if args[0] in advertisements:
@@ -6867,8 +6891,8 @@
                 logger.info('set fastjet to %s' % args[1])
                 self.options[args[0]] = args[1]
 
-        elif args[0] in ['pjfry','golem','samurai','ninja'] and \
-                           not (args[0]=='ninja' and args[1]=='./HEPTools/lib'):
+        elif args[0] in ['pjfry','golem','samurai','ninja','collier'] and \
+             not (args[0] in ['ninja','collier'] and args[1]=='./HEPTools/lib'):
             if args[1] in ['None',"''",'""']:
                 self.options[args[0]] = None
             else:
@@ -7268,6 +7292,8 @@
         if self._export_format == 'madevent':
             calls += self._curr_exporter.export_processes(self._curr_matrix_elements,
                                                  self._curr_fortran_model)
+            self._curr_exporter.write_global_specs(
+                               self._curr_matrix_elements.get_matrix_elements())
             
             # Write the procdef_mg5.dat file with process info
             card_path = pjoin(path, os.path.pardir, 'SubProcesses', \

=== modified file 'madgraph/iolibs/export_fks.py'
--- madgraph/iolibs/export_fks.py	2016-05-12 10:46:19 +0000
+++ madgraph/iolibs/export_fks.py	2016-07-02 20:37:49 +0000
@@ -187,7 +187,8 @@
         writer = writers.FortranWriter(os.path.join(self.dir_path, 
                                              "SubProcesses","MadLoopCommons.f"))
         writer.writelines(MadLoopCommon%{
-                                   'print_banner_commands':self.MadLoop_banner})
+                                   'print_banner_commands':self.MadLoop_banner},
+                                            context={'collier_available':False})
         writer.close()
                                        
         # Write the cts_mpc.h and cts_mprec.h files imported from CutTools
@@ -1573,12 +1574,19 @@
             matrix_element.get('processes')[0].nice_string(print_weighted=False))
         plot.draw()
 
+        # We also need to write the overall maximum quantities for this group
+        # of processes in 'global_specs.inc'. In aMCatNLO, there is always
+        # only one process, so this is trivial
+        self.write_global_specs(matrix_element)
+        open('unique_id.inc','w').write(
+"""      integer UNIQUE_ID
+      parameter(UNIQUE_ID=1)""")
+
         linkfiles = ['coupl.inc', 'mp_coupl.inc', 'mp_coupl_same_name.inc',
                      'cts_mprec.h', 'cts_mpc.h', 'MadLoopParamReader.f',
-                     'MadLoopCommons.f','MadLoopParams.inc']
+                     'MadLoopCommons.f','MadLoopParams.inc','global_specs.inc']
 
         # We should move to MadLoop5_resources directory from the SubProcesses
-
         ln(pjoin(os.path.pardir,os.path.pardir,'MadLoopParams.dat'),
                                               pjoin('..','MadLoop5_resources'))
 
@@ -3090,17 +3098,31 @@
                                               libpath,"lib%s.a"%tir,tir_name=tir)
             setattr(self,tir_dir,libpath)
             if libpath != "":
-                if tir in ['pjfry','ninja','golem', 'samurai']:
+                if tir in ['pjfry','ninja','golem', 'samurai','collier']:
                     # We should link dynamically when possible, so we use the original
                     # location of these libraries.
                     link_tir_libs.append('-L%s/ -l%s'%(libpath,tir))
                     tir_libs.append('%s/lib%s.$(libext)'%(libpath,tir))
-                    # We must add the corresponding includes for golem, samurai
-                    # and ninja
-                    if tir in ['golem','samurai','ninja']:
+                    # For Ninja, we must also link against OneLoop.
+                    if tir in ['ninja']:
+                        if not any(os.path.isfile(pjoin(libpath,'libavh_olo.%s'%ext)) 
+                                              for ext in ['a','dylib','so']):
+                            raise MadGraph5Error(
+"The OneLOop library 'libavh_olo.(a|dylib|so)' could no be found in path '%s'. Please place a symlink to it there."%libpath)
+                        link_tir_libs.append('-L%s/ -l%s'%(libpath,'avh_olo'))
+                        tir_libs.append('%s/lib%s.$(libext)'%(libpath,'avh_olo'))
+                    # We must add the corresponding includes for these TIR
+                    if tir in ['golem','samurai','ninja','collier']:
                         trg_path = pjoin(os.path.dirname(libpath),'include')
-                        to_include = misc.find_includes_path(trg_path,
+                        if os.path.isdir(trg_path):
+                            to_include = misc.find_includes_path(trg_path,
                                                         self.include_names[tir])
+                        else:
+                            to_include = None
+                        # Special possible location for collier
+                        if to_include is None and tir=='collier':
+                            to_include = misc.find_includes_path(
+                               pjoin(libpath,'modules'),self.include_names[tir])
                         if to_include is None:
                             logger.error(
 'Could not find the include directory for %s, looking in %s.\n' % (tir ,str(trg_path))+
@@ -3170,7 +3192,8 @@
         writer = writers.FortranWriter(os.path.join(self.dir_path, 
                                              "SubProcesses","MadLoopCommons.f"))
         writer.writelines(MadLoopCommon%{
-                                   'print_banner_commands':self.MadLoop_banner})
+                                   'print_banner_commands':self.MadLoop_banner},
+               context={'collier_available':self.tir_available_dict['collier']})
         writer.close()
 
         # link the files from the MODEL
@@ -3280,14 +3303,21 @@
                                                           print_weighted=False))
         plot.draw()
 
+        # We also need to write the overall maximum quantities for this group
+        # of processes in 'global_specs.inc'. In aMCatNLO, there is always
+        # only one process, so this is trivial
+        self.write_global_specs(matrix_element)
+        open('unique_id.inc','w').write(
+"""      integer UNIQUE_ID
+      parameter(UNIQUE_ID=1)""")
+
         linkfiles = ['coupl.inc', 'mp_coupl.inc', 'mp_coupl_same_name.inc',
                      'cts_mprec.h', 'cts_mpc.h', 'MadLoopParamReader.f',
-                     'MadLoopParams.inc','MadLoopCommons.f']
+                     'MadLoopParams.inc','MadLoopCommons.f','global_specs.inc']
 
         for file in linkfiles:
             ln('../../%s' % file)
-
-
+                
         os.system("ln -s ../../makefile_loop makefile")
         
 # We should move to MadLoop5_resources directory from the SubProcesses

=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2016-04-16 23:51:09 +0000
+++ madgraph/iolibs/export_v4.py	2016-07-02 20:37:49 +0000
@@ -102,9 +102,11 @@
         
         calls = 0
         if isinstance(matrix_elements, group_subprocs.SubProcessGroupList):
+            unique_id=1
             for (group_number, me_group) in enumerate(matrix_elements):
                 calls = calls + self.generate_subprocess_directory_v4(\
-                                          me_group, fortran_model, group_number)
+                    me_group, fortran_model, group_number, unique_id=unique_id)
+                unique_id += len(me_group.get('matrix_elements'))
         else:
             for me_number, me in enumerate(matrix_elements.get_matrix_elements()):
                 calls = calls + self.generate_subprocess_directory_v4(\
@@ -112,7 +114,16 @@
                         
         return calls    
         
-
+    #===========================================================================
+    #  Generate an include file with global quantities about all ME's output.
+    #===========================================================================
+    def write_global_specs(self, matrix_elements):
+        """ Writes the file global_specs.inc which contains general information
+        about *all* the ME's output in this directory."""
+        
+        # Do nothing here, but daughter classes such as LoopOptimizedExporterFortran
+        # overwrites this.
+        pass
 
     #===========================================================================
     #  create the run_card 
@@ -4462,10 +4473,12 @@
     #===========================================================================
     def generate_subprocess_directory_v4(self, subproc_group,
                                          fortran_model,
-                                         group_number):
+                                         group_number,
+                                         unique_id=None):
         """Generate the Pn directory for a subprocess group in MadEvent,
         including the necessary matrix_N.f files, configs.inc and various
-        other helper files"""
+        other helper files. Unique_id is dummy here, but not in 
+        LoopInducedExporterMEGroup, so it must be kept."""
 
         assert isinstance(subproc_group, group_subprocs.SubProcessGroup), \
                                       "subproc_group object not SubProcessGroup"
@@ -6261,34 +6274,58 @@
     else:
         curr_proc = None
 
-    requires_ninja = opt['loop_optimized_output'] and (not curr_proc is None) and \
-      (curr_proc.get('perturbation_couplings') != [] and \
-      not curr_proc.get('NLO_mode') in [None,'real','tree','LO','LOonly'])
+    requires_reduction_tool = opt['loop_optimized_output'] and \
+                (not curr_proc is None) and \
+                (curr_proc.get('perturbation_couplings') != [] and \
+                not curr_proc.get('NLO_mode') in [None,'real','tree','LO','LOonly'])
     # An installation is required then, but only if the specified path is the
     # default local one and that the Ninja library appears missing.
-    if requires_ninja and (not opt['ninja'] is None) and\
+    if requires_reduction_tool:
+        if (not opt['ninja'] is None) and\
             os.path.abspath(opt['ninja'])==pjoin(MG5DIR,'HEPTools','lib') and\
             not os.path.isfile(pjoin(MG5DIR,'HEPTools','lib','libninja.a')):
-                # Then install Ninja here from the tarballs in the vendor
-                # directory so that it would work offline too.
-                logger.info(
+            # Then install Ninja here from the tarballs in the vendor
+            # directory so that it would work offline too.
+            logger.info(
 """MG5aMC will now install the loop reduction tool 'Ninja' from the local offline installer.
 Use the command 'install ninja' if you want to update to the latest online version.
 This installation can take some time but only needs to be performed once.""",'$MG:color:GREEN')
-                try:
-                    cmd.do_install('ninja',paths={'HEPToolsInstaller':
+            try:
+                cmd.do_install('ninja',paths={'HEPToolsInstaller':
                         pjoin(MG5DIR,'vendor','OfflineHEPToolsInstaller.tar.gz')},
-                  additional_options=[
+                additional_options=[
                   '--ninja_tarball=%s'%pjoin(MG5DIR,'vendor','ninja.tar.gz'),
                   '--oneloop_tarball=%s'%pjoin(MG5DIR,'vendor','oneloop.tar.gz')])
-                except InvalidCmd:
+            except InvalidCmd:
                     logger.warning(
 """The offline installation of Ninja was unsuccessful, and MG5aMC disabled it.
 In the future, if you want to reactivate Ninja, you can do so by re-attempting
 its online installation with the command 'install ninja' or install it on your
 own and set the path to its library in the MG5aMC option 'ninja'.""")
                     cmd.exec_cmd("set ninja ''")
-                    cmd.exec_cmd('save options')  
+                    cmd.exec_cmd('save options')
+        if (not opt['collier'] is None) and\
+            os.path.abspath(opt['collier'])==pjoin(MG5DIR,'HEPTools','lib') and\
+            not os.path.isfile(pjoin(MG5DIR,'HEPTools','lib','libcollier.a')):
+            # Then install Collier here from the tarballs in the vendor
+            # directory so that it would work offline too.
+            logger.info(
+"""MG5aMC will now install the loop reduction tool 'COLLIER' from the local offline installer.
+Use the command 'install collier' if you want to update to the latest online version.
+This installation can take some time but only needs to be performed once.""",'$MG:color:GREEN')
+            try:
+                cmd.do_install('collier',paths={'HEPToolsInstaller':
+                        pjoin(MG5DIR,'vendor','OfflineHEPToolsInstaller.tar.gz')},
+                additional_options=[
+                  '--collier_tarball=%s'%pjoin(MG5DIR,'vendor','collier.tar.gz')])
+            except InvalidCmd:
+                    logger.warning(
+"""The offline installation of COLLIER was unsuccessful, and MG5aMC disabled it.
+In the future, if you want to reactivate COLLIER, you can do so by re-attempting
+its online installation with the command 'install collier' or install it on your
+own and set the path to its library in the MG5aMC option 'collier'.""")
+                    cmd.exec_cmd("set collier ''")
+                    cmd.exec_cmd('save options')
 
 
     # ==========================================================================
@@ -6304,6 +6341,7 @@
       'golem_dir':cmd.options['golem'],
       'samurai_dir':cmd.options['samurai'],
       'ninja_dir':cmd.options['ninja'],
+      'collier_dir':cmd.options['collier'],
       'fortran_compiler':cmd.options['fortran_compiler'],
       'f2py_compiler':cmd.options['f2py_compiler'],
       'output_dependencies':cmd.options['output_dependencies'],
@@ -6431,10 +6469,12 @@
     #===========================================================================
     def generate_subprocess_directory_v4(self, subproc_group,
                                          fortran_model,
-                                         group_number):
+                                         group_number,
+                                         unique_id=None):
         """Generate the Pn directory for a subprocess group in MadEvent,
         including the necessary matrix_N.f files, configs.inc and various
-        other helper files"""
+        other helper files. Unique_id is dummy here, but not in 
+        LoopInducedExporterMEGroup, so it must be kept."""
 
         if not isinstance(subproc_group, group_subprocs.SubProcessGroup):
             raise base_objects.PhysicsObject.PhysicsObjectError,\

=== modified file 'madgraph/iolibs/template_files/loop/check_sa.inc'
--- madgraph/iolibs/template_files/loop/check_sa.inc	2016-02-17 21:24:01 +0000
+++ madgraph/iolibs/template_files/loop/check_sa.inc	2016-07-02 20:37:49 +0000
@@ -262,6 +262,7 @@
           write(*,*) 'Estimated accuracy could not be computed for an unknown reason.'
         endif
       endif
+	  Write (*,'(1x,a23,3x,i3)') 'MadLoop return code   =',RETURNCODE
       write (*,*) "---------------------------------"
 	  IF (NBORNCHOSEN.EQ.0) THEN
         write (*,*) "No Born contribution satisfied the squared order constraints."

=== modified file 'madgraph/iolibs/template_files/loop/check_sa_loop_induced.inc'
--- madgraph/iolibs/template_files/loop/check_sa_loop_induced.inc	2016-02-17 21:24:01 +0000
+++ madgraph/iolibs/template_files/loop/check_sa_loop_induced.inc	2016-07-02 20:37:49 +0000
@@ -241,6 +241,7 @@
           write(*,*) 'Estimated accuracy could not be computed for an unknown reason.'
         endif
       endif
+	  Write (*,'(1x,a23,3x,i3)') 'MadLoop return code   =',RETURNCODE
       write (*,*) "---------------------------------"
       IF (NLOOPCHOSEN.ne.NSQUAREDSO) THEN
         write (*,*) "Selected squared coupling orders combination for the loop summed result below:"

=== added file 'madgraph/iolibs/template_files/loop_optimized/COLLIER_interface.inc'
--- madgraph/iolibs/template_files/loop_optimized/COLLIER_interface.inc	1970-01-01 00:00:00 +0000
+++ madgraph/iolibs/template_files/loop_optimized/COLLIER_interface.inc	2016-07-02 20:37:49 +0000
@@ -0,0 +1,646 @@
+	  SUBROUTINE %(proc_prefix)sCOLLIERLOOP(CTMODE, NLOOPLINE,RANK,PL,PDEN,M2L,TIRCOEFS,TIRCOEFSERRORS)
+C  
+%(info_lines)s
+C 
+C Interface between MG5 and COLLIER.
+C It supports any rank when 1 < NLOOPLINE < 7.
+C  
+%(process_lines)s
+C
+C
+C MODULES
+C
+	  USE COLLIER
+	  IMPLICIT NONE
+C  
+C CONSTANTS 
+C  
+      LOGICAL checkPConservation
+      PARAMETER (checkPConservation=.TRUE.)
+      %(real_dp_format)s PIVALUE
+      PARAMETER(PIVALUE=3.14159265358979323846D0)
+	  integer MAXRANK
+	  parameter (MAXRANK=%(maxrank)d)
+      INTEGER NLOOPGROUPS
+      PARAMETER (NLOOPGROUPS=%(nloop_groups)d)
+      include 'loop_max_coefs.inc'	  
+C  
+C ARGUMENTS 
+C
+      INTEGER NLOOPLINE, RANK, CTMODE
+      %(real_dp_format)s PL(0:3,NLOOPLINE)
+      %(real_dp_format)s PDEN(0:3,NLOOPLINE-1)	  
+      %(mass_dp_format)s M2L(NLOOPLINE)
+      %(complex_dp_format)s RES(3)
+	  %(complex_dp_format)s TIRCOEFS(0:LOOPMAXCOEFS-1,3)
+	  %(complex_dp_format)s TIRCOEFSERRORS(0:LOOPMAXCOEFS-1,3)
+
+C  
+C LOCAL VARIABLES 
+C
+	  INTEGER N, I, J, K, L
+	  INTEGER C0,C1,C2,C3
+	  INTEGER N_CACHES
+	  INTEGER CURR_RANK, SGN
+	  INTEGER CURR_INDEX
+	  INTEGER CURR_MAXCOEF
+	  %(real_dp_format)s RBUFF(0:3)
+
+      INTEGER COEFMAP_ZERO(0:LOOPMAXCOEFS-1)
+      INTEGER COEFMAP_ONE(0:LOOPMAXCOEFS-1)
+      INTEGER COEFMAP_TWO(0:LOOPMAXCOEFS-1)
+      INTEGER COEFMAP_THREE(0:LOOPMAXCOEFS-1)	  
+	  %(collier_coefmap)s
+
+	  %(real_dp_format)s REF_NORMALIZATION
+
+	  double complex M2LCOLLIER(0:NLOOPLINE-1)
+      double complex MomVec(0:3,NLOOPLINE-1)
+	  double complex MomInv((NLOOPLINE*(NLOOPLINE-1))/2)
+	 
+      double complex TNten(0:RANK,0:RANK,0:RANK,0:RANK)
+      double complex TNtenuv(0:RANK,0:RANK,0:RANK,0:RANK)
+      double precision TNtenerr(0:RANK)
+
+	  %(real_dp_format)s MaxCoefForRank(0:RANK)
+
+C     These quantities are for the pole evaluation
+      double complex TNten_UV(0:RANK,0:RANK,0:RANK,0:RANK)
+      double complex TNtenuv_UV(0:RANK,0:RANK,0:RANK,0:RANK)
+      double precision TNtenerr_UV(0:RANK)
+	  double complex TNten_IR1(0:RANK,0:RANK,0:RANK,0:RANK)
+      double complex TNtenuv_IR1(0:RANK,0:RANK,0:RANK,0:RANK)
+      double precision TNtenerr_IR1(0:RANK)
+      double complex TNten_IR2(0:RANK,0:RANK,0:RANK,0:RANK)
+      double complex TNtenuv_IR2(0:RANK,0:RANK,0:RANK,0:RANK)
+      double precision TNtenerr_IR2(0:RANK)
+C
+C GLOBAL VARIABLES
+C
+      include 'coupl.inc'
+      include 'MadLoopParams.inc'
+	  include 'unique_id.inc'
+	  include 'global_specs.inc'
+
+## if(not AmplitudeReduction) {
+      INTEGER ID,SQSOINDEX,R
+	  common/%(proc_prefix)sLOOP/ID,SQSOINDEX,R
+## }else{
+      INTEGER ID,R
+	  common/%(proc_prefix)sLOOP/ID,R
+## }
+
+C
+C     These global variables cache the coefficients already computed by COLLIER
+C     so that they can be reused. In particular, in CTMODE 2, the cached ERROR 
+C     quantities for eahc rank will be used to make use of COLLIER's internal
+C     assessment of the numerical accuracy.
+C
+      logical LOOP_ID_DONE(NLOOPGROUPS)
+	  DATA LOOP_ID_DONE/NLOOPGROUPS*.FALSE./
+	  %(complex_dp_format)s TIR_COEFS_DIRECT_MODE(0:LOOPMAXCOEFS-1,3,NLOOPGROUPS)
+	  %(real_dp_format)s TIR_COEFS_ERROR(0:MAXRANK,NLOOPGROUPS)
+	  common/%(proc_prefix)sCOLLIER_TIR_COEFS/TIR_COEFS_DIRECT_MODE,TIR_COEFS_ERROR,LOOP_ID_DONE
+
+      INTEGER COLLIER_CACHE_ACTIVE, NCALLS_IN_CACHE(4), NCOLLIERCALLS(4)
+      LOGICAL MUST_INIT_EVENT
+	  COMMON/%(proc_prefix)sCOLLIER_CACHE_STATUS/COLLIER_CACHE_ACTIVE,NCALLS_IN_CACHE, NCOLLIERCALLS,MUST_INIT_EVENT
+
+      LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT, COLLIERINIT
+      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT, COLLIERINIT
+
+C ----------
+C BEGIN CODE
+C ----------
+
+	  IF (COLLIERUseCacheForPoles) THEN
+	    N_CACHES = 4
+	  ELSE
+	    N_CACHES =1
+	  ENDIF
+
+C     Initialize COLLIER if needed 
+      IF (COLLIERINIT) THEN
+	    COLLIERINIT=.FALSE.
+	    CALL INITCOLLIER()
+      ENDIF
+
+C     Initialize the event if it is the first time collier is called on this PS point
+      IF(MUST_INIT_EVENT) THEN
+	    MUST_INIT_EVENT = .False.
+        IF (COLLIERGlobalCache.eq.0) THEN
+	      CALL InitEvent_cll()
+	    ELSE
+		  DO I=1,N_CACHES
+C           Record how many events where put in the cache. On the first PS point.
+		    IF (NCALLS_IN_CACHE(I).eq.-1.AND.NCOLLIERCALLS(I).gt.0) THEN
+		      NCALLS_IN_CACHE(I) = NCOLLIERCALLS(I)
+		    ENDIF
+	      ENDDO
+		  DO I=1,N_CACHES		  
+C           Now apply a safety check that our last event had as many calls as the cache is setup for.
+C           The only case for now when it can be half of the calls when we are doing the true loop-direction test with also the computation of a rotated PS point (which is computed for one mode only).
+            IF (NCALLS_IN_CACHE(I).ne.-1.and..NOT. ( NCOLLIERCALLS(I).eq.NCALLS_IN_CACHE(I).or. ( CTModeRun.eq.-1.and..not.COLLIERUseInternalStabilityTest.and.NRotations_DP.gt.0.and.MOD(NCALLS_IN_CACHE(I),2).eq.0.and.(NCALLS_IN_CACHE(I)/2).eq.NCOLLIERCALLS(I) ) ) ) then
+		      WRITE(*,*) 'WARNING: A consistency check in MadLoop failed and, for safety, forced MadLoop to reinitialize the global cache of COLLIER. Report this to MadLoop authors. The problematic cache was number ',I
+              IF (COLLIERGLOBALCACHE.EQ.-1) THEN
+                CALL INITCACHESYSTEM_CLL(N_CACHES*NPROCS,MAXNEXTERNAL)
+              ELSE
+		        CALL INITCACHESYSTEM_CLL(N_CACHES*NPROCS,COLLIERGLOBALCACHE)
+              ENDIF
+C             Make sure all caches are switched off at first.
+              call SwitchOffCacheSystem_cll()
+C             Reset the cache design property
+		      NCALLS_IN_CACHE(:) = -1
+			  NCOLLIERCALLS(:)   = 0
+		      IF (COLLIER_CACHE_ACTIVE.eq.1) THEN			  
+		        CALL SwitchOnCache_cll((unique_id-1)*N_CACHES+1)
+			  ENDIF
+C             No need to check the other caches since we already had to reset here.
+			  EXIT
+		    ENDIF
+		  ENDDO
+	      CALL InitEvent_cll((UNIQUE_ID-1)*N_CACHES+1)
+		  NCOLLIERCALLS(1) = 0		  
+		  IF(COLLIERComputeUVpoles.and.COLLIERUseCacheForPoles) THEN
+	        CALL InitEvent_cll((UNIQUE_ID-1)*N_CACHES+2)
+			NCOLLIERCALLS(2) = 0
+		  ENDIF
+		  IF(COLLIERComputeIRpoles.and.COLLIERUseCacheForPoles) THEN
+	        CALL InitEvent_cll((UNIQUE_ID-1)*N_CACHES+3)
+			NCOLLIERCALLS(3) = 0
+	        CALL InitEvent_cll((UNIQUE_ID-1)*N_CACHES+4)
+			NCOLLIERCALLS(4) = 0
+		  ENDIF
+        ENDIF
+        CALL SetDeltaIR_cll(0.0d0,(PIVALUE**2)/6.0d0)
+        CALL SetDeltaUV_cll(0.0d0)
+        CALL SetMuIR2_cll(mu_r**2)      
+        CALL SetMuUV2_cll(mu_r**2)
+	  ENDIF
+
+C     Now really start the reduction with COLLIER
+
+C     Number of coefficients for the current rank
+      CURR_MAXCOEF = 0
+      DO I=0,RANK
+        CURR_MAXCOEF=CURR_MAXCOEF+(3+I)*(2+I)*(1+I)/6
+      ENDDO
+
+	  IF (CTMODE.ne.1.and.CTMODE.ne.2) THEN
+	    WRITE(*,*) 'ERROR: COLLIER only support the computational mode CTMODE equal to 1 or 2, not',CTMODE
+	    stop 'Incorrect computational mode specified to the COLLIER MG5aMC interface.'
+	  ENDIF
+
+      DO I=0,CURR_MAXCOEF-1
+        DO K=1,3			
+          TIRCOEFSERRORS(I,K)=DCMPLX(0.0d0,0.0d0)
+	    ENDDO
+      ENDDO
+
+	  IF (COLLIERUseInternalStabilityTest) THEN
+C       Use MADLOOP internal cache dedicated to COLLIER that emulates the CTMODE 2
+	    IF (LOOP_ID_DONE(ID)) THEN
+		  DO I=0,CURR_MAXCOEF-1
+		    DO K=1,3
+			  TIRCOEFS(I,K)=TIR_COEFS_DIRECT_MODE(I,K,ID)
+			ENDDO
+	      ENDDO
+		  DO I=0,RANK
+			MaxCoefForRank(I) = 0.0d0
+	      ENDDO
+		  DO I=0,CURR_MAXCOEF-1
+		    CURR_RANK = COEFMAP_ZERO(I)+COEFMAP_ONE(I)+COEFMAP_TWO(I)+COEFMAP_THREE(I)			  
+			MaxCoefForRank(CURR_RANK)=max(MaxCoefForRank(CURR_RANK),abs(TIR_COEFS_DIRECT_MODE(I,1,ID)))
+		  ENDDO
+		  DO I=0,CURR_MAXCOEF-1
+		    CURR_RANK = COEFMAP_ZERO(I)+COEFMAP_ONE(I)+COEFMAP_TWO(I)+COEFMAP_THREE(I)
+            IF (MaxCoefForRank(CURR_RANK).ne.0.0d0) THEN
+              DO K=1,3
+C               The expression below is like taking the absolute value when summing errors linearly
+C                TIRCOEFSERRORS(I,K)=(TIR_COEFS_ERROR(CURR_RANK,ID)/MaxCoefForRank(CURR_RANK))*DCMPLX( ABS(DBLE(TIRCOEFS(I,K))),ABS(DIMAG(TIRCOEFS(I,K))) )
+C               But empirically, I observed that retaining the original complex phase leads to slightly more accurate estimates
+			    TIRCOEFSERRORS(I,K)=(TIR_COEFS_ERROR(CURR_RANK,ID)/MaxCoefForRank(CURR_RANK))*TIRCOEFS(I,K)           
+			  ENDDO
+            ENDIF	
+	      ENDDO    
+		  RETURN
+		ENDIF	    
+      ELSE
+C       Apply the loop-direction switching here.
+        CALL %(proc_prefix)sSWITCH_ORDER(CTMODE,NLOOPLINE,PL,PDEN,M2L)	  
+	  ENDIF
+
+C     Make sure masses are complex 
+      do I=1,NLOOPLINE
+	    M2LCOLLIER(I-1)=DCMPLX(M2L(I))
+	  ENDDO
+
+C     Now convert the loop offset momenta to COLLIER conventions
+      DO i=0,3
+	    DO j=1,NLOOPLINE-1
+	      MomVec(i,j)=DCMPLX(PDEN(i,j),0.0d0)
+	    ENDDO
+      ENDDO
+
+C     This first do loop spans over 'N' in '\foreach_N \foreach_i(k_i+k_{i+1}+..+k_{i+N})^2'
+      CURR_INDEX = 0
+      DO N=0,NLOOPLINE-1
+C       We stop whenever we reached the target number of invariants
+        IF (CURR_INDEX.GE.((NLOOPLINE*(NLOOPLINE-1))/2)) THEN
+		  EXIT
+		ENDIF
+C       This do loop spans over 'i' in the expression above.
+        DO I=1,NLOOPLINE
+          IF (CURR_INDEX.GE.((NLOOPLINE*(NLOOPLINE-1))/2)) THEN
+		    EXIT
+		  ENDIF
+
+          CURR_INDEX = CURR_INDEX+1
+          RBUFF(:) = 0.0D0
+          DO J=I,I+N
+            RBUFF(:) = RBUFF(:) + PL(:,MOD(J-1,NLOOPLINE)+1)
+          ENDDO
+          MomInv(CURR_INDEX) = DCMPLX(RBUFF(0)**2-RBUFF(1)**2-RBUFF(2)**2-RBUFF(3)**2,0.0d0)
+
+C         Now regularize the onshell behavior of the kinematic invarients
+C         All loop masses are tested, but that might be a bit too inclusive.
+	      DO K=1,NLOOPLINE
+		    IF(ABS(M2L(K)).NE.0.0d0) THEN
+              IF(ABS((MomInv(CURR_INDEX)-M2L(K))/M2L(K)).LT.OSTHRES) THEN
+                MomInv(CURR_INDEX)=DCMPLX(M2L(K))
+              ENDIF
+		    ENDIF
+		  ENDDO
+C         For the massless onshell-case, we base the threshold only on the energy component
+          REF_NORMALIZATION=0.0D0
+          DO L=0,0
+            REF_NORMALIZATION = REF_NORMALIZATION + ABS(RBUFF(L))
+          ENDDO
+          REF_NORMALIZATION = (REF_NORMALIZATION/(N+1))**2
+          IF(REF_NORMALIZATION.NE.0.0D0)THEN
+            IF(ABS(MomInv(CURR_INDEX)/REF_NORMALIZATION).LT.OSTHRES)THEN
+              MomInv(CURR_INDEX)=DCMPLX(0.0d0,0.0d0)
+            ENDIF
+          ENDIF
+
+		ENDDO
+	  ENDDO
+
+C     We can now call COLLIER
+	  IF (NLOOPLINE.ne.1) THEN
+        CALL TNten_cll(TNten, TNtenuv, MomVec, MomInv, M2LCOLLIER, NLOOPLINE, RANK, TNtenerr)
+      ELSE
+        CALL TNten_cll(TNten, TNtenuv, M2LCOLLIER, NLOOPLINE, RANK, TNtenerr)
+	  ENDIF
+	  IF (COLLIER_CACHE_ACTIVE.eq.1) THEN
+	    NCOLLIERCALLS(1) = NCOLLIERCALLS(1)+1
+	  ENDIF
+
+C     Now compute the UV poles if asked for
+      IF (COLLIERComputeUVpoles) THEN
+	    IF(COLLIER_CACHE_ACTIVE.eq.1) THEN
+		  CALL SwitchOffCache_cll((unique_id-1)*N_CACHES+1)
+		  IF(COLLIERUseCacheForPoles) THEN
+		    CALL SwitchOnCache_cll((unique_id-1)*N_CACHES+2)
+		  ENDIF
+        ENDIF
+        CALL SetDeltaUV_cll(1.0d0)
+	    IF (NLOOPLINE.ne.1) THEN
+          CALL TNten_cll(TNten_UV, TNtenuv_UV, MomVec, MomInv, M2LCOLLIER, NLOOPLINE, RANK, TNtenerr_UV)
+        ELSE
+          CALL TNten_cll(TNten_UV, TNtenuv_UV, M2LCOLLIER, NLOOPLINE, RANK, TNtenerr_UV)
+	    ENDIF
+	    IF (COLLIER_CACHE_ACTIVE.eq.1) THEN
+	      NCOLLIERCALLS(2) = NCOLLIERCALLS(2)+1
+	    ENDIF
+	    IF(COLLIER_CACHE_ACTIVE.eq.1) THEN
+		  IF(COLLIERUseCacheForPoles) THEN
+		    CALL SwitchOffCache_cll((unique_id-1)*N_CACHES+2)
+		  ENDIF
+		  CALL SwitchOnCache_cll((unique_id-1)*N_CACHES+1)
+        ENDIF
+        CALL SetDeltaUV_cll(0.0d0)
+	  ENDIF
+
+C     Now compute the IR poles if asked for
+      IF (COLLIERComputeIRpoles) THEN
+	    IF(COLLIER_CACHE_ACTIVE.eq.1) THEN
+		  CALL SwitchOffCache_cll((unique_id-1)*N_CACHES+1)
+		  IF(COLLIERUseCacheForPoles) THEN
+		    CALL SwitchOnCache_cll((unique_id-1)*N_CACHES+3)
+		  ENDIF
+        ENDIF
+        CALL SetDeltaIR_cll(1.0d0,(PIVALUE**2)/6.0d0)
+	    IF (NLOOPLINE.ne.1) THEN
+          CALL TNten_cll(TNten_IR1, TNtenuv_IR1, MomVec, MomInv, M2LCOLLIER, NLOOPLINE, RANK, TNtenerr_IR1)
+        ELSE
+          CALL TNten_cll(TNten_IR1, TNtenuv_IR1, M2LCOLLIER, NLOOPLINE, RANK, TNtenerr_IR1)
+	    ENDIF
+	    IF (COLLIER_CACHE_ACTIVE.eq.1) THEN
+	      NCOLLIERCALLS(3) = NCOLLIERCALLS(3)+1
+	    ENDIF
+        IF(COLLIER_CACHE_ACTIVE.eq.1.AND.COLLIERUseCacheForPoles) THEN
+		  CALL SwitchOffCache_cll((unique_id-1)*N_CACHES+3)
+		  CALL SwitchOnCache_cll((unique_id-1)*N_CACHES+4)			
+        ENDIF
+        CALL SetDeltaIR_cll(0.0d0,1.0d0+(PIVALUE**2)/6.0d0)
+	    IF (NLOOPLINE.ne.1) THEN
+          CALL TNten_cll(TNten_IR2, TNtenuv_IR2, MomVec, MomInv, M2LCOLLIER, NLOOPLINE, RANK, TNtenerr_IR2)
+        ELSE
+          CALL TNten_cll(TNten_IR2, TNtenuv_IR2, M2LCOLLIER, NLOOPLINE, RANK, TNtenerr_IR2)
+	    ENDIF
+	    IF (COLLIER_CACHE_ACTIVE.eq.1) THEN
+	      NCOLLIERCALLS(4) = NCOLLIERCALLS(4)+1
+	    ENDIF
+	    IF(COLLIER_CACHE_ACTIVE.eq.1) THEN
+		  IF(COLLIERUseCacheForPoles) THEN
+		    CALL SwitchOffCache_cll((unique_id-1)*N_CACHES+4)
+		  ENDIF
+		  CALL SwitchOnCache_cll((unique_id-1)*N_CACHES+1)		  
+        ENDIF
+		CALL SetDeltaIR_cll(0.0d0,(PIVALUE**2)/6.0d0)
+	  ENDIF
+
+	  DO I=0,(CURR_MAXCOEF-1)
+	    C0 = COEFMAP_ZERO(I)
+		C1 = COEFMAP_ONE(I)
+		C2 = COEFMAP_TWO(I)
+		C3 = COEFMAP_THREE(I)
+	    CURR_RANK = C0+C1+C2+C3
+C       Because we must set q -> -q, we apply a minus sign to coefs of odd rank
+	    IF (MOD(CURR_RANK,2).eq.1) THEN
+		  SGN = -1
+		ELSE
+		  SGN = 1
+		ENDIF
+		TIRCOEFS(I,1) = SGN*TNten(C0,C1,C2,C3)
+        IF (COLLIERComputeUVpoles) THEN
+		  TIRCOEFS(I,2) = SGN*( TNten_UV(C0,C1,C2,C3)-TNten(C0,C1,C2,C3) )
+		ELSE
+		  TIRCOEFS(I,2) = DCMPLX(0.0d0,0.0d0)
+		ENDIF
+        IF (COLLIERComputeIRpoles) THEN
+		  TIRCOEFS(I,2) = TIRCOEFS(I,2) + SGN*( TNten_IR1(C0,C1,C2,C3)-TNten(C0,C1,C2,C3) )
+		  TIRCOEFS(I,3) = SGN*( TNten_IR2(C0,C1,C2,C3)-TNten(C0,C1,C2,C3) )
+		ELSE
+		  TIRCOEFS(I,3) = DCMPLX(0.0d0,0.0d0)		
+		ENDIF
+      ENDDO
+
+	  IF (COLLIERUseInternalStabilityTest) THEN
+c       Finish by caching the coefficients and error just computed
+	    LOOP_ID_DONE(ID) = .TRUE.
+	    DO J=0,CURR_MAXCOEF-1
+	      DO K=1,3
+	        TIR_COEFS_DIRECT_MODE(J,K,ID) = TIRCOEFS(J,K)
+	      ENDDO
+	    ENDDO
+	    DO J=0,RANK
+	      TIR_COEFS_ERROR(J,ID)=TNtenerr(J)
+	    ENDDO
+
+C       Now compute the errors on each coefficient
+        DO I=0,RANK
+          MaxCoefForRank(I) = 0.0d0
+        ENDDO
+        DO I=0,CURR_MAXCOEF-1
+          CURR_RANK = COEFMAP_ZERO(I)+COEFMAP_ONE(I)+COEFMAP_TWO(I)+COEFMAP_THREE(I)			  
+          MaxCoefForRank(CURR_RANK)=max(MaxCoefForRank(CURR_RANK),abs(TIRCOEFS(I,1)))
+        ENDDO
+        DO I=0,CURR_MAXCOEF-1
+          CURR_RANK = COEFMAP_ZERO(I)+COEFMAP_ONE(I)+COEFMAP_TWO(I)+COEFMAP_THREE(I)
+          IF (MaxCoefForRank(CURR_RANK).ne.0.0d0) THEN
+            DO K=1,3
+C               The expression below is like taking the absolute value when summing errors linearly
+C			  TIRCOEFSERRORS(I,K)=(TNtenerr(CURR_RANK)/MaxCoefForRank(CURR_RANK))*DCMPLX( ABS(DBLE(TIRCOEFS(I,K))),ABS(DIMAG(TIRCOEFS(I,K))) )			
+C               But empirically, I observed that retaining the original complex phase leads to slightly more accurate estimates
+			  TIRCOEFSERRORS(I,K)=(TNtenerr(CURR_RANK)/MaxCoefForRank(CURR_RANK))*TIRCOEFS(I,K)
+			ENDDO
+          ENDIF
+        ENDDO
+	    
+	  ENDIF
+
+	  END SUBROUTINE %(proc_prefix)sCOLLIERLOOP 
+
+      SUBROUTINE %(proc_prefix)sCLEAR_COLLIER_CACHE()
+	    
+	    USE COLLIER
+
+        include 'loop_max_coefs.inc'
+        INTEGER NLOOPGROUPS
+        PARAMETER (NLOOPGROUPS=%(nloop_groups)d)
+	    integer MAXRANK
+	    parameter (MAXRANK=%(maxrank)d)
+		
+		integer I,J,K
+
+        include 'MadLoopParams.inc'
+
+        logical LOOP_ID_DONE(NLOOPGROUPS)
+	    %(complex_dp_format)s TIR_COEFS_DIRECT_MODE(0:LOOPMAXCOEFS-1,3,NLOOPGROUPS)
+	    %(complex_dp_format)s TIR_COEFS_ERROR(0:MAXRANK,NLOOPGROUPS)
+	    common/%(proc_prefix)sCOLLIER_TIR_COEFS/TIR_COEFS_DIRECT_MODE,TIR_COEFS_ERROR,LOOP_ID_DONE
+
+        INTEGER COLLIER_CACHE_ACTIVE, NCALLS_IN_CACHE(4), NCOLLIERCALLS(4)
+        LOGICAL MUST_INIT_EVENT
+	    COMMON/%(proc_prefix)sCOLLIER_CACHE_STATUS/COLLIER_CACHE_ACTIVE,NCALLS_IN_CACHE, NCOLLIERCALLS,MUST_INIT_EVENT
+
+C       Make sure that next time the COLLIER Subroutine is called it will call the subroutine initEvent of Collier.
+        MUST_INIT_EVENT = .True.
+
+C       Reinitialize the ML cache for COLLIER
+        DO I=1,NLOOPGROUPS
+		  LOOP_ID_DONE(I) = .FALSE.
+		  DO J=0,LOOPMAXCOEFS-1
+		    DO K=1,3
+			  TIR_COEFS_DIRECT_MODE(J,K,I) = DCMPLX(0.0d0,0.0d0)
+			ENDDO
+		  ENDDO
+		  DO J=0,MAXRANK
+		    TIR_COEFS_ERROR(J,I)=0.0d0
+	      ENDDO
+		ENDDO
+
+	  END
+
+	  SUBROUTINE %(proc_prefix)sSET_COLLIER_GLOBAL_CACHE(ONOFF)
+C
+C       This routine is used by loop_matrix.f to turn on the global
+C       cache of COLLIER when it the main SLOOP subroutine starts and
+C       turn it off when it ends.
+C       However several checks are performed to verify that it is safe
+C       to turn it on and to reinitialize it if necessary.
+C
+C       MODULES
+C
+        USE COLLIER
+		implicit none
+C
+C       ARGUMENTS
+C
+        LOGICAL ONOFF
+C
+C       LOCAL VARIABLES
+C
+        LOGICAL NEED_REINITIALIZATION
+		INTEGER N_CACHES
+C		
+C       GLOBAL VARIABLES
+C
+C
+C       This common blocks saves the relevant ML parameters when activating the
+C       global cache of COLLIER so that we know when we must reinitialize it.
+C       COLLIER_CACHE_ACTIVE is set to -1 when it has never been turned on yet and
+C       to 1 for 'Active' and 0 for 'Inactive'.
+C       The integer NCALLS_IN_CACHE saves how many calls the cache is setup for, for each of the up to four caches.
+C       When it is the first PS points, it is set to -1.
+        INTEGER COLLIER_CACHE_ACTIVE, NCALLS_IN_CACHE(4), NCOLLIERCALLS(4)
+	    DATA COLLIER_CACHE_ACTIVE/-1/
+		DATA NCALLS_IN_CACHE/-1,-1,-1,-1/
+		DATA NCOLLIERCALLS/0,0,0,0/
+C       This is a flag to tell the COLLIER subroutine that it must init the event when called.
+		LOGICAL MUST_INIT_EVENT
+		DATA MUST_INIT_EVENT/.TRUE./
+		COMMON/%(proc_prefix)sCOLLIER_CACHE_STATUS/COLLIER_CACHE_ACTIVE, NCALLS_IN_CACHE, NCOLLIERCALLS,MUST_INIT_EVENT
+
+		LOGICAL COLLIERUseInternalStabilityTest_BU
+		INTEGER USERHEL_BU, SQSO_TARGET_BU, COLLIERMode_BU,CTModeRun_BU
+        COMMON/%(proc_prefix)sCOLLIER_CACHE_RELEVANT_PARAMS/USERHEL_BU,SQSO_TARGET_BU,COLLIERMode_BU,CTModeRun_BU,COLLIERUseInternalStabilityTest_BU
+
+C       The common blocks below are to retrieve the necessary information about
+C       MadLoop running mode and store it in the sCOLLIER_CACHE_RELEVANT_PARAMS common block.
+
+        include 'MadLoopParams.inc'
+		include 'unique_id.inc'
+		include 'global_specs.inc'
+
+        LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT, COLLIERINIT
+        COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT, COLLIERINIT
+
+	    LOGICAL CHECKPHASE
+	    LOGICAL HELDOUBLECHECKED
+        common/%(proc_prefix)sINIT/CHECKPHASE, HELDOUBLECHECKED
+
+	    INTEGER USERHEL
+	    common/%(proc_prefix)sUSERCHOICE/USERHEL
+
+	    INTEGER SQSO_TARGET
+	    common/%(proc_prefix)sSOCHOICE/SQSO_TARGET
+		
+C
+C       BEGIN CODE
+C
+
+		IF (COLLIERUseCacheForPoles) THEN
+		  N_CACHES = 4
+		ELSE
+		  N_CACHES =1
+		ENDIF
+
+C       Do nothing if COLLIER still has to be initialized or if global caches are disabled
+        IF(COLLIERINIT.OR.COLLIERGlobalCache.eq.0) THEN
+		  RETURN
+		ENDIF
+
+C       Never activate anything in the checkphase
+        IF (ONOFF.AND.CHECKPHASE) THEN
+		  RETURN
+		ENDIF
+
+C       Handle the request of turning off the caching
+		IF (.NOT.ONOFF) THEN
+		  IF (COLLIER_CACHE_ACTIVE.eq.1) THEN
+		    CALL SwitchOffCache_cll((unique_id-1)*4+1)
+		    COLLIER_CACHE_ACTIVE = 0
+		  ENDIF
+C         If we were asked to turn the cache off but it was already so, then do nothing		  
+		  RETURN
+		ENDIF
+
+C       Handle the request of turning on the caching
+
+C       If asked to activate it but already active, then do nothing
+		IF (ONOFF.AND.COLLIER_CACHE_ACTIVE.eq.1) THEN
+		  RETURN
+        ENDIF
+
+C       We are now in the position where we are asked to activate the global cache but it was *not* already active.
+
+C       If we activate it for the first time, make sure to store the value of the relevant parameters, activate and return.
+        IF (COLLIER_CACHE_ACTIVE.eq.-1) THEN
+		  USERHEL_BU         = USERHEL 
+		  SQSO_TARGET_BU     = SQSO_TARGET
+		  COLLIERMode_BU     = COLLIERMode
+		  COLLIERUseInternalStabilityTest_BU = COLLIERUseInternalStabilityTest
+		  CTModeRun_BU       = CTModeRun
+		  CALL SwitchOnCache_cll((unique_id-1)*4+1)
+		  COLLIER_CACHE_ACTIVE = 1
+		  RETURN
+		ENDIF
+
+C       Now perform sanity check before the activation to decide if we need to reinitialize the cache system first.
+		NEED_REINITIALIZATION = .FALSE.
+
+        IF (SQSO_TARGET.NE.SQSO_TARGET_BU) THEN
+		  NEED_REINITIALIZATION = .TRUE.		
+		ENDIF
+
+        IF (COLLIERMode.NE.COLLIERMode_BU) THEN
+		  NEED_REINITIALIZATION = .TRUE.		
+		ENDIF
+
+        IF (COLLIERUseInternalStabilityTest.NEQV.COLLIERUseInternalStabilityTest_BU) THEN
+		  NEED_REINITIALIZATION = .TRUE.		
+		ENDIF
+
+		IF (CTModeRun_BU.ne.CTModeRun.and.(.not.COLLIERUseInternalStabilityTest)) THEN
+		  NEED_REINITIALIZATION = .TRUE.			  
+		ENDIF
+
+## if(AmplitudeReduction){
+C       When doing amplitude reduction, the parameter USERHEL impacts the number of COLLIER calls/order if
+C       it involves a different number of helicity configurations or always if the Loop filter is on, which
+C       it should almost never be by now.
+        IF(UseLoopFilter.and.(USERHEL.ne.USERHEL_BU)) THEN
+		  NEED_REINITIALIZATION = .TRUE.
+		ENDIF
+		IF((USERHEL.eq.-1).and.(USERHEL_BU.ne.-1)) THEN
+		  NEED_REINITIALIZATION = .TRUE.
+		ENDIF
+		IF((USERHEL.ne.-1).and.(USERHEL_BU.eq.-1)) THEN
+		  NEED_REINITIALIZATION = .TRUE.
+		ENDIF
+## }else{
+C       When doing amplitude reduction the parameter USERHEL does not impact the number/order of COLLIER calls
+C       except if the LoopFilter is ON which really shouldn't be the case anymore.
+        IF(UseLoopFilter.and.(USERHEL.ne.USERHEL_BU)) THEN
+		  NEED_REINITIALIZATION = .TRUE.
+		ENDIF
+## }
+
+		IF(NEED_REINITIALIZATION) THEN
+C         Log the event because if it happens a lot of time and floods the screen, the user must see it
+C         and either change its usage of MadLoop or turnoff COLLIER cache
+          WRITE(*,*) 'INFO: MadLoop detected that the global cache of COLLIER had to be reset because of a change in your use of MadLoop. This should not happend for each event.'
+		  USERHEL_BU         = USERHEL 
+		  SQSO_TARGET_BU     = SQSO_TARGET
+		  COLLIERMode_BU     = COLLIERMode
+		  COLLIERUseInternalStabilityTest_BU = COLLIERUseInternalStabilityTest
+		  CTModeRun_BU       = CTModeRun		  
+		  IF (COLLIERGLOBALCACHE.EQ.-1) THEN
+            CALL INITCACHESYSTEM_CLL(N_CACHES*NPROCS,MAXNEXTERNAL)
+          ELSE
+		    CALL INITCACHESYSTEM_CLL(N_CACHES*NPROCS,COLLIERGLOBALCACHE)
+          ENDIF
+		  NCOLLIERCALLS(:)   = 0
+		  NCALLS_IN_CACHE(:) = -1 
+C         Make sure all caches are switched off at first.
+          call SwitchOffCacheSystem_cll()
+		ENDIF
+
+C       Now we can finally activate the cache
+		CALL SwitchOnCache_cll((unique_id-1)*4+1)
+		COLLIER_CACHE_ACTIVE = 1
+
+	  END

=== modified file 'madgraph/iolibs/template_files/loop_optimized/CT_interface.inc'
--- madgraph/iolibs/template_files/loop_optimized/CT_interface.inc	2016-04-28 23:52:58 +0000
+++ madgraph/iolibs/template_files/loop_optimized/CT_interface.inc	2016-07-02 20:37:49 +0000
@@ -34,8 +34,8 @@
 C
 	  %(complex_dp_format)s R1, ACC
 	  INTEGER I, J, K
-	  LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
-	  COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT
+	  LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT,COLLIERINIT
+	  COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT,COLLIERINIT
 C
 C EXTERNAL FUNCTIONS
 C
@@ -351,8 +351,8 @@
 	  %(complex_dp_format)s R1, ACC
 	  INTEGER I, J, K
       %(real_dp_format)s PDEN_DUMMY(0:3,NLOOPLINE-1)
-	  LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
-	  COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT
+	  LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT,COLLIERINIT
+	  COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT,COLLIERINIT
 C
 C EXTERNAL FUNCTIONS
 C
@@ -649,8 +649,8 @@
 C
       include 'coupl.inc'
 
-	  LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
-	  COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT
+	  LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT,COLLIERINIT
+	  COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT,COLLIERINIT
 
 	  %(real_dp_format)s LSCALE
 	  INTEGER CTMODE
@@ -875,8 +875,8 @@
 C
       include 'coupl.inc'
 
-	  LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
-	  COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT
+	  LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT,COLLIERINIT
+	  COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT,COLLIERINIT
 
 	  %(real_dp_format)s LSCALE
 	  INTEGER CTMODE

=== modified file 'madgraph/iolibs/template_files/loop_optimized/GOLEM_interface.inc'
--- madgraph/iolibs/template_files/loop_optimized/GOLEM_interface.inc	2016-03-10 00:23:21 +0000
+++ madgraph/iolibs/template_files/loop_optimized/GOLEM_interface.inc	2016-07-02 20:37:49 +0000
@@ -103,8 +103,8 @@
 	  common/%(proc_prefix)sLOOP/ID,R
 ## }
 
-      LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
-      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT
+      LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT,COLLIERINIT
+      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT,COLLIERINIT
 
       INTEGER NLOOPGROUPS
       PARAMETER (NLOOPGROUPS=%(nloop_groups)d)

=== modified file 'madgraph/iolibs/template_files/loop_optimized/TIR_interface.inc'
--- madgraph/iolibs/template_files/loop_optimized/TIR_interface.inc	2016-03-10 00:23:21 +0000
+++ madgraph/iolibs/template_files/loop_optimized/TIR_interface.inc	2016-07-02 20:37:49 +0000
@@ -43,14 +43,14 @@
 C
 	  INTEGER I, J, K
 	  INTEGER NLOOPCOEFS
-	  LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
-	  COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT
+	  LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT,COLLIERINIT
+	  COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT,COLLIERINIT
 
 c This variable will be used to detect changes in the TIR library used so as to force the reset of the TIR filter.
 	  INTEGER LAST_LIB_USED
 	  DATA LAST_LIB_USED/-1/
 
-       %(complex_dp_format)s TIRCOEFS(0:LOOPMAXCOEFS-1,3)
+       %(complex_dp_format)s TIRCOEFS(0:LOOPMAXCOEFS-1,3),TIRCOEFSERRORS(0:LOOPMAXCOEFS-1,3)
        %(complex_dp_format)s PJCOEFS(0:LOOPMAXCOEFS-1,3)
 ## if(TIRCaching){
        integer TIR_CACHE_INDEX 
@@ -72,9 +72,13 @@
 	  %(real_dp_format)s LSCALE
       common/%(proc_prefix)sCT/LSCALE,CTMODE
 
-C     The variables below are just for monitoring purposes. 
+## if(not AmplitudeReduction) {
       INTEGER ID,SQSOINDEX,R
-      COMMON/%(proc_prefix)sLOOP/ID,SQSOINDEX,R
+	  common/%(proc_prefix)sLOOP/ID,SQSOINDEX,R
+## }else{
+      INTEGER ID,R
+	  common/%(proc_prefix)sLOOP/ID,R
+## }
 
 c The argument ILIB is the TIR library to be used for that specific library.
       INTEGER LIBINDEX
@@ -178,9 +182,9 @@
       enddo
 
       do I=0,3
-      	 do J=1,(NLOOPLINE-1)
-	    PDEN(I,J)=PCT(I,J)
-	 enddo
+        do J=1,(NLOOPLINE-1)
+	      PDEN(I,J)=PCT(I,J)
+	    enddo
       enddo
 C NUMBER OF INDEPEDENT LOOPCOEFS FOR RANK=RANK
       NLOOPCOEFS=0
@@ -220,6 +224,23 @@
       WRITE(*,*) 'ERROR:: IREGI is not interfaced.'
       STOP
 ## }
+      CASE(7)
+C     COLLIER
+## if(collier_available){
+      CALL %(proc_prefix)sCOLLIERLOOP(CTMODE,NLOOPLINE,RANK,PL,PDEN,M2L,TIRCOEFS,TIRCOEFSERRORS)
+C     Shift the TIR coefficients by the corresponding COLLIER error if in CTMODE 2.
+	  if (COLLIERUseInternalStabilityTest.and.CTMODE.eq.2) THEN
+C       We add here the numerical inaccuracies linearly to be conservative 
+	    DO I=1,3
+      	  DO J=0,NLOOPCOEFS-1
+		    TIRCOEFS(J,I)=TIRCOEFS(J,I)+TIRCOEFSERRORS(J,I)
+		  ENDDO
+		ENDDO
+	  ENDIF
+## } else {
+      WRITE(*,*) 'ERROR:: COLLIER is not interfaced.'
+      STOP
+## }
       END SELECT
 ## if(TIRCaching){
 C The zero index is dummy and means no caching at all since TIR_CACHE_SIZE=0.
@@ -262,6 +283,8 @@
 C       WRITE(*,*) 'PJFry: Loop ID',ID,' =',RES(1),RES(2),RES(3)
 C     ELSEIF(MLReductionLib(I_LIB).EQ.3) THEN
 C       WRITE(*,*) 'Iregi: Loop ID',ID,' =',RES(1),RES(2),RES(3)
+C     ELSEIF(MLReductionLib(I_LIB).EQ.7) THEN
+C       WRITE(*,*) 'COLLIER: Loop ID',ID,' =',RES(1),RES(2),RES(3)
 C     ENDIF
       END
 
@@ -385,8 +408,8 @@
 C GLOBAL VARIABLES 
 C
 	  include 'MadLoopParams.inc'
-	  LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
-	  COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT
+	  LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT,COLLIERINIT
+	  COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT,NINJAINIT,COLLIERINIT
 	  
 C ----------
 C BEGIN CODE
@@ -428,7 +451,7 @@
 C CONSTANTS
 C
             INTEGER NLOOPLIB
-			PARAMETER (NLOOPLIB=4)
+			PARAMETER (NLOOPLIB=7)
 			INTEGER QP_NLOOPLIB
 ## if(ninja_supports_quad_prec) {
       		PARAMETER (QP_NLOOPLIB=2)
@@ -479,8 +502,8 @@
 		      ENDIF
 		   ENDDO
 		   IF(J_LIB.EQ.0)THEN
-                      SELECT_LIBINDEX=SELECT_LIBINDEX+1		
-                      IF(SELECT_LIBINDEX.GT.NLOOPLIB.OR.MLReductionLib(SELECT_LIBINDEX).EQ.0)SELECT_LIBINDEX=1
+              SELECT_LIBINDEX=SELECT_LIBINDEX+1		
+              IF(SELECT_LIBINDEX.GT.NLOOPLIB.OR.MLReductionLib(SELECT_LIBINDEX).EQ.0)SELECT_LIBINDEX=1
 		   ENDIF
 		ENDDO
 		I=J_LIB

=== modified file 'madgraph/iolibs/template_files/loop_optimized/loop_matrix_standalone.inc'
--- madgraph/iolibs/template_files/loop_optimized/loop_matrix_standalone.inc	2016-04-05 23:57:17 +0000
+++ madgraph/iolibs/template_files/loop_optimized/loop_matrix_standalone.inc	2016-07-02 20:37:49 +0000
@@ -97,9 +97,9 @@
 	  INTEGER NSQUAREDSOP1
 	  PARAMETER (NSQUAREDSOP1=NSQUAREDSO+1)
 C     The total number of loop reduction libraries
-C     At present, there are only CutTools,PJFry++,IREGI,Golem95,Samurai and Ninja
+C     At present, there are only CutTools,PJFry++,IREGI,Golem95,Samurai, Ninja and COLLIER
       INTEGER NLOOPLIB
-      PARAMETER (NLOOPLIB=6)
+      PARAMETER (NLOOPLIB=7)
 C     Only CutTools or possibly Ninja (if installed with qp support) provide QP
       INTEGER QP_NLOOPLIB
 ## if(ninja_supports_quad_prec) {
@@ -235,6 +235,7 @@
 C GLOBAL VARIABLES
 C 
       include 'process_info.inc'
+      include 'unique_id.inc'
       include 'coupl.inc'
 	  include 'mp_coupl.inc'
 	  include 'MadLoopParams.inc'
@@ -267,20 +268,20 @@
 	  DATA MP_DONE/.FALSE./
 	  common/%(proc_prefix)sMP_DONE/MP_DONE
 C     A FLAG TO DENOTE WHETHER THE CORRESPONDING LOOPLIBS ARE AVAILABLE OR NOT
-      LOGICAL LOOPLIBS_AVAILABLE(6)
+      LOGICAL LOOPLIBS_AVAILABLE(NLOOPLIB)
       DATA LOOPLIBS_AVAILABLE/%(data_looplibs_av)s/
 	  common/%(proc_prefix)sLOOPLIBS_AV/ LOOPLIBS_AVAILABLE
 C     A FLAG TO DENOTE WHETHER THE CORRESPONDING DIRECTION TESTS AVAILABLE OR NOT IN THE LOOPLIBS
-      LOGICAL LOOPLIBS_DIRECTEST(6)
-	  DATA LOOPLIBS_DIRECTEST /.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE./
+      LOGICAL LOOPLIBS_DIRECTEST(NLOOPLIB)
+	  DATA LOOPLIBS_DIRECTEST /.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE./
 C     Specifying for which reduction tool quadruple precision is available.
 C     The index 0 is dummy and simply means that the corresponding loop_library is not available
 C     in which case neither is its quadruple precision version.
-      LOGICAL LOOPLIBS_QPAVAILABLE(0:6)
+      LOGICAL LOOPLIBS_QPAVAILABLE(0:7)
 ## if(ninja_supports_quad_prec) {
-	  DATA LOOPLIBS_QPAVAILABLE /.FALSE.,.TRUE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.TRUE./
+	  DATA LOOPLIBS_QPAVAILABLE /.FALSE.,.TRUE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.TRUE.,.FALSE./
 ## } else {
-	  DATA LOOPLIBS_QPAVAILABLE /.FALSE.,.TRUE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE./
+	  DATA LOOPLIBS_QPAVAILABLE /.FALSE.,.TRUE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE./
 ## }
 C     PS CAN POSSIBILY BE PASSED THROUGH IMPROVE_PS BUT IS NOT MODIFIED FOR THE PURPOSE OF THE STABILITY TEST
 C     EVEN THOUGH THEY ARE PUT IN COMMON BLOCK, FOR NOW THEY ARE NOT USED ANYWHERE ELSE
@@ -329,7 +330,7 @@
 	  DATA CT_REQ_SO_DONE/.FALSE./
 	  DATA MP_CT_REQ_SO_DONE/.FALSE./	  
 	  DATA LOOP_REQ_SO_DONE/.FALSE./
-	  DATA MP_LOOP_REQ_SO_DONE/.FALSE./	  
+	  DATA MP_LOOP_REQ_SO_DONE/.FALSE./
 	  DATA CTCALL_REQ_SO_DONE/.FALSE./
 	  DATA FILTER_SO/.FALSE./
 	  common/%(proc_prefix)sSO_REQS/UVCT_REQ_SO_DONE,MP_UVCT_REQ_SO_DONE,CT_REQ_SO_DONE,MP_CT_REQ_SO_DONE,LOOP_REQ_SO_DONE,MP_LOOP_REQ_SO_DONE,CTCALL_REQ_SO_DONE,FILTER_SO
@@ -345,11 +346,10 @@
 	  integer I_LIB
 	  DATA I_LIB/1/
 	  common/%(proc_prefix)sI_LIB/I_LIB
-C UTIL NOW, ONLY CUTTOOLS PROVIDE QP
-       	  LOGICAL QP_TOOLS_AVAILABLE
-       	  DATA QP_TOOLS_AVAILABLE/.FALSE./
-       	  INTEGER INDEX_QP_TOOLS(QP_NLOOPLIB+1)
-       	  common/%(proc_prefix)sLOOP_TOOLS/QP_TOOLS_AVAILABLE,INDEX_QP_TOOLS
+      LOGICAL QP_TOOLS_AVAILABLE
+      DATA QP_TOOLS_AVAILABLE/.FALSE./
+      INTEGER INDEX_QP_TOOLS(QP_NLOOPLIB+1)
+      common/%(proc_prefix)sLOOP_TOOLS/QP_TOOLS_AVAILABLE,INDEX_QP_TOOLS
 
 ## if(not LoopInduced) {
 	  %(complex_dp_format)s AMP(NBORNAMPS)
@@ -420,9 +420,16 @@
 	  DATA MP_DONE_ONCE/.FALSE./
 	  common/%(proc_prefix)sMP_DONE_ONCE/MP_DONE_ONCE
 
-	  character(512) MLPath
+                                               	  character(512) MLPath
       common/MLPATH/MLPath
 
+C     This is just so that if the user disabled the computation of poles by COLLIER
+C     using the MadLoop subroutine, we don't overwrite his choice when reading the parameters
+      LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, COLLIER_IR_POLE_COMPUTATION_CHOICE
+	  DATA  FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION/.FALSE.,.FALSE./
+	  COMMON/%(proc_prefix)sCOLLIERPOLESFORCEDCHOICE/FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION,COLLIER_UV_POLE_COMPUTATION_CHOICE,COLLIER_IR_POLE_COMPUTATION_CHOICE
+
 C     This variable controls the general initialization which is *common* between all MadLoop SubProcesses.
 C     For example setting the MadLoopPath or reading the ML runtime parameters.
 	  LOGICAL ML_INIT
@@ -440,6 +447,7 @@
       LOGICAL FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION
 	  data FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION/.FALSE.,.FALSE./      
 	  COMMON/%(proc_prefix)sFPE_IN_REDUCTION/FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION
+
 C ----------
 C BEGIN CODE
 C ----------
@@ -451,6 +459,12 @@
   CALL setMadLoopPath(TMP)
   CALL JOINPATH(MLPATH,PARAMFNAME,PARAMFN)
   CALL MADLOOPPARAMREADER(PARAMFN,.TRUE.)
+  IF (FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION) THEN
+    COLLIERComputeUVpoles = COLLIER_UV_POLE_COMPUTATION_CHOICE
+  ENDIF
+  IF (FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION) THEN
+    COLLIERComputeIRpoles = COLLIER_IR_POLE_COMPUTATION_CHOICE
+  ENDIF
   IF (FORBID_HEL_DOUBLECHECK) THEN
     DoubleCheckHelicityFilter = .False. 
   ENDIF 
@@ -801,6 +815,11 @@
   CALL %(proc_prefix)sCLEAR_CACHES()
 ENDIF
 
+## if(collier_available){
+C Now make sure to turn on the global COLLIER cache if applicable
+CALL %(proc_prefix)sSET_COLLIER_GLOBAL_CACHE(.TRUE.)
+## }
+
 IF (ImprovePSPoint.ge.0) THEN
 C Make the input PS more precise (exact onshell and energy-momentum conservation)
   CALL %(proc_prefix)sIMPROVE_PS_POINT_PRECISION(PS)
@@ -1812,14 +1831,27 @@
 IF (AUTOMATIC_CACHE_CLEARING) THEN
   CALL %(proc_prefix)sCLEAR_CACHES()
 ENDIF
+
+## if(collier_available){
+C Now make sure to turn off the global COLLIER cache if applicable
+CALL %(proc_prefix)sSET_COLLIER_GLOBAL_CACHE(.FALSE.)
+## }
+
       END
 
       SUBROUTINE %(proc_prefix)sCLEAR_CACHES()
-C Clears all the caches used at some point in MadLoop
+C
+C     This routine can be called directly from the user if
+C     AUTOMATIC_CACHE_CLEARING is set to False. It must then be called after
+C     ech event
+C
         CALL %(proc_prefix)sCLEAR_TIR_CACHE()
 ## if(ninja_available){
         call ninja_clear_integral_cache()
 ## }
+## if(collier_available){
+        call %(proc_prefix)sCLEAR_COLLIER_CACHE()
+## }
       END
 
 C --=========================================--
@@ -1952,6 +1984,13 @@
 	  real*8 accuracies(3)
 	  real*8 list(maxstabilitylength)
 
+C
+C GLOBAL VARIABLES
+C
+	  integer I_LIB
+	  common/%(proc_prefix)sI_LIB/I_LIB
+	  include 'MadLoopParams.inc'
+
 C ----------
 C BEGIN CODE
 C ----------
@@ -1998,6 +2037,13 @@
 		if (avg.ne.0.0d0) then
           acc(k)  = acc(k) / ( ABS(AVG) / 3.0d0)
 		endif
+        
+C       When using COLLIER with the internal stability test, the first evaluation is typically more reliable so we do not want to use the average but rather the first evaluation.
+		IF (MLREDUCTIONLIB(I_LIB).eq.7.and.COLLIERUseInternalStabilityTest) THEN
+		  DO I=1,3
+		    estimate(i,k) = fulllist(i,k,1)
+		  ENDDO
+		ENDIF
 
 	  enddo
 
@@ -2326,6 +2372,46 @@
 C   Definition of additional access routines
 C --=========================================--
 
+	  SUBROUTINE %(proc_prefix)sCOLLIER_COMPUTE_UV_POLES(ONOFF)
+C
+C This function can be called by the MadLoop user so as to chose to have COLLIER
+C compute the UV pole or not (it costs more time).
+C
+      LOGICAL ONOFF
+
+	  include 'MadLoopParams.inc'
+
+      LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, COLLIER_IR_POLE_COMPUTATION_CHOICE
+	  COMMON/%(proc_prefix)sCOLLIERPOLESFORCEDCHOICE/FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION,COLLIER_UV_POLE_COMPUTATION_CHOICE,COLLIER_IR_POLE_COMPUTATION_CHOICE
+
+	  COLLIERComputeUVpoles                        = ONOFF
+C     This is just so that if we read the param again, we don't overwrite the choice made here
+	  FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION = .TRUE.
+      COLLIER_UV_POLE_COMPUTATION_CHOICE           = ONOFF
+
+	  END SUBROUTINE
+
+	  SUBROUTINE %(proc_prefix)sCOLLIER_COMPUTE_IR_POLES(ONOFF)
+C
+C This function can be called by the MadLoop user so as to chose to have COLLIER
+C compute the IR pole or not (it costs more time).
+C
+      LOGICAL ONOFF
+
+	  include 'MadLoopParams.inc'
+
+      LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, COLLIER_IR_POLE_COMPUTATION_CHOICE
+	  COMMON/%(proc_prefix)sCOLLIERPOLESFORCEDCHOICE/FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION,COLLIER_UV_POLE_COMPUTATION_CHOICE,COLLIER_IR_POLE_COMPUTATION_CHOICE
+
+	  COLLIERComputeIRpoles         = ONOFF
+C     This is just so that if we read the param again, we don't overwrite the choice made here
+	  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION = .TRUE.
+      COLLIER_IR_POLE_COMPUTATION_CHOICE           = ONOFF
+
+	  END SUBROUTINE
+
 	  SUBROUTINE %(proc_prefix)sFORCE_STABILITY_CHECK(ONOFF)
 C
 C This function can be called by the MadLoop user so as to always have stability

=== modified file 'madgraph/iolibs/template_files/matrix_loop_induced_madevent.inc'
--- madgraph/iolibs/template_files/matrix_loop_induced_madevent.inc	2016-02-18 03:09:50 +0000
+++ madgraph/iolibs/template_files/matrix_loop_induced_madevent.inc	2016-07-02 20:37:49 +0000
@@ -125,6 +125,9 @@
 
 	IF(.NOT.INIT_MADLOOP) THEN
 	  INIT_MADLOOP = .TRUE.
+C We don't use the poles, so let's not compute them with COLLIER
+      CALL %(ml_prefix)sCOLLIER_COMPUTE_UV_POLES(.False.)
+      CALL %(ml_prefix)sCOLLIER_COMPUTE_IR_POLES(.False.)
 C Unless this is the first iteration, make sure to never double check the helicity filter again
       IF (read_grid_file) THEN
         CALL SET_FORBID_HEL_DOUBLECHECK(.True.)
@@ -392,7 +395,7 @@
 	  return
 	endif
 
-	IF (((OneEps+TwoEps)/finite).gt.(1000d0*PREC).AND.(((OneEps+TwoEps)/finite).gt.1.0d-5)) THEN
+	IF (U.NE.7.AND.((OneEps+TwoEps)/finite).gt.(1000d0*PREC).AND.(((OneEps+TwoEps)/finite).gt.1.0d-5)) THEN
 	  WARNING_COUNTERS(2) = WARNING_COUNTERS(2) + 1
 	  IF (WARNING_COUNTERS(2).le.10) THEN
 	    WRITE(*,*) "WARNING : The residue of the single and double pole of the loop matrix element being integrated does not seem to vanish."

=== modified file 'madgraph/iolibs/template_files/matrix_loop_induced_madevent_group.inc'
--- madgraph/iolibs/template_files/matrix_loop_induced_madevent_group.inc	2016-02-18 03:09:50 +0000
+++ madgraph/iolibs/template_files/matrix_loop_induced_madevent_group.inc	2016-07-02 20:37:49 +0000
@@ -133,6 +133,9 @@
 
 	IF(.NOT.INIT_MADLOOP) THEN
 	  INIT_MADLOOP = .TRUE.
+C We don't use the poles, so let's not compute them with COLLIER
+      CALL %(ml_prefix)sCOLLIER_COMPUTE_UV_POLES(.False.)
+      CALL %(ml_prefix)sCOLLIER_COMPUTE_IR_POLES(.False.)
 C Unless this is the first iteration, make sure to never double check the helicity filter again
       IF (read_grid_file) THEN 
         CALL SET_FORBID_HEL_DOUBLECHECK(.True.)
@@ -405,7 +408,7 @@
 	  return
 	endif
 
-	IF (((OneEps+TwoEps)/finite).gt.(1000d0*PREC).AND.(((OneEps+TwoEps)/finite).gt.1.0d-5)) THEN
+	IF (U.NE.7.and.((OneEps+TwoEps)/finite).gt.(1000d0*PREC).AND.(((OneEps+TwoEps)/finite).gt.1.0d-5)) THEN
 	  WARNING_COUNTERS(2) = WARNING_COUNTERS(2) + 1
 	  IF (WARNING_COUNTERS(2).le.10) THEN
 	    WRITE(*,*) "WARNING : The residue of the single and double pole of the loop matrix element being integrated does not seem to vanish."

=== modified file 'madgraph/loop/loop_exporters.py'
--- madgraph/loop/loop_exporters.py	2016-05-11 21:44:52 +0000
+++ madgraph/loop/loop_exporters.py	2016-07-02 20:37:49 +0000
@@ -90,7 +90,8 @@
 
     include_names    = {'ninja' : 'mninja.mod',
                         'golem' : 'generic_function_1p.mod',
-                        'samurai':'msamurai.mod'}
+                        'samurai':'msamurai.mod',
+                        'collier': 'collier.mod'}
 
     def __init__(self, mgme_dir="", dir_path = "", opt=None):
         """Initiate the LoopExporterFortran with directory information on where
@@ -263,17 +264,18 @@
         calls = self.write_loop_makefile_definitions(
                         writers.MakefileWriter(filePath),link_tir_libs,tir_libs)
 
-            
         # We need minimal editing of MadLoopCommons.f
+        # For the optimized output, this file will be overwritten once the 
+        # availability of COLLIER has been determined.
         MadLoopCommon = open(os.path.join(self.loop_dir,'StandAlone', 
                                     "SubProcesses","MadLoopCommons.inc")).read()
         writer = writers.FortranWriter(os.path.join(self.dir_path, 
                                              "SubProcesses","MadLoopCommons.f"))
         writer.writelines(MadLoopCommon%{
-                                   'print_banner_commands':self.MadLoop_banner})
+          'print_banner_commands':self.MadLoop_banner}, context={
+                'collier_available':False})
         writer.close()
         
-        
         # Copy the whole MadLoop5_resources directory (empty at this stage)
         if not os.path.exists(pjoin(self.dir_path,'SubProcesses',
                                                         'MadLoop5_resources')):
@@ -678,7 +680,7 @@
     # generate_subprocess_directory_v4
     #===========================================================================
     def generate_loop_subprocess(self, matrix_element, fortran_model,
-                          group_number = None, proc_id = None, config_map=None):
+          group_number = None, proc_id = None, config_map=None, unique_id=None):
         """Generate the Pxxxxx directory for a loop subprocess in MG4 standalone,
         including the necessary loop_matrix.f, born_matrix.f and include files.
         Notice that this is too different from generate_subprocess_directory_v4
@@ -705,6 +707,14 @@
 
         logger.info('Creating files in directory %s' % dirpath)
 
+        if unique_id is None:
+            raise MadGraph5Error, 'A unique id must be provided to the function'+\
+                     'generate_loop_subprocess of LoopProcessExporterFortranSA.'
+        # Create an include with the unique consecutive ID assigned
+        open('unique_id.inc','w').write(
+"""      integer UNIQUE_ID
+      parameter(UNIQUE_ID=%d)"""%unique_id)
+        
         # Extract number of external particles
         (nexternal, ninitial) = matrix_element.get_nexternal_ninitial()
 
@@ -778,7 +788,7 @@
                      'cts_mprec.h', 'cts_mpc.h', 'mp_coupl.inc', 
                      'mp_coupl_same_name.inc',
                      'MadLoopParamReader.f','MadLoopCommons.f',
-                     'MadLoopParams.inc']
+                     'MadLoopParams.inc','global_specs.inc']
         
         for file in linkfiles:
             ln('../%s' % file)
@@ -1014,11 +1024,13 @@
 
         writer.writelines(proc_include)
                                 
-    def generate_subprocess_directory_v4(self, matrix_element, fortran_model):
+    def generate_subprocess_directory_v4(self, matrix_element, fortran_model,
+                                                                     unique_id):
         """ To overload the default name for this function such that the correct
         function is used when called from the command interface """
         
-        return self.generate_loop_subprocess(matrix_element,fortran_model)
+        return self.generate_loop_subprocess(matrix_element,fortran_model,
+                                                            unique_id=unique_id)
 
     def write_check_sa(self, writer, matrix_element):
         """Writes out the steering code check_sa. In the optimized output mode,
@@ -1576,7 +1588,7 @@
     # make sure that ninja appears first in the list of -L because 
     # it is the tool for which the user is most susceptible of 
     # using a standalone verison independent of gosam_contrib
-    all_tir=['pjfry','iregi','ninja','golem','samurai']
+    all_tir=['pjfry','iregi','ninja','golem','samurai','collier']
     
     def __init__(self, mgme_dir="", dir_path = "", opt=None):
         """Initiate the LoopProcessOptimizedExporterFortranSA with directory 
@@ -1588,7 +1600,7 @@
 
         # TIR available ones
         self.tir_available_dict={'pjfry':True,'iregi':True,'golem':True,
-                                 'samurai':True,'ninja':True}
+                                 'samurai':True,'ninja':True,'collier':True}
 
         for tir in self.all_tir:
             tir_dir="%s_dir"%tir
@@ -1627,7 +1639,7 @@
         for tir in self.all_tir:
             context['%s_available'%tir]=self.tir_available_dict[tir]
             # safety check
-            if tir not in ['golem','pjfry','iregi','samurai','ninja']:
+            if tir not in ['golem','pjfry','iregi','samurai','ninja','collier']:
                 raise MadGraph5Error,"%s was not a TIR currently interfaced."%tir_name
 
         return context
@@ -1648,16 +1660,30 @@
             tir_name=tir
             libpath = self.link_TIR(os.path.join(self.dir_path, 'lib'),
                                               libpath,libname,tir_name=tir_name)
-            setattr(self,tir_dir,libpath)
             if libpath != "":
-                if tir in ['ninja','pjfry','golem','samurai']:
+                if tir in ['ninja','pjfry','golem','samurai','collier']:
                     # It is cleaner to use the original location of the libraries
                     link_tir_libs.append('-L%s/ -l%s'%(libpath,tir))
                     tir_libs.append('%s/lib%s.$(libext)'%(libpath,tir))
-                    if tir in ['ninja','golem', 'samurai']:
+                    # For Ninja, we must also link against OneLoop.
+                    if tir in ['ninja']:
+                        if not any(os.path.isfile(pjoin(libpath,'libavh_olo.%s'%ext)) 
+                                              for ext in ['a','dylib','so']):
+                            raise MadGraph5Error(
+"The OneLOop library 'libavh_olo.(a|dylib|so)' could no be found in path '%s'. Please place a symlink to it there."%libpath)
+                        link_tir_libs.append('-L%s/ -l%s'%(libpath,'avh_olo'))
+                        tir_libs.append('%s/lib%s.$(libext)'%(libpath,'avh_olo'))
+                    if tir in ['ninja','golem', 'samurai','collier']:
                         trgt_path = pjoin(os.path.dirname(libpath),'include')
-                        to_include = misc.find_includes_path(trgt_path,
+                        if os.path.isdir(trgt_path):
+                            to_include = misc.find_includes_path(trgt_path,
                                                         self.include_names[tir])
+                        else:
+                            to_include = None
+                        # Special possible location for collier
+                        if to_include is None and tir=='collier':
+                            to_include = misc.find_includes_path(
+                               pjoin(libpath,'modules'),self.include_names[tir])
                         if to_include is None:
                             logger.error(
 'Could not find the include directory for %s, looking in %s.\n' % (tir, str(trgt_path))+
@@ -1669,7 +1695,7 @@
                         # (such as what is done with the Sherpa interface), we
                         # place here an easy handle on the golem includes
                         name_map = {'golem':'golem95','samurai':'samurai',
-                                    'ninja':'ninja'}
+                                    'ninja':'ninja','collier':'collier'}
                         ln(to_include, starting_dir=pjoin(self.dir_path,'lib'),
                                    name='%s_include'%name_map[tir],abspath=True)
                         ln(libpath, starting_dir=pjoin(self.dir_path,'lib'),
@@ -1687,6 +1713,17 @@
                         writers.MakefileWriter(MadLoop_makefile_definitions),
                                 link_tir_libs,tir_libs, tir_include=tir_include)
 
+        # Finally overwrite MadLoopCommons.f now that we know the availibility of
+        # COLLIER.
+        MadLoopCommon = open(os.path.join(self.loop_dir,'StandAlone', 
+                                    "SubProcesses","MadLoopCommons.inc")).read()
+        writer = writers.FortranWriter(os.path.join(self.dir_path, 
+                                             "SubProcesses","MadLoopCommons.f"))
+        writer.writelines(MadLoopCommon%{
+          'print_banner_commands':self.MadLoop_banner}, context={
+                'collier_available':self.tir_available_dict['collier']})
+        writer.close()
+
     def link_files_from_Subprocesses(self,proc_name):
         """ Does the same as the mother routine except that it also links
         coef_specs.inc in the HELAS folder."""
@@ -1704,7 +1741,7 @@
         """Link the TIR source directory inside the target path given
         in argument"""
         
-        if tir_name in ['pjfry','golem','samurai','ninja']:
+        if tir_name in ['pjfry','golem','samurai','ninja','collier']:
             # not self-contained libraries
             if (not isinstance(libpath,str)) or (not os.path.exists(libpath)) \
             or (not os.path.isfile(pjoin(libpath,libname))):
@@ -1742,7 +1779,7 @@
                 return ""
        
         if self.dependencies=='internal':
-            if tir_name in ['pjfry','golem','samurai','ninja']:
+            if tir_name in ['pjfry','golem','samurai','ninja','collier']:
                 self.tir_available_dict[tir_name]=False
                 logger.info("When using the 'output_dependencies=internal' "+\
 " MG5_aMC option, the (optional) reduction library %s cannot be employed because"%tir_name+\
@@ -1793,7 +1830,7 @@
                     self.tir_available_dict[tir_name]=False
                     return ""
             # We link the tools below directly to the lib directory of the output 
-            if not tir_name in ['pjfry','golem','samurai','ninja']:
+            if not tir_name in ['pjfry','golem','samurai','ninja','collier']:
                 ln(os.path.join(libpath,libname),targetPath,abspath=True)
 
         elif self.dependencies=='environment_paths':
@@ -1804,7 +1841,7 @@
                 logger.info('MG5_aMC is using %s installation found at %s.'%\
                                                           (tir_name,newlibpath)) 
                 # We link the tools below directly to directly where the library is detected
-                if not tir_name in ['pjfry','golem','samurai','ninja']:
+                if not tir_name in ['pjfry','golem','samurai','ninja','collier']:
                     ln(newlibpath,targetPath,abspath=True)
                 self.tir_available_dict[tir_name]=True
                 return os.path.dirname(newlibpath)
@@ -1911,6 +1948,11 @@
             self.write_GOLEM_interface(writers.FortranWriter(filename),
                                        matrix_element)
 
+        if 'collier' in self.tir_available_dict and self.tir_available_dict['collier']:
+            filename = 'COLLIER_interface.f'
+            self.write_COLLIER_interface(writers.FortranWriter(filename),
+                                       matrix_element)
+
         filename = 'loop_num.f'
         self.write_loop_num(writers.FortranWriter(filename),\
                                            matrix_element,OptimizedFortranModel)
@@ -2042,6 +2084,40 @@
         else:
             return file
 
+    def write_COLLIER_interface(self, writer, matrix_element):
+        """ Create the file COLLIER_interface.f"""
+
+        # First write GOLEM_interface which interfaces MG5 with TIR.
+        replace_dict=copy.copy(matrix_element.rep_dict)
+            
+        file = open(os.path.join(self.template_dir,'COLLIER_interface.inc')).read()
+ 
+        FPR = q_polynomial.FortranPolynomialRoutines(replace_dict['maxrank'],\
+                                                    coef_format=replace_dict['complex_dp_format'],\
+                                                    sub_prefix=replace_dict['proc_prefix'])
+        map_definition = []
+        collier_map = FPR.get_COLLIER_mapping()
+        
+        chunk_size = 10
+        for map_name, indices_list in \
+                           [('COEFMAP_ZERO',[c[0] for c in collier_map]),
+                            ('COEFMAP_ONE',[c[1] for c in collier_map]),
+                            ('COEFMAP_TWO',[c[2] for c in collier_map]),
+                            ('COEFMAP_THREE',[c[3] for c in collier_map])]:
+            for k in xrange(0, len(indices_list), chunk_size):
+                map_definition.append("DATA (%s(I),I=%3r,%3r) /%s/" % \
+                (map_name,k, min(k + chunk_size, len(indices_list))-1,
+                 ','.join('%2r'%ind for ind in indices_list[k:k + chunk_size])))
+
+        replace_dict['collier_coefmap'] = '\n'.join(map_definition) 
+ 
+        file = file % replace_dict
+        
+        if writer:
+            writer.writelines(file,context=self.get_context(matrix_element))
+        else:
+            return file
+
     def write_GOLEM_interface(self, writer, matrix_element):
         """ Create the file GOLEM_interface.f which does NOT contain the subroutine
          defining the loop HELAS-like calls along with the general interfacing 
@@ -2371,6 +2447,27 @@
 
         writer.writelines(file,context=self.get_context(matrix_element))
     
+    def write_global_specs(self, matrix_element_list):
+        """ From the list of matrix element, or the single matrix element, derive
+        the global quantities to write in global_coef_specs.inc"""
+        
+        if isinstance(matrix_element_list, list):
+            me_list = matrix_element_list
+        else:
+            me_list = [matrix_element_list]    
+
+        open(pjoin(self.dir_path,'SubProcesses','global_specs.inc'),'w').write(
+"""      integer MAXNEXTERNAL
+      parameter(MAXNEXTERNAL=%d)
+      integer OVERALLMAXRANK
+      parameter(OVERALLMAXRANK=%d)
+      integer NPROCS
+      parameter(NPROCS=%d)"""%(
+         max(me.get_nexternal_ninitial()[0] for me in me_list),
+         max(me.get_max_loop_rank() for me in me_list),
+         len(me_list)))    
+
+    
     def fix_coef_specs(self, overall_max_lwf_spin, overall_max_loop_vert_rank):
         """ If processes with different maximum loop wavefunction size or
         different maximum loop vertex rank have to be output together, then
@@ -2572,7 +2669,7 @@
         looplibs_av=['.TRUE.']
         # one should be careful about the order in the following as it must match
         # the ordering in MadLoopParamsCard.
-        for tir_lib in ['pjfry','iregi','golem','samurai','ninja']:
+        for tir_lib in ['pjfry','iregi','golem','samurai','ninja','collier']:
             looplibs_av.append('.TRUE.' if tir_lib in self.all_tir and \
                                 self.tir_available_dict[tir_lib] else '.FALSE.')
         replace_dict['data_looplibs_av']=','.join(looplibs_av)
@@ -2944,11 +3041,15 @@
         LoopInducedExporterME.finalize_v4_directory(self,*args,**opts)
         
     def generate_subprocess_directory_v4(self, subproc_group,
-                                                    fortran_model,group_number):
+                            fortran_model,group_number, unique_id=None):
         """Generate the Pn directory for a subprocess group in MadEvent,
         including the necessary matrix_N.f files, configs.inc and various
         other helper files"""
-            
+        
+        if unique_id is None:
+            raise MadGraph5Error, 'A unique_id must be provided to '+\
+                'generate_subprocess_directory_v4 in LoopInducedExporterMEGroup.'
+
         # Generate the MadLoop files
         calls = 0
         matrix_elements = subproc_group.get('matrix_elements')
@@ -2956,7 +3057,8 @@
             calls += self.generate_loop_subprocess(matrix_element,fortran_model,
           group_number = group_number, proc_id = str(ime+1),
 #          group_number = str(subproc_group.get('number')), proc_id = str(ime+1),
-          config_map = subproc_group.get('diagram_maps')[ime])
+          config_map = subproc_group.get('diagram_maps')[ime],
+          unique_id=unique_id+ime)
         
         # Then generate the MadEvent files
         export_v4.ProcessExporterFortranMEGroup.generate_subprocess_directory_v4(
@@ -3107,14 +3209,15 @@
         # care of MadLoop virtuals initialization
         LoopInducedExporterME.finalize_v4_directory(self,*args,**opts)
 
-    def generate_subprocess_directory_v4(self, matrix_element, fortran_model, me_number):
+    def generate_subprocess_directory_v4(self, matrix_element, fortran_model,
+                                                                     me_number):
         """Generate the Pn directory for a subprocess group in MadEvent,
         including the necessary matrix_N.f files, configs.inc and various
         other helper files"""
     
         # Then generate the MadLoop files
         calls = self.generate_loop_subprocess(matrix_element,fortran_model,                           
-                                                       group_number = me_number)
+                                  group_number = me_number, unique_id=me_number)
         
         
         # First generate the MadEvent files

=== modified file 'madgraph/madevent/sum_html.py'
--- madgraph/madevent/sum_html.py	2016-06-28 15:36:47 +0000
+++ madgraph/madevent/sum_html.py	2016-07-02 20:37:49 +0000
@@ -57,6 +57,7 @@
           'Samurai_usage'      : 0,
           'Ninja_usage'        : 0,
           'Ninja_QP_usage'     : 0,
+          'COLLIER_usage'      : 0,
           'max_precision'      : 1.0e99,
           'min_precision'      : 0.0,
           'averaged_timing'    : 0.0,
@@ -119,6 +120,7 @@
         self['Golem_usage']       = u_codes[4]
         self['Samurai_usage']     = u_codes[5]
         self['Ninja_usage']       = u_codes[6]
+        self['COLLIER_usage']     = u_codes[7]        
         self['Ninja_QP_usage']    = u_codes[8]
         self['CutTools_QP_usage'] = u_codes[9]
         t_return_code = xml_node.getElementsByTagName('t_return_code')
@@ -177,6 +179,7 @@
           ('Golem',float(self['Golem_usage'])/self['n_madloop_calls']),
           ('IREGI',float(self['IREGI_usage'])/self['n_madloop_calls']),
           ('Samurai',float(self['Samurai_usage'])/self['n_madloop_calls']),
+          ('COLLIER',float(self['COLLIER_usage'])/self['n_madloop_calls']),          
           ('Ninja_DP',float(self['Ninja_usage'])/self['n_madloop_calls']),
           ('Ninja_QP',float(self['Ninja_QP_usage'])/self['n_madloop_calls'])]
 

=== modified file 'madgraph/various/banner.py'
--- madgraph/various/banner.py	2016-06-28 12:22:14 +0000
+++ madgraph/various/banner.py	2016-07-02 20:37:49 +0000
@@ -2189,17 +2189,23 @@
     contains a parser to read it, facilities to write a new file,...
     """
     
-
-            
+    _ID_reduction_tool_map = {1:'CutTools',
+                             2:'PJFry++',
+                             3:'IREGI',
+                             4:'Golem95',
+                             5:'Samurai',
+                             6:'Ninja',
+                             7:'COLLIER'}
+    
     def default_setup(self):
         """initialize the directory to the default value"""
         
-        self.add_param("MLReductionLib", "6|1|2")
+        self.add_param("MLReductionLib", "6|7|1")
         self.add_param("IREGIMODE", 2)
         self.add_param("IREGIRECY", True)
         self.add_param("CTModeRun", -1)
         self.add_param("MLStabThres", 1e-3)
-        self.add_param("NRotations_DP", 1)
+        self.add_param("NRotations_DP", 0)
         self.add_param("NRotations_QP", 0)
         self.add_param("ImprovePSPoint", 2)
         self.add_param("CTLoopLibrary", 2)
@@ -2215,8 +2221,16 @@
         self.add_param("HelicityFilterLevel", 2)
         self.add_param("LoopInitStartOver", False)
         self.add_param("HelInitStartOver", False)
-        self.add_param("UseQPIntegrandForNinja", False)        
+        self.add_param("UseQPIntegrandForNinja", True)        
         self.add_param("UseQPIntegrandForCutTools", True)
+        self.add_param("COLLIERMode", 1)
+        self.add_param("COLLIERComputeUVpoles", True)
+        self.add_param("COLLIERComputeIRpoles", True)
+        self.add_param("COLLIERRequiredAccuracy", 1.0e-8)
+        self.add_param("COLLIERCanOutput",False)
+        self.add_param("COLLIERGlobalCache",-1)
+        self.add_param("COLLIERUseCacheForPoles",False)
+        self.add_param("COLLIERUseInternalStabilityTest",True)
 
     def read(self, finput):
         """Read the input file, this can be a path to a file, 

=== modified file 'madgraph/various/misc.py'
--- madgraph/various/misc.py	2016-05-18 14:36:57 +0000
+++ madgraph/various/misc.py	2016-07-02 20:37:49 +0000
@@ -308,7 +308,7 @@
             log(msg)
     
 
-    if dependency in ['pjfry','golem','samurai','ninja']:
+    if dependency in ['pjfry','golem','samurai','ninja','collier']:
         if cmd.options[dependency] not in ['None',None,'']:
             tell("Deactivating MG5_aMC dependency '%s'"%dependency)
             cmd.options[dependency] = None
@@ -350,6 +350,13 @@
             tell("Installing ninja...")
             cmd.do_install('ninja')
  
+    if dependency=='collier':
+        if cmd.options['collier'] in ['None',None,''] or\
+         (cmd.options['collier'] == 'auto' and which_lib('libcollier.a') is None) or\
+         which_lib(pjoin(cmd.options['collier'],'libcollier.a')) is None:
+            tell("Installing COLLIER...")
+            cmd.do_install('collier')
+
 #===============================================================================
 # find a library
 #===============================================================================

=== modified file 'madgraph/various/process_checks.py'
--- madgraph/various/process_checks.py	2016-04-15 11:14:25 +0000
+++ madgraph/various/process_checks.py	2016-07-02 20:37:49 +0000
@@ -613,7 +613,8 @@
                 self.mg_root, export_dir, MLoptions)
             FortranModel = helas_call_writers.FortranUFOHelasCallWriter(model)
             FortranExporter.copy_v4template(modelname=model.get('name'))
-            FortranExporter.generate_subprocess_directory_v4(matrix_element, FortranModel)
+            FortranExporter.generate_subprocess_directory_v4(matrix_element, FortranModel, 1)
+            FortranExporter.write_global_specs(matrix_element)
             wanted_lorentz = list(set(matrix_element.get_used_lorentz()))
             wanted_couplings = list(set([c for l in matrix_element.get_used_couplings() \
                                                                     for c in l]))
@@ -1155,7 +1156,8 @@
             FortranExporter = exporter_class(self.mg_root, export_dir, MLoptions)
             FortranModel = helas_call_writers.FortranUFOHelasCallWriter(model)
             FortranExporter.copy_v4template(modelname=model.get('name'))
-            FortranExporter.generate_subprocess_directory_v4(matrix_element, FortranModel)
+            FortranExporter.generate_subprocess_directory_v4(matrix_element, FortranModel, 1)
+            FortranExporter.write_global_specs(matrix_element)
             wanted_lorentz = list(set(matrix_element.get_used_lorentz()))
             wanted_couplings = list(set([c for l in matrix_element.get_used_couplings() \
                                                                 for c in l]))
@@ -1470,8 +1472,8 @@
             tools=list(set(tools)) # remove the duplication ones
             
         # not self-contained tir libraries
-        tool_var={'pjfry':2,'golem':4,'samurai':5,'ninja':6}
-        for tool in ['pjfry','golem','samurai','ninja']:
+        tool_var={'pjfry':2,'golem':4,'samurai':5,'ninja':6,'collier':7}
+        for tool in ['pjfry','golem','samurai','ninja','collier']:
             tool_dir='%s_dir'%tool
             if not tool_dir in self.tir_dir:
                 continue
@@ -1493,8 +1495,8 @@
         export_dir=pjoin(self.mg_root,("SAVED" if keep_folder else "")+\
                                                 temp_dir_prefix+"_%s"%proc_name)
         
-        tools_name={1:'CutTools',2:'PJFry++',3:'IREGI',4:'Golem95',5:'Samurai',
-                    6:'Ninja'}
+        tools_name=bannermod.MadLoopParam._ID_reduction_tool_map
+        
         return_dict={}
         return_dict['Stability']={}
         infos_save={'Process_output': None,
@@ -1517,7 +1519,7 @@
             # The exceptional PS points are those which stay unstable in quad prec.
             Exceptional_PS_points = []
         
-            MLoptions={}
+            MLoptions=MLOptions
             MLoptions["MLReductionLib"]=tool
             clean = (tool==tools[0]) and not nPoints==0
             if infos_IN==None or (tool_name not in infos_IN):
@@ -1606,7 +1608,7 @@
                             str(datetime.timedelta(seconds=sec_needed)),\
                             datetime.datetime.now().strftime("%d-%m-%Y %H:%M")))
             if logger.getEffectiveLevel()<logging.WARNING and \
-                (sec_needed>5 or (reusing and infos['Initialization'] == None)):
+                              (sec_needed>5 or infos['Initialization'] == None):
                 widgets = ['Stability check:', pbar.Percentage(), ' ', 
                                             pbar.Bar(),' ', pbar.ETA(), ' ']
                 progress_bar = pbar.ProgressBar(widgets=widgets, maxval=nPoints, 
@@ -2125,8 +2127,10 @@
 
     timing2 = myProfiler.time_matrix_element(matrix_element, reusing, 
                             param_card, keep_folder=keep_folder,options=options,
-                            MLOptions = MLoptions)
+                            MLOptions = MLoptions) 
     
+    timing2['reduction_tool'] = MLoptions['MLReductionLib'][0]
+
     if timing2 == None:
         return None, None
 
@@ -2172,6 +2176,22 @@
         MLoptions = {}
     else:
         MLoptions = MLOptions
+        # Make sure that the poles computation is disabled for COLLIER
+        if 'COLLIERComputeUVpoles' not in MLoptions:
+            MLoptions['COLLIERComputeUVpoles']=False
+        if 'COLLIERComputeIRpoles' not in MLoptions:
+            MLoptions['COLLIERComputeIRpoles']=False
+        # Use high required accuracy in COLLIER's requirement if not specified
+        if 'COLLIERRequiredAccuracy' not in MLoptions:
+            MLoptions['COLLIERRequiredAccuracy']=1e-13
+        # Use loop-direction switching as stability test if not specifed (more reliable)
+        if 'COLLIERUseInternalStabilityTest' not in MLoptions:
+            MLoptions['COLLIERUseInternalStabilityTest']=False
+        # Finally we *must* forbid the use of COLLIER global cache here, because it
+        # does not work with the way we call independently CTMODERun 1 and 2
+        # with the StabilityChecker.
+        MLoptions['COLLIERGlobalCache'] = 0
+
         if "MLReductionLib" not in MLOptions:
             MLoptions["MLReductionLib"] = []
             if cuttools:
@@ -2186,6 +2206,8 @@
                 MLoptions["MLReductionLib"].extend([5])
             if "ninja_dir" in tir:
                 MLoptions["MLReductionLib"].extend([6])
+            if "collier_dir" in tir:
+                MLoptions["MLReductionLib"].extend([7])
 
     stability = myStabilityChecker.check_matrix_element_stability(matrix_element, 
                         options=options,param_card=param_card, 
@@ -2223,6 +2245,19 @@
         MLoptions = {}
     else:
         MLoptions = MLOptions
+        # Make sure that the poles computation is disabled for COLLIER
+        if 'COLLIERComputeUVpoles' not in MLoptions:
+            MLoptions['COLLIERComputeUVpoles']=False
+        if 'COLLIERComputeIRpoles' not in MLoptions:
+            MLoptions['COLLIERComputeIRpoles']=False
+        # And the COLLIER global cache is active, if not specified
+        if 'COLLIERGlobalCache' not in MLoptions:
+            MLoptions['COLLIERGlobalCache']=-1
+        # And time NINJA by default if not specified:
+        if 'MLReductionLib' not in MLoptions or \
+                                            len(MLoptions['MLReductionLib'])==0:
+            MLoptions['MLReductionLib'] = [6]
+            
     timing2 = myTimer.time_matrix_element(matrix_element, reusing, param_card,
                                      keep_folder = keep_folder, options=options,
                                      MLOptions = MLoptions)
@@ -2233,6 +2268,7 @@
         # Return the merged two dictionaries
         res = dict(timing1.items()+timing2.items())
         res['loop_optimized_output']=myTimer.loop_optimized_output
+        res['reduction_tool'] = MLoptions['MLReductionLib'][0]
         return res
 
 #===============================================================================
@@ -2678,7 +2714,8 @@
         # Use nicer name for the XML tag in the log file
         xml_toolname = {'GOLEM95':'GOLEM','IREGI':'IREGI',
                         'CUTTOOLS':'CUTTOOLS','PJFRY++':'PJFRY',
-                        'NINJA':'NINJA','SAMURAI':'SAMURAI'}[toolname.upper()]
+                        'NINJA':'NINJA','SAMURAI':'SAMURAI',
+                        'COLLIER':'COLLIER'}[toolname.upper()]
         if len(UPS)>0:
             res_str_i = "\nDetails of the %d/%d UPS encountered by %s\n"\
                                                         %(len(UPS),nPS,toolname)
@@ -2819,7 +2856,7 @@
 
     try:
         import matplotlib.pyplot as plt
-        colorlist=['b','r','g','y','m','c']
+        colorlist=['b','r','g','y','m','c','k']
         for i,key in enumerate(data_plot_dict.keys()):
             color=colorlist[i]
             data_plot=data_plot_dict[key]
@@ -2860,6 +2897,8 @@
     # Define shortcut
     f = format_output
     loop_optimized_output = timings['loop_optimized_output']
+    reduction_tool        = bannermod.MadLoopParam._ID_reduction_tool_map[
+                                                      timings['reduction_tool']]
     
     res_str = "%s \n"%process.nice_string()
     try:
@@ -2885,6 +2924,7 @@
     res_str += "|= Initialization............ %s\n"\
                                             %f(timings['Initialization'],'%.3gs')
 
+    res_str += "\n= Reduction tool tested...... %s\n"%reduction_tool
     res_str += "\n= Helicity sum time / PSpoint ========== %.3gms\n"\
                                     %(timings['run_unpolarized_total']*1000.0)
     if loop_optimized_output:
@@ -2894,7 +2934,7 @@
         total=coef_time+loop_time
         res_str += "|= Coefs. computation time... %.3gms (%d%%)\n"\
                                   %(coef_time,int(round(100.0*coef_time/total)))
-        res_str += "|= Loop evaluation (OPP) time %.3gms (%d%%)\n"\
+        res_str += "|= Loop evaluation time...... %.3gms (%d%%)\n"\
                                   %(loop_time,int(round(100.0*loop_time/total)))
     res_str += "\n= One helicity time / PSpoint ========== %.3gms\n"\
                                     %(timings['run_polarized_total']*1000.0)
@@ -2905,7 +2945,7 @@
         total=coef_time+loop_time        
         res_str += "|= Coefs. computation time... %.3gms (%d%%)\n"\
                                   %(coef_time,int(round(100.0*coef_time/total)))
-        res_str += "|= Loop evaluation (OPP) time %.3gms (%d%%)\n"\
+        res_str += "|= Loop evaluation time...... %.3gms (%d%%)\n"\
                                   %(loop_time,int(round(100.0*loop_time/total)))
     res_str += "\n= Miscellaneous ========================\n"
     res_str += "|= Number of hel. computed... %s/%s\n"\

=== modified file 'madgraph/various/q_polynomial.py'
--- madgraph/various/q_polynomial.py	2016-06-09 00:59:23 +0000
+++ madgraph/various/q_polynomial.py	2016-07-02 20:37:49 +0000
@@ -362,6 +362,20 @@
         
         return '\n'.join(lines)
     
+    def get_COLLIER_mapping(self):
+        """ Returns a list of tuples of the form:
+          [ (COLLIER_ind0, COLLIER_ind1, COLLIER_ind2, COLLIER_ind3), ]
+          where the position in the list is the coef_ID in MadLoop ordering.
+        """        
+        res = []
+        for coef_pos in range(0,get_number_of_coefs_for_rank(self.pq.rank)):
+            indices_list = self.pq.get_coef_at_position(coef_pos)
+            res.append((indices_list.count(0),
+                        indices_list.count(1),
+                        indices_list.count(2),
+                        indices_list.count(3)))
+        return res
+    
     def write_golem95_mapping(self):
         """ Returns a fortran subroutine which fills in the array of tensorial
         coefficients following golem95 standards using MadLoop coefficients."""

=== modified file 'tests/IOTests.py'
--- tests/IOTests.py	2015-10-01 16:00:08 +0000
+++ tests/IOTests.py	2016-07-02 20:37:49 +0000
@@ -155,7 +155,7 @@
                                         myloopamp,optimized_output=isOptimized)
 
         self.exporter.copy_v4template(model.get('name'))
-        self.exporter.generate_loop_subprocess(hel_amp, self.helasModel)
+        self.exporter.generate_loop_subprocess(hel_amp, self.helasModel,unique_id=1)
         wanted_lorentz = hel_amp.get_used_lorentz()
         wanted_couplings = list(set(sum(hel_amp.get_used_couplings(),[])))
         self.exporter.convert_model_to_mg4(model,wanted_lorentz,wanted_couplings)

=== modified file 'tests/acceptance_tests/test_cmd.py'
--- tests/acceptance_tests/test_cmd.py	2016-05-11 21:07:38 +0000
+++ tests/acceptance_tests/test_cmd.py	2016-07-02 20:37:49 +0000
@@ -189,6 +189,7 @@
                     'f2py_compiler':None,
                     'cluster_retry_wait': 300,
                     'syscalc_path':'./SysCalc',
+                    'collier':'./HEPTools/lib',
                     'hepmc_path': './hepmc',
                     'hwpp_path': './herwigPP',
                     'thepeg_path': './thepeg',

=== modified file 'tests/acceptance_tests/test_histograms.py'
--- tests/acceptance_tests/test_histograms.py	2016-06-09 00:59:23 +0000
+++ tests/acceptance_tests/test_histograms.py	2016-07-02 20:37:49 +0000
@@ -47,6 +47,5 @@
           histo_list.output(pjoin('%s','HistoOut'), format = 'gnuplot');
         ''' % (sys.path, _file_path, self.IOpath)
 
-        import os
         os.system('echo "%s" | python' % line) 
 

=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%COLLIER_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%COLLIER_interface.f	1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%COLLIER_interface.f	2016-07-02 20:37:49 +0000
@@ -0,0 +1,732 @@
+      SUBROUTINE COLLIERLOOP(CTMODE, NLOOPLINE,RANK,PL,PDEN,M2L
+     $ ,TIRCOEFS,TIRCOEFSERRORS)
+C     
+C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
+C     By the MadGraph5_aMC@NLO Development Team
+C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
+C     
+C     Interface between MG5 and COLLIER.
+C     It supports any rank when 1 < NLOOPLINE < 7.
+C     
+C     Process: d~ u > w+ WEIGHTED<=2 QED<=1 [ all = QCD ]
+C     Process: s~ c > w+ WEIGHTED<=2 QED<=1 [ all = QCD ]
+C     
+C     
+C     MODULES
+C     
+      USE COLLIER
+      IMPLICIT NONE
+C     
+C     CONSTANTS 
+C     
+      LOGICAL CHECKPCONSERVATION
+      PARAMETER (CHECKPCONSERVATION=.TRUE.)
+      REAL*8 PIVALUE
+      PARAMETER(PIVALUE=3.14159265358979323846D0)
+      INTEGER MAXRANK
+      PARAMETER (MAXRANK=2)
+      INTEGER NLOOPGROUPS
+      PARAMETER (NLOOPGROUPS=1)
+      INCLUDE 'loop_max_coefs.inc'
+C     
+C     ARGUMENTS 
+C     
+      INTEGER NLOOPLINE, RANK, CTMODE
+      REAL*8 PL(0:3,NLOOPLINE)
+      REAL*8 PDEN(0:3,NLOOPLINE-1)
+      COMPLEX*16 M2L(NLOOPLINE)
+      COMPLEX*16 RES(3)
+      COMPLEX*16 TIRCOEFS(0:LOOPMAXCOEFS-1,3)
+      COMPLEX*16 TIRCOEFSERRORS(0:LOOPMAXCOEFS-1,3)
+
+C     
+C     LOCAL VARIABLES 
+C     
+      INTEGER N, I, J, K, L
+      INTEGER C0,C1,C2,C3
+      INTEGER N_CACHES
+      INTEGER CURR_RANK, SGN
+      INTEGER CURR_INDEX
+      INTEGER CURR_MAXCOEF
+      REAL*8 RBUFF(0:3)
+
+      INTEGER COEFMAP_ZERO(0:LOOPMAXCOEFS-1)
+      INTEGER COEFMAP_ONE(0:LOOPMAXCOEFS-1)
+      INTEGER COEFMAP_TWO(0:LOOPMAXCOEFS-1)
+      INTEGER COEFMAP_THREE(0:LOOPMAXCOEFS-1)
+      DATA (COEFMAP_ZERO(I),I=  0,  9) / 0, 1, 0, 0, 0, 2, 1, 0, 1, 0/
+      DATA (COEFMAP_ZERO(I),I= 10, 14) / 0, 1, 0, 0, 0/
+      DATA (COEFMAP_ONE(I),I=  0,  9) / 0, 0, 1, 0, 0, 0, 1, 2, 0, 1/
+      DATA (COEFMAP_ONE(I),I= 10, 14) / 0, 0, 1, 0, 0/
+      DATA (COEFMAP_TWO(I),I=  0,  9) / 0, 0, 0, 1, 0, 0, 0, 0, 1, 1/
+      DATA (COEFMAP_TWO(I),I= 10, 14) / 2, 0, 0, 1, 0/
+      DATA (COEFMAP_THREE(I),I=  0,  9) / 0, 0, 0, 0, 1, 0, 0, 0, 0, 0/
+      DATA (COEFMAP_THREE(I),I= 10, 14) / 0, 1, 1, 1, 2/
+
+      REAL*8 REF_NORMALIZATION
+
+      DOUBLE COMPLEX M2LCOLLIER(0:NLOOPLINE-1)
+      DOUBLE COMPLEX MOMVEC(0:3,NLOOPLINE-1)
+      DOUBLE COMPLEX MOMINV((NLOOPLINE*(NLOOPLINE-1))/2)
+
+      DOUBLE COMPLEX TNTEN(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE COMPLEX TNTENUV(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE PRECISION TNTENERR(0:RANK)
+
+      REAL*8 MAXCOEFFORRANK(0:RANK)
+
+C     These quantities are for the pole evaluation
+      DOUBLE COMPLEX TNTEN_UV(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE COMPLEX TNTENUV_UV(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE PRECISION TNTENERR_UV(0:RANK)
+      DOUBLE COMPLEX TNTEN_IR1(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE COMPLEX TNTENUV_IR1(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE PRECISION TNTENERR_IR1(0:RANK)
+      DOUBLE COMPLEX TNTEN_IR2(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE COMPLEX TNTENUV_IR2(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE PRECISION TNTENERR_IR2(0:RANK)
+C     
+C     GLOBAL VARIABLES
+C     
+      INCLUDE 'coupl.inc'
+      INCLUDE 'MadLoopParams.inc'
+      INCLUDE 'unique_id.inc'
+      INCLUDE 'global_specs.inc'
+
+      INTEGER ID,SQSOINDEX,R
+      COMMON/LOOP/ID,SQSOINDEX,R
+
+C     
+C     These global variables cache the coefficients already computed
+C      by COLLIER
+C     so that they can be reused. In particular, in CTMODE 2, the
+C      cached ERROR 
+C     quantities for eahc rank will be used to make use of COLLIER's
+C      internal
+C     assessment of the numerical accuracy.
+C     
+      LOGICAL LOOP_ID_DONE(NLOOPGROUPS)
+      DATA LOOP_ID_DONE/NLOOPGROUPS*.FALSE./
+      COMPLEX*16 TIR_COEFS_DIRECT_MODE(0:LOOPMAXCOEFS-1,3,NLOOPGROUPS)
+      REAL*8 TIR_COEFS_ERROR(0:MAXRANK,NLOOPGROUPS)
+      COMMON/COLLIER_TIR_COEFS/TIR_COEFS_DIRECT_MODE,TIR_COEFS_ERROR
+     $ ,LOOP_ID_DONE
+
+      INTEGER COLLIER_CACHE_ACTIVE, NCALLS_IN_CACHE(4),
+     $  NCOLLIERCALLS(4)
+      LOGICAL MUST_INIT_EVENT
+      COMMON/COLLIER_CACHE_STATUS/COLLIER_CACHE_ACTIVE,NCALLS_IN_CACHE
+     $ , NCOLLIERCALLS,MUST_INIT_EVENT
+
+      LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT,
+     $  COLLIERINIT
+      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT,GOLEMINIT,SAMURAIINIT
+     $ ,NINJAINIT, COLLIERINIT
+
+C     ----------
+C     BEGIN CODE
+C     ----------
+
+      IF (COLLIERUSECACHEFORPOLES) THEN
+        N_CACHES = 4
+      ELSE
+        N_CACHES =1
+      ENDIF
+
+C     Initialize COLLIER if needed 
+      IF (COLLIERINIT) THEN
+        COLLIERINIT=.FALSE.
+        CALL INITCOLLIER()
+      ENDIF
+
+C     Initialize the event if it is the first time collier is called
+C      on this PS point
+      IF(MUST_INIT_EVENT) THEN
+        MUST_INIT_EVENT = .FALSE.
+        IF (COLLIERGLOBALCACHE.EQ.0) THEN
+          CALL INITEVENT_CLL()
+        ELSE
+          DO I=1,N_CACHES
+C           Record how many events where put in the cache. On the
+C            first PS point.
+            IF (NCALLS_IN_CACHE(I).EQ.-1.AND.NCOLLIERCALLS(I).GT.0)
+     $        THEN
+              NCALLS_IN_CACHE(I) = NCOLLIERCALLS(I)
+            ENDIF
+          ENDDO
+          DO I=1,N_CACHES
+C           Now apply a safety check that our last event had as many
+C            calls as the cache is setup for.
+C           The only case for now when it can be half of the calls
+C            when we are doing the true loop-direction test with also
+C            the computation of a rotated PS point (which is computed
+C            for one mode only).
+            IF (NCALLS_IN_CACHE(I).NE.-1.AND..NOT. ( NCOLLIERCALLS(I)
+     $       .EQ.NCALLS_IN_CACHE(I).OR. ( CTMODERUN.EQ.
+     $       -1.AND..NOT.COLLIERUSEINTERNALSTABILITYTEST.AND.NROTATIONS
+     $_DP.GT.0.AND.MOD(NCALLS_IN_CACHE(I),2).EQ.0.AND.(NCALLS_IN_CACHE(
+     $I)/2).EQ.NCOLLIERCALLS(I) ) ) ) THEN
+              WRITE(*,*) 'WARNING: A consistency check in MadLoop'
+     $         //' failed and, for safety, forced MadLoop to'
+     $         //' reinitialize the global cache of COLLIER. Report'
+     $         //' this to MadLoop authors. The problematic cache was'
+     $         //' number ',I
+              IF (COLLIERGLOBALCACHE.EQ.-1) THEN
+                CALL INITCACHESYSTEM_CLL(N_CACHES*NPROCS,MAXNEXTERNAL)
+              ELSE
+                CALL INITCACHESYSTEM_CLL(N_CACHES*NPROCS
+     $           ,COLLIERGLOBALCACHE)
+              ENDIF
+C             Make sure all caches are switched off at first.
+              CALL SWITCHOFFCACHESYSTEM_CLL()
+C             Reset the cache design property
+              NCALLS_IN_CACHE(:) = -1
+              NCOLLIERCALLS(:)   = 0
+              IF (COLLIER_CACHE_ACTIVE.EQ.1) THEN
+                CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*N_CACHES+1)
+              ENDIF
+C             No need to check the other caches since we already had
+C              to reset here.
+              EXIT
+            ENDIF
+          ENDDO
+          CALL INITEVENT_CLL((UNIQUE_ID-1)*N_CACHES+1)
+          NCOLLIERCALLS(1) = 0
+          IF(COLLIERCOMPUTEUVPOLES.AND.COLLIERUSECACHEFORPOLES) THEN
+            CALL INITEVENT_CLL((UNIQUE_ID-1)*N_CACHES+2)
+            NCOLLIERCALLS(2) = 0
+          ENDIF
+          IF(COLLIERCOMPUTEIRPOLES.AND.COLLIERUSECACHEFORPOLES) THEN
+            CALL INITEVENT_CLL((UNIQUE_ID-1)*N_CACHES+3)
+            NCOLLIERCALLS(3) = 0
+            CALL INITEVENT_CLL((UNIQUE_ID-1)*N_CACHES+4)
+            NCOLLIERCALLS(4) = 0
+          ENDIF
+        ENDIF
+        CALL SETDELTAIR_CLL(0.0D0,(PIVALUE**2)/6.0D0)
+        CALL SETDELTAUV_CLL(0.0D0)
+        CALL SETMUIR2_CLL(MU_R**2)
+        CALL SETMUUV2_CLL(MU_R**2)
+      ENDIF
+
+C     Now really start the reduction with COLLIER
+
+C     Number of coefficients for the current rank
+      CURR_MAXCOEF = 0
+      DO I=0,RANK
+        CURR_MAXCOEF=CURR_MAXCOEF+(3+I)*(2+I)*(1+I)/6
+      ENDDO
+
+      IF (CTMODE.NE.1.AND.CTMODE.NE.2) THEN
+        WRITE(*,*) 'ERROR: COLLIER only support the computational mode'
+     $   //' CTMODE equal to 1 or 2, not',CTMODE
+        STOP 'Incorrect computational mode specified to the COLLIER'
+     $   //' MG5aMC interface.'
+      ENDIF
+
+      DO I=0,CURR_MAXCOEF-1
+        DO K=1,3
+          TIRCOEFSERRORS(I,K)=DCMPLX(0.0D0,0.0D0)
+        ENDDO
+      ENDDO
+
+      IF (COLLIERUSEINTERNALSTABILITYTEST) THEN
+C       Use MADLOOP internal cache dedicated to COLLIER that emulates
+C        the CTMODE 2
+        IF (LOOP_ID_DONE(ID)) THEN
+          DO I=0,CURR_MAXCOEF-1
+            DO K=1,3
+              TIRCOEFS(I,K)=TIR_COEFS_DIRECT_MODE(I,K,ID)
+            ENDDO
+          ENDDO
+          DO I=0,RANK
+            MAXCOEFFORRANK(I) = 0.0D0
+          ENDDO
+          DO I=0,CURR_MAXCOEF-1
+            CURR_RANK = COEFMAP_ZERO(I)+COEFMAP_ONE(I)+COEFMAP_TWO(I)
+     $       +COEFMAP_THREE(I)
+            MAXCOEFFORRANK(CURR_RANK)=MAX(MAXCOEFFORRANK(CURR_RANK)
+     $       ,ABS(TIR_COEFS_DIRECT_MODE(I,1,ID)))
+          ENDDO
+          DO I=0,CURR_MAXCOEF-1
+            CURR_RANK = COEFMAP_ZERO(I)+COEFMAP_ONE(I)+COEFMAP_TWO(I)
+     $       +COEFMAP_THREE(I)
+            IF (MAXCOEFFORRANK(CURR_RANK).NE.0.0D0) THEN
+              DO K=1,3
+C               The expression below is like taking the absolute value
+C                when summing errors linearly
+C               TIRCOEFSERRORS(I,K)=(TIR_COEFS_ERROR(CURR_RANK,ID)/MaxC
+C               oefForRank(CURR_RANK))*DCMPLX( ABS(DBLE(TIRCOEFS(I,K)))
+C               ,ABS(DIMAG(TIRCOEFS(I,K))) )
+C               But empirically, I observed that retaining the
+C                original complex phase leads to slightly more
+C                accurate estimates
+                TIRCOEFSERRORS(I,K)=(TIR_COEFS_ERROR(CURR_RANK,ID)
+     $           /MAXCOEFFORRANK(CURR_RANK))*TIRCOEFS(I,K)
+              ENDDO
+            ENDIF
+          ENDDO
+          RETURN
+        ENDIF
+      ELSE
+C       Apply the loop-direction switching here.
+        CALL SWITCH_ORDER(CTMODE,NLOOPLINE,PL,PDEN,M2L)
+      ENDIF
+
+C     Make sure masses are complex 
+      DO I=1,NLOOPLINE
+        M2LCOLLIER(I-1)=DCMPLX(M2L(I))
+      ENDDO
+
+C     Now convert the loop offset momenta to COLLIER conventions
+      DO I=0,3
+        DO J=1,NLOOPLINE-1
+          MOMVEC(I,J)=DCMPLX(PDEN(I,J),0.0D0)
+        ENDDO
+      ENDDO
+
+C     This first do loop spans over 'N' in '\foreach_N'
+C     //' \foreach_i(k_i+k_{i+1}+..+k_{i+N})^2'
+      CURR_INDEX = 0
+      DO N=0,NLOOPLINE-1
+C       We stop whenever we reached the target number of invariants
+        IF (CURR_INDEX.GE.((NLOOPLINE*(NLOOPLINE-1))/2)) THEN
+          EXIT
+        ENDIF
+C       This do loop spans over 'i' in the expression above.
+        DO I=1,NLOOPLINE
+          IF (CURR_INDEX.GE.((NLOOPLINE*(NLOOPLINE-1))/2)) THEN
+            EXIT
+          ENDIF
+
+          CURR_INDEX = CURR_INDEX+1
+          RBUFF(:) = 0.0D0
+          DO J=I,I+N
+            RBUFF(:) = RBUFF(:) + PL(:,MOD(J-1,NLOOPLINE)+1)
+          ENDDO
+          MOMINV(CURR_INDEX) = DCMPLX(RBUFF(0)**2-RBUFF(1)**2-RBUFF(2)
+     $     **2-RBUFF(3)**2,0.0D0)
+
+C         Now regularize the onshell behavior of the kinematic
+C          invarients
+C         All loop masses are tested, but that might be a bit too
+C          inclusive.
+          DO K=1,NLOOPLINE
+            IF(ABS(M2L(K)).NE.0.0D0) THEN
+              IF(ABS((MOMINV(CURR_INDEX)-M2L(K))/M2L(K)).LT.OSTHRES)
+     $          THEN
+                MOMINV(CURR_INDEX)=DCMPLX(M2L(K))
+              ENDIF
+            ENDIF
+          ENDDO
+C         For the massless onshell-case, we base the threshold only on
+C          the energy component
+          REF_NORMALIZATION=0.0D0
+          DO L=0,0
+            REF_NORMALIZATION = REF_NORMALIZATION + ABS(RBUFF(L))
+          ENDDO
+          REF_NORMALIZATION = (REF_NORMALIZATION/(N+1))**2
+          IF(REF_NORMALIZATION.NE.0.0D0)THEN
+            IF(ABS(MOMINV(CURR_INDEX)/REF_NORMALIZATION).LT.OSTHRES)
+     $       THEN
+              MOMINV(CURR_INDEX)=DCMPLX(0.0D0,0.0D0)
+            ENDIF
+          ENDIF
+
+        ENDDO
+      ENDDO
+
+C     We can now call COLLIER
+      IF (NLOOPLINE.NE.1) THEN
+        CALL TNTEN_CLL(TNTEN, TNTENUV, MOMVEC, MOMINV, M2LCOLLIER,
+     $    NLOOPLINE, RANK, TNTENERR)
+      ELSE
+        CALL TNTEN_CLL(TNTEN, TNTENUV, M2LCOLLIER, NLOOPLINE, RANK,
+     $    TNTENERR)
+      ENDIF
+      IF (COLLIER_CACHE_ACTIVE.EQ.1) THEN
+        NCOLLIERCALLS(1) = NCOLLIERCALLS(1)+1
+      ENDIF
+
+C     Now compute the UV poles if asked for
+      IF (COLLIERCOMPUTEUVPOLES) THEN
+        IF(COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          CALL SWITCHOFFCACHE_CLL((UNIQUE_ID-1)*N_CACHES+1)
+          IF(COLLIERUSECACHEFORPOLES) THEN
+            CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*N_CACHES+2)
+          ENDIF
+        ENDIF
+        CALL SETDELTAUV_CLL(1.0D0)
+        IF (NLOOPLINE.NE.1) THEN
+          CALL TNTEN_CLL(TNTEN_UV, TNTENUV_UV, MOMVEC, MOMINV,
+     $      M2LCOLLIER, NLOOPLINE, RANK, TNTENERR_UV)
+        ELSE
+          CALL TNTEN_CLL(TNTEN_UV, TNTENUV_UV, M2LCOLLIER, NLOOPLINE,
+     $      RANK, TNTENERR_UV)
+        ENDIF
+        IF (COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          NCOLLIERCALLS(2) = NCOLLIERCALLS(2)+1
+        ENDIF
+        IF(COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          IF(COLLIERUSECACHEFORPOLES) THEN
+            CALL SWITCHOFFCACHE_CLL((UNIQUE_ID-1)*N_CACHES+2)
+          ENDIF
+          CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*N_CACHES+1)
+        ENDIF
+        CALL SETDELTAUV_CLL(0.0D0)
+      ENDIF
+
+C     Now compute the IR poles if asked for
+      IF (COLLIERCOMPUTEIRPOLES) THEN
+        IF(COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          CALL SWITCHOFFCACHE_CLL((UNIQUE_ID-1)*N_CACHES+1)
+          IF(COLLIERUSECACHEFORPOLES) THEN
+            CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*N_CACHES+3)
+          ENDIF
+        ENDIF
+        CALL SETDELTAIR_CLL(1.0D0,(PIVALUE**2)/6.0D0)
+        IF (NLOOPLINE.NE.1) THEN
+          CALL TNTEN_CLL(TNTEN_IR1, TNTENUV_IR1, MOMVEC, MOMINV,
+     $      M2LCOLLIER, NLOOPLINE, RANK, TNTENERR_IR1)
+        ELSE
+          CALL TNTEN_CLL(TNTEN_IR1, TNTENUV_IR1, M2LCOLLIER, NLOOPLINE
+     $     , RANK, TNTENERR_IR1)
+        ENDIF
+        IF (COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          NCOLLIERCALLS(3) = NCOLLIERCALLS(3)+1
+        ENDIF
+        IF(COLLIER_CACHE_ACTIVE.EQ.1.AND.COLLIERUSECACHEFORPOLES) THEN
+          CALL SWITCHOFFCACHE_CLL((UNIQUE_ID-1)*N_CACHES+3)
+          CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*N_CACHES+4)
+        ENDIF
+        CALL SETDELTAIR_CLL(0.0D0,1.0D0+(PIVALUE**2)/6.0D0)
+        IF (NLOOPLINE.NE.1) THEN
+          CALL TNTEN_CLL(TNTEN_IR2, TNTENUV_IR2, MOMVEC, MOMINV,
+     $      M2LCOLLIER, NLOOPLINE, RANK, TNTENERR_IR2)
+        ELSE
+          CALL TNTEN_CLL(TNTEN_IR2, TNTENUV_IR2, M2LCOLLIER, NLOOPLINE
+     $     , RANK, TNTENERR_IR2)
+        ENDIF
+        IF (COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          NCOLLIERCALLS(4) = NCOLLIERCALLS(4)+1
+        ENDIF
+        IF(COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          IF(COLLIERUSECACHEFORPOLES) THEN
+            CALL SWITCHOFFCACHE_CLL((UNIQUE_ID-1)*N_CACHES+4)
+          ENDIF
+          CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*N_CACHES+1)
+        ENDIF
+        CALL SETDELTAIR_CLL(0.0D0,(PIVALUE**2)/6.0D0)
+      ENDIF
+
+      DO I=0,(CURR_MAXCOEF-1)
+        C0 = COEFMAP_ZERO(I)
+        C1 = COEFMAP_ONE(I)
+        C2 = COEFMAP_TWO(I)
+        C3 = COEFMAP_THREE(I)
+        CURR_RANK = C0+C1+C2+C3
+C       Because we must set q -> -q, we apply a minus sign to coefs of
+C        odd rank
+        IF (MOD(CURR_RANK,2).EQ.1) THEN
+          SGN = -1
+        ELSE
+          SGN = 1
+        ENDIF
+        TIRCOEFS(I,1) = SGN*TNTEN(C0,C1,C2,C3)
+        IF (COLLIERCOMPUTEUVPOLES) THEN
+          TIRCOEFS(I,2) = SGN*( TNTEN_UV(C0,C1,C2,C3)-TNTEN(C0,C1,C2
+     $     ,C3) )
+        ELSE
+          TIRCOEFS(I,2) = DCMPLX(0.0D0,0.0D0)
+        ENDIF
+        IF (COLLIERCOMPUTEIRPOLES) THEN
+          TIRCOEFS(I,2) = TIRCOEFS(I,2) + SGN*( TNTEN_IR1(C0,C1,C2,C3)
+     $     -TNTEN(C0,C1,C2,C3) )
+          TIRCOEFS(I,3) = SGN*( TNTEN_IR2(C0,C1,C2,C3)-TNTEN(C0,C1,C2
+     $     ,C3) )
+        ELSE
+          TIRCOEFS(I,3) = DCMPLX(0.0D0,0.0D0)
+        ENDIF
+      ENDDO
+
+      IF (COLLIERUSEINTERNALSTABILITYTEST) THEN
+C       Finish by caching the coefficients and error just computed
+        LOOP_ID_DONE(ID) = .TRUE.
+        DO J=0,CURR_MAXCOEF-1
+          DO K=1,3
+            TIR_COEFS_DIRECT_MODE(J,K,ID) = TIRCOEFS(J,K)
+          ENDDO
+        ENDDO
+        DO J=0,RANK
+          TIR_COEFS_ERROR(J,ID)=TNTENERR(J)
+        ENDDO
+
+C       Now compute the errors on each coefficient
+        DO I=0,RANK
+          MAXCOEFFORRANK(I) = 0.0D0
+        ENDDO
+        DO I=0,CURR_MAXCOEF-1
+          CURR_RANK = COEFMAP_ZERO(I)+COEFMAP_ONE(I)+COEFMAP_TWO(I)
+     $     +COEFMAP_THREE(I)
+          MAXCOEFFORRANK(CURR_RANK)=MAX(MAXCOEFFORRANK(CURR_RANK)
+     $     ,ABS(TIRCOEFS(I,1)))
+        ENDDO
+        DO I=0,CURR_MAXCOEF-1
+          CURR_RANK = COEFMAP_ZERO(I)+COEFMAP_ONE(I)+COEFMAP_TWO(I)
+     $     +COEFMAP_THREE(I)
+          IF (MAXCOEFFORRANK(CURR_RANK).NE.0.0D0) THEN
+            DO K=1,3
+C             The expression below is like taking the absolute value
+C              when summing errors linearly
+C             TIRCOEFSERRORS(I,K)=(TNtenerr(CURR_RANK)/MaxCoefForRank(C
+C             URR_RANK))*DCMPLX( ABS(DBLE(TIRCOEFS(I,K))),ABS(DIMAG(TIR
+C             COEFS(I,K))) )			
+C             But empirically, I observed that retaining the original
+C              complex phase leads to slightly more accurate estimates
+              TIRCOEFSERRORS(I,K)=(TNTENERR(CURR_RANK)
+     $         /MAXCOEFFORRANK(CURR_RANK))*TIRCOEFS(I,K)
+            ENDDO
+          ENDIF
+        ENDDO
+
+      ENDIF
+
+      END SUBROUTINE COLLIERLOOP
+
+      SUBROUTINE CLEAR_COLLIER_CACHE()
+
+      USE COLLIER
+
+      INCLUDE 'loop_max_coefs.inc'
+      INTEGER NLOOPGROUPS
+      PARAMETER (NLOOPGROUPS=1)
+      INTEGER MAXRANK
+      PARAMETER (MAXRANK=2)
+
+      INTEGER I,J,K
+
+      INCLUDE 'MadLoopParams.inc'
+
+      LOGICAL LOOP_ID_DONE(NLOOPGROUPS)
+      COMPLEX*16 TIR_COEFS_DIRECT_MODE(0:LOOPMAXCOEFS-1,3,NLOOPGROUPS)
+      COMPLEX*16 TIR_COEFS_ERROR(0:MAXRANK,NLOOPGROUPS)
+      COMMON/COLLIER_TIR_COEFS/TIR_COEFS_DIRECT_MODE,TIR_COEFS_ERROR
+     $ ,LOOP_ID_DONE
+
+      INTEGER COLLIER_CACHE_ACTIVE, NCALLS_IN_CACHE(4),
+     $  NCOLLIERCALLS(4)
+      LOGICAL MUST_INIT_EVENT
+      COMMON/COLLIER_CACHE_STATUS/COLLIER_CACHE_ACTIVE,NCALLS_IN_CACHE
+     $ , NCOLLIERCALLS,MUST_INIT_EVENT
+
+C     Make sure that next time the COLLIER Subroutine is called it
+C      will call the subroutine initEvent of Collier.
+      MUST_INIT_EVENT = .TRUE.
+
+C     Reinitialize the ML cache for COLLIER
+      DO I=1,NLOOPGROUPS
+        LOOP_ID_DONE(I) = .FALSE.
+        DO J=0,LOOPMAXCOEFS-1
+          DO K=1,3
+            TIR_COEFS_DIRECT_MODE(J,K,I) = DCMPLX(0.0D0,0.0D0)
+          ENDDO
+        ENDDO
+        DO J=0,MAXRANK
+          TIR_COEFS_ERROR(J,I)=0.0D0
+        ENDDO
+      ENDDO
+
+      END
+
+      SUBROUTINE SET_COLLIER_GLOBAL_CACHE(ONOFF)
+C     
+C     This routine is used by loop_matrix.f to turn on the global
+C     cache of COLLIER when it the main SLOOP subroutine starts and
+C     turn it off when it ends.
+C     However several checks are performed to verify that it is safe
+C     to turn it on and to reinitialize it if necessary.
+C     
+C     MODULES
+C     
+      USE COLLIER
+      IMPLICIT NONE
+C     
+C     ARGUMENTS
+C     
+      LOGICAL ONOFF
+C     
+C     LOCAL VARIABLES
+C     
+      LOGICAL NEED_REINITIALIZATION
+      INTEGER N_CACHES
+C     
+C     GLOBAL VARIABLES
+C     
+C     
+C     This common blocks saves the relevant ML parameters when
+C      activating the
+C     global cache of COLLIER so that we know when we must
+C      reinitialize it.
+C     COLLIER_CACHE_ACTIVE is set to -1 when it has never been turned
+C      on yet and
+C     to 1 for 'Active' and 0 for 'Inactive'.
+C     The integer NCALLS_IN_CACHE saves how many calls the cache is
+C      setup for, for each of the up to four caches.
+C     When it is the first PS points, it is set to -1.
+      INTEGER COLLIER_CACHE_ACTIVE, NCALLS_IN_CACHE(4),
+     $  NCOLLIERCALLS(4)
+      DATA COLLIER_CACHE_ACTIVE/-1/
+      DATA NCALLS_IN_CACHE/-1,-1,-1,-1/
+      DATA NCOLLIERCALLS/0,0,0,0/
+C     This is a flag to tell the COLLIER subroutine that it must init
+C      the event when called.
+      LOGICAL MUST_INIT_EVENT
+      DATA MUST_INIT_EVENT/.TRUE./
+      COMMON/COLLIER_CACHE_STATUS/COLLIER_CACHE_ACTIVE,
+     $  NCALLS_IN_CACHE, NCOLLIERCALLS,MUST_INIT_EVENT
+
+      LOGICAL COLLIERUSEINTERNALSTABILITYTEST_BU
+      INTEGER USERHEL_BU, SQSO_TARGET_BU, COLLIERMODE_BU,CTMODERUN_BU
+      COMMON/COLLIER_CACHE_RELEVANT_PARAMS/USERHEL_BU,SQSO_TARGET_BU
+     $ ,COLLIERMODE_BU,CTMODERUN_BU,COLLIERUSEINTERNALSTABILITYTEST_BU
+
+C     The common blocks below are to retrieve the necessary
+C      information about
+C     MadLoop running mode and store it in the sCOLLIER_CACHE_RELEVANT_
+C     PARAMS common block.
+
+      INCLUDE 'MadLoopParams.inc'
+      INCLUDE 'unique_id.inc'
+      INCLUDE 'global_specs.inc'
+
+      LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT,
+     $  COLLIERINIT
+      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT,GOLEMINIT,SAMURAIINIT
+     $ ,NINJAINIT, COLLIERINIT
+
+      LOGICAL CHECKPHASE
+      LOGICAL HELDOUBLECHECKED
+      COMMON/INIT/CHECKPHASE, HELDOUBLECHECKED
+
+      INTEGER USERHEL
+      COMMON/USERCHOICE/USERHEL
+
+      INTEGER SQSO_TARGET
+      COMMON/SOCHOICE/SQSO_TARGET
+
+C     
+C     BEGIN CODE
+C     
+
+      IF (COLLIERUSECACHEFORPOLES) THEN
+        N_CACHES = 4
+      ELSE
+        N_CACHES =1
+      ENDIF
+
+C     Do nothing if COLLIER still has to be initialized or if global
+C      caches are disabled
+      IF(COLLIERINIT.OR.COLLIERGLOBALCACHE.EQ.0) THEN
+        RETURN
+      ENDIF
+
+C     Never activate anything in the checkphase
+      IF (ONOFF.AND.CHECKPHASE) THEN
+        RETURN
+      ENDIF
+
+C     Handle the request of turning off the caching
+      IF (.NOT.ONOFF) THEN
+        IF (COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          CALL SWITCHOFFCACHE_CLL((UNIQUE_ID-1)*4+1)
+          COLLIER_CACHE_ACTIVE = 0
+        ENDIF
+C       If we were asked to turn the cache off but it was already so,
+C        then do nothing		  
+        RETURN
+      ENDIF
+
+C     Handle the request of turning on the caching
+
+C     If asked to activate it but already active, then do nothing
+      IF (ONOFF.AND.COLLIER_CACHE_ACTIVE.EQ.1) THEN
+        RETURN
+      ENDIF
+
+C     We are now in the position where we are asked to activate the
+C      global cache but it was *not* already active.
+
+C     If we activate it for the first time, make sure to store the
+C      value of the relevant parameters, activate and return.
+      IF (COLLIER_CACHE_ACTIVE.EQ.-1) THEN
+        USERHEL_BU         = USERHEL
+        SQSO_TARGET_BU     = SQSO_TARGET
+        COLLIERMODE_BU     = COLLIERMODE
+        COLLIERUSEINTERNALSTABILITYTEST_BU = COLLIERUSEINTERNALSTABILIT
+     $YTEST
+        CTMODERUN_BU       = CTMODERUN
+        CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*4+1)
+        COLLIER_CACHE_ACTIVE = 1
+        RETURN
+      ENDIF
+
+C     Now perform sanity check before the activation to decide if we
+C      need to reinitialize the cache system first.
+      NEED_REINITIALIZATION = .FALSE.
+
+      IF (SQSO_TARGET.NE.SQSO_TARGET_BU) THEN
+        NEED_REINITIALIZATION = .TRUE.
+      ENDIF
+
+      IF (COLLIERMODE.NE.COLLIERMODE_BU) THEN
+        NEED_REINITIALIZATION = .TRUE.
+      ENDIF
+
+      IF (COLLIERUSEINTERNALSTABILITYTEST.NEQV.COLLIERUSEINTERNALSTABIL
+     $ITYTEST_BU) THEN
+        NEED_REINITIALIZATION = .TRUE.
+      ENDIF
+
+      IF (CTMODERUN_BU.NE.CTMODERUN.AND.(.NOT.COLLIERUSEINTERNALSTABILI
+     $TYTEST)) THEN
+        NEED_REINITIALIZATION = .TRUE.
+      ENDIF
+
+C     When doing amplitude reduction the parameter USERHEL does not
+C      impact the number/order of COLLIER calls
+C     except if the LoopFilter is ON which really shouldn't be the
+C      case anymore.
+      IF(USELOOPFILTER.AND.(USERHEL.NE.USERHEL_BU)) THEN
+        NEED_REINITIALIZATION = .TRUE.
+      ENDIF
+
+      IF(NEED_REINITIALIZATION) THEN
+C       Log the event because if it happens a lot of time and floods
+C        the screen, the user must see it
+C       and either change its usage of MadLoop or turnoff COLLIER cache
+        WRITE(*,*) 'INFO: MadLoop detected that the global cache of'
+     $   //' COLLIER had to be reset because of a change in your use'
+     $   //' of MadLoop. This should not happend for each event.'
+        USERHEL_BU         = USERHEL
+        SQSO_TARGET_BU     = SQSO_TARGET
+        COLLIERMODE_BU     = COLLIERMODE
+        COLLIERUSEINTERNALSTABILITYTEST_BU = COLLIERUSEINTERNALSTABILIT
+     $YTEST
+        CTMODERUN_BU       = CTMODERUN
+        IF (COLLIERGLOBALCACHE.EQ.-1) THEN
+          CALL INITCACHESYSTEM_CLL(N_CACHES*NPROCS,MAXNEXTERNAL)
+        ELSE
+          CALL INITCACHESYSTEM_CLL(N_CACHES*NPROCS,COLLIERGLOBALCACHE)
+        ENDIF
+        NCOLLIERCALLS(:)   = 0
+        NCALLS_IN_CACHE(:) = -1
+C       Make sure all caches are switched off at first.
+        CALL SWITCHOFFCACHESYSTEM_CLL()
+      ENDIF
+
+C     Now we can finally activate the cache
+      CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*4+1)
+      COLLIER_CACHE_ACTIVE = 1
+
+      END
+

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%CT_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%CT_interface.f	2016-05-06 09:08:11 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%CT_interface.f	2016-07-02 20:37:49 +0000
@@ -39,8 +39,9 @@
       COMPLEX*16 R1, ACC
       INTEGER I, J, K
       LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
+     $ ,COLLIERINIT
       COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT
-     $ ,NINJAINIT
+     $ ,NINJAINIT,COLLIERINIT
 C     
 C     EXTERNAL FUNCTIONS
 C     
@@ -392,8 +393,9 @@
       INCLUDE 'coupl.inc'
 
       LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
+     $ ,COLLIERINIT
       COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT
-     $ ,NINJAINIT
+     $ ,NINJAINIT,COLLIERINIT
 
       REAL*8 LSCALE
       INTEGER CTMODE
@@ -542,267 +544,9 @@
       END
 
 C     
-C     Quadruple precision version of loop_ninja
-C     
-      SUBROUTINE MP_NINJA_LOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
-C     
-C     Module used
-C     
-      USE MNINJA
-C     
-C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
-C     By the MadGraph5_aMC@NLO Development Team
-C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
-C     
-C     Interface between MG5 and Ninja.
-C     
-C     Process: d~ u > w+ WEIGHTED<=2 QED<=1 [ all = QCD ]
-C     Process: s~ c > w+ WEIGHTED<=2 QED<=1 [ all = QCD ]
-C     
-C     
-C     CONSTANTS 
-C     
-      INTEGER    NEXTERNAL
-      PARAMETER (NEXTERNAL=3)
-      LOGICAL CHECKPCONSERVATION
-      PARAMETER (CHECKPCONSERVATION=.TRUE.)
-      REAL*8 NORMALIZATION
-      PARAMETER (NORMALIZATION = 1.D0/(16.D0*3.14159265358979323846D0*
-     $ *2))
-      INTEGER NLOOPGROUPS
-      PARAMETER (NLOOPGROUPS=1)
-C     These are constants related to the split orders
-      INTEGER NSQUAREDSO
-      PARAMETER (NSQUAREDSO=1)
-      INCLUDE 'loop_max_coefs.inc'
-C     
-C     ARGUMENTS 
-C     
-      INTEGER NLOOPLINE, RANK
-      REAL*16 PL(0:3,NLOOPLINE)
-      COMPLEX*16 M2L(NLOOPLINE)
-      COMPLEX*16 RES(3)
-      LOGICAL STABLE
-C     
-C     LOCAL VARIABLES 
-C     
-      REAL*16 NINJA_SCALE
-      REAL*16 P_TMP(0:3,0:NLOOPLINE-1), ABSP_TMP(0:3)
-      REAL*8 REF_P
-      REAL(KI_QNIN) MP_P_NINJA(0:3,NLOOPLINE)
-      REAL*16 MP_P(0:3,NLOOPLINE)
-      REAL*16 P_S_MAT(NLOOPLINE,0:3)
-      COMPLEX*32 MP_M2L(NLOOPLINE)
-      COMPLEX(KI_QNIN) MP_M2L_NINJA(NLOOPLINE)
-      COMPLEX(KI_QNIN) NINJA_RES(0:2)
-      COMPLEX(KI_QNIN) NINJA_R1
-      COMPLEX*16 R1
-      COMPLEX*16 DP_RES(0:2)
-      INTEGER NINJA_STATUS
-      INTEGER I, J, K
-      REAL*16 PDEN_DUMMY(0:3,NLOOPLINE-1)
-
-      COMPLEX*32 MP_S_MAT(NLOOPLINE,NLOOPLINE)
-      REAL*16 MP_REAL_S_MAT(NLOOPLINE,NLOOPLINE)
-      REAL(KI_QNIN) MP_REAL_S_MAT_NINJA(NLOOPLINE,NLOOPLINE)
-
-      INTEGER CURR_MAXCOEF
-      COMPLEX*32, ALLOCATABLE :: MP_TENSORCOEFS(:)
-      COMPLEX(KI_QNIN), ALLOCATABLE :: MP_NINJA_TENSORCOEFS(:)
-
-C     
-C     GLOBAL VARIABLES
-C     
-      INCLUDE 'coupl.inc'
-
-      LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
-      COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT
-     $ ,NINJAINIT
-
-      REAL*8 LSCALE
-      INTEGER CTMODE
-      COMMON/CT/LSCALE,CTMODE
-
-      INTEGER ID,SQSOINDEX,R
-      COMMON/LOOP/ID,SQSOINDEX,R
-      COMPLEX*32 MP_LOOPCOEFS(0:LOOPMAXCOEFS-1,NSQUAREDSO,NLOOPGROUPS)
-      COMMON/MP_LCOEFS/MP_LOOPCOEFS
-
-      LOGICAL FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION
-      COMMON/FPE_IN_REDUCTION/FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION
-
-C     ----------
-C     BEGIN CODE
-C     ----------
-
-C     Cast the masses in complex quadruple precision
-      DO I=1,NLOOPLINE
-        MP_M2L(I) = CMPLX(M2L(I),KIND=16)
-      ENDDO
-
-C     For the direction test, we must switch the direction in which
-C      the loop is read for CTMode equal to 2 or 4.
-      CALL MP_SWITCH_ORDER(CTMODE,NLOOPLINE,PL,PDEN_DUMMY,MP_M2L)
-
-C     The CT initialization is also performed here if not done already
-C      because it calls MPINIT of OneLOop which is necessary on some
-C      system
-      IF (CTINIT) THEN
-        CTINIT=.FALSE.
-        CALL INITCT()
-      ENDIF
-
-C     INITIALIZE NINJA IF NEEDED
-      IF (NINJAINIT) THEN
-        NINJAINIT=.FALSE.
-        CALL INITNINJA()
-      ENDIF
-
-C     CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO NINJA
-C      CONVENTIONS
-      DO I=0,3
-        ABSP_TMP(I)=0.E0+0_16
-        DO J=0,(NLOOPLINE-1)
-          P_TMP(I,J)=0.E0+0_16
-        ENDDO
-      ENDDO
-      DO I=0,3
-        DO J=1,NLOOPLINE
-          P_TMP(I,0)=P_TMP(I,0)+PL(I,J)
-          ABSP_TMP(I)=ABSP_TMP(I)+ABS(PL(I,J))
-        ENDDO
-      ENDDO
-      REF_P = MAX(ABSP_TMP(0), ABSP_TMP(1),ABSP_TMP(2),ABSP_TMP(3))
-      DO I=0,3
-        ABSP_TMP(I) = MAX(REF_P*1E-6, ABSP_TMP(I))
-      ENDDO
-      IF (CHECKPCONSERVATION.AND.REF_P.GT.1.E-8_16) THEN
-        IF ((P_TMP(0,0)/ABSP_TMP(0)).GT.1.E-6_16) THEN
-          WRITE(*,*) 'energy is not conserved (flag:CT968)'
-     $     ,DBLE(P_TMP(0,0))
-          STOP 'energy is not conserved (flag:CT968)'
-        ELSEIF ((P_TMP(1,0)/ABSP_TMP(1)).GT.1.E-6_16) THEN
-          WRITE(*,*) 'px is not conserved (flag:CT968)',DBLE(P_TMP(1,0)
-     $     )
-          STOP 'px is not conserved (flag:CT968)'
-        ELSEIF ((P_TMP(2,0)/ABSP_TMP(2)).GT.1.E-6_16) THEN
-          WRITE(*,*) 'py is not conserved (flag:CT968)',DBLE(P_TMP(2,0)
-     $     )
-          STOP 'py is not conserved (flag:CT968)'
-        ELSEIF ((P_TMP(3,0)/ABSP_TMP(3)).GT.1.E-6_16) THEN
-          WRITE(*,*) 'pz is not conserved (flag:CT968)',DBLE(P_TMP(3,0)
-     $     )
-          STOP 'pz is not conserved (flag:CT968)'
-        ENDIF
-      ENDIF
-      DO I=0,3
-        DO J=1,(NLOOPLINE-1)
-          DO K=1,J
-            P_TMP(I,J)=P_TMP(I,J)+PL(I,K)
-          ENDDO
-        ENDDO
-      ENDDO
-C     In Ninja, the loop line index starts at 1
-      DO I=0,NLOOPLINE-1
-        MP_P(0,I+1) = P_TMP(0,I)
-        MP_P(1,I+1) = P_TMP(1,I)
-        MP_P(2,I+1) = P_TMP(2,I)
-        MP_P(3,I+1) = P_TMP(3,I)
-      ENDDO
-
-C     Number of coefficients for the current rank
-      CURR_MAXCOEF = 0
-      DO I=0,RANK
-        CURR_MAXCOEF=CURR_MAXCOEF+(3+I)*(2+I)*(1+I)/6
-      ENDDO
-C     Now write the tensor coefficients for Ninja
-C     It should never be allocated at this stage
-      IF (.NOT. ALLOCATED(MP_TENSORCOEFS)) THEN
-        ALLOCATE(MP_TENSORCOEFS(0:CURR_MAXCOEF-1))
-      ENDIF
-      IF (.NOT. ALLOCATED(MP_NINJA_TENSORCOEFS)) THEN
-        ALLOCATE(MP_NINJA_TENSORCOEFS(0:CURR_MAXCOEF-1))
-      ENDIF
-      DO I=0,CURR_MAXCOEF-1
-        MP_TENSORCOEFS(I) = MP_LOOPCOEFS(I,SQSOINDEX,ID)
-      ENDDO
-C     The loop momentum is in fact q_loop -> -q_loop, so that the
-C     coefficients must be changed accordingly
-      CALL MP_INVERT_MOMENTA_IN_POLYNOMIAL(CURR_MAXCOEF,MP_TENSORCOEFS)
-
-C     Compute the kinematic matrix
-      DO J=1,NLOOPLINE
-        DO I=0,3
-          P_S_MAT(J,I)=MP_P(I,J)
-        ENDDO
-      ENDDO
-      CALL MP_BUILD_KINEMATIC_MATRIX(NLOOPLINE,P_S_MAT,MP_M2L,MP_S_MAT)
-
-      DO I=1,NLOOPLINE
-        DO J=1,NLOOPLINE
-          MP_REAL_S_MAT(I,J) = REAL(MP_S_MAT(I,J)+MP_M2L(I)+MP_M2L(J)
-     $     ,KIND=16)
-        ENDDO
-      ENDDO
-
-C     Now typecast to Ninja's quadruple precision format
-      DO I=0,CURR_MAXCOEF-1
-        MP_NINJA_TENSORCOEFS(I)=CMPLX(MP_TENSORCOEFS(I),KIND=KI_QNIN)
-      ENDDO
-      DO I=1,NLOOPLINE
-        DO J=1,NLOOPLINE
-          MP_REAL_S_MAT_NINJA(I,J) = REAL(MP_REAL_S_MAT(I,J)
-     $     ,KIND=KI_QNIN)
-        ENDDO
-      ENDDO
-      DO I=1,NLOOPLINE
-        MP_M2L_NINJA(I)=CMPLX(MP_M2L(I),KIND=KI_QNIN)
-      ENDDO
-      DO I=1,NLOOPLINE
-        MP_P_NINJA(0,I) = REAL(MP_P(0,I),KIND=KI_QNIN)
-        MP_P_NINJA(1,I) = REAL(MP_P(1,I),KIND=KI_QNIN)
-        MP_P_NINJA(2,I) = REAL(MP_P(2,I),KIND=KI_QNIN)
-        MP_P_NINJA(3,I) = REAL(MP_P(3,I),KIND=KI_QNIN)
-      ENDDO
-      NINJA_SCALE = REAL(MU_R**2,KIND=KI_QNIN)
-
-
-C     Below is the call specifying the kinematic matrix
-      CALL NINJA_TENSOR_EVALUATE(MP_NINJA_TENSORCOEFS,NLOOPLINE,RANK
-     $ ,MP_REAL_S_MAT_NINJA,MP_P_NINJA,MP_M2L_NINJA,NINJA_SCALE
-     $ ,NINJA_RES,NINJA_R1,NINJA_STATUS)
-C     Below is the call without specification of the kinematic matrix
-C     call ninja_tensor_evaluate(MP_NINJA_TENSORCOEFS,NLOOPLINE,RANK,MP
-C     _P_NINJA,MP_M2L_NINJA,NINJA_SCALE,NINJA_RES,NINJA_R1,NINJA_STATUS)
-C     
-
-C     If a floating point exception was found in Ninja (e.g. exactly
-C      zero gram. det.)
-C     Then warn loop_matrix.f so that it will flag this kinematic
-C      point as unstable no matter what.
-      IF (NINJA_STATUS.EQ.NINJA_UNSTABLE_KINEMATICS) THEN
-        FPE_IN_QP_REDUCTION = .TRUE.
-      ENDIF
-
-C     Typecast the result back
-      R1 = DCMPLX(R1)
-      DO I=0,2
-        DP_RES(I)=DCMPLX(NINJA_RES(I))
-      ENDDO
-
-C     Make sure to deallocate the tensor of coefficients
-      IF (ALLOCATED(MP_TENSORCOEFS)) THEN
-        DEALLOCATE(MP_TENSORCOEFS)
-      ENDIF
-      IF (ALLOCATED(MP_NINJA_TENSORCOEFS)) THEN
-        DEALLOCATE(MP_NINJA_TENSORCOEFS)
-      ENDIF
-
-      RES(1)=NORMALIZATION*2.0D0*DBLE(DP_RES(0))
-      RES(2)=NORMALIZATION*2.0D0*DBLE(DP_RES(1))
-      RES(3)=NORMALIZATION*2.0D0*DBLE(DP_RES(2))
-C     WRITE(*,*) 'QP Ninja:  Loop ID',ID,' =',RES(1),RES(2),RES(3)
-      END
+C     The Ninja version installed does not support quadruple precision
+C     so that the corresponding subroutines are not output.
+C     
 
       SUBROUTINE INITNINJA()
 C     
@@ -975,8 +719,10 @@
             CALL NINJA_LOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1
      $       ,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
           ELSE
-            CALL MP_NINJA_LOOP(NLOOPLINE,MP_PL,M2L,RANK,LOOPRES(1
-     $       ,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
+            WRITE(*,*) 'ERROR: Ninja should not be called in quadruple'
+     $       //' precision since the installed version considered does'
+     $       //' not support it.'
+            STOP 9
           ENDIF
         ELSE
 C         Tensor Integral Reduction is used 

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%GOLEM_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%GOLEM_interface.f	2016-03-31 10:53:12 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%GOLEM_interface.f	2016-07-02 20:37:49 +0000
@@ -108,8 +108,9 @@
       COMMON/LOOP/ID,SQSOINDEX,R
 
       LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
+     $ ,COLLIERINIT
       COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT,GOLEMINIT,SAMURAIINIT
-     $ ,NINJAINIT
+     $ ,NINJAINIT,COLLIERINIT
 
       INTEGER NLOOPGROUPS
       PARAMETER (NLOOPGROUPS=1)

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f	2016-05-06 09:08:11 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f	2016-07-02 20:37:49 +0000
@@ -45,15 +45,17 @@
       INTEGER I, J, K
       INTEGER NLOOPCOEFS
       LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
+     $ ,COLLIERINIT
       COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT
-     $ ,NINJAINIT
+     $ ,NINJAINIT,COLLIERINIT
 
 C     This variable will be used to detect changes in the TIR library
 C      used so as to force the reset of the TIR filter.
       INTEGER LAST_LIB_USED
       DATA LAST_LIB_USED/-1/
 
-      COMPLEX*16 TIRCOEFS(0:LOOPMAXCOEFS-1,3)
+      COMPLEX*16 TIRCOEFS(0:LOOPMAXCOEFS-1,3),TIRCOEFSERRORS(0:LOOPMAXC
+     $OEFS-1,3)
       COMPLEX*16 PJCOEFS(0:LOOPMAXCOEFS-1,3)
 C     
 C     EXTERNAL FUNCTIONS
@@ -67,7 +69,6 @@
       REAL*8 LSCALE
       COMMON/CT/LSCALE,CTMODE
 
-C     The variables below are just for monitoring purposes. 
       INTEGER ID,SQSOINDEX,R
       COMMON/LOOP/ID,SQSOINDEX,R
 
@@ -184,6 +185,21 @@
      $ ,M2L,MU_R,PJCOEFS,STABLE)
 C     CONVERT TO MADLOOP CONVENTION
       CALL CONVERT_IREGI_COEFFS(RANK,PJCOEFS,TIRCOEFS)
+      CASE(7)
+C     COLLIER
+      CALL COLLIERLOOP(CTMODE,NLOOPLINE,RANK,PL,PDEN,M2L,TIRCOEFS
+     $ ,TIRCOEFSERRORS)
+C     Shift the TIR coefficients by the corresponding COLLIER error if
+C      in CTMODE 2.
+      IF (COLLIERUSEINTERNALSTABILITYTEST.AND.CTMODE.EQ.2) THEN
+C       We add here the numerical inaccuracies linearly to be
+C        conservative 
+        DO I=1,3
+          DO J=0,NLOOPCOEFS-1
+            TIRCOEFS(J,I)=TIRCOEFS(J,I)+TIRCOEFSERRORS(J,I)
+          ENDDO
+        ENDDO
+      ENDIF
       END SELECT
       DO I=1,3
         RES(I)=(0.0D0,0.0D0)
@@ -198,6 +214,8 @@
 C     WRITE(*,*) 'PJFry: Loop ID',ID,' =',RES(1),RES(2),RES(3)
 C     ELSEIF(MLReductionLib(I_LIB).EQ.3) THEN
 C     WRITE(*,*) 'Iregi: Loop ID',ID,' =',RES(1),RES(2),RES(3)
+C     ELSEIF(MLReductionLib(I_LIB).EQ.7) THEN
+C     WRITE(*,*) 'COLLIER: Loop ID',ID,' =',RES(1),RES(2),RES(3)
 C     ENDIF
       END
 
@@ -322,8 +340,9 @@
 C     
       INCLUDE 'MadLoopParams.inc'
       LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
+     $ ,COLLIERINIT
       COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT,SAMURAIINIT
-     $ ,NINJAINIT
+     $ ,NINJAINIT,COLLIERINIT
 
 C     ----------
 C     BEGIN CODE
@@ -366,9 +385,9 @@
 C     CONSTANTS
 C     
       INTEGER NLOOPLIB
-      PARAMETER (NLOOPLIB=4)
+      PARAMETER (NLOOPLIB=7)
       INTEGER QP_NLOOPLIB
-      PARAMETER (QP_NLOOPLIB=2)
+      PARAMETER (QP_NLOOPLIB=1)
       INTEGER NLOOPGROUPS
       PARAMETER (NLOOPGROUPS=1)
 C     

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa.f	2016-03-31 10:53:12 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa.f	2016-07-02 20:37:49 +0000
@@ -294,6 +294,8 @@
      $         //' an unknown reason.'
             ENDIF
           ENDIF
+          WRITE (*,'(1x,a23,3x,i3)') 'MadLoop return code   ='
+     $     ,RETURNCODE
           WRITE (*,*) '---------------------------------'
           IF (NBORNCHOSEN.EQ.0) THEN
             WRITE (*,*) 'No Born contribution satisfied the squared'

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f	2016-05-06 09:08:11 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f	2016-07-02 20:37:49 +0000
@@ -77,13 +77,13 @@
       PARAMETER (NSQUAREDSOP1=NSQUAREDSO+1)
 C     The total number of loop reduction libraries
 C     At present, there are only CutTools,PJFry++,IREGI,Golem95,Samurai
-C      and Ninja
+C     , Ninja and COLLIER
       INTEGER NLOOPLIB
-      PARAMETER (NLOOPLIB=6)
+      PARAMETER (NLOOPLIB=7)
 C     Only CutTools or possibly Ninja (if installed with qp support)
 C      provide QP
       INTEGER QP_NLOOPLIB
-      PARAMETER (QP_NLOOPLIB=2)
+      PARAMETER (QP_NLOOPLIB=1)
       INTEGER MAXSTABILITYLENGTH
       DATA MAXSTABILITYLENGTH/20/
       COMMON/STABILITY_TESTS/MAXSTABILITYLENGTH
@@ -207,6 +207,7 @@
 C     GLOBAL VARIABLES
 C     
       INCLUDE 'process_info.inc'
+      INCLUDE 'unique_id.inc'
       INCLUDE 'coupl.inc'
       INCLUDE 'mp_coupl.inc'
       INCLUDE 'MadLoopParams.inc'
@@ -236,23 +237,23 @@
       COMMON/MP_DONE/MP_DONE
 C     A FLAG TO DENOTE WHETHER THE CORRESPONDING LOOPLIBS ARE
 C      AVAILABLE OR NOT
-      LOGICAL LOOPLIBS_AVAILABLE(6)
+      LOGICAL LOOPLIBS_AVAILABLE(NLOOPLIB)
       DATA LOOPLIBS_AVAILABLE/.TRUE.,.TRUE.,.TRUE.,.TRUE.,.FALSE.
-     $ ,.TRUE./
+     $ ,.TRUE.,.TRUE./
       COMMON/LOOPLIBS_AV/ LOOPLIBS_AVAILABLE
 C     A FLAG TO DENOTE WHETHER THE CORRESPONDING DIRECTION TESTS
 C      AVAILABLE OR NOT IN THE LOOPLIBS
-      LOGICAL LOOPLIBS_DIRECTEST(6)
+      LOGICAL LOOPLIBS_DIRECTEST(NLOOPLIB)
       DATA LOOPLIBS_DIRECTEST /.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.
-     $ ,.TRUE./
+     $ ,.TRUE.,.TRUE./
 C     Specifying for which reduction tool quadruple precision is
 C      available.
 C     The index 0 is dummy and simply means that the corresponding
 C      loop_library is not available
 C     in which case neither is its quadruple precision version.
-      LOGICAL LOOPLIBS_QPAVAILABLE(0:6)
+      LOGICAL LOOPLIBS_QPAVAILABLE(0:7)
       DATA LOOPLIBS_QPAVAILABLE /.FALSE.,.TRUE.,.FALSE.,.FALSE.
-     $ ,.FALSE.,.FALSE.,.TRUE./
+     $ ,.FALSE.,.FALSE.,.FALSE.,.FALSE./
 C     PS CAN POSSIBILY BE PASSED THROUGH IMPROVE_PS BUT IS NOT
 C      MODIFIED FOR THE PURPOSE OF THE STABILITY TEST
 C     EVEN THOUGH THEY ARE PUT IN COMMON BLOCK, FOR NOW THEY ARE NOT
@@ -333,7 +334,6 @@
       INTEGER I_LIB
       DATA I_LIB/1/
       COMMON/I_LIB/I_LIB
-C     UTIL NOW, ONLY CUTTOOLS PROVIDE QP
       LOGICAL QP_TOOLS_AVAILABLE
       DATA QP_TOOLS_AVAILABLE/.FALSE./
       INTEGER INDEX_QP_TOOLS(QP_NLOOPLIB+1)
@@ -392,6 +392,21 @@
       CHARACTER(512) MLPATH
       COMMON/MLPATH/MLPATH
 
+C     This is just so that if the user disabled the computation of
+C      poles by COLLIER
+C     using the MadLoop subroutine, we don't overwrite his choice when
+C      reading the parameters
+      LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
+     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, COLLIER_IR_POLE_COMPU
+     $TATION_CHOICE
+      DATA  FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION
+     $ ,FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION/.FALSE.,.FALSE./
+      COMMON/COLLIERPOLESFORCEDCHOICE/FORCED_CHOICE_OF_COLLIER_UV_POLE_
+     $COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE,COLLIER_IR_POLE_COMPUTATION_
+     $CHOICE
+
 C     This variable controls the general initialization which is
 C      *common* between all MadLoop SubProcesses.
 C     For example setting the MadLoopPath or reading the ML runtime
@@ -413,6 +428,7 @@
       LOGICAL FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION
       DATA FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION/.FALSE.,.FALSE./
       COMMON/FPE_IN_REDUCTION/FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION
+
 C     ----------
 C     BEGIN CODE
 C     ----------
@@ -424,6 +440,12 @@
         CALL SETMADLOOPPATH(TMP)
         CALL JOINPATH(MLPATH,PARAMFNAME,PARAMFN)
         CALL MADLOOPPARAMREADER(PARAMFN,.TRUE.)
+        IF (FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION) THEN
+          COLLIERCOMPUTEUVPOLES = COLLIER_UV_POLE_COMPUTATION_CHOICE
+        ENDIF
+        IF (FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION) THEN
+          COLLIERCOMPUTEIRPOLES = COLLIER_IR_POLE_COMPUTATION_CHOICE
+        ENDIF
         IF (FORBID_HEL_DOUBLECHECK) THEN
           DOUBLECHECKHELICITYFILTER = .FALSE.
         ENDIF
@@ -766,6 +788,9 @@
         CALL CLEAR_CACHES()
       ENDIF
 
+C     Now make sure to turn on the global COLLIER cache if applicable
+      CALL SET_COLLIER_GLOBAL_CACHE(.TRUE.)
+
       IF (IMPROVEPSPOINT.GE.0) THEN
 C       Make the input PS more precise (exact onshell and
 C        energy-momentum conservation)
@@ -1570,12 +1595,22 @@
       IF (AUTOMATIC_CACHE_CLEARING) THEN
         CALL CLEAR_CACHES()
       ENDIF
+
+C     Now make sure to turn off the global COLLIER cache if applicable
+      CALL SET_COLLIER_GLOBAL_CACHE(.FALSE.)
+
       END
 
       SUBROUTINE CLEAR_CACHES()
-C     Clears all the caches used at some point in MadLoop
+C     
+C     This routine can be called directly from the user if
+C     AUTOMATIC_CACHE_CLEARING is set to False. It must then be called
+C      after
+C     ech event
+C     
       CALL CLEAR_TIR_CACHE()
       CALL NINJA_CLEAR_INTEGRAL_CACHE()
+      CALL CLEAR_COLLIER_CACHE()
       END
 
 C     --=========================================--
@@ -1716,6 +1751,13 @@
       REAL*8 ACCURACIES(3)
       REAL*8 LIST(MAXSTABILITYLENGTH)
 
+C     
+C     GLOBAL VARIABLES
+C     
+      INTEGER I_LIB
+      COMMON/I_LIB/I_LIB
+      INCLUDE 'MadLoopParams.inc'
+
 C     ----------
 C     BEGIN CODE
 C     ----------
@@ -1765,6 +1807,16 @@
           ACC(K)  = ACC(K) / ( ABS(AVG) / 3.0D0)
         ENDIF
 
+C       When using COLLIER with the internal stability test, the first
+C        evaluation is typically more reliable so we do not want to
+C        use the average but rather the first evaluation.
+        IF (MLREDUCTIONLIB(I_LIB).EQ.7.AND.COLLIERUSEINTERNALSTABILITYT
+     $EST) THEN
+          DO I=1,3
+            ESTIMATE(I,K) = FULLLIST(I,K,1)
+          ENDDO
+        ENDIF
+
       ENDDO
 
       END
@@ -2119,6 +2171,60 @@
 C     Definition of additional access routines
 C     --=========================================--
 
+      SUBROUTINE COLLIER_COMPUTE_UV_POLES(ONOFF)
+C     
+C     This function can be called by the MadLoop user so as to chose
+C      to have COLLIER
+C     compute the UV pole or not (it costs more time).
+C     
+      LOGICAL ONOFF
+
+      INCLUDE 'MadLoopParams.inc'
+
+      LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
+     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, COLLIER_IR_POLE_COMPU
+     $TATION_CHOICE
+      COMMON/COLLIERPOLESFORCEDCHOICE/FORCED_CHOICE_OF_COLLIER_UV_POLE_
+     $COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE,COLLIER_IR_POLE_COMPUTATION_
+     $CHOICE
+
+      COLLIERCOMPUTEUVPOLES                        = ONOFF
+C     This is just so that if we read the param again, we don't
+C      overwrite the choice made here
+      FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION = .TRUE.
+      COLLIER_UV_POLE_COMPUTATION_CHOICE           = ONOFF
+
+      END SUBROUTINE
+
+      SUBROUTINE COLLIER_COMPUTE_IR_POLES(ONOFF)
+C     
+C     This function can be called by the MadLoop user so as to chose
+C      to have COLLIER
+C     compute the IR pole or not (it costs more time).
+C     
+      LOGICAL ONOFF
+
+      INCLUDE 'MadLoopParams.inc'
+
+      LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION,
+     $  FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+      LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, COLLIER_IR_POLE_COMPU
+     $TATION_CHOICE
+      COMMON/COLLIERPOLESFORCEDCHOICE/FORCED_CHOICE_OF_COLLIER_UV_POLE_
+     $COMPUTATION, FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION
+     $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE,COLLIER_IR_POLE_COMPUTATION_
+     $CHOICE
+
+      COLLIERCOMPUTEIRPOLES         = ONOFF
+C     This is just so that if we read the param again, we don't
+C      overwrite the choice made here
+      FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION = .TRUE.
+      COLLIER_IR_POLE_COMPUTATION_CHOICE           = ONOFF
+
+      END SUBROUTINE
+
       SUBROUTINE FORCE_STABILITY_CHECK(ONOFF)
 C     
 C     This function can be called by the MadLoop user so as to always

=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%mp_compute_loop_coefs.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%mp_compute_loop_coefs.f	2016-04-06 18:45:52 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%mp_compute_loop_coefs.f	2016-07-02 20:37:49 +0000
@@ -250,9 +250,15 @@
             NHEL(I)=HELC(I,H)
           ENDDO
 
-          MP_UVCT_REQ_SO_DONE=.FALSE.
-          MP_CT_REQ_SO_DONE=.FALSE.
-          MP_LOOP_REQ_SO_DONE=.FALSE.
+          IF (COMPUTE_INTEGRAND_IN_QP) THEN
+            MP_UVCT_REQ_SO_DONE=.FALSE.
+            MP_CT_REQ_SO_DONE=.FALSE.
+            MP_LOOP_REQ_SO_DONE=.FALSE.
+          ELSE
+            UVCT_REQ_SO_DONE=.FALSE.
+            CT_REQ_SO_DONE=.FALSE.
+            LOOP_REQ_SO_DONE=.FALSE.
+          ENDIF
 
           IF (.NOT.CHECKPHASE.AND.HELDOUBLECHECKED.AND.HELPICKED.EQ.-1)
      $      THEN
@@ -264,8 +270,10 @@
 
           IF (COMPUTE_INTEGRAND_IN_QP) THEN
             CALL MP_HELAS_CALLS_AMPB_1(MP_P,NHEL,H,IC)
+            CONTINUE
           ELSE
             CALL HELAS_CALLS_AMPB_1(P,NHEL,H,IC)
+            CONTINUE
           ENDIF
 
  2000     CONTINUE
@@ -273,8 +281,10 @@
 
           IF (COMPUTE_INTEGRAND_IN_QP) THEN
 
+            CONTINUE
           ELSE
 
+            CONTINUE
           ENDIF
 
           IF (.NOT.COMPUTE_INTEGRAND_IN_QP) THEN

=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%unique_id.inc'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%unique_id.inc	1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%unique_id.inc	2016-07-02 20:37:49 +0000
@@ -0,0 +1,2 @@
+      integer UNIQUE_ID
+      parameter(UNIQUE_ID=1)
\ No newline at end of file

=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%COLLIER_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%COLLIER_interface.f	1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%COLLIER_interface.f	2016-07-02 20:37:49 +0000
@@ -0,0 +1,732 @@
+      SUBROUTINE COLLIERLOOP(CTMODE, NLOOPLINE,RANK,PL,PDEN,M2L
+     $ ,TIRCOEFS,TIRCOEFSERRORS)
+C     
+C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
+C     By the MadGraph5_aMC@NLO Development Team
+C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
+C     
+C     Interface between MG5 and COLLIER.
+C     It supports any rank when 1 < NLOOPLINE < 7.
+C     
+C     Process: u d~ > w+ WEIGHTED<=2 QED<=1 [ all = QCD ]
+C     Process: c s~ > w+ WEIGHTED<=2 QED<=1 [ all = QCD ]
+C     
+C     
+C     MODULES
+C     
+      USE COLLIER
+      IMPLICIT NONE
+C     
+C     CONSTANTS 
+C     
+      LOGICAL CHECKPCONSERVATION
+      PARAMETER (CHECKPCONSERVATION=.TRUE.)
+      REAL*8 PIVALUE
+      PARAMETER(PIVALUE=3.14159265358979323846D0)
+      INTEGER MAXRANK
+      PARAMETER (MAXRANK=2)
+      INTEGER NLOOPGROUPS
+      PARAMETER (NLOOPGROUPS=1)
+      INCLUDE 'loop_max_coefs.inc'
+C     
+C     ARGUMENTS 
+C     
+      INTEGER NLOOPLINE, RANK, CTMODE
+      REAL*8 PL(0:3,NLOOPLINE)
+      REAL*8 PDEN(0:3,NLOOPLINE-1)
+      COMPLEX*16 M2L(NLOOPLINE)
+      COMPLEX*16 RES(3)
+      COMPLEX*16 TIRCOEFS(0:LOOPMAXCOEFS-1,3)
+      COMPLEX*16 TIRCOEFSERRORS(0:LOOPMAXCOEFS-1,3)
+
+C     
+C     LOCAL VARIABLES 
+C     
+      INTEGER N, I, J, K, L
+      INTEGER C0,C1,C2,C3
+      INTEGER N_CACHES
+      INTEGER CURR_RANK, SGN
+      INTEGER CURR_INDEX
+      INTEGER CURR_MAXCOEF
+      REAL*8 RBUFF(0:3)
+
+      INTEGER COEFMAP_ZERO(0:LOOPMAXCOEFS-1)
+      INTEGER COEFMAP_ONE(0:LOOPMAXCOEFS-1)
+      INTEGER COEFMAP_TWO(0:LOOPMAXCOEFS-1)
+      INTEGER COEFMAP_THREE(0:LOOPMAXCOEFS-1)
+      DATA (COEFMAP_ZERO(I),I=  0,  9) / 0, 1, 0, 0, 0, 2, 1, 0, 1, 0/
+      DATA (COEFMAP_ZERO(I),I= 10, 14) / 0, 1, 0, 0, 0/
+      DATA (COEFMAP_ONE(I),I=  0,  9) / 0, 0, 1, 0, 0, 0, 1, 2, 0, 1/
+      DATA (COEFMAP_ONE(I),I= 10, 14) / 0, 0, 1, 0, 0/
+      DATA (COEFMAP_TWO(I),I=  0,  9) / 0, 0, 0, 1, 0, 0, 0, 0, 1, 1/
+      DATA (COEFMAP_TWO(I),I= 10, 14) / 2, 0, 0, 1, 0/
+      DATA (COEFMAP_THREE(I),I=  0,  9) / 0, 0, 0, 0, 1, 0, 0, 0, 0, 0/
+      DATA (COEFMAP_THREE(I),I= 10, 14) / 0, 1, 1, 1, 2/
+
+      REAL*8 REF_NORMALIZATION
+
+      DOUBLE COMPLEX M2LCOLLIER(0:NLOOPLINE-1)
+      DOUBLE COMPLEX MOMVEC(0:3,NLOOPLINE-1)
+      DOUBLE COMPLEX MOMINV((NLOOPLINE*(NLOOPLINE-1))/2)
+
+      DOUBLE COMPLEX TNTEN(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE COMPLEX TNTENUV(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE PRECISION TNTENERR(0:RANK)
+
+      REAL*8 MAXCOEFFORRANK(0:RANK)
+
+C     These quantities are for the pole evaluation
+      DOUBLE COMPLEX TNTEN_UV(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE COMPLEX TNTENUV_UV(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE PRECISION TNTENERR_UV(0:RANK)
+      DOUBLE COMPLEX TNTEN_IR1(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE COMPLEX TNTENUV_IR1(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE PRECISION TNTENERR_IR1(0:RANK)
+      DOUBLE COMPLEX TNTEN_IR2(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE COMPLEX TNTENUV_IR2(0:RANK,0:RANK,0:RANK,0:RANK)
+      DOUBLE PRECISION TNTENERR_IR2(0:RANK)
+C     
+C     GLOBAL VARIABLES
+C     
+      INCLUDE 'coupl.inc'
+      INCLUDE 'MadLoopParams.inc'
+      INCLUDE 'unique_id.inc'
+      INCLUDE 'global_specs.inc'
+
+      INTEGER ID,SQSOINDEX,R
+      COMMON/LOOP/ID,SQSOINDEX,R
+
+C     
+C     These global variables cache the coefficients already computed
+C      by COLLIER
+C     so that they can be reused. In particular, in CTMODE 2, the
+C      cached ERROR 
+C     quantities for eahc rank will be used to make use of COLLIER's
+C      internal
+C     assessment of the numerical accuracy.
+C     
+      LOGICAL LOOP_ID_DONE(NLOOPGROUPS)
+      DATA LOOP_ID_DONE/NLOOPGROUPS*.FALSE./
+      COMPLEX*16 TIR_COEFS_DIRECT_MODE(0:LOOPMAXCOEFS-1,3,NLOOPGROUPS)
+      REAL*8 TIR_COEFS_ERROR(0:MAXRANK,NLOOPGROUPS)
+      COMMON/COLLIER_TIR_COEFS/TIR_COEFS_DIRECT_MODE,TIR_COEFS_ERROR
+     $ ,LOOP_ID_DONE
+
+      INTEGER COLLIER_CACHE_ACTIVE, NCALLS_IN_CACHE(4),
+     $  NCOLLIERCALLS(4)
+      LOGICAL MUST_INIT_EVENT
+      COMMON/COLLIER_CACHE_STATUS/COLLIER_CACHE_ACTIVE,NCALLS_IN_CACHE
+     $ , NCOLLIERCALLS,MUST_INIT_EVENT
+
+      LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT,
+     $  COLLIERINIT
+      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT,GOLEMINIT,SAMURAIINIT
+     $ ,NINJAINIT, COLLIERINIT
+
+C     ----------
+C     BEGIN CODE
+C     ----------
+
+      IF (COLLIERUSECACHEFORPOLES) THEN
+        N_CACHES = 4
+      ELSE
+        N_CACHES =1
+      ENDIF
+
+C     Initialize COLLIER if needed 
+      IF (COLLIERINIT) THEN
+        COLLIERINIT=.FALSE.
+        CALL INITCOLLIER()
+      ENDIF
+
+C     Initialize the event if it is the first time collier is called
+C      on this PS point
+      IF(MUST_INIT_EVENT) THEN
+        MUST_INIT_EVENT = .FALSE.
+        IF (COLLIERGLOBALCACHE.EQ.0) THEN
+          CALL INITEVENT_CLL()
+        ELSE
+          DO I=1,N_CACHES
+C           Record how many events where put in the cache. On the
+C            first PS point.
+            IF (NCALLS_IN_CACHE(I).EQ.-1.AND.NCOLLIERCALLS(I).GT.0)
+     $        THEN
+              NCALLS_IN_CACHE(I) = NCOLLIERCALLS(I)
+            ENDIF
+          ENDDO
+          DO I=1,N_CACHES
+C           Now apply a safety check that our last event had as many
+C            calls as the cache is setup for.
+C           The only case for now when it can be half of the calls
+C            when we are doing the true loop-direction test with also
+C            the computation of a rotated PS point (which is computed
+C            for one mode only).
+            IF (NCALLS_IN_CACHE(I).NE.-1.AND..NOT. ( NCOLLIERCALLS(I)
+     $       .EQ.NCALLS_IN_CACHE(I).OR. ( CTMODERUN.EQ.
+     $       -1.AND..NOT.COLLIERUSEINTERNALSTABILITYTEST.AND.NROTATIONS
+     $_DP.GT.0.AND.MOD(NCALLS_IN_CACHE(I),2).EQ.0.AND.(NCALLS_IN_CACHE(
+     $I)/2).EQ.NCOLLIERCALLS(I) ) ) ) THEN
+              WRITE(*,*) 'WARNING: A consistency check in MadLoop'
+     $         //' failed and, for safety, forced MadLoop to'
+     $         //' reinitialize the global cache of COLLIER. Report'
+     $         //' this to MadLoop authors. The problematic cache was'
+     $         //' number ',I
+              IF (COLLIERGLOBALCACHE.EQ.-1) THEN
+                CALL INITCACHESYSTEM_CLL(N_CACHES*NPROCS,MAXNEXTERNAL)
+              ELSE
+                CALL INITCACHESYSTEM_CLL(N_CACHES*NPROCS
+     $           ,COLLIERGLOBALCACHE)
+              ENDIF
+C             Make sure all caches are switched off at first.
+              CALL SWITCHOFFCACHESYSTEM_CLL()
+C             Reset the cache design property
+              NCALLS_IN_CACHE(:) = -1
+              NCOLLIERCALLS(:)   = 0
+              IF (COLLIER_CACHE_ACTIVE.EQ.1) THEN
+                CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*N_CACHES+1)
+              ENDIF
+C             No need to check the other caches since we already had
+C              to reset here.
+              EXIT
+            ENDIF
+          ENDDO
+          CALL INITEVENT_CLL((UNIQUE_ID-1)*N_CACHES+1)
+          NCOLLIERCALLS(1) = 0
+          IF(COLLIERCOMPUTEUVPOLES.AND.COLLIERUSECACHEFORPOLES) THEN
+            CALL INITEVENT_CLL((UNIQUE_ID-1)*N_CACHES+2)
+            NCOLLIERCALLS(2) = 0
+          ENDIF
+          IF(COLLIERCOMPUTEIRPOLES.AND.COLLIERUSECACHEFORPOLES) THEN
+            CALL INITEVENT_CLL((UNIQUE_ID-1)*N_CACHES+3)
+            NCOLLIERCALLS(3) = 0
+            CALL INITEVENT_CLL((UNIQUE_ID-1)*N_CACHES+4)
+            NCOLLIERCALLS(4) = 0
+          ENDIF
+        ENDIF
+        CALL SETDELTAIR_CLL(0.0D0,(PIVALUE**2)/6.0D0)
+        CALL SETDELTAUV_CLL(0.0D0)
+        CALL SETMUIR2_CLL(MU_R**2)
+        CALL SETMUUV2_CLL(MU_R**2)
+      ENDIF
+
+C     Now really start the reduction with COLLIER
+
+C     Number of coefficients for the current rank
+      CURR_MAXCOEF = 0
+      DO I=0,RANK
+        CURR_MAXCOEF=CURR_MAXCOEF+(3+I)*(2+I)*(1+I)/6
+      ENDDO
+
+      IF (CTMODE.NE.1.AND.CTMODE.NE.2) THEN
+        WRITE(*,*) 'ERROR: COLLIER only support the computational mode'
+     $   //' CTMODE equal to 1 or 2, not',CTMODE
+        STOP 'Incorrect computational mode specified to the COLLIER'
+     $   //' MG5aMC interface.'
+      ENDIF
+
+      DO I=0,CURR_MAXCOEF-1
+        DO K=1,3
+          TIRCOEFSERRORS(I,K)=DCMPLX(0.0D0,0.0D0)
+        ENDDO
+      ENDDO
+
+      IF (COLLIERUSEINTERNALSTABILITYTEST) THEN
+C       Use MADLOOP internal cache dedicated to COLLIER that emulates
+C        the CTMODE 2
+        IF (LOOP_ID_DONE(ID)) THEN
+          DO I=0,CURR_MAXCOEF-1
+            DO K=1,3
+              TIRCOEFS(I,K)=TIR_COEFS_DIRECT_MODE(I,K,ID)
+            ENDDO
+          ENDDO
+          DO I=0,RANK
+            MAXCOEFFORRANK(I) = 0.0D0
+          ENDDO
+          DO I=0,CURR_MAXCOEF-1
+            CURR_RANK = COEFMAP_ZERO(I)+COEFMAP_ONE(I)+COEFMAP_TWO(I)
+     $       +COEFMAP_THREE(I)
+            MAXCOEFFORRANK(CURR_RANK)=MAX(MAXCOEFFORRANK(CURR_RANK)
+     $       ,ABS(TIR_COEFS_DIRECT_MODE(I,1,ID)))
+          ENDDO
+          DO I=0,CURR_MAXCOEF-1
+            CURR_RANK = COEFMAP_ZERO(I)+COEFMAP_ONE(I)+COEFMAP_TWO(I)
+     $       +COEFMAP_THREE(I)
+            IF (MAXCOEFFORRANK(CURR_RANK).NE.0.0D0) THEN
+              DO K=1,3
+C               The expression below is like taking the absolute value
+C                when summing errors linearly
+C               TIRCOEFSERRORS(I,K)=(TIR_COEFS_ERROR(CURR_RANK,ID)/MaxC
+C               oefForRank(CURR_RANK))*DCMPLX( ABS(DBLE(TIRCOEFS(I,K)))
+C               ,ABS(DIMAG(TIRCOEFS(I,K))) )
+C               But empirically, I observed that retaining the
+C                original complex phase leads to slightly more
+C                accurate estimates
+                TIRCOEFSERRORS(I,K)=(TIR_COEFS_ERROR(CURR_RANK,ID)
+     $           /MAXCOEFFORRANK(CURR_RANK))*TIRCOEFS(I,K)
+              ENDDO
+            ENDIF
+          ENDDO
+          RETURN
+        ENDIF
+      ELSE
+C       Apply the loop-direction switching here.
+        CALL SWITCH_ORDER(CTMODE,NLOOPLINE,PL,PDEN,M2L)
+      ENDIF
+
+C     Make sure masses are complex 
+      DO I=1,NLOOPLINE
+        M2LCOLLIER(I-1)=DCMPLX(M2L(I))
+      ENDDO
+
+C     Now convert the loop offset momenta to COLLIER conventions
+      DO I=0,3
+        DO J=1,NLOOPLINE-1
+          MOMVEC(I,J)=DCMPLX(PDEN(I,J),0.0D0)
+        ENDDO
+      ENDDO
+
+C     This first do loop spans over 'N' in '\foreach_N'
+C     //' \foreach_i(k_i+k_{i+1}+..+k_{i+N})^2'
+      CURR_INDEX = 0
+      DO N=0,NLOOPLINE-1
+C       We stop whenever we reached the target number of invariants
+        IF (CURR_INDEX.GE.((NLOOPLINE*(NLOOPLINE-1))/2)) THEN
+          EXIT
+        ENDIF
+C       This do loop spans over 'i' in the expression above.
+        DO I=1,NLOOPLINE
+          IF (CURR_INDEX.GE.((NLOOPLINE*(NLOOPLINE-1))/2)) THEN
+            EXIT
+          ENDIF
+
+          CURR_INDEX = CURR_INDEX+1
+          RBUFF(:) = 0.0D0
+          DO J=I,I+N
+            RBUFF(:) = RBUFF(:) + PL(:,MOD(J-1,NLOOPLINE)+1)
+          ENDDO
+          MOMINV(CURR_INDEX) = DCMPLX(RBUFF(0)**2-RBUFF(1)**2-RBUFF(2)
+     $     **2-RBUFF(3)**2,0.0D0)
+
+C         Now regularize the onshell behavior of the kinematic
+C          invarients
+C         All loop masses are tested, but that might be a bit too
+C          inclusive.
+          DO K=1,NLOOPLINE
+            IF(ABS(M2L(K)).NE.0.0D0) THEN
+              IF(ABS((MOMINV(CURR_INDEX)-M2L(K))/M2L(K)).LT.OSTHRES)
+     $          THEN
+                MOMINV(CURR_INDEX)=DCMPLX(M2L(K))
+              ENDIF
+            ENDIF
+          ENDDO
+C         For the massless onshell-case, we base the threshold only on
+C          the energy component
+          REF_NORMALIZATION=0.0D0
+          DO L=0,0
+            REF_NORMALIZATION = REF_NORMALIZATION + ABS(RBUFF(L))
+          ENDDO
+          REF_NORMALIZATION = (REF_NORMALIZATION/(N+1))**2
+          IF(REF_NORMALIZATION.NE.0.0D0)THEN
+            IF(ABS(MOMINV(CURR_INDEX)/REF_NORMALIZATION).LT.OSTHRES)
+     $       THEN
+              MOMINV(CURR_INDEX)=DCMPLX(0.0D0,0.0D0)
+            ENDIF
+          ENDIF
+
+        ENDDO
+      ENDDO
+
+C     We can now call COLLIER
+      IF (NLOOPLINE.NE.1) THEN
+        CALL TNTEN_CLL(TNTEN, TNTENUV, MOMVEC, MOMINV, M2LCOLLIER,
+     $    NLOOPLINE, RANK, TNTENERR)
+      ELSE
+        CALL TNTEN_CLL(TNTEN, TNTENUV, M2LCOLLIER, NLOOPLINE, RANK,
+     $    TNTENERR)
+      ENDIF
+      IF (COLLIER_CACHE_ACTIVE.EQ.1) THEN
+        NCOLLIERCALLS(1) = NCOLLIERCALLS(1)+1
+      ENDIF
+
+C     Now compute the UV poles if asked for
+      IF (COLLIERCOMPUTEUVPOLES) THEN
+        IF(COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          CALL SWITCHOFFCACHE_CLL((UNIQUE_ID-1)*N_CACHES+1)
+          IF(COLLIERUSECACHEFORPOLES) THEN
+            CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*N_CACHES+2)
+          ENDIF
+        ENDIF
+        CALL SETDELTAUV_CLL(1.0D0)
+        IF (NLOOPLINE.NE.1) THEN
+          CALL TNTEN_CLL(TNTEN_UV, TNTENUV_UV, MOMVEC, MOMINV,
+     $      M2LCOLLIER, NLOOPLINE, RANK, TNTENERR_UV)
+        ELSE
+          CALL TNTEN_CLL(TNTEN_UV, TNTENUV_UV, M2LCOLLIER, NLOOPLINE,
+     $      RANK, TNTENERR_UV)
+        ENDIF
+        IF (COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          NCOLLIERCALLS(2) = NCOLLIERCALLS(2)+1
+        ENDIF
+        IF(COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          IF(COLLIERUSECACHEFORPOLES) THEN
+            CALL SWITCHOFFCACHE_CLL((UNIQUE_ID-1)*N_CACHES+2)
+          ENDIF
+          CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*N_CACHES+1)
+        ENDIF
+        CALL SETDELTAUV_CLL(0.0D0)
+      ENDIF
+
+C     Now compute the IR poles if asked for
+      IF (COLLIERCOMPUTEIRPOLES) THEN
+        IF(COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          CALL SWITCHOFFCACHE_CLL((UNIQUE_ID-1)*N_CACHES+1)
+          IF(COLLIERUSECACHEFORPOLES) THEN
+            CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*N_CACHES+3)
+          ENDIF
+        ENDIF
+        CALL SETDELTAIR_CLL(1.0D0,(PIVALUE**2)/6.0D0)
+        IF (NLOOPLINE.NE.1) THEN
+          CALL TNTEN_CLL(TNTEN_IR1, TNTENUV_IR1, MOMVEC, MOMINV,
+     $      M2LCOLLIER, NLOOPLINE, RANK, TNTENERR_IR1)
+        ELSE
+          CALL TNTEN_CLL(TNTEN_IR1, TNTENUV_IR1, M2LCOLLIER, NLOOPLINE
+     $     , RANK, TNTENERR_IR1)
+        ENDIF
+        IF (COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          NCOLLIERCALLS(3) = NCOLLIERCALLS(3)+1
+        ENDIF
+        IF(COLLIER_CACHE_ACTIVE.EQ.1.AND.COLLIERUSECACHEFORPOLES) THEN
+          CALL SWITCHOFFCACHE_CLL((UNIQUE_ID-1)*N_CACHES+3)
+          CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*N_CACHES+4)
+        ENDIF
+        CALL SETDELTAIR_CLL(0.0D0,1.0D0+(PIVALUE**2)/6.0D0)
+        IF (NLOOPLINE.NE.1) THEN
+          CALL TNTEN_CLL(TNTEN_IR2, TNTENUV_IR2, MOMVEC, MOMINV,
+     $      M2LCOLLIER, NLOOPLINE, RANK, TNTENERR_IR2)
+        ELSE
+          CALL TNTEN_CLL(TNTEN_IR2, TNTENUV_IR2, M2LCOLLIER, NLOOPLINE
+     $     , RANK, TNTENERR_IR2)
+        ENDIF
+        IF (COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          NCOLLIERCALLS(4) = NCOLLIERCALLS(4)+1
+        ENDIF
+        IF(COLLIER_CACHE_ACTIVE.EQ.1) THEN
+          IF(COLLIERUSECACHEFORPOLES) THEN
+            CALL SWITCHOFFCACHE_CLL((UNIQUE_ID-1)*N_CACHES+4)
+          ENDIF
+          CALL SWITCHONCACHE_CLL((UNIQUE_ID-1)*N_CACHES+1)
+        ENDIF
+        CALL SETDELTAIR_CLL(0.0D0,(PIVALUE**2)/6.0D0)
+      ENDIF
+
+      DO I=0,(CURR_MAXCOEF-1)
+        C0 = COEFMAP_ZERO(I)
+        C1 = COEFMAP_ONE(I)
+        C2 = COEFMAP_TWO(I)
+        C3 = COEFMAP_THREE(I)
+        CURR_RANK = C0+C1+C2+C3
+C       Because we must set q -> -q, we apply a minus sign to coefs of
+C        odd rank
+        IF (MOD(CURR_RANK,2).EQ.1) THEN
+          SGN = -1
+        ELSE
+          SGN = 1
+        ENDIF
+        TIRCOEFS(I,1) = SGN*TNTEN(C0,C1,C2,C3)
+        IF (COLLIERCOMPUTEUVPOLES) THEN
+          TIRCOEFS(I,2) = SGN*( TNTEN_UV(C0,C1,C2,C3)-TNTEN(C0,C1,C2
+     $     ,C3) )
+        ELSE
+          TIRCOEFS(I,2) = DCMPLX(0.0D0,0.0D0)
+        ENDIF
+        IF (COLLIERCOMPUTEIRPOLES) THEN
+          TIRCOEFS(I,2) = TIRCOEFS(I,2) + SGN*( TNTEN_IR1(C0,C1,C2,C3)
+     $     -TNTEN(C0,C1,C2,C3) )
+          TIRCOEFS(I,3) = SGN*( TNTEN_IR2(C0,C1,C2,C3)-TNTEN(C0,C1,C2
+     $     ,C3) )
+        ELSE
+          TIRCOEFS(I,3) = DCMPLX(0.0D0,0.0D0)
+        ENDIF
+      ENDDO
+
+      IF (COLLIERUSEINTERNALSTABILITYTEST) THEN
+C       Finish by caching the coefficients and error just computed
+        LOOP_ID_DONE(ID) = .TRUE.
+        DO J=0,CURR_MAXCOEF-1
+          DO K=1,3
+            TIR_COEFS_DIRECT_MODE(J,K,ID) = TIRCOEFS(J,K)
+          ENDDO
+        ENDDO
+        DO J=0,RANK
+          TIR_COEFS_ERROR(J,ID)=TNTENERR(J)
+        ENDDO
+
+C       Now compute the errors on each coefficient
+        DO I=0,RANK
+          MAXCOEFFORRANK(I) = 0.0D0
+        ENDDO
+        DO I=0,CURR_MAXCOEF-1
+          CURR_RANK = COEFMAP_ZERO(I)+COEFMAP_ONE(I)+COEFMAP_TWO(I)
+     $     +COEFMAP_THREE(I)
+          MAXCOEFFORRANK(CURR_RANK)=MAX(MAXCOEFFORRANK(CURR_RANK)
+     $     ,ABS(TIRCOEFS(I,1)))
+        ENDDO
+        DO I=0,CURR_MAXCOEF-1
+          CURR_RANK = COEFMAP_ZERO(I)+COEFMAP_ONE(I)+COEFMAP_TWO(I)
+     $     +COEFMAP_THREE(I)
+          IF (MAXCOEFFORRANK(CURR_RANK).NE.0.0D0) THEN
+            DO K=1,3
+C             The expression below is like taking the absolute value
+C              when summing errors linearly
+C             TIRCOEFSERRORS(I,K)=(TNtenerr(CURR_RANK)/MaxCoefForRank(C
+C             URR_RANK))*DCMPLX( ABS(DBLE(TIRCOEFS(I,K))),ABS(DIMAG(TIR
+C             COEFS(I,K))) )			
+C             But empirically, I observed that retaining the original
+C              complex phase leads to slightly more accurate estimates
+              TIRCOEFSERRORS(I,K)=(TNTENERR(CURR_RANK)
+     $         /MAXCOEFFORRANK(CURR_RANK))*TIRCOEFS(I,K)
+            ENDDO
+          ENDIF
+        ENDDO
+
+      ENDIF
+
+      END SUBROUTINE COLLIERLOOP
+
+      SUBROUTINE CLEAR_COLLIER_CACHE()
+
+      USE COLLIER
+
+      INCLUDE 'loop_max_coefs.inc'
+      INTEGER NLOOPGROUPS
+      PARAMETER (NLOOPGROUPS=1)
+      INTEGER MAXRANK
+      PARAMETER (MAXRANK=2)
+
+      INTEGER I,J,K
+
+      INCLUDE 'MadLoopParams.inc'
+
+      LOGICAL LOOP_ID_DONE(NLOOPGROUPS)
+      COMPLEX*16 TIR_COEFS_DIRECT_MODE(0:LOOPMAXCOEFS-1,3,NLOOPGROUPS)
+      COMPLEX*16 TIR_COEFS_ERROR(0:MAXRANK,NLOOPGROUPS)
+      COMMON/COLLIER_TIR_COEFS/TIR_COEFS_DIRECT_MODE,TIR_COEFS_ERROR
+     $ ,LOOP_ID_DONE
+
+      INTEGER COLLIER_CACHE_ACTIVE, NCALLS_IN_CACHE(4),
+     $  NCOLLIERCALLS(4)
+      LOGICAL MUST_INIT_EVENT
+      COMMON/COLLIER_CACHE_STATUS/COLLIER_CACHE_ACTIVE,NCALLS_IN_CACHE
+     $ , NCOLLIERCALLS,MUST_INIT_EVENT
+
+C     Make sure that next time the COLLIER Subroutine is called it
+C      will call the subroutine initEvent of Collier.
+      MUST_INIT_EVENT = .TRUE.
+
+C     Reinitialize the ML cache for COLLIER
+      DO I=1,NLOOPGROUPS
+        LOOP_ID_DONE(I) = .FALSE.
+        DO J=0,LOOPMAXCOEFS-1
+          DO K=1,3
+            TIR_COEFS_DIRECT_MODE(J,K,I) = DCMPLX(0.0D0,0.0D0)
+          ENDDO
+        ENDDO
+        DO J=0,MAXRANK
+          TIR_COEFS_ERROR(J,I)=0.0D0
+        ENDDO
+      ENDDO
+
+      END
+
+      SUBROUTINE SET_COLLIER_GLOBAL_CACHE(ONOFF)
+C     
+C     This routine is used by loop_matrix.f to turn on the global
+C     cache of COLLIER when it the main SLOOP subroutine starts and
+C     turn it off when it ends.
+C     However several checks are performed to verify that it is safe
+C     to turn it on and to reinitialize it if necessary.
+C     
+C     MODULES
+C     
+      USE COLLIER
+      IMPLICIT NONE
+C     
+C     ARGUMENTS
+C     
+      LOGICAL ONOFF
+C     
+C     LOCAL VARIABLES
+C     
+      LOGICAL NEED_REINITIALIZATION
+      INTEGER N_CACHES
+C     
+C     GLOBAL VARIABLES
+C     
+C     
+C     This common blocks saves the relevant ML parameters when
+C      activating the
+C     global cache of COLLIER so that we know when we must
+C      reinitialize it.
+C     COLLIER_CACHE_ACTIVE is set to -1 when it has never been turned
+C      on yet and
+C     to 1 for 'Active' and 0 for 'Inactive'.
+C     The integer NCALLS_IN_CACHE saves how many calls the cache is
+C      setup for, for each of the up to four caches.
+C     When it is the first PS points, it is set to -1.
+      INTEGER COLLIER_CACHE_ACTIVE, NCALLS_IN_CACHE(4),
+     $  NCOLLIERCALLS(4)
+      DATA COLLIER_CACHE_ACTIVE/-1/
+      DATA NCALLS_IN_CACHE/-1,-1,-1,-1/
+      DATA NCOLLIERCALLS/0,0,0,0/
+C     This is a flag to tell the COLLIER subroutine that it must init
+C      the event when called.
+      LOGICAL MUST_INIT_EVENT
+      DATA MUST_INIT_EVENT/.TRUE./
+      COMMON/COLLIER_CACHE_STATUS/COLLIER_CACHE_ACTIVE,
+     $  NCALLS_IN_CACHE, NCOLLIERCALLS,MUST_INIT_EVENT
+
+      LOGICAL COLLIERUSEINTERNALSTABILITYTEST_BU
+      INTEGER USERHEL_BU, SQSO_TARGET_BU, COLLIERMODE_BU,CTMODERUN_BU
+      COMMON/COLLIER_CACHE_RELEVANT_PARAMS/USERHEL_BU,SQSO_TARGET_BU
+     $ ,COLLIERMODE_BU,CTMODERUN_BU,COLLIERUSEINTERNALSTABILITYTEST_BU
+
+C     The common blocks below are to retrieve the necessary
+C      information about
+C     MadLoop running mode and store it in the sCOLLIER_CACHE_RELEVANT_