Merge lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4

[marco zaro]

A new generation mode at NLO has been implemented, in order to spread the generation on all available CPUs and to use less memory (writing partial results to disk). The new generation mode can be triggered with
set new_nlo_generation True from the interface.
It can really be amazingly faster and more memory-efficient than before (Rikkert generated ttjj@NLO in 10 -ten- minutes on a 20-core cluster node), so that we can try generating higher multiplicities and look what the new bottlenecks are.
 The output should be identical (I have added an unit test to check this, yet for a very simple process : pp>w+)

As many threads start together, the logging is not as nice as the old one.

I have to thank Josh Bendavid for having provided the code for this.

Thanks for having a look and send comments...
Cheers,

Marco

[marco zaro (marco-zaro)]

just to be clear, we are speaking about code generation (not event generation, unfortunately...)

[Rikkert Frederix (frederix)]

Hi Marco,

I would rename the option "new_nlo_generation" to something like "low_memory_multicore_nlo_generation" (or something shorter/better).

It seems like Ctrl-C doesn't work correctly during the 'output' of the process.

Cheers,
Rik

review: Needs Fixing

[marco zaro (marco-zaro)]

Hi Rik
On 18 Feb 2016, at 09:29, Rikkert Frederix <email address hidden> wrote:

> Review: Needs Fixing
>
> Hi Marco,
>
> I would rename the option "new_nlo_generation" to something like "low_memory_multicore_nlo_generation" (or something shorter/better).
>
I did not think much about the name

> It seems like Ctrl-C doesn't work correctly during the 'output' of the process.
>
argh, that was something we have already seen when we used the same multicore threading for compiling…
I will have a look and let you know…
Cheers,

Marco

> Cheers,
> Rik
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.

[lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory]

311. By marco zaro on 2016-02-18

ctrl+c should be caught now

[marco zaro (marco-zaro)]

Ciao Rik,
catching ctrl+c should be ok now.
Another thing i have noticed, is that by default all available cores are used… The nb_core option should be used instead for that, i will fix also this issue
Ciao,

Marco
On 18 Feb 2016, at 10:33, marco zaro <email address hidden> wrote:

> Hi Rik
> On 18 Feb 2016, at 09:29, Rikkert Frederix <email address hidden> wrote:
>
>> Review: Needs Fixing
>>
>> Hi Marco,
>>
>> I would rename the option "new_nlo_generation" to something like "low_memory_multicore_nlo_generation" (or something shorter/better).
>>
> I did not think much about the name
>
>> It seems like Ctrl-C doesn't work correctly during the 'output' of the process.
>>
> argh, that was something we have already seen when we used the same multicore threading for compiling…
> I will have a look and let you know…
> Cheers,
>
> Marco
>
>> Cheers,
>> Rik
>>
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
>> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.

[lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory]

312. By marco zaro on 2016-02-18

changed the name of the option (new_nlo_generation -> low_mem_multicore_nlo_generation).
Now if nb_core is set, the number of used cores is determined by that

[marco zaro (marco-zaro)]

Hi,
I have changed the name of the option (low_mem_multicore_nlo_generation). Also the number of cores is dictated by nb_cores (if set)
cheers,

Marco
On 18 Feb 2016, at 12:24, marco zaro <email address hidden> wrote:

> Ciao Rik,
> catching ctrl+c should be ok now.
> Another thing i have noticed, is that by default all available cores are used… The nb_core option should be used instead for that, i will fix also this issue
> Ciao,
>
> Marco
> On 18 Feb 2016, at 10:33, marco zaro <email address hidden> wrote:
>
>> Hi Rik
>> On 18 Feb 2016, at 09:29, Rikkert Frederix <email address hidden> wrote:
>>
>>> Review: Needs Fixing
>>>
>>> Hi Marco,
>>>
>>> I would rename the option "new_nlo_generation" to something like "low_memory_multicore_nlo_generation" (or something shorter/better).
>>>
>> I did not think much about the name
>>
>>> It seems like Ctrl-C doesn't work correctly during the 'output' of the process.
>>>
>> argh, that was something we have already seen when we used the same multicore threading for compiling…
>> I will have a look and let you know…
>> Cheers,
>>
>> Marco
>>
>>> Cheers,
>>> Rik
>>>
>>> --
>>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
>>> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.
>>
>>
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
>> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.

[Rikkert Frederix (frederix)]

review: Approve

[Olivier Mattelaer (olivier-mattelaer)]

Hi Marco,

That's so cool, I missed the fact that a library was allowing to evade the GIL limitation.
This is so much nicer than using the Threading library (like the multi-core option of MG).
I guess that we should think to change the multi-core options of MG to move to this multiprocessing library.

1) Is there any reason why this mode is not set as default?
2) It should be indicated in the UpdateNotes how to use it.
3) concerning the ncores options, would not be better to pass to the class the self.options of the interface? Such that we do not have to multiply the optional argument when we to have access to any of the options define in that dictionary (this should have been done for OLP)
4)can you check why you set self['has_loops'] = False (line 466 of the diff) around line 220 of madgraph/fks/fks_helas_objects.py
6) Would not be better to not write the data on a file (via pickle), but add them in a multiprocessing.queue object that can be used by the main process? Or do you see a problem with that?
7) in madgraph/interface/amcatnlo_interface.py you define glob_directories_map
Defining a non static global is always annoying and creates troubles very often ( I have to apologize to have done it in aloha). Is it really necessary? and is it thread safe?

8) Is it possible to have a printout when the last process is generated and that you collect the info (I guess that you do that) since it hangs on one core on a line of type:
INFO: Generating born process: s~ s > t t t~ t~ [ all = QCD ]
I guess he is doing something else actually.

So this is very very nice. Nice work!!

Cheers,

Olivier

review: Needs Information

[lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory]

313. By Valentin Hirschi on 2016-02-18

1. Fixed function 'fix_coef_specs'

[Valentin Hirschi (valentin-hirschi)]

Hi Marco,

Awesome job, really.

On the testing side, I noticed the following:

a) Doing a full comparison of p p > t t~ j [QCD] generated in single and multicore, I noticed that the file 'SubProcesses/proc_characteristics' was incorrect in mutlicore, since it contains:
 ninitial = 0
 nexternal = 0

b) This implementation broke process generation with another OLP than MadLoop (i.e. Gosam).
When selecting Gosam as OLP, the generation crashes at:

FKSProcessError : Wrong type of argument for matrix_elements in function write_lh_order.

from

[...]
  File "/Users/valentin/Documents/Work/MG5/2.3.4-low-memory/madgraph/iolibs/export_fks.py", line 1355, in generate_virtuals_from_OLP
    self.write_lh_order(filename, FKSHMultiproc.get('matrix_elements'),OLP)
  File "/Users/valentin/Documents/Work/MG5/2.3.4-low-memory/madgraph/iolibs/export_fks.py", line 1623, in write_lh_order
    'matrix_elements in function write_lh_order.')

Of course we don't want to have any sort of parallelization on behalf of GoSam, but normally here the only thing that is done is the writing of the leshouches order. The actual generation with GoSam should take place at the very end, basically after that the process is generated and the parallelization is done. So if this is not too difficult to fix, I think it would be worth doing so. Otherwise, explicitly forbid this mode when OLP is not MadLoop, but that would be weird because there is no reason for doing so.
Note that you don't need to have GoSam installed to test this I think, because it crashes when writing the leshouches order file, so before actually calling GoSam.

c) This is a problem which was not introduced by this branch but that I realized now and needs to be fixed. It showed up because I saw that you had to play around with the overall 'max_loop_coefs' value for different FKSHelasMatrixElements.
So imagine you do the following:

MG5_aMC>import model HC_NLO_X0_UFO-heft
MG5_aMC>generate g g > x0 [virt=QCD]
MG5_aMC>add process g g > z [virt=QCD]
MG5_aMC>output MyTest

Then you will see the message:
"""
ML5 has just output processes which do not share the same maximum loop wavefunction size or the same maximum loop vertex rank. This is potentially dangerous. Please prefer to output them separately.
"""

This is because if you look at the different coef_specs.inc, you will find:

> cat MyTest/SubProcesses/P0_gg_x0/coef_specs.inc
      INTEGER MAXLWFSIZE
      PARAMETER (MAXLWFSIZE=4)
      INTEGER LOOP_MAXCOEFS
      PARAMETER (LOOP_MAXCOEFS=70)
      INTEGER VERTEXMAXCOEFS
      PARAMETER (VERTEXMAXCOEFS=15)

> cat MyTest/SubProcesses/P1_gg_z/coef_specs.inc
      INTEGER MAXLWFSIZE
      PARAMETER (MAXLWFSIZE=4)
      INTEGER LOOP_MAXCOEFS
      PARAMETER (LOOP_MAXCOEFS=35)
      INTEGER VERTEXMAXCOEFS
      PARAMETER (VERTEXMAXCOEFS=5)

But then you see that in Source/DHELAS you have:

> cat MyTest/Source/DHELAS/coef_specs.inc
      INTEGER MAXLWFSIZE
      PARAMETER (MAXLWFSIZE=4)
      INTEGER VERTEXMAXCOEFS
      PARAMETER (VERTEXMAXCOEFS=15)

The tricky thing here is that DHELAS, which is used by both subprocesses, depends on VERTEXMAXCOEFS which is different for both subprocesses.
It is actually fine but *...

review: Needs Fixing

[marco zaro (marco-zaro)]

Ciao Olivier,
thanks for the review
On 18 Feb 2016, at 20:44, Olivier Mattelaer <email address hidden> wrote:

> Review: Needs Information
>
> Hi Marco,
>
> That's so cool, I missed the fact that a library was allowing to evade the GIL limitation.
> This is so much nicer than using the Threading library (like the multi-core option of MG).
> I guess that we should think to change the multi-core options of MG to move to this multiprocessing library.
>
>
> 1) Is there any reason why this mode is not set as default?
for simple processes, the old mode works OK, and has a more detailed logging (number of diagrams, etc).
> 2) It should be indicated in the UpdateNotes how to use it.
done
> 3) concerning the ncores options, would not be better to pass to the class the self.options of the interface? Such that we do not have to multiply the optional argument when we to have access to any of the options define in that dictionary (this should have been done for OLP)
i don’t know… many of the interface options do not make sense for the amplitude generation, and will therefore have to be removed…

> 4)can you check why you set self['has_loops'] = False (line 466 of the diff) around line 220 of madgraph/fks/fks_helas_objects.py
it starts with false, but then, after the calls to
async_generate_born
it is updated
                has_loops = bornout[2]
                self['has_loops'] = self['has_loops'] or has_loops

> 6) Would not be better to not write the data on a file (via pickle), but add them in a multiprocessing.queue object that can be used by the main process? Or do you see a problem with that?
 The main reason to use pkl files is to keep the RAM as free as possible
> 7) in madgraph/interface/amcatnlo_interface.py you define glob_directories_map
> Defining a non static global is always annoying and creates troubles very often ( I have to apologize to have done it in aloha). Is it really necessary? and is it thread safe?
It can be put as global. Josh said that this is related to memory/cpu performance. I think it is thread safe
> 8) Is it possible to have a printout when the last process is generated and that you collect the info (I guess that you do that) since it hangs on one core on a line of type:
> INFO: Generating born process: s~ s > t t t~ t~ [ all = QCD ]
> I guess he is doing something else actually.
i have added some logging
>
> So this is very very nice. Nice work!!
>
I have pushed the changes, now I’ll look at valentin’s comments…

Cheers,

Marco

> Cheers,
>
> Olivier
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.

[lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory]

314. By marco zaro on 2016-02-19

update with olivier's comments

[Olivier Mattelaer (olivier-mattelaer)]

Thanks,

>> 3) concerning the ncores options, would not be better to pass to the class the self.options of the interface? Such that we do not have to multiply the optional argument when we to have access to any of the options define in that dictionary (this should have been done for OLP)
> i don’t know… many of the interface options do not make sense for the amplitude generation, and will therefore have to be removed…

No just keep self.options and do not touch to it
such that this is just a pointer to the one of the interface
This will not use any memory (just a pointer in term of memory), actually less that what you do right now. (since you store two additional value if we count OLP)
So yes you will have access to pointless information but I do not see the problem.

Thanks a lot for the rest.

Cheers,

Olivier

> On Feb 19, 2016, at 09:33, marco zaro <email address hidden> wrote:
>
> Ciao Olivier,
> thanks for the review
> On 18 Feb 2016, at 20:44, Olivier Mattelaer <email address hidden> wrote:
>
>> Review: Needs Information
>>
>> Hi Marco,
>>
>> That's so cool, I missed the fact that a library was allowing to evade the GIL limitation.
>> This is so much nicer than using the Threading library (like the multi-core option of MG).
>> I guess that we should think to change the multi-core options of MG to move to this multiprocessing library.
>>
>>
>> 1) Is there any reason why this mode is not set as default?
> for simple processes, the old mode works OK, and has a more detailed logging (number of diagrams, etc).
>> 2) It should be indicated in the UpdateNotes how to use it.
> done
>> 3) concerning the ncores options, would not be better to pass to the class the self.options of the interface? Such that we do not have to multiply the optional argument when we to have access to any of the options define in that dictionary (this should have been done for OLP)
> i don’t know… many of the interface options do not make sense for the amplitude generation, and will therefore have to be removed…
>
>> 4)can you check why you set self['has_loops'] = False (line 466 of the diff) around line 220 of madgraph/fks/fks_helas_objects.py
> it starts with false, but then, after the calls to
> async_generate_born
> it is updated
> has_loops = bornout[2]
> self['has_loops'] = self['has_loops'] or has_loops
>
>> 6) Would not be better to not write the data on a file (via pickle), but add them in a multiprocessing.queue object that can be used by the main process? Or do you see a problem with that?
> The main reason to use pkl files is to keep the RAM as free as possible
>> 7) in madgraph/interface/amcatnlo_interface.py you define glob_directories_map
>> Defining a non static global is always annoying and creates troubles very often ( I have to apologize to have done it in aloha). Is it really necessary? and is it thread safe?
> It can be put as global. Josh said that this is related to memory/cpu performance. I think it is thread safe
>> 8) Is it possible to have a printout when the last process is generated and that you collect the info (I guess that you do that) sinc...

[lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory]

315. By marco zaro on 2016-02-19

nexternal and ninitial are correctly written in proc_characteristics

316. By marco zaro on 2016-02-19

fixed bug occurring with write_lh_order for OLP!=MadLoop

[marco zaro (marco-zaro)]

Ciao Valentin,
thanks for your review, so:

On 18 Feb 2016, at 23:48, Valentin Hirschi <email address hidden> wrote:

> Review: Needs Fixing
>
> Hi Marco,
>
> Awesome job, really.
>
> On the testing side, I noticed the following:
>
> a) Doing a full comparison of p p > t t~ j [QCD] generated in single and multicore, I noticed that the file 'SubProcesses/proc_characteristics' was incorrect in mutlicore, since it contains:
> ninitial = 0
> nexternal = 0
hmm right, i was not checking this file with the unit test. I have added it to the test, it should be fixed now (rev 315)
>
> b) This implementation broke process generation with another OLP than MadLoop (i.e. Gosam).
> When selecting Gosam as OLP, the generation crashes at:
>
> FKSProcessError : Wrong type of argument for matrix_elements in function write_lh_order.
>
> from
>
> [...]
> File "/Users/valentin/Documents/Work/MG5/2.3.4-low-memory/madgraph/iolibs/export_fks.py", line 1355, in generate_virtuals_from_OLP
> self.write_lh_order(filename, FKSHMultiproc.get('matrix_elements'),OLP)
> File "/Users/valentin/Documents/Work/MG5/2.3.4-low-memory/madgraph/iolibs/export_fks.py", line 1623, in write_lh_order
> 'matrix_elements in function write_lh_order.’)
>
> Of course we don't want to have any sort of parallelization on behalf of GoSam, but normally here the only thing that is done is the writing of the leshouches order. The actual generation with GoSam should take place at the very end, basically after that the process is generated and the parallelization is done. So if this is not too difficult to fix, I think it would be worth doing so. Otherwise, explicitly forbid this mode when OLP is not MadLoop, but that would be weird because there is no reason for doing so.
> Note that you don't need to have GoSam installed to test this I think, because it crashes when writing the leshouches order file, so before actually calling GoSam.
>
I have added a unit test also for this and reproduced your error. It should be fixed now (rev 316), please check that i did not screw anything, as some modifications to the code were needed (the unit tests pass after trival fixes)

> c) This is a problem which was not introduced by this branch but that I realized now and needs to be fixed. It showed up because I saw that you had to play around with the overall 'max_loop_coefs' value for different FKSHelasMatrixElements.
> So imagine you do the following:
>
> MG5_aMC>import model HC_NLO_X0_UFO-heft
> MG5_aMC>generate g g > x0 [virt=QCD]
> MG5_aMC>add process g g > z [virt=QCD]
> MG5_aMC>output MyTest
>
> Then you will see the message:
> """
> ML5 has just output processes which do not share the same maximum loop wavefunction size or the same maximum loop vertex rank. This is potentially dangerous. Please prefer to output them separately.
> """
>
> This is because if you look at the different coef_specs.inc, you will find:
>
>> cat MyTest/SubProcesses/P0_gg_x0/coef_specs.inc
> INTEGER MAXLWFSIZE
> PARAMETER (MAXLWFSIZE=4)
> INTEGER LOOP_MAXCOEFS
> PARAMETER (LOOP_MAXCOEFS=70)
> INTEGER VERTEXMAXCOEFS
> PARAMETER (VERTEXMAXCOEFS=15)
>
>> cat MyTest/SubProcess...

[lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory]

317. By Valentin Hirschi on 2016-02-19

1. Fixed an issue of the FKS output when the MG option loop_color_flow is set to True.

[Valentin Hirschi (valentin-hirschi)]

Hi Marco,

One more thing, please check how you defined the possible values for the option 'low_mem_multicore_nlo_generation' that you have introduced because the auto completion to 'True False default' does not work.

> Ciao Valentin,
> thanks for your review, so:

> On 18 Feb 2016, at 23:48, Valentin Hirschi <email address hidden> wrote:
>
> > Review: Needs Fixing
> >
> > Hi Marco,
> >
> > Awesome job, really.
> >
> > On the testing side, I noticed the following:
> >
> > a) Doing a full comparison of p p > t t~ j [QCD] generated in single and
> multicore, I noticed that the file 'SubProcesses/proc_characteristics' was
> incorrect in mutlicore, since it contains:
> > ninitial = 0
> > nexternal = 0
> hmm right, i was not checking this file with the unit test. I have added it to
> the test, it should be fixed now (rev 315)
> >
> > b) This implementation broke process generation with another OLP than
> MadLoop (i.e. Gosam).
> > When selecting Gosam as OLP, the generation crashes at:
> >
> > FKSProcessError : Wrong type of argument for matrix_elements in function
> write_lh_order.
> >
> > from
> >
> > [...]
> > File "/Users/valentin/Documents/Work/MG5/2.3.4-low-
> memory/madgraph/iolibs/export_fks.py", line 1355, in
> generate_virtuals_from_OLP
> > self.write_lh_order(filename, FKSHMultiproc.get('matrix_elements'),OLP)
> > File "/Users/valentin/Documents/Work/MG5/2.3.4-low-
> memory/madgraph/iolibs/export_fks.py", line 1623, in write_lh_order
> > 'matrix_elements in function write_lh_order.’)
> >
> > Of course we don't want to have any sort of parallelization on behalf of
> GoSam, but normally here the only thing that is done is the writing of the
> leshouches order. The actual generation with GoSam should take place at the
> very end, basically after that the process is generated and the
> parallelization is done. So if this is not too difficult to fix, I think it
> would be worth doing so. Otherwise, explicitly forbid this mode when OLP is
> not MadLoop, but that would be weird because there is no reason for doing so.
> > Note that you don't need to have GoSam installed to test this I think,
> because it crashes when writing the leshouches order file, so before actually
> calling GoSam.
> >
> I have added a unit test also for this and reproduced your error. It should be
> fixed now (rev 316), please check that i did not screw anything, as some
> modifications to the code were needed (the unit tests pass after trival fixes)

So when not using the low-mem option, now OLP = GoSam crashes with:

AttributeError : 'MasterCmd' object has no attribute 'born_processes_for_olp'

at

[...]
  File "/Users/valentin/Documents/Work/MG5/2.3.4-low-memory/madgraph/interface/amcatnlo_interface.py", line 557, in do_output
    self.born_processes_for_olp,self._export_dir,self.options['OLP'])
AttributeError: 'MasterCmd' object has no attribute 'born_processes_for_olp'

And when doing it with the new option turned on, it goes through but crashes at the compilation because apparently you still linking the version of 'BinothLHA.f' for MadLoop and not the one for GoSam as it is done when not using the new option.

>
> > c) This is a problem which was n...

[lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory]

318. By marco zaro on 2016-02-23

fixed autocompletion, behaviour with gosam and [real=QCD]

[marco zaro (marco-zaro)]

Ciao Valentin,
Thanks for the review…
I should have addressed all points now, see below…
cheers,

Marco
On 19 Feb 2016, at 21:04, Valentin Hirschi <email address hidden> wrote:

> Hi Marco,
>
> One more thing, please check how you defined the possible values for the option 'low_mem_multicore_nlo_generation' that you have introduced because the auto completion to 'True False default' does not work.

fixed
>
>> Ciao Valentin,
>> thanks for your review, so:
>
>> On 18 Feb 2016, at 23:48, Valentin Hirschi <email address hidden> wrote:
>>
>>> Review: Needs Fixing
>>>
>>> Hi Marco,
>>>
>>> Awesome job, really.
>>>
>>> On the testing side, I noticed the following:
>>>
>>> a) Doing a full comparison of p p > t t~ j [QCD] generated in single and
>> multicore, I noticed that the file 'SubProcesses/proc_characteristics' was
>> incorrect in mutlicore, since it contains:
>>> ninitial = 0
>>> nexternal = 0
>> hmm right, i was not checking this file with the unit test. I have added it to
>> the test, it should be fixed now (rev 315)
>>>
>>> b) This implementation broke process generation with another OLP than
>> MadLoop (i.e. Gosam).
>>> When selecting Gosam as OLP, the generation crashes at:
>>>
>>> FKSProcessError : Wrong type of argument for matrix_elements in function
>> write_lh_order.
>>>
>>> from
>>>
>>> [...]
>>> File "/Users/valentin/Documents/Work/MG5/2.3.4-low-
>> memory/madgraph/iolibs/export_fks.py", line 1355, in
>> generate_virtuals_from_OLP
>>> self.write_lh_order(filename, FKSHMultiproc.get('matrix_elements'),OLP)
>>> File "/Users/valentin/Documents/Work/MG5/2.3.4-low-
>> memory/madgraph/iolibs/export_fks.py", line 1623, in write_lh_order
>>> 'matrix_elements in function write_lh_order.’)
>>>
>>> Of course we don't want to have any sort of parallelization on behalf of
>> GoSam, but normally here the only thing that is done is the writing of the
>> leshouches order. The actual generation with GoSam should take place at the
>> very end, basically after that the process is generated and the
>> parallelization is done. So if this is not too difficult to fix, I think it
>> would be worth doing so. Otherwise, explicitly forbid this mode when OLP is
>> not MadLoop, but that would be weird because there is no reason for doing so.
>>> Note that you don't need to have GoSam installed to test this I think,
>> because it crashes when writing the leshouches order file, so before actually
>> calling GoSam.
>>>
>> I have added a unit test also for this and reproduced your error. It should be
>> fixed now (rev 316), please check that i did not screw anything, as some
>> modifications to the code were needed (the unit tests pass after trival fixes)
>
> So when not using the low-mem option, now OLP = GoSam crashes with:
>
> AttributeError : 'MasterCmd' object has no attribute ‘born_processes_for_olp'

>
> at
>
> [...]
> File "/Users/valentin/Documents/Work/MG5/2.3.4-low-memory/madgraph/interface/amcatnlo_interface.py", line 557, in do_output
> self.born_processes_for_olp,self._export_dir,self.options['OLP'])
> AttributeError: 'MasterCmd' object has no attribute ‘born_processes_for_olp’
that should be ok now, please chec...

[Olivier Mattelaer (olivier-mattelaer)]

Hi Marco,

I will accept the merge even without the fix associate to this:

>> 3) concerning the ncores options, would not be better to pass to the class the self.options of the interface? Such that we do not have to multiply the optional argument when we to have access to any of the options define in that dictionary (this should have been done for OLP)

But please do like that next time, since it is easier if we have access to the full user options rather that duplicate the information.

Cheers,

Olivier

review: Approve

[lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory]

319. By marco zaro on 2016-02-24

changed the arguments passed to FKSMultiProcess to options

[marco zaro (marco-zaro)]

Hi Olivier,
I have changed the arguments to options, please have a look. I have also removed a chunk of code from fks_base (10 lines starting with #check limitations of FKS) which did not make any sense and was never used.

Let me know if this is fine…
Cheers,

Marco
On 24 Feb 2016, at 13:31, Olivier Mattelaer <email address hidden> wrote:

> Review: Approve
>
> Hi Marco,
>
> I will accept the merge even without the fix associate to this:
>
>>> 3) concerning the ncores options, would not be better to pass to the class the self.options of the interface? Such that we do not have to multiply the optional argument when we to have access to any of the options define in that dictionary (this should have been done for OLP)
>
> But please do like that next time, since it is easier if we have access to the full user options rather that duplicate the information.
>
> Cheers,
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.

[lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory]

320. By marco zaro on 2016-02-24

imerged with latest 2.3.4

[Valentin Hirschi (valentin-hirschi)]

Hi Marco,

I tested the OLP=GoSam behavior and a couple other things and everything went smooth.

So please just confirm that you checked that there couldn't be any border effect on the modification of the function 'RunCardNLO.create_default_for_process(' from other parts of the code that might use it and also the fact that no deepcopy at all is used in the new low-mem mode.

Anyway, I assume that the above is indeed the case so that I already give my green light.

Great work,

Cheers

review: Approve

[Olivier Mattelaer (olivier-mattelaer)]

> Hi Olivier,
> I have changed the arguments to options, please have a look. I have also
> removed a chunk of code from fks_base (10 lines starting with #check
> limitations of FKS) which did not make any sense and was never used.

That's already much nicer. Thanks.

Actually my argument was to avoid to use the full options of MadGraph:
and therefore have instead of
self._fks_multi_proc = fks_base.FKSMultiProcess(myprocdef,fks_options)
the following line:
self._fks_multi_proc = fks_base.FKSMultiProcess(myprocdef,self.options)

Obviously, this force to have a slightly different behavior for
ncores_for_proc_gen
since that one is not part of the self.options dictionary.
The two advantages: all the configuration options are available and less memory consumption (since you have one less dictionary in memory).

But this is already so much nicer that I'm more than happy with your nice change.
So this comment was just to clarify my point and not to encourage you to perform any changes.

Thanks,

Olivier

Olivier

[marco zaro (marco-zaro)]

Ciao Olivier,
thanks for the comment
On 29 Feb 2016, at 11:50, Olivier Mattelaer <email address hidden> wrote:

>> Hi Olivier,
>> I have changed the arguments to options, please have a look. I have also
>> removed a chunk of code from fks_base (10 lines starting with #check
>> limitations of FKS) which did not make any sense and was never used.
>
> That's already much nicer. Thanks.
>
> Actually my argument was to avoid to use the full options of MadGraph:
> and therefore have instead of
> self._fks_multi_proc = fks_base.FKSMultiProcess(myprocdef,fks_options)
> the following line:
> self._fks_multi_proc = fks_base.FKSMultiProcess(myprocdef,self.options)
anyway the same options will have to be passed to the MultiProcess Instance, and many of them are not valid.
So sooner or later one has to create a new dictionary…
I will proceed with the merge

Cheers,

Marco
>
> Obviously, this force to have a slightly different behavior for
> ncores_for_proc_gen
> since that one is not part of the self.options dictionary.
> The two advantages: all the configuration options are available and less memory consumption (since you have one less dictionary in memory).
>
> But this is already so much nicer that I'm more than happy with your nice change.
> So this comment was just to clarify my point and not to encourage you to perform any changes.
>
> Thanks,
>
> Olivier
>
> Olivier
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2016-02-23 22:55:27 +0000
+++ UpdateNotes.txt	2016-02-24 13:54:50 +0000
@@ -1,6 +1,12 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
 2.4.0(XX/XX/XX)
+        MZ: new NLO generation mode. It is more efficient from the memory and CPU point of 
+            view, in particular for high-multiplicity processes. 
+            Many thanks to Josh Bendavid for his fundamental contribution for this.
+            The mode can be enabled with
+            > set low_mem_multicore_nlo_generation True
+            before generating the process.
         VH: Interfaced MadLoop to Samurai and Ninja (the latter is now the default)
         OM: Adding the possibility to use new syntax for tree-level processes:
             QED==2 and QCD>2: The first allows to select exactly a power of the coupling (at amplitude level

=== modified file 'madgraph/fks/fks_base.py'
--- madgraph/fks/fks_base.py	2015-10-20 10:56:04 +0000
+++ madgraph/fks/fks_base.py	2016-02-24 13:54:50 +0000
@@ -44,15 +44,18 @@
     def default_setup(self):
         """Default values for all properties"""
         super(FKSMultiProcess, self).default_setup()
+        self['real_amplitudes'] = diagram_generation.AmplitudeList()
+        self['pdgs'] = []
         self['born_processes'] = FKSProcessList()
         if not 'OLP' in self.keys():
             self['OLP'] = 'MadLoop'
+            self['ncores_for_proc_gen'] = 0
     
     def get_sorted_keys(self):
         """Return particle property names as a nicely sorted list."""
         keys = super(FKSMultiProcess, self).get_sorted_keys()
         keys += ['born_processes', 'real_amplitudes', 'real_pdgs', 'has_isr', 
-                 'has_fsr', 'OLP']
+                 'has_fsr', 'OLP', 'ncores_for_proc_gen']
         return keys
 
     def filter(self, name, value):
@@ -77,10 +80,15 @@
             if not isinstance(value,str):
                 raise self.PhysicsObjectError, \
                     "%s is not a valid string for OLP " % str(value)
+
+        if name == 'ncores_for_proc_gen':
+            if not isinstance(value,int):
+                raise self.PhysicsObjectError, \
+                    "%s is not a valid value for ncores_for_proc_gen " % str(value)
                                                      
         return super(FKSMultiProcess,self).filter(name, value)
     
-    def __init__(self, *arguments, **options):
+    def __init__(self, procdef=None, options={}):
         """Initializes the original multiprocess, then generates the amps for the 
         borns, then generate the born processes and the reals.
         Real amplitudes are stored in real_amplitudes according on the pdgs of their
@@ -94,17 +102,25 @@
         for logg in loggers_off:
             logg.setLevel(logging.WARNING)
         
-        self['real_amplitudes'] = diagram_generation.AmplitudeList()
-        self['pdgs'] = []
-        
+        # OLP option
+        olp='MadLoop'
         if 'OLP' in options.keys():
-            self['OLP']=options['OLP']
+            olp = options['OLP']
             del options['OLP']
 
+        ncores_for_proc_gen = 0
+        # ncores_for_proc_gen has the following meaning
+        #   0 : do things the old way
+        #   > 0 use ncores_for_proc_gen
+        #   -1 : use all cores
+        if 'ncores_for_proc_gen' in options.keys():
+            ncores_for_proc_gen = options['ncores_for_proc_gen']
+            del options['ncores_for_proc_gen']
 
         try:
             # Now generating the borns for the first time.
-            super(FKSMultiProcess, self).__init__(*arguments,**options)
+            super(FKSMultiProcess, self).__init__(procdef, **options)
+
         except diagram_generation.NoDiagramException as error:
             # If no born, then this process most likely does not have any.
             raise NoBornException, "Born diagrams could not be generated for the "+\
@@ -113,17 +129,9 @@
                " processes yet, but you can still use MadLoop if you want to "+\
                "only generate them."+\
                " For this, use the 'virt=' mode, without multiparticle labels."
-        
-        #check limitation of FKS
-        if arguments and isinstance(arguments, MG.Process):
-            myprocdef = arguments[0]
-            misc.sprint( myprocdef.keys())
-            if myprocdef['perturbation_couplings']!=['QCD']:
-                raise InvalidCmd("FKS for reals only available in QCD for now, you asked %s" \
-                            % ', '.join(myprocdef['perturbation_couplings']))
-            elif myprocdef.get_ninitial()==1:
-                raise InvalidCmd("At this stage aMC@NLO cannot handle decay process.\n"+\
-                 "   Only Leading Order (loop-induced and tree level) decay are supported.") 
+
+        self['OLP'] = olp
+        self['ncores_for_proc_gen'] = ncores_for_proc_gen
             
         #check process definition(s):
         # a process such as g g > g g will lead to real emissions 
@@ -164,58 +172,63 @@
                                                                  'Process', ''),
                  i + 1, len(amps)))
 
-            born = FKSProcess(amp)
+            born = FKSProcess(amp, ncores_for_proc_gen = self['ncores_for_proc_gen'])
             self['born_processes'].append(born)
             born.generate_reals(self['pdgs'], self['real_amplitudes'])
 
-        born_pdg_list = [[l['id'] for l in born.born_proc['legs']] \
-                for born in self['born_processes'] ]
-
-        for born in self['born_processes']:
-            for real in born.real_amps:
-                real.find_fks_j_from_i(born_pdg_list)
-
-        if amps:
-            if self['process_definitions'][0].get('NLO_mode') == 'all':
-                self.generate_virtuals()
-            
-            elif not self['process_definitions'][0].get('NLO_mode') in ['all', 'real','LOonly']:
-                raise fks_common.FKSProcessError(\
-                   "Not a valid NLO_mode for a FKSMultiProcess: %s" % \
-                   self['process_definitions'][0].get('NLO_mode'))
-
-            # now get the total number of diagrams
-            n_diag_born = sum([len(amp.get('diagrams')) 
-                     for amp in self.get_born_amplitudes()])
-            n_diag_real = sum([len(amp.get('diagrams')) 
-                     for amp in self.get_real_amplitudes()])
-            n_diag_virt = sum([len(amp.get('loop_diagrams')) 
-                     for amp in self.get_virt_amplitudes()])
-
-            if n_diag_virt == 0 and n_diag_real ==0 and \
-                    not self['process_definitions'][0].get('NLO_mode') == 'LOonly':
-                raise fks_common.FKSProcessError(
-                        'This process does not have any correction up to NLO in %s'\
-                        %','.join(perturbation))
-
-            logger.info(('Generated %d subprocesses with %d real emission diagrams, ' + \
-                        '%d born diagrams and %d virtual diagrams') % \
-                                (len(self['born_processes']), n_diag_real, n_diag_born, n_diag_virt))
-
-        for i, logg in enumerate(loggers_off):
-            logg.setLevel(old_levels[i])
+        if not self['ncores_for_proc_gen']:
+            # old generation mode 
+
+            born_pdg_list = [[l['id'] for l in born.born_proc['legs']] \
+                    for born in self['born_processes'] ]
+
+            for born in self['born_processes']:
+                for real in born.real_amps:
+                    real.find_fks_j_from_i(born_pdg_list)
+            if amps:
+                if self['process_definitions'][0].get('NLO_mode') == 'all':
+                    self.generate_virtuals()
+                
+                elif not self['process_definitions'][0].get('NLO_mode') in ['all', 'real','LOonly']:
+                    raise fks_common.FKSProcessError(\
+                       "Not a valid NLO_mode for a FKSMultiProcess: %s" % \
+                       self['process_definitions'][0].get('NLO_mode'))
+
+                # now get the total number of diagrams
+                n_diag_born = sum([len(amp.get('diagrams')) 
+                         for amp in self.get_born_amplitudes()])
+                n_diag_real = sum([len(amp.get('diagrams')) 
+                         for amp in self.get_real_amplitudes()])
+                n_diag_virt = sum([len(amp.get('loop_diagrams')) 
+                         for amp in self.get_virt_amplitudes()])
+
+                if n_diag_virt == 0 and n_diag_real ==0 and \
+                        not self['process_definitions'][0].get('NLO_mode') == 'LOonly':
+                    raise fks_common.FKSProcessError(
+                            'This process does not have any correction up to NLO in %s'\
+                            %','.join(perturbation))
+
+                logger.info(('Generated %d subprocesses with %d real emission diagrams, ' + \
+                            '%d born diagrams and %d virtual diagrams') % \
+                                    (len(self['born_processes']), n_diag_real, n_diag_born, n_diag_virt))
+
+            for i, logg in enumerate(loggers_off):
+                logg.setLevel(old_levels[i])
 
         self['has_isr'] = any([proc.isr for proc in self['born_processes']])
         self['has_fsr'] = any([proc.fsr for proc in self['born_processes']])
 
     def add(self, other):
         """combines self and other, extending the lists of born/real amplitudes"""
+        self['process_definitions'].extend(other['process_definitions'])
+        self['amplitudes'].extend(other['amplitudes'])
         self['born_processes'].extend(other['born_processes'])
         self['real_amplitudes'].extend(other['real_amplitudes'])
         self['pdgs'].extend(other['pdgs'])
         self['has_isr'] = self['has_isr'] or other['has_isr']
         self['has_fsr'] = self['has_fsr'] or other['has_fsr']
         self['OLP'] = other['OLP']
+        self['ncores_for_proc_gen'] = other['ncores_for_proc_gen']
 
     def get_born_amplitudes(self):
         """return an amplitudelist with the born amplitudes"""
@@ -394,11 +407,16 @@
     """The class for a FKS process. Starts from the born process and finds
     all the possible splittings."""  
     
-    def __init__(self, start_proc = None, remove_reals = True):
+    def __init__(self, start_proc = None, remove_reals = True, ncores_for_proc_gen=0):
         """initialization: starts either from an amplitude or a process,
         then init the needed variables.
         remove_borns tells if the borns not needed for integration will be removed
-        from the born list (mainly used for testing)"""
+        from the born list (mainly used for testing)
+        ncores_for_proc_gen has the following meaning
+           0 : do things the old way
+           > 0 use ncores_for_proc_gen
+           -1 : use all cores
+        """
                 
         self.splittings = {}
         self.reals = []
@@ -416,6 +434,7 @@
         self.nincoming = 0
         self.virt_amp = None
         self.perturbation = 'QCD'
+        self.ncores_for_proc_gen = ncores_for_proc_gen
 
         if not remove_reals in [True, False]:
             raise fks_common.FKSProcessError(\
@@ -509,7 +528,6 @@
         self.real_amps = real_amps
 
 
-        
     def generate_reals(self, pdg_list, real_amp_list, combine=True): #test written
         """For all the possible splittings, creates an FKSRealProcess.
         It removes double counted configorations from the ones to integrates and
@@ -534,8 +552,9 @@
         self.find_reals_to_integrate()
         if combine:
             self.combine_real_amplitudes()
-        self.generate_real_amplitudes(pdg_list, real_amp_list)
-        self.link_born_reals()
+        if not self.ncores_for_proc_gen:
+            self.generate_real_amplitudes(pdg_list, real_amp_list)
+            self.link_born_reals()
 
 
     def link_born_reals(self):

=== modified file 'madgraph/fks/fks_helas_objects.py'
--- madgraph/fks/fks_helas_objects.py	2015-10-01 16:00:08 +0000
+++ madgraph/fks/fks_helas_objects.py	2016-02-24 13:54:50 +0000
@@ -25,20 +25,183 @@
 import madgraph.fks.fks_base as fks_base
 import madgraph.fks.fks_common as fks_common
 import madgraph.loop.loop_helas_objects as loop_helas_objects
+import madgraph.loop.loop_diagram_generation as loop_diagram_generation
 import copy
 import logging
 import array
+import multiprocessing
+import signal
+import tempfile
+import cPickle
+import itertools
+import os
 
 logger = logging.getLogger('madgraph.fks_helas_objects')
 
 
+#functions to be used in the ncores_for_proc_gen mode
+def async_generate_real(args):
+    i = args[0]
+    real_amp = args[1]
+
+    #amplitude generation
+    amplitude = real_amp.generate_real_amplitude()
+    helasreal = helas_objects.HelasMatrixElement(amplitude)
+    logger.info('Generating real %s' % \
+            real_amp.process.nice_string(print_weighted=False).replace('Process', 'process'))
+
+    # Keep track of already generated color objects, to reuse as
+    # much as possible
+    list_colorize = []
+    list_color_basis = []
+    list_color_matrices = []
+    
+    # Now this keeps track of the color matrices created from the loop-born
+    # color basis. Keys are 2-tuple with the index of the loop and born basis
+    # in the list above and the value is the resulting matrix.
+    dict_loopborn_matrices = {}
+    # The dictionary below is simply a container for convenience to be 
+    # passed to the function process_color.
+    color_information = { 'list_colorize' : list_colorize,
+                          'list_color_basis' : list_color_basis,
+                          'list_color_matrices' : list_color_matrices,
+                          'dict_loopborn_matrices' : dict_loopborn_matrices}
+
+    helas_objects.HelasMultiProcess.process_color(helasreal,color_information)
+
+    outdata = [amplitude,helasreal]
+
+    output = tempfile.NamedTemporaryFile(delete = False)   
+    cPickle.dump(outdata,output,protocol=2)
+    output.close()
+    
+    return [output.name,helasreal.get_num_configs(),helasreal.get_nexternal_ninitial()[0]]
+
+
+def async_generate_born(args):
+    i = args[0]
+    born = args[1]
+    born_pdg_list = args[2]
+    loop_orders = args[3]
+    pdg_list = args[4]
+    loop_optimized = args[5]
+    OLP = args[6]
+    realmapout = args[7]
+
+    logger.info('Generating born %s' % \
+            born.born_proc.nice_string(print_weighted=False).replace('Process', 'process'))
+    
+    #load informations on reals from temp files
+    helasreal_list = []
+    for amp in born.real_amps:
+        idx = pdg_list.index(amp.pdgs)
+        infilename = realmapout[idx]
+        infile = open(infilename,'rb')
+        realdata = cPickle.load(infile)
+        infile.close()
+        amp.amplitude = realdata[0]
+        helasreal_list.append(realdata[1])
+        
+    born.link_born_reals()
+        
+    for amp in born.real_amps:
+        amp.find_fks_j_from_i(born_pdg_list)        
+    
+    # generate the virtuals if needed
+    has_loops = False
+    if born.born_proc.get('NLO_mode') == 'all' and OLP == 'MadLoop':
+        myproc = copy.copy(born.born_proc)
+        # take the orders that are actually used by the matrix element
+        myproc['orders'] = loop_orders
+        myproc['legs'] = fks_common.to_legs(copy.copy(myproc['legs']))
+        myamp = loop_diagram_generation.LoopAmplitude(myproc)
+        if myamp.get('diagrams'):
+            has_loops = True
+            born.virt_amp = myamp
+        
+    helasfull = FKSHelasProcess(born, helasreal_list,
+                                loop_optimized = loop_optimized,
+                                decay_ids=[],
+                                gen_color=False)
+
+    processes = helasfull.born_matrix_element.get('processes')
+    
+    metag = helas_objects.IdentifyMETag.create_tag(helasfull.born_matrix_element.get('base_amplitude'))
+    
+    outdata = helasfull
+    
+    output = tempfile.NamedTemporaryFile(delete = False)   
+    cPickle.dump(outdata,output,protocol=2)
+    output.close()
+    
+    return [output.name,metag,has_loops,processes]
+
+
+def async_finalize_matrix_elements(args):
+
+    i = args[0]
+    mefile = args[1]
+    duplist = args[2]
+    
+    infile = open(mefile,'rb')
+    me = cPickle.load(infile)
+    infile.close()    
+
+    #set unique id based on position in unique me list
+    me.get('processes')[0].set('uid', i)
+
+    # Always create an empty color basis, and the
+    # list of raw colorize objects (before
+    # simplification) associated with amplitude
+    col_basis = color_amp.ColorBasis()
+    new_amp = me.born_matrix_element.get_base_amplitude()
+    me.born_matrix_element.set('base_amplitude', new_amp)
+    colorize_obj = col_basis.create_color_dict_list(new_amp)
+
+    col_basis.build()
+    col_matrix = color_amp.ColorMatrix(col_basis)
+
+    me.born_matrix_element.set('color_basis',col_basis)
+    me.born_matrix_element.set('color_matrix',col_matrix)
+    
+    for iother,othermefile in enumerate(duplist):
+        infileother = open(othermefile,'rb')
+        otherme = cPickle.load(infileother)
+        infileother.close()
+        me.add_process(otherme)
+        
+    me.set_color_links()    
+        
+    initial_states=[]
+    for fksreal in me.real_processes:
+        # Pick out all initial state particles for the two beams
+            initial_states.append(sorted(list(set((p.get_initial_pdg(1),p.get_initial_pdg(2)) for \
+                                              p in fksreal.matrix_element.get('processes')))))
+    
+    if me.virt_matrix_element:
+        has_virtual = True
+    else:
+        has_virtual = False
+     
+    #data to write to file
+    outdata = me
+
+    output = tempfile.NamedTemporaryFile(delete = False)   
+    cPickle.dump(outdata,output,protocol=2)
+    output.close()
+    
+    #data to be returned to parent process (filename plus small objects only)
+    return [output.name,initial_states,me.get_used_lorentz(),me.get_used_couplings(),has_virtual]
+
+
 class FKSHelasMultiProcess(helas_objects.HelasMultiProcess):
     """class to generate the helas calls for a FKSMultiProcess"""
 
     def get_sorted_keys(self):
         """Return particle property names as a nicely sorted list."""
         keys = super(FKSHelasMultiProcess, self).get_sorted_keys()
-        keys += ['real_matrix_elements', ['has_isr'], ['has_fsr']]
+        keys += ['real_matrix_elements', ['has_isr'], ['has_fsr'], 
+                 'used_lorentz', 'used_couplings', 'max_configs', 'max_particles', 'processes']
         return keys
 
     def filter(self, name, value):
@@ -61,23 +224,190 @@
 
         self.loop_optimized = loop_optimized
 
-        # generate the real ME's if they are needed.
-        # note that it may not be always the case, e.g. it the NLO_mode is LOonly
-        if fksmulti['real_amplitudes']:
-            logger.info('Generating real emission matrix-elements...')
-            self['real_matrix_elements'] = self.generate_matrix_elements(
-                    copy.copy(fksmulti['real_amplitudes']), combine_matrix_elements = False)
-        else:
-            self['real_matrix_elements'] = helas_objects.HelasMatrixElementList()
-
-        self['matrix_elements'] = self.generate_matrix_elements_fks(
-                                fksmulti, 
-                                gen_color, decay_ids)
-        self['initial_states']=[]
+        self['used_lorentz'] = []
+        self['used_couplings'] = []
+        self['processes'] = []
+
+        self['max_particles'] = -1
+        self['max_configs'] = -1
+
+        if not fksmulti['ncores_for_proc_gen']:
+            # generate the real ME's if they are needed.
+            # note that it may not be always the case, e.g. it the NLO_mode is LOonly
+            if fksmulti['real_amplitudes']:
+                logger.info('Generating real emission matrix-elements...')
+                self['real_matrix_elements'] = self.generate_matrix_elements(
+                        copy.copy(fksmulti['real_amplitudes']), combine_matrix_elements = False)
+            else:
+                self['real_matrix_elements'] = helas_objects.HelasMatrixElementList()
+
+            self['matrix_elements'] = self.generate_matrix_elements_fks(
+                                    fksmulti, 
+                                    gen_color, decay_ids)
+            self['initial_states']=[]
+            self['has_loops'] = len(self.get_virt_matrix_elements()) > 0 
+
+        else: 
+            self['has_loops'] = False
+            #more efficient generation
+            born_procs = fksmulti.get('born_processes')
+            born_pdg_list = [[l['id'] for l in born.born_proc['legs']] \
+            for born in born_procs ]
+            loop_orders = {}
+            for  born in born_procs:
+                for coup, val in fks_common.find_orders(born.born_amp).items():
+                    try:
+                        loop_orders[coup] = max([loop_orders[coup], val])
+                    except KeyError:
+                        loop_orders[coup] = val        
+            pdg_list = []        
+            real_amp_list = []
+            for born in born_procs:
+                for amp in born.real_amps:
+                    if not pdg_list.count(amp.pdgs):
+                        pdg_list.append(amp.pdgs)
+                        real_amp_list.append(amp)
+                        
+            #generating and store in tmp files all output corresponding to each real_amplitude
+            real_out_list = []
+            realmapin = []
+            for i,real_amp in enumerate(real_amp_list):
+                realmapin.append([i,real_amp])
+
+            # start the pool instance with a signal instance to catch ctr+c
+            original_sigint_handler = signal.signal(signal.SIGINT, signal.SIG_IGN)
+            if fksmulti['ncores_for_proc_gen'] < 0: # use all cores
+                pool = multiprocessing.Pool(maxtasksperchild=1)
+            else:
+                pool = multiprocessing.Pool(processes=fksmulti['ncores_for_proc_gen'],maxtasksperchild=1)
+            signal.signal(signal.SIGINT, original_sigint_handler)
+
+            logger.info('Generating real matrix elements...')
+            try:
+                # the very large timeout passed to get is to be able to catch
+                # KeyboardInterrupts
+                realmapout = pool.map_async(async_generate_real,realmapin).get(9999999)
+            except KeyboardInterrupt:
+                pool.terminate()
+                raise KeyboardInterrupt 
+
+            realmapfiles = []
+            for realout in realmapout:
+                realmapfiles.append(realout[0])
+
+            logger.info('Generating born and virtual matrix elements...')
+            #now loop over born and consume reals, generate virtuals
+            bornmapin = []
+            OLP=fksmulti['OLP']
+            for i,born in enumerate(born_procs):
+                bornmapin.append([i,born,born_pdg_list,loop_orders,pdg_list,loop_optimized,OLP,realmapfiles])
+
+            try:
+                bornmapout = pool.map_async(async_generate_born,bornmapin).get(9999999)
+            except KeyboardInterrupt:
+                pool.terminate()
+                raise KeyboardInterrupt 
+
+            #remove real temp files
+            for realtmp in realmapout:
+                os.remove(realtmp[0])
+                
+            logger.info('Collecting infos and finalizing matrix elements...')
+            unique_me_list = []
+            duplicate_me_lists = []
+            for bornout in bornmapout:
+                mefile = bornout[0]
+                metag = bornout[1]
+                has_loops = bornout[2]
+                self['has_loops'] = self['has_loops'] or has_loops
+                processes = bornout[3]
+                self['processes'].extend(processes)
+                unique = True
+                for ime2,bornout2 in enumerate(unique_me_list):
+                    mefile2 = bornout2[0]
+                    metag2 = bornout2[1]
+                    if metag==metag2:
+                        duplicate_me_lists[ime2].append(mefile)
+                        unique = False
+                        break;
+                if unique:
+                    unique_me_list.append(bornout)
+                    duplicate_me_lists.append([])
+            
+            memapin = []
+            for i,bornout in enumerate(unique_me_list):
+                mefile = bornout[0]
+                memapin.append([i,mefile, duplicate_me_lists[i]])
+
+            try:
+                memapout = pool.map_async(async_finalize_matrix_elements,memapin).get(9999999)
+            except KeyboardInterrupt:
+                pool.terminate()
+                raise KeyboardInterrupt 
+
+            #remove born+virtual temp files
+            for bornout in bornmapout:
+                mefile = bornout[0]
+                os.remove(mefile)
+
+            pool.close()
+            pool.join()
+
+            #set final list of matrix elements (paths to temp files)
+            matrix_elements = []
+            for meout in memapout:
+                matrix_elements.append(meout[0])
+  
+            self['matrix_elements']=matrix_elements
+  
+            #cache information needed for output which will not be available from
+            #the matrix elements later
+            initial_states = []
+            for meout in memapout:
+                me_initial_states = meout[1]
+                for state in me_initial_states:
+                    initial_states.append(state)
+                              
+            # remove doubles from the list
+            checked = []
+            for e in initial_states:
+                if e not in checked:
+                    checked.append(e)
+            initial_states=checked
+
+            self['initial_states']=initial_states
+            
+            helas_list = []
+            for meout in memapout:
+                helas_list.extend(meout[2])
+            self['used_lorentz']=list(set(helas_list))        
+            
+            coupling_list = []
+            for meout in memapout:
+                coupling_list.extend([c for l in meout[3] for c in l])
+            self['used_couplings'] = list(set(coupling_list)) 
+            
+            has_virtuals = False
+            for meout in memapout:
+                if meout[4]:
+                    has_virtuals = True
+                    break
+            self['has_virtuals'] = has_virtuals
+            
+            configs_list = []
+            for meout in realmapout:
+                configs_list.append(meout[1])
+            self['max_configs'] = max(configs_list)
+            
+            nparticles_list = []
+            for meout in realmapout:
+                nparticles_list.append(meout[2])
+            self['max_particles'] = max(nparticles_list)        
 
         self['has_isr'] = fksmulti['has_isr']
         self['has_fsr'] = fksmulti['has_fsr']
-        self['has_loops'] = len(self.get_virt_matrix_elements()) > 0 
+
+        logger.info('... Done')
 
         for i, logg in enumerate(loggers_off):
             logg.setLevel(old_levels[i])
@@ -85,19 +415,66 @@
     def get_used_lorentz(self):
         """Return a list of (lorentz_name, conjugate, outgoing) with
         all lorentz structures used by this HelasMultiProcess."""
-        helas_list = []
-        for me in self.get('matrix_elements'):
-            helas_list.extend(me.get_used_lorentz())
-        return list(set(helas_list))
+
+        if not self['used_lorentz']:
+            helas_list = []
+            for me in self.get('matrix_elements'):
+                helas_list.extend(me.get_used_lorentz())
+            self['used_lorentz'] = list(set(helas_list))
+
+        return self['used_lorentz']
+
 
     def get_used_couplings(self):
         """Return a list with all couplings used by this
         HelasMatrixElement."""
-        coupling_list = []
-        for me in self.get('matrix_elements'):
-            coupling_list.extend([c for l in me.get_used_couplings() for c in l])
-        return list(set(coupling_list))
+
+        if not self['used_couplings']:
+            coupling_list = []
+            for me in self.get('matrix_elements'):
+                coupling_list.extend([c for l in me.get_used_couplings() for c in l])
+            self['used_couplings'] = list(set(coupling_list))
+
+        return self['used_couplings']
+
+
+    def get_processes(self):
+        """Return a list with all couplings used by this
+        HelasMatrixElement."""
+
+        if not self['processes']:
+            process_list = []
+            for me in self.get('matrix_elements'):
+                process_list.extend(me.born_matrix_element.get('processes'))
+            self['processes'] = process_list
+
+        return self['processes']
+
+
+    def get_max_configs(self):
+        """Return max_configs"""
+
+        if self['max_configs'] < 0:
+            try:
+                self['max_configs'] = max([me.get_num_configs() \
+                                  for me in self['real_matrix_elements']])
+            except ValueError:
+                self['max_configs'] = max([me.born_matrix_element.get_num_configs() \
+                                  for me in self['matrix_elements']])
+
+        return self['max_configs']
+
+
+    def get_max_particles(self):
+        """Return max_paricles"""
+
+        if self['max_particles'] < 0:
+            self['max_particles'] = max([me.get_nexternal_ninitial()[0] \
+                                for me in self['matrix_elements']])
+
+        return self['max_particles']
     
+
     def get_matrix_elements(self):
         """Extract the list of matrix elements"""
         return self.get('matrix_elements')        
@@ -240,17 +617,32 @@
             self.perturbation = fksproc.perturbation
             real_amps_new = []
             # combine for example u u~ > t t~ and d d~ > t t~
-            for proc in fksproc.real_amps:
-                fksreal_me = FKSHelasRealProcess(proc, real_me_list, real_amp_list, **opts)
-                try:
-                    other = self.real_processes[self.real_processes.index(fksreal_me)]
-                    other.matrix_element.get('processes').extend(\
-                            fksreal_me.matrix_element.get('processes') )
-                except ValueError:
-                    if fksreal_me.matrix_element.get('processes') and \
-                            fksreal_me.matrix_element.get('diagrams'):
-                        self.real_processes.append(fksreal_me)
-                        real_amps_new.append(proc)
+            if fksproc.ncores_for_proc_gen:
+                # new NLO (multicore) generation mode 
+                for real_me, proc in itertools.izip(real_me_list,fksproc.real_amps):
+                    fksreal_me = FKSHelasRealProcess(proc, real_me, **opts)
+                    try:
+                        other = self.real_processes[self.real_processes.index(fksreal_me)]
+                        other.matrix_element.get('processes').extend(\
+                                fksreal_me.matrix_element.get('processes') )
+                    except ValueError:
+                        if fksreal_me.matrix_element.get('processes') and \
+                                fksreal_me.matrix_element.get('diagrams'):
+                            self.real_processes.append(fksreal_me)
+                            real_amps_new.append(proc)
+            else:
+                #old mode
+                for proc in fksproc.real_amps:
+                    fksreal_me = FKSHelasRealProcess(proc, real_me_list, real_amp_list, **opts)
+                    try:
+                        other = self.real_processes[self.real_processes.index(fksreal_me)]
+                        other.matrix_element.get('processes').extend(\
+                                fksreal_me.matrix_element.get('processes') )
+                    except ValueError:
+                        if fksreal_me.matrix_element.get('processes') and \
+                                fksreal_me.matrix_element.get('diagrams'):
+                            self.real_processes.append(fksreal_me)
+                            real_amps_new.append(proc)
             fksproc.real_amps = real_amps_new
             if fksproc.virt_amp:
                 self.virt_matrix_element = \
@@ -410,23 +802,35 @@
             self.fks_infos = fksrealproc.fks_infos
             self.is_to_integrate = fksrealproc.is_to_integrate
 
-            if len(real_me_list) != len(real_amp_list):
+            # real_me_list is a list in the old NLO generation mode;
+            # in the new one it is a matrix element
+            if type(real_me_list) == list and len(real_me_list) != len(real_amp_list):
                 raise fks_common.FKSProcessError(
                         'not same number of amplitudes and matrix elements: %d, %d' % \
                                 (len(real_amp_list), len(real_me_list)))
-            if real_me_list and real_amp_list:
+            if type(real_me_list) == list and real_me_list and real_amp_list:
                 self.matrix_element = copy.deepcopy(real_me_list[real_amp_list.index(fksrealproc.amplitude)])
                 self.matrix_element['processes'] = copy.deepcopy(self.matrix_element['processes'])
+
+            elif type(real_me_list) == helas_objects.HelasMatrixElement: 
+                #new NLO generation mode
+                self.matrix_element = real_me_list
+
             else:
-                logger.info('generating matrix element...')
-                self.matrix_element = helas_objects.HelasMatrixElement(
-                                                  fksrealproc.amplitude, **opts)
-                #generate the color for the real
-                self.matrix_element.get('color_basis').build(
-                                    self.matrix_element.get('base_amplitude'))
-                self.matrix_element.set('color_matrix',
-                                 color_amp.ColorMatrix(
-                                    self.matrix_element.get('color_basis')))
+
+                if real_me_list and real_amp_list:
+                    self.matrix_element = copy.deepcopy(real_me_list[real_amp_list.index(fksrealproc.amplitude)])
+                    self.matrix_element['processes'] = copy.deepcopy(self.matrix_element['processes'])
+                else:
+                    logger.info('generating matrix element...')
+                    self.matrix_element = helas_objects.HelasMatrixElement(
+                                                      fksrealproc.amplitude, **opts)
+                    #generate the color for the real
+                    self.matrix_element.get('color_basis').build(
+                                        self.matrix_element.get('base_amplitude'))
+                    self.matrix_element.set('color_matrix',
+                                     color_amp.ColorMatrix(
+                                        self.matrix_element.get('color_basis')))
             #self.fks_j_from_i = fksrealproc.find_fks_j_from_i()
             self.fks_j_from_i = fksrealproc.fks_j_from_i
 

=== modified file 'madgraph/interface/amcatnlo_interface.py'
--- madgraph/interface/amcatnlo_interface.py	2015-10-20 14:39:50 +0000
+++ madgraph/interface/amcatnlo_interface.py	2016-02-24 13:54:50 +0000
@@ -24,6 +24,13 @@
 import optparse
 import subprocess
 import shutil
+import multiprocessing
+import signal
+import tempfile
+import itertools
+import os
+import cPickle
+
 
 import madgraph
 from madgraph import MG4DIR, MG5DIR, MadGraph5Error
@@ -52,6 +59,37 @@
 logger = logging.getLogger('cmdprint') # -> stdout
 logger_stderr = logging.getLogger('fatalerror') # ->stderr
 
+# a new function for the improved NLO generation
+glob_directories_map = []
+def generate_directories_fks_async(i):
+        
+    arglist = glob_directories_map[i]
+    
+    curr_exporter = arglist[0]
+    mefile = arglist[1]
+    curr_fortran_model = arglist[2]
+    ime = arglist[3]
+    nme = arglist[4]
+    path = arglist[5]
+    olpopts = arglist[6]
+    
+    infile = open(mefile,'rb')
+    me = cPickle.load(infile)
+    infile.close()      
+    
+    calls = curr_exporter.generate_directories_fks(me, curr_fortran_model, ime, nme, path, olpopts)
+    nexternal = curr_exporter.proc_characteristic['nexternal']
+    ninitial = curr_exporter.proc_characteristic['ninitial']
+    processes = me.born_matrix_element.get('processes')
+    
+    #only available after export has been done, so has to be returned from here
+    max_loop_vertex_rank = -99
+    if me.virt_matrix_element:
+        max_loop_vertex_rank = me.virt_matrix_element.get_max_loop_vertex_rank()  
+    
+    return [calls, curr_exporter.fksdirs, max_loop_vertex_rank, ninitial, nexternal, processes]
+
+
 class CheckFKS(mg_interface.CheckValidForCmd):
 
 
@@ -431,16 +469,30 @@
 #                
 #                raise self.InvalidCmd("FKS for reals only available in QCD for now, you asked %s" \
 #                        % ', '.join(myprocdef['perturbation_couplings']))
+        ##
+
+        # if the new nlo process generation mode is enabled, the number of cores to be
+        # used has to be passed
+        # ncores_for_proc_gen has the following meaning
+        #   0 : do things the old way
+        #   > 0 use ncores_for_proc_gen
+        #   -1 : use all cores
+        if self.options['low_mem_multicore_nlo_generation']:
+            if self.options['nb_core']:
+                self.ncores_for_proc_gen = int(self.options['nb_core'])
+            else:
+                self.ncores_for_proc_gen = -1
+        else:
+            self.ncores_for_proc_gen = 0
+
+        # this is the options dictionary to pass to the FKSMultiProcess
+        fks_options = {'OLP': self.options['OLP'],
+                       'ignore_six_quark_processes': self.options['ignore_six_quark_processes'],
+                       'ncores_for_proc_gen': self.ncores_for_proc_gen}
         try:
-            self._fks_multi_proc.add(fks_base.FKSMultiProcess(myprocdef,
-                                   collect_mirror_procs,
-                                   ignore_six_quark_processes,
-                                   OLP=self.options['OLP']))
+            self._fks_multi_proc.add(fks_base.FKSMultiProcess(myprocdef,fks_options))
         except AttributeError: 
-            self._fks_multi_proc = fks_base.FKSMultiProcess(myprocdef,
-                                   collect_mirror_procs,
-                                   ignore_six_quark_processes,
-                                   OLP=self.options['OLP'])
+            self._fks_multi_proc = fks_base.FKSMultiProcess(myprocdef,fks_options)
 
 
     def do_output(self, line):
@@ -498,7 +550,7 @@
         # Generate the virtuals if from OLP
         if self.options['OLP']!='MadLoop':
             self._curr_exporter.generate_virtuals_from_OLP(
-              self._curr_matrix_elements,self._export_dir,self.options['OLP'])
+              self.born_processes_for_olp,self._export_dir,self.options['OLP'])
                 
         # Remember that we have done export
         self._done_export = (self._export_dir, self._export_format)
@@ -531,41 +583,50 @@
                                 self._fks_multi_proc, 
                                 loop_optimized= self.options['loop_optimized_output'])
                     
-                    ndiags = sum([len(me.get('diagrams')) for \
-                                  me in self._curr_matrix_elements.\
-                                  get_matrix_elements()])
-                    # assign a unique id number to all process and
-                    # generate a list of possible PDF combinations
-                    uid = 0 
-                    initial_states=[]
-                    for me in self._curr_matrix_elements.get_matrix_elements():
-                        uid += 1 # update the identification number
-                        me.get('processes')[0].set('uid', uid)
-                        try:
-                            initial_states.append(sorted(list(set((p.get_initial_pdg(1),p.get_initial_pdg(2)) for \
-                                                                  p in me.born_matrix_element.get('processes')))))
-                        except IndexError:
-                            initial_states.append(sorted(list(set((p.get_initial_pdg(1)) for \
-                                                                  p in me.born_matrix_element.get('processes')))))
-                    
-                        for fksreal in me.real_processes:
-                        # Pick out all initial state particles for the two beams
+                    if not self.options['low_mem_multicore_nlo_generation']: 
+                        # generate the code the old way
+                        ndiags = sum([len(me.get('diagrams')) for \
+                                      me in self._curr_matrix_elements.\
+                                      get_matrix_elements()])
+                        # assign a unique id number to all process and
+                        # generate a list of possible PDF combinations
+                        uid = 0 
+                        initial_states=[]
+                        for me in self._curr_matrix_elements.get_matrix_elements():
+                            uid += 1 # update the identification number
+                            me.get('processes')[0].set('uid', uid)
                             try:
                                 initial_states.append(sorted(list(set((p.get_initial_pdg(1),p.get_initial_pdg(2)) for \
-                                                             p in fksreal.matrix_element.get('processes')))))
+                                                                      p in me.born_matrix_element.get('processes')))))
                             except IndexError:
                                 initial_states.append(sorted(list(set((p.get_initial_pdg(1)) for \
-                                                             p in fksreal.matrix_element.get('processes')))))
-                                
+                                                                      p in me.born_matrix_element.get('processes')))))
                         
-                    # remove doubles from the list
-                    checked = []
-                    for e in initial_states:
-                        if e not in checked:
-                            checked.append(e)
-                    initial_states=checked
-
-                    self._curr_matrix_elements.set('initial_states',initial_states)
+                            for fksreal in me.real_processes:
+                            # Pick out all initial state particles for the two beams
+                                try:
+                                    initial_states.append(sorted(list(set((p.get_initial_pdg(1),p.get_initial_pdg(2)) for \
+                                                                 p in fksreal.matrix_element.get('processes')))))
+                                except IndexError:
+                                    initial_states.append(sorted(list(set((p.get_initial_pdg(1)) for \
+                                                                 p in fksreal.matrix_element.get('processes')))))
+                                    
+                            
+                        # remove doubles from the list
+                        checked = []
+                        for e in initial_states:
+                            if e not in checked:
+                                checked.append(e)
+                        initial_states=checked
+
+                        self._curr_matrix_elements.set('initial_states',initial_states)
+
+                    else:
+                        #new NLO generation
+                        if self._curr_matrix_elements['has_loops']:
+                            self._curr_exporter.opt['mp'] = True
+                        self._curr_exporter.model = self._curr_model
+                        ndiags = 0
 
             cpu_time2 = time.time()
             return ndiags, cpu_time2 - cpu_time1
@@ -586,23 +647,84 @@
             for charac in ['has_isr', 'has_fsr', 'has_loops']:
                 proc_charac[charac] = self._curr_matrix_elements[charac]
 
+            # prepare for the generation
+            # glob_directories_map is for the new NLO generation
+            global glob_directories_map
+            glob_directories_map = []
 
+            self.born_processes_for_olp = []
             for ime, me in \
                 enumerate(self._curr_matrix_elements.get('matrix_elements')):
-                #me is a FKSHelasProcessFromReals
-                calls = calls + \
-                        self._curr_exporter.generate_directories_fks(me, 
-                        self._curr_fortran_model, 
-                        ime, len(self._curr_matrix_elements.get('matrix_elements')), 
-                        path,self.options['OLP'])
-                self._fks_directories.extend(self._curr_exporter.fksdirs)
+                if not self.options['low_mem_multicore_nlo_generation']:
+                    #me is a FKSHelasProcessFromReals
+                    calls = calls + \
+                            self._curr_exporter.generate_directories_fks(me, 
+                            self._curr_fortran_model, 
+                            ime, len(self._curr_matrix_elements.get('matrix_elements')), 
+                            path,self.options['OLP'])
+                    self._fks_directories.extend(self._curr_exporter.fksdirs)
+                    self.born_processes_for_olp.append(me.born_matrix_element.get('processes')[0])
+                else:
+                    glob_directories_map.append(\
+                            [self._curr_exporter, me, self._curr_fortran_model, 
+                             ime, len(self._curr_matrix_elements.get('matrix_elements')), 
+                             path, self.options['OLP']])
+
+            if self.options['low_mem_multicore_nlo_generation']:
+                # start the pool instance with a signal instance to catch ctr+c
+                logger.info('Writing directories...')
+                original_sigint_handler = signal.signal(signal.SIGINT, signal.SIG_IGN)
+                if self.ncores_for_proc_gen < 0: # use all cores
+                    pool = multiprocessing.Pool(maxtasksperchild=1)
+                else:
+                    pool = multiprocessing.Pool(processes=self.ncores_for_proc_gen,maxtasksperchild=1)
+                signal.signal(signal.SIGINT, original_sigint_handler)
+                try:
+                    # the very large timeout passed to get is to be able to catch
+                    # KeyboardInterrupts
+                    diroutputmap = pool.map_async(generate_directories_fks_async,
+                                                  range(len(glob_directories_map))).get(9999999)
+                except KeyboardInterrupt:
+                    pool.terminate()
+                    raise KeyboardInterrupt 
+    
+                pool.close()
+                pool.join()
+                
+                #clean up tmp files containing final matrix elements
+                for mefile in self._curr_matrix_elements.get('matrix_elements'):
+                    os.remove(mefile)
+
+                for charac in ['nexternal', 'ninitial']:
+                    proc_charac[charac] = self._curr_exporter.proc_characteristic[charac]
+                # ninitial and nexternal
+                proc_charac['nexternal'] = max([diroutput[4] for diroutput in diroutputmap])
+                ninitial_set = set([diroutput[3] for diroutput in diroutputmap])
+                if len(ninitial_set) != 1:
+                    raise MadGraph5Error, ("Invalid ninitial values: %s" % ' ,'.join(list(ninitial_set)))    
+                proc_charac['ninitial'] = list(ninitial_set)[0]
+
+                self.born_processes = []
+                self.born_processes_for_olp = []
+                max_loop_vertex_ranks = []
+                
+                for diroutput in diroutputmap:
+                    calls = calls + diroutput[0]
+                    self._fks_directories.extend(diroutput[1])
+                    max_loop_vertex_ranks.append(diroutput[2])
+                    self.born_processes.extend(diroutput[5])
+                    self.born_processes_for_olp.append(diroutput[5][0])
+
+            else:
+                max_loop_vertex_ranks = [me.get_max_loop_vertex_rank() for \
+                                         me in self._curr_matrix_elements.get_virt_matrix_elements()]
+
             card_path = os.path.join(path, os.path.pardir, 'SubProcesses', \
                                      'procdef_mg5.dat')
             
             if self.options['loop_optimized_output'] and \
-                    len(self._curr_matrix_elements.get_virt_matrix_elements()) > 0:
-                self._curr_exporter.write_coef_specs_file(\
-                        self._curr_matrix_elements.get_virt_matrix_elements())
+               len(max_loop_vertex_ranks) > 0:
+                self._curr_exporter.write_coef_specs_file(max_loop_vertex_ranks)
             if self._generate_info:
                 self._curr_exporter.write_procdef_mg5(card_path, #
                                 self._curr_model['name'],

=== modified file 'madgraph/interface/madgraph_interface.py'
--- madgraph/interface/madgraph_interface.py	2016-02-23 22:55:27 +0000
+++ madgraph/interface/madgraph_interface.py	2016-02-24 13:54:50 +0000
@@ -1332,7 +1332,8 @@
 
         if len(args) == 1 and args[0] in ['complex_mass_scheme',\
                                           'loop_optimized_output',\
-                                          'loop_color_flows']:
+                                          'loop_color_flows',\
+                                          'low_mem_multicore_nlo_generation']:
             args.append('True')
 
         if len(args) > 2 and '=' == args[1]:
@@ -1381,7 +1382,7 @@
             if not args[1].isdigit():
                 raise self.InvalidCmd('%s values should be a integer' % args[0])
 
-        if args[0] in ['loop_optimized_output', 'loop_color_flows']:
+        if args[0] in ['loop_optimized_output', 'loop_color_flows', 'low_mem_multicore_nlo_generation']:
             try:
                 args[1] = banner_module.ConfigFile.format_variable(args[1], bool, args[0])
             except Exception:
@@ -2381,7 +2382,8 @@
 
         if len(args) == 2:
             if args[1] in ['group_subprocesses', 'complex_mass_scheme',\
-                           'loop_optimized_output', 'loop_color_flows']:
+                           'loop_optimized_output', 'loop_color_flows',\
+                           'low_mem_multicore_nlo_generation']:
                 return self.list_completion(text, ['False', 'True', 'default'])
             elif args[1] in ['ignore_six_quark_processes']:
                 return self.list_completion(text, self._multiparticles.keys())
@@ -2698,6 +2700,7 @@
 
     options_madgraph= {'group_subprocesses': 'Auto',
                           'ignore_six_quark_processes': False,
+                          'low_mem_multicore_nlo_generation': False,
                           'complex_mass_scheme': False,
                           'gauge':'unitary',
                           'stdout_level':None,

=== modified file 'madgraph/iolibs/export_fks.py'
--- madgraph/iolibs/export_fks.py	2016-02-23 19:44:10 +0000
+++ madgraph/iolibs/export_fks.py	2016-02-24 13:54:50 +0000
@@ -55,6 +55,20 @@
 _file_path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0] + '/'
 logger = logging.getLogger('madgraph.export_fks')
 
+
+def make_jpeg_async(args):
+    Pdir = args[0]
+    old_pos = args[1]
+    dir_path = args[2]
+  
+    devnull = os.open(os.devnull, os.O_RDWR)
+  
+    os.chdir(Pdir)
+    subprocess.call([os.path.join(old_pos, dir_path, 'bin', 'internal', 'gen_jpeg-pl')],
+                    stdout = devnull)
+    os.chdir(os.path.pardir)  
+
+
 #=================================================================================
 # Class for used of the (non-optimized) Loop process
 #=================================================================================
@@ -294,12 +308,9 @@
     #===========================================================================
     # write_maxparticles_file
     #===========================================================================
-    def write_maxparticles_file(self, writer, matrix_elements):
+    def write_maxparticles_file(self, writer, maxparticles):
         """Write the maxparticles.inc file for MadEvent"""
 
-        maxparticles = max([me.get_nexternal_ninitial()[0] \
-                              for me in matrix_elements['matrix_elements']])
-
         lines = "integer max_particles, max_branch\n"
         lines += "parameter (max_particles=%d) \n" % maxparticles
         lines += "parameter (max_branch=max_particles-1)"
@@ -313,16 +324,9 @@
     #===========================================================================
     # write_maxconfigs_file
     #===========================================================================
-    def write_maxconfigs_file(self, writer, matrix_elements):
+    def write_maxconfigs_file(self, writer, maxconfigs):
         """Write the maxconfigs.inc file for MadEvent"""
 
-        try:
-            maxconfigs = max([me.get_num_configs() \
-                            for me in matrix_elements['real_matrix_elements']])
-        except ValueError:
-            maxconfigs = max([me.born_matrix_element.get_num_configs() \
-                            for me in matrix_elements['matrix_elements']])
-
         lines = "integer lmaxconfigs\n"
         lines += "parameter (lmaxconfigs=%d)" % maxconfigs
 
@@ -455,7 +459,7 @@
         # likely also generate it for each subproc.
         if OLP=='NJET':
             filename = 'OLE_order.lh'
-            self.write_lh_order(filename, matrix_element, OLP)
+            self.write_lh_order(filename, [matrix_element.born_matrix_element.get('processes')[0]], OLP)
         
         if matrix_element.virt_matrix_element:
                     calls += self.generate_virt_directory( \
@@ -683,14 +687,11 @@
     #===========================================================================
     #  create the run_card 
     #===========================================================================
-    def create_run_card(self, matrix_elements, history):
+    def create_run_card(self, processes, history):
         """ """
  
         run_card = banner_mod.RunCardNLO()
         
-        processes = [me.get('processes') 
-                                 for me in matrix_elements['matrix_elements']]
-        
         run_card.create_default_for_process(self.proc_characteristic, 
                                             history,
                                             processes)
@@ -709,7 +710,7 @@
         self.proc_characteristic['grouped_matrix'] = False
         self.create_proc_charac()
 
-        self.create_run_card(matrix_elements, history)
+        self.create_run_card(matrix_elements.get_processes(), history)
 #        modelname = self.model.get('name')
 #        if modelname == 'mssm' or modelname.startswith('mssm-'):
 #            param_card = os.path.join(self.dir_path, 'Cards','param_card.dat')
@@ -723,14 +724,15 @@
         self.write_get_mass_width_file(writers.FortranWriter(filename), makeinc, self.model)
 
 #        # Write maxconfigs.inc based on max of ME's/subprocess groups
+
         filename = os.path.join(self.dir_path,'Source','maxconfigs.inc')
         self.write_maxconfigs_file(writers.FortranWriter(filename),
-                                   matrix_elements)
+                                   matrix_elements.get_max_configs())
         
 #        # Write maxparticles.inc based on max of ME's/subprocess groups
         filename = os.path.join(self.dir_path,'Source','maxparticles.inc')
         self.write_maxparticles_file(writers.FortranWriter(filename),
-                                     matrix_elements)
+                                     matrix_elements.get_max_particles())
 
         # Touch "done" file
         os.system('touch %s/done' % os.path.join(self.dir_path,'SubProcesses'))
@@ -1341,7 +1343,8 @@
                          matrix_element, i,
                          fortran_model)
 
-    def generate_virtuals_from_OLP(self,FKSHMultiproc,export_path, OLP):
+
+    def generate_virtuals_from_OLP(self,process_list,export_path, OLP):
         """Generates the library for computing the loop matrix elements
         necessary for this process using the OLP specified."""
         
@@ -1350,7 +1353,7 @@
         if not os.path.exists(virtual_path):
             os.makedirs(virtual_path)
         filename = os.path.join(virtual_path,'OLE_order.lh')
-        self.write_lh_order(filename, FKSHMultiproc.get('matrix_elements'),OLP)
+        self.write_lh_order(filename, process_list, OLP)
 
         fail_msg='Generation of the virtuals with %s failed.\n'%OLP+\
             'Please check the virt_generation.log file in %s.'\
@@ -1419,20 +1422,19 @@
         proc_to_label = self.parse_contract_file(
                                             pjoin(virtual_path,'OLE_order.olc'))
 
-        self.write_BinothLHA_inc(FKSHMultiproc,proc_to_label,\
+        self.write_BinothLHA_inc(process_list,proc_to_label,\
                                               pjoin(export_path,'SubProcesses'))
         
         # Link the contract file to within the SubProcess directory
         ln(pjoin(virtual_path,'OLE_order.olc'),pjoin(export_path,'SubProcesses'))
         
-    def write_BinothLHA_inc(self, FKSHMultiproc, proc_to_label, SubProcPath):
+    def write_BinothLHA_inc(self, processes, proc_to_label, SubProcPath):
         """ Write the file Binoth_proc.inc in each SubProcess directory so as 
         to provide the right process_label to use in the OLP call to get the
         loop matrix element evaluation. The proc_to_label is the dictionary of
         the format of the one returned by the function parse_contract_file."""
         
-        for matrix_element in FKSHMultiproc.get('matrix_elements'):
-            proc = matrix_element.get('processes')[0]
+        for proc in processes:
             name = "P%s"%proc.shell_string()
             proc_pdgs=(tuple([leg.get('id') for leg in proc.get('legs') if \
                                                          not leg.get('state')]),
@@ -1608,26 +1610,19 @@
     # write_lh_order
     #===============================================================================
     #test written
-    def write_lh_order(self, filename, matrix_elements, OLP='MadLoop'):
+    def write_lh_order(self, filename, process_list, OLP='MadLoop'):
         """Creates the OLE_order.lh file. This function should be edited according
         to the OLP which is used. For now it is generic."""
         
-        if isinstance(matrix_elements,fks_helas_objects.FKSHelasProcess):
-            fksborns=fks_helas_objects.FKSHelasProcessList([matrix_elements])
-        elif isinstance(matrix_elements,fks_helas_objects.FKSHelasProcessList):
-            fksborns= matrix_elements
-        else:
-            raise fks_common.FKSProcessError('Wrong type of argument for '+\
-                                  'matrix_elements in function write_lh_order.')
         
-        if len(fksborns)==0:
+        if len(process_list)==0:
             raise fks_common.FKSProcessError('No matrix elements provided to '+\
                                                  'the function write_lh_order.')
             return
         
         # We assume the orders to be common to all Subprocesses
         
-        orders = fksborns[0].orders 
+        orders = process_list[0].get('orders') 
         if 'QED' in orders.keys() and 'QCD' in orders.keys():
             QED=orders['QED']
             QCD=orders['QCD']
@@ -1639,12 +1634,12 @@
             QCD=orders['QCD']
         else:
             QED, QCD = self.get_qed_qcd_orders_from_weighted(\
-                    fksborns[0].get_nexternal_ninitial()[0]-1, # -1 is because the function returns nexternal of the real emission
+                    len(process_list[0].get('legs')),
                     orders['WEIGHTED'])
 
         replace_dict = {}
         replace_dict['mesq'] = 'CHaveraged'
-        replace_dict['corr'] = ' '.join(matrix_elements[0].get('processes')[0].\
+        replace_dict['corr'] = ' '.join(process_list[0].\
                                                   get('perturbation_couplings'))
         replace_dict['irreg'] = 'CDR'
         replace_dict['aspow'] = QCD
@@ -1652,8 +1647,10 @@
         replace_dict['modelfile'] = './param_card.dat'
         replace_dict['params'] = 'alpha_s'
         proc_lines=[]
-        for fksborn in fksborns:
-            proc_lines.append(fksborn.get_lh_pdg_string())
+        for proc in process_list:
+            proc_lines.append('%s -> %s' % \
+                    (' '.join(str(l['id']) for l in proc['legs'] if not l['state']),
+                     ' '.join(str(l['id']) for l in proc['legs'] if l['state'])))
         replace_dict['pdgs'] = '\n'.join(proc_lines)
         replace_dict['symfin'] = 'Yes'
         content = \
@@ -3300,7 +3297,7 @@
     #===============================================================================
     # write_coef_specs
     #===============================================================================
-    def write_coef_specs_file(self, virt_me_list):
+    def write_coef_specs_file(self, max_loop_vertex_ranks):
         """ writes the coef_specs.inc in the DHELAS folder. Should not be called in the 
         non-optimized mode"""
         filename = os.path.join(self.dir_path, 'Source', 'DHELAS', 'coef_specs.inc')
@@ -3308,7 +3305,6 @@
         replace_dict = {}
         replace_dict['max_lwf_size'] = 4 
 
-        max_loop_vertex_ranks = [me.get_max_loop_vertex_rank() for me in virt_me_list]
         replace_dict['vertex_max_coefs'] = max(\
                 [q_polynomial.get_number_of_coefs_for_rank(n) 
                     for n in max_loop_vertex_ranks])

=== modified file 'madgraph/iolibs/template_files/loop_optimized/compute_color_flows.inc'
--- madgraph/iolibs/template_files/loop_optimized/compute_color_flows.inc	2015-07-03 00:47:03 +0000
+++ madgraph/iolibs/template_files/loop_optimized/compute_color_flows.inc	2016-02-24 13:54:50 +0000
@@ -357,7 +357,7 @@
 C  
 C LOCAL VARIABLES 
 C
-	  INTEGER I, J, M, N, K, SQSOINDEX, DOUBLEFACT
+	  INTEGER I, J, M, N, K, ISQSO, DOUBLEFACT
 	  %(complex_dp_format)s COLOR_COEF
 C  
 C FUNCTIONS
@@ -401,8 +401,8 @@
   IF (LoopColorFlowMatrix(I,I)%%Denom.LT.0) COLOR_COEF=COLOR_COEF*IMAG1
   DO M=1,NLOOPAMPSO
     DO N=M,NLOOPAMPSO
-      SQSOINDEX = %(proc_prefix)sML5SQSOINDEX(M,N)
-      IF(CHOSEN_SO_CONFIGS(SQSOINDEX)) THEN
+      ISQSO = %(proc_prefix)sML5SQSOINDEX(M,N)
+      IF(CHOSEN_SO_CONFIGS(ISQSO)) THEN
         IF(M.ne.N) THEN
           DOUBLEFACT=2 
         ELSE 
@@ -452,7 +452,7 @@
 C  
 C LOCAL VARIABLES 
 C
-	  INTEGER I, J, M, N, K, SQSOINDEX, DOUBLEFACT
+	  INTEGER I, J, M, N, K, ISQSO, DOUBLEFACT
 ## if(LoopInduced) {
 	  INTEGER LOWERBOUND
 ## }
@@ -550,15 +550,15 @@
       	CALL %(proc_prefix)sML5GET_ORDERS_FOR_AMPSOINDEX(M,ORDERS_A)
       	CALL %(proc_prefix)sML5GET_ORDERS_FOR_AMPSOINDEX(N,ORDERS_B)
 C Now figure out to which SQSOINDEX these orders together correspond to *in the Born subroutine* (Notice the absence of the ML5 prefix in the functions used then)
-        SQSOINDEX = %(proc_prefix)sSQSOINDEX(%(proc_prefix)sSOINDEX_FOR_AMPORDERS(ORDERS_A),%(proc_prefix)sSOINDEX_FOR_AMPORDERS(ORDERS_B))
+        ISQSO = %(proc_prefix)sSQSOINDEX(%(proc_prefix)sSOINDEX_FOR_AMPORDERS(ORDERS_A),%(proc_prefix)sSOINDEX_FOR_AMPORDERS(ORDERS_B))
         IF(J.ne.I) THEN
           DOUBLEFACT=2 
         ELSE 
           DOUBLEFACT=1
         ENDIF
         TEMP(1) = DOUBLEFACT*HEL_MULT*DBLE(COLOR_COEF*JAMPB(I,M)*DCONJG(JAMPB(J,N)))
-        RES(0,SQSOINDEX) = RES(0,SQSOINDEX) + TEMP(1)
-        IF((.not.FILTER_SO).or.SQSO_TARGET.eq.-1.or.SQSO_TARGET.eq.SQSOINDEX) THEN
+        RES(0,ISQSO) = RES(0,ISQSO) + TEMP(1)
+        IF((.not.FILTER_SO).or.SQSO_TARGET.eq.-1.or.SQSO_TARGET.eq.ISQSO) THEN
           RES(0,0) = RES(0,0) + TEMP(1)
         ENDIF
       ENDDO
@@ -591,7 +591,7 @@
 	  ENDIF
       DO N=LOWERBOUND,NLOOPAMPSO
 ## }
-        SQSOINDEX = %(proc_prefix)sML5SQSOINDEX(M,N)
+        ISQSO = %(proc_prefix)sML5SQSOINDEX(M,N)
 ## if(LoopInduced){
         IF(J.ne.I) THEN
           DOUBLEFACT=2 
@@ -626,14 +626,14 @@
 ## if(not LoopInduced){      
           TEMP(K) = 2.0d0*HEL_MULT*DBLE(COLOR_COEF*JAMPL(K,I,M)*DCONJG(JAMPB(J,N)))
 ## }
-          RES(K,SQSOINDEX) = RES(K,SQSOINDEX) + TEMP(K)
+          RES(K,ISQSO) = RES(K,ISQSO) + TEMP(K)
 ## if(LoopInduced and MadEventOutput){
 		  DO config_i=1,nconfigs
-            AMP2_ALL(K,config_i,SQSOINDEX) = AMP2_ALL(K,config_i,SQSOINDEX) + TEMP_AMP2(K,config_i)
+            AMP2_ALL(K,config_i,ISQSO) = AMP2_ALL(K,config_i,ISQSO) + TEMP_AMP2(K,config_i)
 		  ENDDO
 ## }
         ENDDO
-        IF((.NOT.FILTER_SO).OR.SQSO_TARGET.eq.-1.or.SQSO_TARGET.eq.SQSOINDEX) THEN
+        IF((.NOT.FILTER_SO).OR.SQSO_TARGET.eq.-1.or.SQSO_TARGET.eq.ISQSO) THEN
           DO K=1,3
             RES(K,0) = RES(K,0) + TEMP(K)
 ## if(LoopInduced and MadEventOutput){

=== modified file 'madgraph/loop/loop_exporters.py'
--- madgraph/loop/loop_exporters.py	2016-02-23 22:55:27 +0000
+++ madgraph/loop/loop_exporters.py	2016-02-24 13:54:50 +0000
@@ -45,6 +45,7 @@
 import madgraph.various.diagram_symmetry as diagram_symmetry
 import madgraph.various.process_checks as process_checks
 import madgraph.various.progressbar as pbar
+import madgraph.various.q_polynomial as q_polynomial
 import madgraph.core.color_amp as color_amp
 import madgraph.iolibs.helas_call_writers as helas_call_writers
 import models.check_param_card as check_param_card
@@ -2314,7 +2315,7 @@
 
         writer.writelines(file,context=self.get_context(matrix_element))
     
-    def fix_coef_specs(self, overall_max_lwf_size, overall_max_loop_vert_rank):
+    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
         the file 'coef.inc' in the HELAS Source folder must contain the overall
@@ -2325,7 +2326,11 @@
         coef_specs_path=os.path.join(self.dir_path,'Source','DHELAS',\
                                                                'coef_specs.inc')
         os.remove(coef_specs_path)
-        
+       
+        spin_to_wf_size = {1:4,2:4,3:4,4:16,5:16}
+        overall_max_lwf_size = spin_to_wf_size[overall_max_lwf_spin]
+        overall_max_loop_vert_coefs = q_polynomial.get_number_of_coefs_for_rank(
+                                                     overall_max_loop_vert_rank)
         # Replace it by the appropriate value
         IncWriter=writers.FortranWriter(coef_specs_path,'w')
         IncWriter.writelines("""INTEGER MAXLWFSIZE
@@ -2333,7 +2338,7 @@
                            INTEGER VERTEXMAXCOEFS
                            PARAMETER (VERTEXMAXCOEFS=%(vertex_max_coefs)d)"""\
                            %{'max_lwf_size':overall_max_lwf_size,
-                             'vertex_max_coefs':overall_max_loop_vert_rank})
+                             'vertex_max_coefs':overall_max_loop_vert_coefs})
         IncWriter.close()
 
     def setup_check_sa_replacement_dictionary(self, matrix_element, \

=== modified file 'madgraph/various/banner.py'
--- madgraph/various/banner.py	2016-02-18 11:42:56 +0000
+++ madgraph/various/banner.py	2016-02-24 13:54:50 +0000
@@ -1990,7 +1990,7 @@
         # check for beam_id
         beam_id = set()
         for proc in proc_def:
-            for leg in proc[0]['legs']:
+            for leg in proc['legs']:
                 if not leg['state']:
                     beam_id.add(leg['id'])
         if any(i in beam_id for i in [1,-1,2,-2,3,-3,4,-4,5,-5,21,22]):

=== added directory 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall'
=== added 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	1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%CT_interface.f	2016-02-24 13:54:50 +0000
@@ -0,0 +1,278 @@
+      SUBROUTINE CTLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
+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 CutTools.
+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))
+C     
+C     ARGUMENTS 
+C     
+      INTEGER NLOOPLINE, RANK
+      REAL*8 PL(0:3,NLOOPLINE)
+      REAL*8 PCT(0:3,0:NLOOPLINE-1)
+      COMPLEX*16 M2L(NLOOPLINE)
+      COMPLEX*16 M2LCT(0:NLOOPLINE-1)
+      COMPLEX*16 RES(3)
+      LOGICAL STABLE
+C     
+C     LOCAL VARIABLES 
+C     
+      COMPLEX*16 R1, ACC
+      INTEGER I, J, K
+      LOGICAL CTINIT, TIRINIT, GOLEMINIT
+      COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT
+C     
+C     EXTERNAL FUNCTIONS
+C     
+      EXTERNAL LOOPNUM
+      EXTERNAL MPLOOPNUM
+C     
+C     GLOBAL VARIABLES
+C     
+      INCLUDE 'coupl.inc'
+      INTEGER CTMODE
+      REAL*8 LSCALE
+      COMMON/CT/LSCALE,CTMODE
+
+      INTEGER ID,SQSOINDEX,R
+      COMMON/LOOP/ID,SQSOINDEX,R
+
+C     ----------
+C     BEGIN CODE
+C     ----------
+
+C     INITIALIZE CUTTOOLS IF NEEDED
+      IF (CTINIT) THEN
+        CTINIT=.FALSE.
+        CALL INITCT()
+      ENDIF
+
+C     YOU CAN FIND THE DETAILS ABOUT THE DIFFERENT CTMODE AT THE
+C      BEGINNING OF THE FILE CTS_CUTS.F90 IN THE CUTTOOLS DISTRIBUTION
+
+C     CONVERT THE MASSES TO BE COMPLEX
+      DO I=1,NLOOPLINE
+        M2LCT(I-1)=M2L(I)
+      ENDDO
+
+C     CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO CT CONVENTIONS
+      DO I=0,3
+        DO J=0,(NLOOPLINE-1)
+          PCT(I,J)=0.D0
+        ENDDO
+      ENDDO
+      DO I=0,3
+        DO J=1,NLOOPLINE
+          PCT(I,0)=PCT(I,0)+PL(I,J)
+        ENDDO
+      ENDDO
+      IF (CHECKPCONSERVATION) THEN
+        IF (PCT(0,0).GT.1.D-6) THEN
+          WRITE(*,*) 'energy is not conserved ',PCT(0,0)
+          STOP 'energy is not conserved'
+        ELSEIF (PCT(1,0).GT.1.D-6) THEN
+          WRITE(*,*) 'px is not conserved ',PCT(1,0)
+          STOP 'px is not conserved'
+        ELSEIF (PCT(2,0).GT.1.D-6) THEN
+          WRITE(*,*) 'py is not conserved ',PCT(2,0)
+          STOP 'py is not conserved'
+        ELSEIF (PCT(3,0).GT.1.D-6) THEN
+          WRITE(*,*) 'pz is not conserved ',PCT(3,0)
+          STOP 'pz is not conserved'
+        ENDIF
+      ENDIF
+      DO I=0,3
+        DO J=1,(NLOOPLINE-1)
+          DO K=1,J
+            PCT(I,J)=PCT(I,J)+PL(I,K)
+          ENDDO
+        ENDDO
+      ENDDO
+
+      CALL CTSXCUT(CTMODE,LSCALE,MU_R,NLOOPLINE,LOOPNUM,MPLOOPNUM,RANK
+     $ ,PCT,M2LCT,RES,ACC,R1,STABLE)
+      RES(1)=NORMALIZATION*2.0D0*DBLE(RES(1))
+      RES(2)=NORMALIZATION*2.0D0*DBLE(RES(2))
+      RES(3)=NORMALIZATION*2.0D0*DBLE(RES(3))
+C     WRITE(*,*) 'Loop ID',ID,' =',RES(1),RES(2),RES(3)
+      END
+
+      SUBROUTINE INITCT()
+C     
+C     INITIALISATION OF CUTTOOLS
+C     
+C     LOCAL VARIABLES 
+C     
+      REAL*8 THRS
+      LOGICAL EXT_NUM_FOR_R1
+C     
+C     GLOBAL VARIABLES 
+C     
+      INCLUDE 'MadLoopParams.inc'
+C     ----------
+C     BEGIN CODE
+C     ----------
+
+C     DEFAULT PARAMETERS FOR CUTTOOLS
+C     -------------------------------  
+C     THRS1 IS THE PRECISION LIMIT BELOW WHICH THE MP ROUTINES
+C      ACTIVATES
+      THRS=CTSTABTHRES
+C     LOOPLIB SET WHAT LIBRARY CT USES
+C     1 -> LOOPTOOLS
+C     2 -> AVH
+C     3 -> QCDLOOP
+      LOOPLIB=CTLOOPLIBRARY
+C     MADLOOP'S NUMERATOR IN THE OPEN LOOP IS MUCH FASTER THAN THE
+C      RECONSTRUCTED ONE IN CT. SO WE BETTER USE MADLOOP ONE IN THIS
+C      CASE.
+      EXT_NUM_FOR_R1=.TRUE.
+C     -------------------------------	  
+
+C     The initialization below is for CT v1.8.+
+      CALL CTSINIT(THRS,LOOPLIB,EXT_NUM_FOR_R1)
+C     The initialization below is for the older stable CT v1.7, still
+C      used for now in the beta release.
+C     CALL CTSINIT(THRS,LOOPLIB)
+
+      END
+
+      SUBROUTINE LOOP_3(W1, W2, W3, M1, M2, M3,  RANK, SQUAREDSOINDEX
+     $ , LOOPNUM)
+      INTEGER    NEXTERNAL
+      PARAMETER (NEXTERNAL=3)
+      INTEGER    NLOOPLINE
+      PARAMETER (NLOOPLINE=3)
+      INTEGER    NWAVEFUNCS
+      PARAMETER (NWAVEFUNCS=3)
+      INTEGER    NLOOPGROUPS
+      PARAMETER (NLOOPGROUPS=1)
+      INTEGER    NCOMB
+      PARAMETER (NCOMB=12)
+C     These are constants related to the split orders
+      INTEGER    NSQUAREDSO
+      PARAMETER (NSQUAREDSO=1)
+C     
+C     ARGUMENTS 
+C     
+      INTEGER W1, W2, W3
+      COMPLEX*16 M1, M2, M3
+
+      INTEGER RANK, LSYMFACT
+      INTEGER LOOPNUM, SQUAREDSOINDEX
+C     
+C     LOCAL VARIABLES 
+C     
+      REAL*8 PL(0:3,NLOOPLINE)
+      COMPLEX*16 M2L(NLOOPLINE)
+      INTEGER PAIRING(NLOOPLINE),WE(3)
+      INTEGER I, J, K, TEMP,I_LIB
+      LOGICAL COMPLEX_MASS,DOING_QP
+C     
+C     GLOBAL VARIABLES
+C     
+      INCLUDE 'MadLoopParams.inc'
+      INTEGER ID,SQSOINDEX,R
+      COMMON/LOOP/ID,SQSOINDEX,R
+
+      LOGICAL CHECKPHASE, HELDOUBLECHECKED
+      COMMON/INIT/CHECKPHASE, HELDOUBLECHECKED
+
+      INTEGER HELOFFSET
+      INTEGER GOODHEL(NCOMB)
+      LOGICAL GOODAMP(NSQUAREDSO,NLOOPGROUPS)
+      COMMON/FILTERS/GOODAMP,GOODHEL,HELOFFSET
+
+      COMPLEX*16 LOOPRES(3,NSQUAREDSO,NLOOPGROUPS)
+      LOGICAL S(NSQUAREDSO,NLOOPGROUPS)
+      COMMON/LOOPRES/LOOPRES,S
+
+
+      COMPLEX*16 W(20,NWAVEFUNCS)
+      COMMON/W/W
+      REAL*8 LSCALE
+      INTEGER CTMODE
+      COMMON/CT/LSCALE,CTMODE
+      INTEGER LIBINDEX
+      COMMON/I_LIB/LIBINDEX
+
+C     ----------
+C     BEGIN CODE
+C     ----------
+
+      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.GOODAMP(SQUAREDSOIND
+     $ EX,LOOPNUM)) THEN
+        WE(1)=W1
+        WE(2)=W2
+        WE(3)=W3
+        M2L(1)=M3**2
+        M2L(2)=M1**2
+        M2L(3)=M2**2
+        DO I=1,NLOOPLINE
+          PAIRING(I)=1
+        ENDDO
+
+        R=RANK
+        ID=LOOPNUM
+        SQSOINDEX=SQUAREDSOINDEX
+        DO I=0,3
+          TEMP=1
+          DO J=1,NLOOPLINE
+            PL(I,J)=0.D0
+            DO K=TEMP,(TEMP+PAIRING(J)-1)
+              PL(I,J)=PL(I,J)-DBLE(W(1+I,WE(K)))
+            ENDDO
+            TEMP=TEMP+PAIRING(J)
+          ENDDO
+        ENDDO
+C       Determine whether the integral is with complex masses or not
+C       since some reduction libraries, e.g.PJFry++ and IREGI, are
+C        still
+C       not able to deal with complex masses
+        COMPLEX_MASS=.FALSE.
+        DO I=1,NLOOPLINE
+          IF(DIMAG(M2L(I)).EQ.0D0)CYCLE
+          IF(ABS(DIMAG(M2L(I)))/MAX(ABS(M2L(I)),1D-2).GT.1D-15)THEN
+            COMPLEX_MASS=.TRUE.
+            EXIT
+          ENDIF
+        ENDDO
+C       Determine it uses qp or not
+        DOING_QP=.FALSE.
+        IF(CTMODE.GE.4)DOING_QP=.TRUE.
+C       Choose the correct loop library
+        CALL CHOOSE_LOOPLIB(LIBINDEX,NLOOPLINE,RANK,COMPLEX_MASS
+     $   ,DOING_QP,I_LIB)
+        IF(MLREDUCTIONLIB(I_LIB).EQ.1)THEN
+C         USE CUTTOOLS
+          CALL CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX
+     $     ,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
+        ELSE
+C         USE TIR
+          CALL TIRLOOP(SQUAREDSOINDEX,LOOPNUM,I_LIB,NLOOPLINE,PL,M2L
+     $     ,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX
+     $     ,LOOPNUM))
+        ENDIF
+      ELSE
+        LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
+        LOOPRES(2,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
+        LOOPRES(3,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
+        S(SQUAREDSOINDEX,LOOPNUM)=.TRUE.
+      ENDIF
+      END
+

=== added 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	1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%GOLEM_interface.f	2016-02-24 13:54:50 +0000
@@ -0,0 +1,760 @@
+      SUBROUTINE GOLEMLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
+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 Golem95.
+C     The Golem95 version should be higher than 1.3.0.
+C     It supports RANK = NLOOPLINE + 1 tensor integrals when 1 <
+C      NLOOPLINE < 6.
+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 MATRICE_S
+      USE FORM_FACTOR_TYPE, ONLY: FORM_FACTOR
+      USE PRECISION_GOLEM, ONLY: KI
+      USE TENS_COMB
+      USE TENS_REC
+      USE FORM_FACTOR_1P, ONLY: A10
+      USE FORM_FACTOR_2P, ONLY: A20
+      USE FORM_FACTOR_3P, ONLY: A30
+      USE FORM_FACTOR_4P, ONLY: A40
+      USE FORM_FACTOR_5P, ONLY: A50
+      USE FORM_FACTOR_6P, ONLY: A60
+      IMPLICIT NONE
+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))
+      REAL(KI),DIMENSION(0:3),PARAMETER::NULL_VEC = (/0.0_KI,0.0_KI
+     $ ,0.0_KI,0.0_KI/)
+C     GOLEM_RUN_MODE = 1: Use directly MadLoop tensorial coefficients
+C     GOLEM_RUN_MODE = 2: Reconstruct the tensorial coefficeints
+C      directly from 
+C     numerator using golem internal reconstruction routine
+C     GOLEM_RUN_MODE = 3: Cross-checked reconstructed coefficients
+C      against
+C     MadLoop internal ones.
+      INTEGER GOLEM_RUN_MODE
+      PARAMETER (GOLEM_RUN_MODE=1)
+C     The following is the acceptance threshold used for GOLEM_RUN_MODE
+C      = 3
+      REAL*8 COEF_CHECK_THRS
+      DATA COEF_CHECK_THRS/1.0D-13/
+      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
+
+      LOGICAL PASS_COEF_CHECK
+C     
+C     ARGUMENTS 
+C     
+      INTEGER NLOOPLINE, RANK
+      REAL*8 PL(0:3,NLOOPLINE)
+      REAL*8 PCT(0:3,0:NLOOPLINE-1)
+      REAL(KI) PGOLEM(NLOOPLINE,0:3)
+      COMPLEX*16 M2L(NLOOPLINE)
+      COMPLEX(KI) M2LGOLEM(NLOOPLINE)
+      COMPLEX*16 RES(3)
+      LOGICAL STABLE
+C     
+C     LOCAL VARIABLES 
+C     
+      INTEGER I, J, K
+      TYPE(FORM_FACTOR)::RES_GOLEM
+
+      COMPLEX(KI)::COEFFS0,COEFFS0_REC
+      TYPE(COEFF_TYPE_1)::COEFFS1,COEFFS1_REC
+      TYPE(COEFF_TYPE_2)::COEFFS2,COEFFS2_REC
+      TYPE(COEFF_TYPE_3)::COEFFS3,COEFFS3_REC
+      TYPE(COEFF_TYPE_4)::COEFFS4,COEFFS4_REC
+      TYPE(COEFF_TYPE_5)::COEFFS5,COEFFS5_REC
+      TYPE(COEFF_TYPE_6)::COEFFS6,COEFFS6_REC
+
+C     The pinch propagator optimization is not used, so for now it is
+C     always 0.
+      INTEGER PINCH
+C     
+C     EXTERNAL FUNCTIONS
+C     
+      COMPLEX(KI) GOLEM_LOOPNUM
+      EXTERNAL GOLEM_LOOPNUM
+      LOGICAL COMPARE_COEFS_0
+      LOGICAL COMPARE_COEFS_1
+      LOGICAL COMPARE_COEFS_2
+      LOGICAL COMPARE_COEFS_3
+      LOGICAL COMPARE_COEFS_4
+      LOGICAL COMPARE_COEFS_5
+      LOGICAL COMPARE_COEFS_6
+C     
+C     GLOBAL VARIABLES
+C     
+      INCLUDE 'coupl.inc'
+      INTEGER CTMODE
+      REAL*8 LSCALE
+      COMMON/CT/LSCALE,CTMODE
+
+      INTEGER ID,SQSOINDEX,R
+      COMMON/LOOP/ID,SQSOINDEX,R
+
+      LOGICAL CTINIT, TIRINIT, GOLEMINIT
+      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT,GOLEMINIT
+
+      INTEGER NLOOPGROUPS
+      PARAMETER (NLOOPGROUPS=1)
+      INTEGER LOOPMAXCOEFS
+      PARAMETER (LOOPMAXCOEFS=15)
+      INTEGER NSQUAREDSO
+      PARAMETER (NSQUAREDSO=1)
+
+      COMPLEX*16 LOOPCOEFS(0:LOOPMAXCOEFS-1,NSQUAREDSO,NLOOPGROUPS)
+      COMMON/LCOEFS/LOOPCOEFS
+C     ----------
+C     BEGIN CODE
+C     ----------
+
+C     INITIALIZE CUTTOOLS IF NEEDED
+      IF (GOLEMINIT) THEN
+        GOLEMINIT=.FALSE.
+        CALL INITGOLEM()
+      ENDIF
+
+C     No stability test intrisic to Golem95 now
+      STABLE=.TRUE.
+
+C     This initialization must be done for each reduction because we
+C     have not setup anyoptimization using pinched propagators yet.
+      CALL INITGOLEM95(NLOOPLINE)
+      PINCH = 0
+
+C     YOU CAN FIND THE DETAILS ABOUT THE DIFFERENT CTMODE AT THE
+C      BEGINNING OF THE FILE CTS_CUTS.F90 IN THE CUTTOOLS DISTRIBUTION
+
+C     CONVERT THE MASSES TO BE COMPLEX
+      DO I=1,NLOOPLINE
+        M2LGOLEM(I)=M2L(I)
+      ENDDO
+
+C     CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO CT CONVENTIONS
+      DO I=0,3
+        DO J=0,(NLOOPLINE-1)
+          PCT(I,J)=0.D0
+        ENDDO
+      ENDDO
+      DO I=0,3
+        DO J=1,NLOOPLINE
+          PCT(I,0)=PCT(I,0)+PL(I,J)
+        ENDDO
+      ENDDO
+      IF (CHECKPCONSERVATION) THEN
+        IF (PCT(0,0).GT.1.D-6) THEN
+          WRITE(*,*) 'energy is not conserved ',PCT(0,0)
+          STOP 'energy is not conserved'
+        ELSEIF (PCT(1,0).GT.1.D-6) THEN
+          WRITE(*,*) 'px is not conserved ',PCT(1,0)
+          STOP 'px is not conserved'
+        ELSEIF (PCT(2,0).GT.1.D-6) THEN
+          WRITE(*,*) 'py is not conserved ',PCT(2,0)
+          STOP 'py is not conserved'
+        ELSEIF (PCT(3,0).GT.1.D-6) THEN
+          WRITE(*,*) 'pz is not conserved ',PCT(3,0)
+          STOP 'pz is not conserved'
+        ENDIF
+      ENDIF
+      DO I=0,3
+        DO J=1,(NLOOPLINE-1)
+          DO K=1,J
+            PCT(I,J)=PCT(I,J)+PL(I,K)
+          ENDDO
+        ENDDO
+      ENDDO
+
+C     Now convert the loop momenta to Golem95 conventions
+      DO I=0,3
+        PGOLEM(1,I)=0.0E0_KI
+        DO J=2,NLOOPLINE
+          PGOLEM(J,I)=PCT(I,J-1)
+        ENDDO
+      ENDDO
+
+C     Fill in the kinematic s-matrix while taking care of on-shell
+C      limits.
+      CALL SETUP_KIN_MATRIX(NLOOPLINE,PGOLEM,M2LGOLEM)
+C     Construct the golem internal matrices derived from the kinetic
+C      one.
+      CALL PREPARESMATRIX()
+
+C     Fill in the golem coefficents and compute the loop
+      IF(GOLEM_RUN_MODE.EQ.2)THEN
+        RES_GOLEM = EVALUATE_B(GOLEM_LOOPNUM,PGOLEM,0,RANK)
+      ELSE
+        PASS_COEF_CHECK=.TRUE.
+        SELECT CASE(RANK)
+        CASE(0)
+        CALL FILL_GOLEM_COEFFS_0(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS0)
+        IF(GOLEM_RUN_MODE.EQ.3)THEN
+          COEFFS0_REC = GOLEM_LOOPNUM(NULL_VEC,0.0_KI)
+          PASS_COEF_CHECK=COMPARE_COEFS_0(COEFFS0,COEFFS0_REC)
+        ENDIF
+        CASE(1)
+        CALL FILL_GOLEM_COEFFS_1(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS1)
+        IF(GOLEM_RUN_MODE.EQ.3)THEN
+          CALL RECONSTRUCT1(GOLEM_LOOPNUM,COEFFS1_REC)
+          PASS_COEF_CHECK=COMPARE_COEFS_1(COEFFS1,COEFFS1_REC)
+        ENDIF
+        CASE(2)
+        CALL FILL_GOLEM_COEFFS_2(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS2)
+        IF(GOLEM_RUN_MODE.EQ.3)THEN
+          CALL RECONSTRUCT2(GOLEM_LOOPNUM,COEFFS2_REC)
+          PASS_COEF_CHECK=COMPARE_COEFS_2(COEFFS2,COEFFS2_REC)
+        ENDIF
+        CASE(3)
+        CALL FILL_GOLEM_COEFFS_3(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS3)
+        IF(GOLEM_RUN_MODE.EQ.3)THEN
+          CALL RECONSTRUCT3(GOLEM_LOOPNUM,COEFFS3_REC)
+          PASS_COEF_CHECK=COMPARE_COEFS_3(COEFFS3,COEFFS3_REC)
+        ENDIF
+        CASE(4)
+        CALL FILL_GOLEM_COEFFS_4(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS4)
+        IF(GOLEM_RUN_MODE.EQ.3)THEN
+          CALL RECONSTRUCT4(GOLEM_LOOPNUM,COEFFS4_REC)
+          PASS_COEF_CHECK=COMPARE_COEFS_4(COEFFS4,COEFFS4_REC)
+        ENDIF
+        CASE(5)
+        CALL FILL_GOLEM_COEFFS_5(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS5)
+        IF(GOLEM_RUN_MODE.EQ.3)THEN
+          CALL RECONSTRUCT5(GOLEM_LOOPNUM,COEFFS5_REC)
+          PASS_COEF_CHECK=COMPARE_COEFS_5(COEFFS5,COEFFS5_REC)
+        ENDIF
+        CASE(6)
+        CALL FILL_GOLEM_COEFFS_6(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS6)
+        IF(GOLEM_RUN_MODE.EQ.3)THEN
+          CALL RECONSTRUCT6(GOLEM_LOOPNUM,COEFFS6_REC)
+          PASS_COEF_CHECK=COMPARE_COEFS_6(COEFFS6,COEFFS6_REC)
+        ENDIF
+        CASE DEFAULT
+        WRITE(*,*)'Not yet implemented in Golem95 for rank= ',RANK
+        STOP
+        END SELECT
+
+        IF(.NOT.PASS_COEF_CHECK)THEN
+          WRITE(*,*)'Coefs mismatch for ID ',ID,' and rank ',RANK
+          WRITE(*,*)'Coefs form MadLoop5:'
+          SELECT CASE(RANK)
+          CASE(0)
+          WRITE(*,*)'Constant coef = ',COEFFS0
+          CASE(1)
+          CALL PRINT_COEFFS(COEFFS1)
+          CASE(2)
+          CALL PRINT_COEFFS(COEFFS2)
+          CASE(3)
+          CALL PRINT_COEFFS(COEFFS3)
+          CASE(4)
+          CALL PRINT_COEFFS(COEFFS4)
+          CASE(5)
+          CALL PRINT_COEFFS(COEFFS5)
+          CASE(6)
+          CALL PRINT_COEFFS(COEFFS6)
+          END SELECT
+          WRITE(*,*)'Coefs reconstructed by Golem95:'
+          SELECT CASE(RANK)
+          CASE(0)
+          WRITE(*,*)'Constant coef = ',COEFFS0_REC
+          CASE(1)
+          CALL PRINT_COEFFS(COEFFS1_REC)
+          CASE(2)
+          CALL PRINT_COEFFS(COEFFS2_REC)
+          CASE(3)
+          CALL PRINT_COEFFS(COEFFS3_REC)
+          CASE(4)
+          CALL PRINT_COEFFS(COEFFS4_REC)
+          CASE(5)
+          CALL PRINT_COEFFS(COEFFS5_REC)
+          CASE(6)
+          CALL PRINT_COEFFS(COEFFS6_REC)
+          END SELECT
+          STOP
+        ENDIF
+
+        SELECT CASE(NLOOPLINE)
+        CASE(1)
+        WRITE(*,*)'Golem95 cannot handle with tadpole yet'
+        STOP
+        CASE(2)
+        SELECT CASE(RANK)
+        CASE(0)
+        RES_GOLEM = COEFFS0*A20(PINCH)
+        CASE(1)
+        RES_GOLEM = CONTRACT2_1(COEFFS1,PGOLEM,PINCH)
+        CASE(2)
+        RES_GOLEM = CONTRACT2_2(COEFFS2,PGOLEM,PINCH)
+        CASE(3)
+        RES_GOLEM = CONTRACT2_3(COEFFS3,PGOLEM,PINCH)
+        CASE DEFAULT
+        WRITE(*,*)'Golem95 cannot handle with: N,r = ',2,RANK
+        STOP
+        END SELECT
+        CASE(3)
+        SELECT CASE(RANK)
+        CASE(0)
+        RES_GOLEM = COEFFS0*A30(PINCH)
+        CASE(1)
+        RES_GOLEM = CONTRACT3_1(COEFFS1,PGOLEM,PINCH)
+        CASE(2)
+        RES_GOLEM = CONTRACT3_2(COEFFS2,PGOLEM,PINCH)
+        CASE(3)
+        RES_GOLEM = CONTRACT3_3(COEFFS3,PGOLEM,PINCH)
+        CASE(4)
+        RES_GOLEM = CONTRACT3_4(COEFFS4,PGOLEM,PINCH)
+        CASE DEFAULT
+        WRITE(*,*)'Golem95 cannot handle with: N,r = ',3,RANK
+        STOP
+        END SELECT
+        CASE(4)
+        SELECT CASE(RANK)
+        CASE(0)
+        RES_GOLEM = COEFFS0*A40(PINCH)
+        CASE(1)
+        RES_GOLEM = CONTRACT4_1(COEFFS1,PGOLEM,PINCH)
+        CASE(2)
+        RES_GOLEM = CONTRACT4_2(COEFFS2,PGOLEM,PINCH)
+        CASE(3)
+        RES_GOLEM = CONTRACT4_3(COEFFS3,PGOLEM,PINCH)
+        CASE(4)
+        RES_GOLEM = CONTRACT4_4(COEFFS4,PGOLEM,PINCH)
+        CASE(5)
+        RES_GOLEM = CONTRACT4_5(COEFFS5,PGOLEM,PINCH)
+        CASE DEFAULT
+        WRITE(*,*)'Golem95 cannot handle with: N,r = ',4,RANK
+        STOP
+        END SELECT
+        CASE(5)
+        SELECT CASE(RANK)
+        CASE(0)
+        RES_GOLEM = COEFFS0*A50(PINCH)
+        CASE(1)
+        RES_GOLEM = CONTRACT5_1(COEFFS1,PGOLEM,PINCH)
+        CASE(2)
+        RES_GOLEM = CONTRACT5_2(COEFFS2,PGOLEM,PINCH)
+        CASE(3)
+        RES_GOLEM = CONTRACT5_3(COEFFS3,PGOLEM,PINCH)
+        CASE(4)
+        RES_GOLEM = CONTRACT5_4(COEFFS4,PGOLEM,PINCH)
+        CASE(5)
+        RES_GOLEM = CONTRACT5_5(COEFFS5,PGOLEM,PINCH)
+        CASE(6)
+        RES_GOLEM = CONTRACT5_6(COEFFS6,PGOLEM,PINCH)
+        CASE DEFAULT
+        WRITE(*,*)'Golem95 cannot handle with: N,r = ',5,RANK
+        STOP
+        END SELECT
+        CASE(6)
+        SELECT CASE(RANK)
+        CASE(0)
+        RES_GOLEM = COEFFS0*A60(PINCH)
+        CASE(1)
+        RES_GOLEM = CONTRACT6_1(COEFFS1,PGOLEM,PINCH)
+        CASE(2)
+        RES_GOLEM = CONTRACT6_2(COEFFS2,PGOLEM,PINCH)
+        CASE(3)
+        RES_GOLEM = CONTRACT6_3(COEFFS3,PGOLEM,PINCH)
+        CASE(4)
+        RES_GOLEM = CONTRACT6_4(COEFFS4,PGOLEM,PINCH)
+        CASE(5)
+        RES_GOLEM = CONTRACT6_5(COEFFS5,PGOLEM,PINCH)
+        CASE(6)
+        RES_GOLEM = CONTRACT6_6(COEFFS6,PGOLEM,PINCH)
+        CASE DEFAULT
+        WRITE(*,*)'Golem95 cannot handle with: N,r = ',6,RANK
+        STOP
+        END SELECT
+        CASE DEFAULT
+        WRITE(*,*)'Golem95 cannot handle with: N = ',NLOOPLINE
+        STOP
+        END SELECT
+      ENDIF
+
+      RES(1)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%C+2.0*LOG(MU_R)
+     $ *RES_GOLEM%%B+2.0*LOG(MU_R)**2*RES_GOLEM%%A)
+      RES(2)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%B+2.0*LOG(MU_R)
+     $ *RES_GOLEM%%A)
+      RES(3)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%A)
+C     WRITE(*,*) 'Loop ID',ID,' =',RES(1),RES(2),RES(3)
+
+C     Finally free golem memory and cache
+      CALL EXITGOLEM95()
+
+      END
+
+      FUNCTION GOLEM_LOOPNUM(Q,MU2)
+      USE PRECISION_GOLEM, ONLY: KI
+      REAL(KI),DIMENSION(0:3),INTENT(IN)::Q
+      REAL(KI),INTENT(IN)::MU2
+      COMPLEX(KI)::GOLEM_LOOPNUM
+
+      COMPLEX*16 QQ(0:3),NUM
+      INTEGER I
+
+      DO I=0,3
+        QQ(I)=CMPLX(Q(I),0.0D0,KIND=16)
+      ENDDO
+
+      CALL LOOPNUM(QQ,NUM)
+      GOLEM_LOOPNUM=NUM
+      RETURN
+      END FUNCTION
+
+      SUBROUTINE INITGOLEM()
+C     
+C     INITIALISATION OF GOLEM
+C     
+C     
+C     MODULE
+C     
+      USE PARAMETRE
+C     
+C     LOCAL VARIABLES 
+C     
+      REAL*8 THRS
+      LOGICAL EXT_NUM_FOR_R1
+C     
+C     GLOBAL VARIABLES 
+C     
+      INCLUDE 'MadLoopParams.inc'
+C     ----------
+C     BEGIN CODE
+C     ----------
+
+C     DEFAULT PARAMETERS FOR GOLEM
+C     -------------------------------  
+C     One can chose here to have either just the rational R1 piece 
+C     or everything but the R2
+      RAT_OR_TOT_PAR = TOT
+
+      END
+
+      SUBROUTINE SETUP_KIN_MATRIX(NLOOPLINE,PGOLEM,M2L)
+C     
+C     MODULE
+C     
+      USE MATRICE_S
+      USE PRECISION_GOLEM, ONLY: KI
+C     
+C     ARGUMENTS
+C     
+      INTEGER NLOOPLINE
+      REAL(KI) PGOLEM(NLOOPLINE,0:3)
+      COMPLEX(KI) M2L(NLOOPLINE)
+
+C     
+C     GLOBAL VARIABLES
+C     
+      INCLUDE 'MadLoopParams.inc'
+      INTEGER CTMODE
+      REAL*8 LSCALE
+      COMMON/CT/LSCALE,CTMODE
+      REAL*8 REF_NORMALIZATION
+C     
+C     LOCAL VARIABLES
+C     
+      INTEGER I,J,K
+      COMPLEX(KI) DIFFSQ
+C     ----------
+C     BEGIN CODE
+C     ----------
+
+      DO I=1,NLOOPLINE
+        DO J=1,NLOOPLINE
+          IF(I.EQ.J)THEN
+            S_MAT(I,J)=-(M2L(I)+M2L(J))
+          ELSE
+            DIFFSQ = (CMPLX(PGOLEM(I,0),0.0_KI,KIND=KI)-CMPLX(PGOLEM(J
+     $       ,0),0.0_KI,KIND=KI))**2
+            DO K=1,3
+              DIFFSQ = DIFFSQ-(CMPLX(PGOLEM(I,K),0.0_KI,KIND=KI)
+     $         -CMPLX(PGOLEM(J,K),0.0_KI,KIND=KI))**2
+            ENDDO
+            S_MAT(I,J)=DIFFSQ-M2L(I)-M2L(J)
+            IF(M2L(I).NE.0.0D0)THEN
+              IF(ABS((DIFFSQ-M2L(I))/M2L(I)).LT.OSTHRES)THEN
+                S_MAT(I,J)=-M2L(J)
+              ENDIF
+            ENDIF
+            IF(M2L(J).NE.0.0D0)THEN
+              IF(ABS((DIFFSQ-M2L(J))/M2L(J)).LT.OSTHRES)THEN
+                S_MAT(I,J)=-M2L(I)
+              ENDIF
+            ENDIF
+C           Chose what seems the most appropriate way to compare
+C           massless onshellness. (here we pick the energy component)
+            REF_NORMALIZATION=(PGOLEM(I,0)+PGOLEM(J,0))**2
+            IF(REF_NORMALIZATION.NE.0.0D0)THEN
+              IF(ABS(DIFFSQ/REF_NORMALIZATION).LT.OSTHRES)THEN
+                S_MAT(I,J)=-(M2L(I)+M2L(J))
+              ENDIF
+            ENDIF
+          ENDIF
+        ENDDO
+      ENDDO
+
+      END
+
+      FUNCTION COMPARE_COEFS_0(COEFS_A,COEFS_B)
+
+      USE PRECISION_GOLEM, ONLY: KI
+      COMPLEX(KI) COEFS_A,COEFS_B
+      REAL*8 COEF_CHECK_THRS
+      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
+      REAL*8 DENOM,NUM
+      LOGICAL COMPARE_COEFS_0
+
+      NUM = ABS(COEFS_A-COEFS_B)
+      DENOM = ABS(COEFS_A+COEFS_B)
+      IF(DENOM.GT.0D0)THEN
+        COMPARE_COEFS_0=((NUM/DENOM).LT.COEF_CHECK_THRS)
+      ELSE
+        COMPARE_COEFS_0=(NUM.LT.COEF_CHECK_THRS)
+      ENDIF
+
+      END
+
+      FUNCTION COMPARE_COEFS_1(COEFS_A,COEFS_B)
+
+      USE TENS_REC
+      TYPE(COEFF_TYPE_1)COEFS_A,COEFS_B
+      REAL*8 COEF_CHECK_THRS
+      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
+      REAL*8 DENOM,NUM
+      LOGICAL COMPARE_COEFS_1
+
+      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
+      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1
+     $ ))
+
+      IF(DENOM.GT.0D0)THEN
+        COMPARE_COEFS_1=((NUM/DENOM).LT.COEF_CHECK_THRS)
+      ELSE
+        COMPARE_COEFS_1=(NUM.LT.COEF_CHECK_THRS)
+      ENDIF
+
+      END
+
+      FUNCTION COMPARE_COEFS_2(COEFS_A,COEFS_B)
+
+      USE TENS_REC
+      TYPE(COEFF_TYPE_2) COEFS_A,COEFS_B
+      REAL*8 COEF_CHECK_THRS
+      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
+      REAL*8 DENOM,NUM
+      LOGICAL COMPARE_COEFS_2
+
+      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
+     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))
+      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1
+     $ ))+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))
+      IF(DENOM.GT.0D0)THEN
+        COMPARE_COEFS_2=((NUM/DENOM).LT.COEF_CHECK_THRS)
+      ELSE
+        COMPARE_COEFS_2=(NUM.LT.COEF_CHECK_THRS)
+      ENDIF
+
+      END
+
+      FUNCTION COMPARE_COEFS_3(COEFS_A,COEFS_B)
+
+      USE TENS_REC
+      TYPE(COEFF_TYPE_3) COEFS_A, COEFS_B
+      REAL*8 COEF_CHECK_THRS
+      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
+      REAL*8 DENOM,NUM
+      LOGICAL COMPARE_COEFS_3
+
+      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
+     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
+      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1
+     $ ))+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3
+     $ ))
+      IF(DENOM.GT.0D0)THEN
+        COMPARE_COEFS_3=((NUM/DENOM).LT.COEF_CHECK_THRS)
+      ELSE
+        COMPARE_COEFS_3=(NUM.LT.COEF_CHECK_THRS)
+      ENDIF
+
+      END
+
+      FUNCTION COMPARE_COEFS_4(COEFS_A,COEFS_B)
+
+      USE TENS_REC
+      TYPE(COEFF_TYPE_4) COEFS_A, COEFS_B
+      REAL*8 COEF_CHECK_THRS
+      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
+      REAL*8 DENOM,NUM
+      LOGICAL COMPARE_COEFS_4
+
+      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
+     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3)
+     $ )+SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
+      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1
+     $ ))+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3
+     $ ))+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
+      IF(DENOM.GT.0D0)THEN
+        COMPARE_COEFS_4=((NUM/DENOM).LT.COEF_CHECK_THRS)
+      ELSE
+        COMPARE_COEFS_4=(NUM.LT.COEF_CHECK_THRS)
+      ENDIF
+
+      END
+
+      FUNCTION COMPARE_COEFS_5(COEFS_A,COEFS_B)
+
+      USE TENS_REC
+      TYPE(COEFF_TYPE_5) COEFS_A,COEFS_B
+      REAL*8 COEF_CHECK_THRS
+      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
+      REAL*8 DENOM,NUM
+      LOGICAL COMPARE_COEFS_5
+
+      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
+     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3)
+     $ )+SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
+      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1
+     $ ))+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3
+     $ ))+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
+      IF(DENOM.GT.0D0)THEN
+        COMPARE_COEFS_5=((NUM/DENOM).LT.COEF_CHECK_THRS)
+      ELSE
+        COMPARE_COEFS_5=(NUM.LT.COEF_CHECK_THRS)
+      ENDIF
+
+      END
+
+      FUNCTION COMPARE_COEFS_6(COEFS_A,COEFS_B)
+
+      USE TENS_REC
+      TYPE(COEFF_TYPE_6) COEFS_A,COEFS_B
+      REAL*8 COEF_CHECK_THRS
+      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
+      REAL*8 DENOM,NUM
+      LOGICAL COMPARE_COEFS_6
+
+      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
+     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3)
+     $ )+SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
+      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1
+     $ ))+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3
+     $ ))+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
+      IF(DENOM.GT.0D0)THEN
+        COMPARE_COEFS_6=((NUM/DENOM).LT.COEF_CHECK_THRS)
+      ELSE
+        COMPARE_COEFS_6=(NUM.LT.COEF_CHECK_THRS)
+      ENDIF
+
+      END
+
+
+      SUBROUTINE FILL_GOLEM_COEFFS_0(ML_COEFS,GOLEM_COEFS)
+      USE PRECISION_GOLEM, ONLY: KI
+      INCLUDE 'coef_specs.inc'
+      COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
+      COMPLEX(KI) GOLEM_COEFS
+      GOLEM_COEFS=ML_COEFS(0)
+      END
+
+      SUBROUTINE FILL_GOLEM_COEFFS_1(ML_COEFS,GOLEM_COEFS)
+      USE TENS_REC, ONLY: COEFF_TYPE_1
+      INCLUDE 'coef_specs.inc'
+      COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
+      TYPE(COEFF_TYPE_1) GOLEM_COEFS
+C     Constant coefficient 
+      GOLEM_COEFS%%C0=ML_COEFS(0)
+C     Coefficient q(0)
+      GOLEM_COEFS%%C1(1,1)=-ML_COEFS(1)
+C     Coefficient q(1)
+      GOLEM_COEFS%%C1(2,1)=-ML_COEFS(2)
+C     Coefficient q(2)
+      GOLEM_COEFS%%C1(3,1)=-ML_COEFS(3)
+C     Coefficient q(3)
+      GOLEM_COEFS%%C1(4,1)=-ML_COEFS(4)
+      END
+
+      SUBROUTINE FILL_GOLEM_COEFFS_2(ML_COEFS,GOLEM_COEFS)
+      USE TENS_REC, ONLY: COEFF_TYPE_2
+      INCLUDE 'coef_specs.inc'
+      COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
+      TYPE(COEFF_TYPE_2) GOLEM_COEFS
+C     Constant coefficient 
+      GOLEM_COEFS%%C0=ML_COEFS(0)
+C     Coefficient q(0)
+      GOLEM_COEFS%%C1(1,1)=-ML_COEFS(1)
+C     Coefficient q(0)^2
+      GOLEM_COEFS%%C1(1,2)= ML_COEFS(5)
+C     Coefficient q(1)
+      GOLEM_COEFS%%C1(2,1)=-ML_COEFS(2)
+C     Coefficient q(1)^2
+      GOLEM_COEFS%%C1(2,2)= ML_COEFS(7)
+C     Coefficient q(2)
+      GOLEM_COEFS%%C1(3,1)=-ML_COEFS(3)
+C     Coefficient q(2)^2
+      GOLEM_COEFS%%C1(3,2)= ML_COEFS(10)
+C     Coefficient q(3)
+      GOLEM_COEFS%%C1(4,1)=-ML_COEFS(4)
+C     Coefficient q(3)^2
+      GOLEM_COEFS%%C1(4,2)= ML_COEFS(14)
+C     Coefficient q(0)*q(1)
+      GOLEM_COEFS%%C2(1,1)= ML_COEFS(6)
+C     Coefficient q(0)*q(2)
+      GOLEM_COEFS%%C2(2,1)= ML_COEFS(8)
+C     Coefficient q(0)*q(3)
+      GOLEM_COEFS%%C2(3,1)= ML_COEFS(11)
+C     Coefficient q(1)*q(2)
+      GOLEM_COEFS%%C2(4,1)= ML_COEFS(9)
+C     Coefficient q(1)*q(3)
+      GOLEM_COEFS%%C2(5,1)= ML_COEFS(12)
+C     Coefficient q(2)*q(3)
+      GOLEM_COEFS%%C2(6,1)= ML_COEFS(13)
+      END
+
+      SUBROUTINE FILL_GOLEM_COEFFS_3(ML_COEFS,GOLEM_COEFS)
+      USE TENS_REC, ONLY: COEFF_TYPE_3
+      INCLUDE 'coef_specs.inc'
+      COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
+      TYPE(COEFF_TYPE_3) GOLEM_COEFS
+C     Dummy routine for FILL_GOLEM_COEFS_3
+      STOP 'ERROR: 3 > 2'
+      END
+
+      SUBROUTINE FILL_GOLEM_COEFFS_4(ML_COEFS,GOLEM_COEFS)
+      USE TENS_REC, ONLY: COEFF_TYPE_4
+      INCLUDE 'coef_specs.inc'
+      COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
+      TYPE(COEFF_TYPE_4) GOLEM_COEFS
+C     Dummy routine for FILL_GOLEM_COEFS_4
+      STOP 'ERROR: 4 > 2'
+      END
+
+      SUBROUTINE FILL_GOLEM_COEFFS_5(ML_COEFS,GOLEM_COEFS)
+      USE TENS_REC, ONLY: COEFF_TYPE_5
+      INCLUDE 'coef_specs.inc'
+      COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
+      TYPE(COEFF_TYPE_5) GOLEM_COEFS
+C     Dummy routine for FILL_GOLEM_COEFS_5
+      STOP 'ERROR: 5 > 2'
+      END
+
+      SUBROUTINE FILL_GOLEM_COEFFS_6(ML_COEFS,GOLEM_COEFS)
+      USE TENS_REC, ONLY: COEFF_TYPE_6
+      INCLUDE 'coef_specs.inc'
+      COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
+      TYPE(COEFF_TYPE_6) GOLEM_COEFS
+C     Dummy routine for FILL_GOLEM_COEFS_6
+      STOP 'ERROR: 6 > 2'
+      END

=== added 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	1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f	2016-02-24 13:54:50 +0000
@@ -0,0 +1,423 @@
+      SUBROUTINE TIRLOOP(I_SQSO,I_LOOPGROUP,I_LIB,NLOOPLINE,PL,M2L
+     $ ,RANK,RES,STABLE)
+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 TIR.
+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 NLOOPGROUPS
+      PARAMETER (NLOOPGROUPS=1)
+C     These are constants related to the split orders
+      INTEGER NSQUAREDSO
+      PARAMETER (NSQUAREDSO=1)
+      INTEGER LOOPMAXCOEFS
+      PARAMETER (LOOPMAXCOEFS=15)
+      INTEGER    NEXTERNAL
+      PARAMETER (NEXTERNAL=3)
+      LOGICAL CHECKPCONSERVATION
+      PARAMETER (CHECKPCONSERVATION=.TRUE.)
+      REAL*8 NORMALIZATION
+      PARAMETER (NORMALIZATION = 1.D0/(16.D0*3.14159265358979323846D0*
+     $ *2))
+C     
+C     ARGUMENTS 
+C     
+      INTEGER I_SQSO,I_LOOPGROUP,I_LIB
+      INTEGER NLOOPLINE, RANK
+      REAL*8 PL(0:3,NLOOPLINE)
+      REAL*8 PCT(0:3,0:NLOOPLINE-1)
+      REAL*8 PDEN(0:3,NLOOPLINE-1)
+      COMPLEX*16 M2L(NLOOPLINE)
+      COMPLEX*16 M2LCT(0:NLOOPLINE-1)
+      COMPLEX*16 RES(3)
+      LOGICAL STABLE
+C     
+C     LOCAL VARIABLES 
+C     
+      INTEGER I, J, K
+      INTEGER NLOOPCOEFS
+      LOGICAL CTINIT, TIRINIT, GOLEMINIT
+      COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT
+
+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 PJCOEFS(0:LOOPMAXCOEFS-1,3)
+C     
+C     EXTERNAL FUNCTIONS
+C     
+C     
+C     GLOBAL VARIABLES
+C     
+      INCLUDE 'MadLoopParams.inc'
+      INCLUDE 'coupl.inc'
+      INTEGER CTMODE
+      REAL*8 LSCALE
+      COMMON/CT/LSCALE,CTMODE
+
+C     The argument ILIB is the TIR library to be used for that
+C      specific library.
+      INTEGER LIBINDEX
+      COMMON/I_LIB/LIBINDEX
+
+      COMPLEX*16 LOOPCOEFS(0:LOOPMAXCOEFS-1,NSQUAREDSO,NLOOPGROUPS)
+      COMMON/LCOEFS/LOOPCOEFS
+C     ----------
+C     BEGIN CODE
+C     ----------
+
+C     Initialize for the very first time ML is called LAST_ILIB with
+C      the first ILIB used.
+      IF(LAST_LIB_USED.EQ.-1) THEN
+        LAST_LIB_USED = MLREDUCTIONLIB(LIBINDEX)
+      ELSE
+C       We changed the TIR library so we must refresh the cache.
+        IF(MLREDUCTIONLIB(LIBINDEX).NE.LAST_LIB_USED) THEN
+          LAST_LIB_USED = MLREDUCTIONLIB(LIBINDEX)
+          CALL CLEAR_TIR_CACHE()
+        ENDIF
+      ENDIF
+
+      IF (MLREDUCTIONLIB(I_LIB).EQ.4) THEN
+C       Using Golem95
+C       PDEN is dummy for Golem95 so we just initialize it to zero
+C        here so as to use it for the function SWITCHORDER
+        DO I=0,3
+          DO J=1,NLOOPLINE-1
+            PDEN(I,J)=0.0D0
+          ENDDO
+        ENDDO
+        CALL SWITCH_ORDER(CTMODE,NLOOPLINE,PL,PDEN,M2L)
+        CALL GOLEMLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
+        RETURN
+      ENDIF
+
+C     INITIALIZE TIR IF NEEDED
+      IF (TIRINIT) THEN
+        TIRINIT=.FALSE.
+        CALL INITTIR()
+      ENDIF
+
+C     CONVERT THE MASSES TO BE COMPLEX
+      DO I=1,NLOOPLINE
+        M2LCT(I-1)=M2L(I)
+      ENDDO
+
+C     CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO CT CONVENTIONS
+      DO I=0,3
+        DO J=0,(NLOOPLINE-1)
+          PCT(I,J)=0.D0
+        ENDDO
+      ENDDO
+      DO I=0,3
+        DO J=1,NLOOPLINE
+          PCT(I,0)=PCT(I,0)+PL(I,J)
+        ENDDO
+      ENDDO
+      IF (CHECKPCONSERVATION) THEN
+        IF (PCT(0,0).GT.1.D-6) THEN
+          WRITE(*,*) 'energy is not conserved ',PCT(0,0)
+          STOP 'energy is not conserved'
+        ELSEIF (PCT(1,0).GT.1.D-6) THEN
+          WRITE(*,*) 'px is not conserved ',PCT(1,0)
+          STOP 'px is not conserved'
+        ELSEIF (PCT(2,0).GT.1.D-6) THEN
+          WRITE(*,*) 'py is not conserved ',PCT(2,0)
+          STOP 'py is not conserved'
+        ELSEIF (PCT(3,0).GT.1.D-6) THEN
+          WRITE(*,*) 'pz is not conserved ',PCT(3,0)
+          STOP 'pz is not conserved'
+        ENDIF
+      ENDIF
+      DO I=0,3
+        DO J=1,(NLOOPLINE-1)
+          DO K=1,J
+            PCT(I,J)=PCT(I,J)+PL(I,K)
+          ENDDO
+        ENDDO
+      ENDDO
+
+      DO I=0,3
+        DO J=1,(NLOOPLINE-1)
+          PDEN(I,J)=PCT(I,J)
+        ENDDO
+      ENDDO
+C     NUMBER OF INDEPEDENT LOOPCOEFS FOR RANK=RANK
+      NLOOPCOEFS=0
+      DO I=0,RANK
+        NLOOPCOEFS=NLOOPCOEFS+(3+I)*(2+I)*(1+I)/6
+      ENDDO
+      SELECT CASE(MLREDUCTIONLIB(I_LIB))
+      CASE(2)
+C     PJFry++
+      WRITE(*,*) 'ERROR:: PJFRY++ is not interfaced.'
+      STOP
+      CASE(3)
+C     IREGI
+      CALL IMLOOP(CTMODE,IREGIMODE,NLOOPLINE,LOOPMAXCOEFS,RANK,PDEN
+     $ ,M2L,MU_R,PJCOEFS,STABLE)
+C     CONVERT TO MADLOOP CONVENTION
+      CALL CONVERT_IREGI_COEFFS(RANK,PJCOEFS,TIRCOEFS)
+      END SELECT
+      DO I=1,3
+        RES(I)=(0.0D0,0.0D0)
+        DO J=0,NLOOPCOEFS-1
+          RES(I)=RES(I)+LOOPCOEFS(J,I_SQSO,I_LOOPGROUP)*TIRCOEFS(J,I)
+        ENDDO
+      ENDDO
+      RES(1)=NORMALIZATION*2.0D0*DBLE(RES(1))
+      RES(2)=NORMALIZATION*2.0D0*DBLE(RES(2))
+      RES(3)=NORMALIZATION*2.0D0*DBLE(RES(3))
+C     WRITE(*,*) 'Loop ID',ID,' =',RES(1),RES(2),RES(3)
+      END
+
+      SUBROUTINE SWITCH_ORDER(CTMODE,NLOOPLINE,PL,PDEN,M2L)
+      IMPLICIT NONE
+
+      INTEGER CTMODE,NLOOPLINE
+
+      REAL*8 PL(0:3,NLOOPLINE)
+      REAL*8 PDEN(0:3,NLOOPLINE-1)
+      COMPLEX*16 M2L(NLOOPLINE)
+      REAL*8 NEW_PL(0:3,NLOOPLINE)
+      REAL*8 NEW_PDEN(0:3,NLOOPLINE-1)
+      COMPLEX*16 NEW_M2L(NLOOPLINE)
+
+      INTEGER I,J,K
+
+      IF (CTMODE.NE.2.AND.CTMODE.NE.4) THEN
+        RETURN
+      ENDIF
+
+      IF (NLOOPLINE.LE.2) THEN
+        RETURN
+      ENDIF
+
+      DO I=1,NLOOPLINE-1
+        DO J=0,3
+          NEW_PDEN(J,NLOOPLINE-I) = PDEN(J,I)
+        ENDDO
+      ENDDO
+      DO I=1,NLOOPLINE-1
+        DO J=0,3
+          PDEN(J,I) = NEW_PDEN(J,I)
+        ENDDO
+      ENDDO
+
+      DO I=2,NLOOPLINE
+        NEW_M2L(I)=M2L(NLOOPLINE-I+2)
+      ENDDO
+      DO I=2,NLOOPLINE
+        M2L(I)=NEW_M2L(I)
+      ENDDO
+
+
+      DO I=1,NLOOPLINE
+        DO J=0,3
+          NEW_PL(J,I) = -PL(J,NLOOPLINE+1-I)
+        ENDDO
+      ENDDO
+      DO I=1,NLOOPLINE
+        DO J=0,3
+          PL(J,I) = NEW_PL(J,I)
+        ENDDO
+      ENDDO
+
+      END
+
+      SUBROUTINE INITTIR()
+C     
+C     INITIALISATION OF TIR
+C     
+C     LOCAL VARIABLES 
+C     
+      REAL*8 THRS
+      LOGICAL EXT_NUM_FOR_R1
+C     
+C     GLOBAL VARIABLES 
+C     
+      INCLUDE 'MadLoopParams.inc'
+      LOGICAL CTINIT, TIRINIT, GOLEMINIT
+      COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT
+
+C     ----------
+C     BEGIN CODE
+C     ----------
+
+C     DEFAULT PARAMETERS FOR TIR
+C     -------------------------------  
+C     THRS1 IS THE PRECISION LIMIT BELOW WHICH THE MP ROUTINES
+C      ACTIVATES
+C     USE THE SAME MADLOOP PARAMETER IN CUTTOOLS AND TIR
+C     IT IS NECESSARY TO INITIALIZE CT BECAUSE IREGI USES THE VERSION
+C      OF ONELOOP
+C     FROM CUTTOOLS LIBRARY
+      THRS=CTSTABTHRES
+C     LOOPLIB SET WHAT LIBRARY CT USES
+C     1 -> LOOPTOOLS
+C     2 -> AVH
+C     3 -> QCDLOOP
+      LOOPLIB=CTLOOPLIBRARY
+      CALL INITIREGI(IREGIRECY,LOOPLIB,1D-6)
+C     The initialization below is for CT v1.9.2+
+      IF (CTINIT) THEN
+        CTINIT=.FALSE.
+        CALL INITCT()
+      ENDIF
+      END
+
+      SUBROUTINE CHOOSE_LOOPLIB(LIBINDEX,NLOOPLINE,RANK,COMPLEX_MASS
+     $ ,DOING_QP,I_LIB)
+C     
+C     CHOOSE THE CORRECT LOOP LIB
+C     Example:
+C     MLReductionLib=3|2|1 and LIBINDEX=3
+C     IF THE LOOP IS BEYOND THE SCOPE OF LOOP LIB MLReductionLib(3)=1
+C     USE LIBINDEX=1, and LIBINDEX=2 ...
+C     IF IT IS STILL NOT GOOD,STOP
+C     
+      IMPLICIT NONE
+C     
+C     CONSTANTS
+C     
+      INTEGER NLOOPLIB,QP_NLOOPLIB
+      PARAMETER (NLOOPLIB=4,QP_NLOOPLIB=1)
+C     
+C     ARGUMENTS
+C     
+      INTEGER LIBINDEX,NLOOPLINE,RANK,I_LIB
+      LOGICAL COMPLEX_MASS,DOING_QP
+C     
+C     LOCAL VARIABLES
+C     
+      INTEGER I,J_LIB,LIBNUM,SELECT_LIBINDEX
+      LOGICAL LPASS
+C     
+C     GLOBAL VARIABLES
+C     
+      INCLUDE 'MadLoopParams.inc'
+C     TILL NOW, ONLY CUTTOOLS PROVIDE QP
+      LOGICAL QP_TOOLS_AVAILABLE
+      INTEGER INDEX_QP_TOOLS(QP_NLOOPLIB)
+      COMMON/LOOP_TOOLS/QP_TOOLS_AVAILABLE,INDEX_QP_TOOLS
+C     ----------
+C     BEGIN CODE
+C     ----------
+
+      IF(DOING_QP)THEN
+C       QP EVALUATION, ONLY CUTTOOLS
+        IF(.NOT.QP_TOOLS_AVAILABLE)THEN
+          STOP 'No qp tools available, please make sure MLReductionLi'
+     $     //'b is correct'
+        ENDIF
+        J_LIB=0
+        SELECT_LIBINDEX=LIBINDEX
+        DO WHILE(J_LIB.EQ.0)
+          DO I=1,QP_NLOOPLIB
+            IF(INDEX_QP_TOOLS(I).EQ.SELECT_LIBINDEX)THEN
+              J_LIB=I
+              EXIT
+            ENDIF
+          ENDDO
+          IF(J_LIB.EQ.0)THEN
+            SELECT_LIBINDEX=SELECT_LIBINDEX+1
+            IF(SELECT_LIBINDEX.GT.NLOOPLIB.OR.MLREDUCTIONLIB(SELECT_LIB
+     $       INDEX).EQ.0)SELECT_LIBINDEX=1
+          ENDIF
+        ENDDO
+        I=J_LIB
+        I_LIB=SELECT_LIBINDEX
+        LIBNUM=MLREDUCTIONLIB(I_LIB)
+        DO
+        CALL DETECT_LOOPLIB(LIBNUM,NLOOPLINE,RANK,COMPLEX_MASS,LPASS)
+        IF(LPASS)EXIT
+        I=I+1
+        IF(I.GT.QP_NLOOPLIB.AND.INDEX_QP_TOOLS(I).EQ.0)THEN
+          I=1
+        ENDIF
+        IF(I.EQ.J_LIB)THEN
+          STOP 'No qp loop library can deal with this integral'
+        ENDIF
+        I_LIB=INDEX_QP_TOOLS(I)
+        LIBNUM=MLREDUCTIONLIB(I_LIB)
+        ENDDO
+      ELSE
+C       DP EVALUATION
+        I_LIB=LIBINDEX
+        LIBNUM=MLREDUCTIONLIB(I_LIB)
+        DO
+        CALL DETECT_LOOPLIB(LIBNUM,NLOOPLINE,RANK,COMPLEX_MASS,LPASS)
+        IF(LPASS)EXIT
+        I_LIB=I_LIB+1
+        IF(I_LIB.GT.NLOOPLIB.OR.MLREDUCTIONLIB(I_LIB).EQ.0)THEN
+          I_LIB=1
+        ENDIF
+        IF(I_LIB.EQ.LIBINDEX)THEN
+          STOP 'No dp loop library can deal with this integral'
+        ENDIF
+        LIBNUM=MLREDUCTIONLIB(I_LIB)
+        ENDDO
+      ENDIF
+      RETURN
+      END
+
+      SUBROUTINE CLEAR_TIR_CACHE()
+C     No TIR caching implemented, this is dummy. (The subroutine is
+C      kept as it might be called by the MC).
+      CONTINUE
+      END SUBROUTINE
+
+
+
+      SUBROUTINE CONVERT_IREGI_COEFFS(RANK,IREGICOEFS,TIRCOEFS)
+C     GLOABLE VARIABLES
+      INCLUDE 'coef_specs.inc'
+C     ARGUMENTS
+      INTEGER RANK
+      COMPLEX*16 IREGICOEFS(0:LOOP_MAXCOEFS-1,3)
+      COMPLEX*16 TIRCOEFS(0:LOOP_MAXCOEFS-1,3)
+C     Reduction Coefficient 1
+      TIRCOEFS(0,1:3)=IREGICOEFS(0,1:3)
+      IF(RANK.LE.0)RETURN
+C     Reduction Coefficient q(0)
+      TIRCOEFS(1,1:3)=IREGICOEFS(1,1:3)
+C     Reduction Coefficient q(1)
+      TIRCOEFS(2,1:3)=IREGICOEFS(2,1:3)
+C     Reduction Coefficient q(2)
+      TIRCOEFS(3,1:3)=IREGICOEFS(3,1:3)
+C     Reduction Coefficient q(3)
+      TIRCOEFS(4,1:3)=IREGICOEFS(4,1:3)
+      IF(RANK.LE.1)RETURN
+C     Reduction Coefficient q(0)^2
+      TIRCOEFS(5,1:3)=IREGICOEFS(5,1:3)
+C     Reduction Coefficient q(0)*q(1)
+      TIRCOEFS(6,1:3)=IREGICOEFS(6,1:3)
+C     Reduction Coefficient q(1)^2
+      TIRCOEFS(7,1:3)=IREGICOEFS(9,1:3)
+C     Reduction Coefficient q(0)*q(2)
+      TIRCOEFS(8,1:3)=IREGICOEFS(7,1:3)
+C     Reduction Coefficient q(1)*q(2)
+      TIRCOEFS(9,1:3)=IREGICOEFS(10,1:3)
+C     Reduction Coefficient q(2)^2
+      TIRCOEFS(10,1:3)=IREGICOEFS(12,1:3)
+C     Reduction Coefficient q(0)*q(3)
+      TIRCOEFS(11,1:3)=IREGICOEFS(8,1:3)
+C     Reduction Coefficient q(1)*q(3)
+      TIRCOEFS(12,1:3)=IREGICOEFS(11,1:3)
+C     Reduction Coefficient q(2)*q(3)
+      TIRCOEFS(13,1:3)=IREGICOEFS(13,1:3)
+C     Reduction Coefficient q(3)^2
+      TIRCOEFS(14,1:3)=IREGICOEFS(14,1:3)
+      IF(RANK.LE.2)RETURN
+      END

=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f	1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f	2016-02-24 13:54:50 +0000
@@ -0,0 +1,540 @@
+      SUBROUTINE SMATRIX(P,ANS_SUMMED)
+C     
+C     Simple routine wrapper to provide the same interface for
+C     backward compatibility for usage without split orders.
+C     
+C     
+C     CONSTANTS
+C     
+      INTEGER    NEXTERNAL
+      PARAMETER (NEXTERNAL=3)
+      INTEGER NSQAMPSO
+      PARAMETER (NSQAMPSO=1)
+C     
+C     ARGUMENTS 
+C     
+      REAL*8 P(0:3,NEXTERNAL), ANS_SUMMED
+C     
+C     VARIABLES
+C     
+      INTEGER I
+      REAL*8 ANS(0:NSQAMPSO)
+C     
+C     BEGIN CODE
+C     
+      CALL SMATRIX_SPLITORDERS(P,ANS)
+      ANS_SUMMED=ANS(0)
+
+      END
+
+      SUBROUTINE SMATRIXHEL(P,HEL,ANS)
+      IMPLICIT NONE
+C     
+C     CONSTANT
+C     
+      INTEGER    NEXTERNAL
+      PARAMETER (NEXTERNAL=3)
+      INTEGER                 NCOMB
+      PARAMETER (             NCOMB=12)
+C     
+C     ARGUMENTS 
+C     
+      REAL*8 P(0:3,NEXTERNAL),ANS
+      INTEGER HEL
+C     
+C     GLOBAL VARIABLES
+C     
+      INTEGER USERHEL
+      COMMON/HELUSERCHOICE/USERHEL
+C     ----------
+C     BEGIN CODE
+C     ----------
+      USERHEL=HEL
+      CALL SMATRIX(P,ANS)
+      USERHEL=-1
+
+      END
+
+      SUBROUTINE SMATRIX_SPLITORDERS(P,ANS)
+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     MadGraph StandAlone Version
+C     
+C     Returns amplitude squared summed/avg over colors
+C     and helicities
+C     for the point in phase space P(0:3,NEXTERNAL)
+C     
+C     Process: d~ u > w+ WEIGHTED=2 QED=1 [ all = QCD ]
+C     Process: s~ c > w+ WEIGHTED=2 QED=1 [ all = QCD ]
+C     
+      IMPLICIT NONE
+C     
+C     CONSTANTS
+C     
+      INTEGER    NEXTERNAL
+      PARAMETER (NEXTERNAL=3)
+      INTEGER                 NCOMB
+      PARAMETER (             NCOMB=12)
+      INTEGER NSQAMPSO
+      PARAMETER (NSQAMPSO=1)
+      INTEGER HELAVGFACTOR
+      PARAMETER (HELAVGFACTOR=4)
+      LOGICAL CHOSEN_SO_CONFIGS(NSQAMPSO)
+      DATA CHOSEN_SO_CONFIGS/.TRUE./
+      COMMON/CHOSEN_BORN_SQSO/CHOSEN_SO_CONFIGS
+C     
+C     ARGUMENTS 
+C     
+      REAL*8 P(0:3,NEXTERNAL),ANS(0:NSQAMPSO)
+C     
+C     LOCAL VARIABLES 
+C     
+      INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
+      REAL*8 T(NSQAMPSO), BUFF
+      INTEGER IHEL,IDEN, I
+      INTEGER JC(NEXTERNAL)
+      LOGICAL GOODHEL(NCOMB)
+      DATA NTRY/0/
+      DATA GOODHEL/NCOMB*.FALSE./
+      DATA (NHEL(I,   1),I=1,3) /-1, 1,-1/
+      DATA (NHEL(I,   2),I=1,3) /-1, 1, 0/
+      DATA (NHEL(I,   3),I=1,3) /-1, 1, 1/
+      DATA (NHEL(I,   4),I=1,3) /-1,-1,-1/
+      DATA (NHEL(I,   5),I=1,3) /-1,-1, 0/
+      DATA (NHEL(I,   6),I=1,3) /-1,-1, 1/
+      DATA (NHEL(I,   7),I=1,3) / 1, 1,-1/
+      DATA (NHEL(I,   8),I=1,3) / 1, 1, 0/
+      DATA (NHEL(I,   9),I=1,3) / 1, 1, 1/
+      DATA (NHEL(I,  10),I=1,3) / 1,-1,-1/
+      DATA (NHEL(I,  11),I=1,3) / 1,-1, 0/
+      DATA (NHEL(I,  12),I=1,3) / 1,-1, 1/
+      DATA IDEN/36/
+C     
+C     GLOBAL VARIABLES
+C     
+      INTEGER USERHEL
+      DATA USERHEL/-1/
+      COMMON/HELUSERCHOICE/USERHEL
+
+C     ----------
+C     BEGIN CODE
+C     ----------
+      NTRY=NTRY+1
+      DO IHEL=1,NEXTERNAL
+        JC(IHEL) = +1
+      ENDDO
+      DO I=1,NSQAMPSO
+        ANS(I) = 0D0
+      ENDDO
+      DO IHEL=1,NCOMB
+        IF (USERHEL.EQ.-1.OR.USERHEL.EQ.IHEL) THEN
+          IF (GOODHEL(IHEL) .OR. NTRY .LT. 2 .OR.USERHEL.NE.-1) THEN
+            CALL MATRIX(P ,NHEL(1,IHEL),JC(1), T)
+            BUFF=0D0
+            DO I=1,NSQAMPSO
+              ANS(I)=ANS(I)+T(I)
+              BUFF=BUFF+T(I)
+            ENDDO
+            IF (BUFF .NE. 0D0 .AND. .NOT.    GOODHEL(IHEL)) THEN
+              GOODHEL(IHEL)=.TRUE.
+            ENDIF
+          ENDIF
+        ENDIF
+      ENDDO
+      ANS(0)=0.0D0
+      DO I=1,NSQAMPSO
+        ANS(I)=ANS(I)/DBLE(IDEN)
+        IF (CHOSEN_SO_CONFIGS(I)) THEN
+          ANS(0)=ANS(0)+ANS(I)
+        ENDIF
+      ENDDO
+      IF(USERHEL.NE.-1) THEN
+        ANS(0)=ANS(0)*HELAVGFACTOR
+        DO I=1,NSQAMPSO
+          ANS(I)=ANS(I)*HELAVGFACTOR
+        ENDDO
+      ENDIF
+      END
+
+      SUBROUTINE SMATRIXHEL_SPLITORDERS(P,HEL,ANS)
+      IMPLICIT NONE
+C     
+C     CONSTANT
+C     
+      INTEGER    NEXTERNAL
+      PARAMETER (NEXTERNAL=3)
+      INTEGER                 NCOMB
+      PARAMETER (             NCOMB=12)
+      INTEGER NSQAMPSO
+      PARAMETER (NSQAMPSO=1)
+C     
+C     ARGUMENTS 
+C     
+      REAL*8 P(0:3,NEXTERNAL),ANS(0:NSQAMPSO)
+      INTEGER HEL
+C     
+C     GLOBAL VARIABLES
+C     
+      INTEGER USERHEL
+      COMMON/HELUSERCHOICE/USERHEL
+C     ----------
+C     BEGIN CODE
+C     ----------
+      USERHEL=HEL
+      CALL SMATRIX_SPLITORDERS(P,ANS)
+      USERHEL=-1
+
+      END
+
+      SUBROUTINE MATRIX(P,NHEL,IC,RES)
+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     Returns amplitude squared summed/avg over colors
+C     for the point with external lines W(0:6,NEXTERNAL)
+C     
+C     Process: d~ u > w+ WEIGHTED=2 QED=1 [ all = QCD ]
+C     Process: s~ c > w+ WEIGHTED=2 QED=1 [ all = QCD ]
+C     
+      IMPLICIT NONE
+C     
+C     CONSTANTS
+C     
+      INTEGER    NGRAPHS
+      PARAMETER (NGRAPHS=1)
+      INTEGER    NEXTERNAL
+      PARAMETER (NEXTERNAL=3)
+      INTEGER    NWAVEFUNCS, NCOLOR
+      PARAMETER (NWAVEFUNCS=3, NCOLOR=1)
+      INTEGER NAMPSO, NSQAMPSO
+      PARAMETER (NAMPSO=1, NSQAMPSO=1)
+      REAL*8     ZERO
+      PARAMETER (ZERO=0D0)
+      COMPLEX*16 IMAG1
+      PARAMETER (IMAG1=(0D0,1D0))
+C     
+C     ARGUMENTS 
+C     
+      REAL*8 P(0:3,NEXTERNAL)
+      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
+      REAL*8 RES(NSQAMPSO)
+C     
+C     LOCAL VARIABLES 
+C     
+      INTEGER I,J,M,N
+      COMPLEX*16 ZTEMP
+      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
+      COMPLEX*16 AMP(NGRAPHS)
+      COMPLEX*16 JAMP(NCOLOR,NAMPSO)
+      COMPLEX*16 W(20,NWAVEFUNCS)
+      COMPLEX*16 DUM0,DUM1
+      DATA DUM0, DUM1/(0D0, 0D0), (1D0, 0D0)/
+C     
+C     FUNCTION
+C     
+      INTEGER SQSOINDEX
+C     
+C     GLOBAL VARIABLES
+C     
+      INCLUDE 'coupl.inc'
+C     
+C     COLOR DATA
+C     
+      DATA DENOM(1)/1/
+      DATA (CF(I,  1),I=  1,  1) /    3/
+C     1 T(1,2)
+C     ----------
+C     BEGIN CODE
+C     ----------
+      CALL OXXXXX(P(0,1),ZERO,NHEL(1),-1*IC(1),W(1,1))
+      CALL IXXXXX(P(0,2),ZERO,NHEL(2),+1*IC(2),W(1,2))
+      CALL VXXXXX(P(0,3),MDL_MW,NHEL(3),+1*IC(3),W(1,3))
+C     Amplitude(s) for diagram number 1
+      CALL FFV2_0(W(1,2),W(1,1),W(1,3),GC_47,AMP(1))
+C     JAMPs contributing to orders QCD=0
+      JAMP(1,1)=-AMP(1)
+
+      RES = 0.D0
+      DO M = 1, NAMPSO
+        DO I = 1, NCOLOR
+          ZTEMP = (0.D0,0.D0)
+          DO J = 1, NCOLOR
+            ZTEMP = ZTEMP + CF(J,I)*JAMP(J,M)
+          ENDDO
+          DO N = 1, NAMPSO
+            RES(SQSOINDEX(M,N)) = RES(SQSOINDEX(M,N)) + ZTEMP
+     $       *DCONJG(JAMP(I,N))/DENOM(I)
+          ENDDO
+        ENDDO
+      ENDDO
+
+      END
+
+      SUBROUTINE GET_ME(P, ALPHAS, NHEL ,ANS)
+      IMPLICIT NONE
+C     
+C     CONSTANT
+C     
+      INTEGER    NEXTERNAL
+      PARAMETER (NEXTERNAL=3)
+C     
+C     ARGUMENTS 
+C     
+      REAL*8 P(0:3,NEXTERNAL),ANS
+      INTEGER NHEL
+      DOUBLE PRECISION ALPHAS
+      REAL*8 PI
+CF2PY INTENT(OUT) :: ANS
+CF2PY INTENT(IN) :: NHEL
+CF2PY INTENT(IN) :: P(0:3,NEXTERNAL)
+CF2PY INTENT(IN) :: ALPHAS
+C     ROUTINE FOR F2PY to read the benchmark point.    
+C     the include file with the values of the parameters and masses 
+      INCLUDE 'coupl.inc'
+
+      PI = 3.141592653589793D0
+      G = 2* DSQRT(ALPHAS*PI)
+      CALL UPDATE_AS_PARAM()
+      IF (NHEL.NE.0) THEN
+        CALL SMATRIXHEL(P, NHEL, ANS)
+      ELSE
+        CALL SMATRIX(P, ANS)
+      ENDIF
+      RETURN
+      END
+
+      SUBROUTINE INITIALISE(PATH)
+C     ROUTINE FOR F2PY to read the benchmark point.    
+      IMPLICIT NONE
+      CHARACTER*180 PATH
+CF2PY INTENT(IN) :: PATH
+      CALL SETPARA(PATH)  !first call to setup the paramaters    
+      RETURN
+      END
+
+
+C     Set of functions to handle the array indices of the split orders
+
+
+      INTEGER FUNCTION SQSOINDEX(ORDERINDEXA, ORDERINDEXB)
+C     
+C     This functions plays the role of the interference matrix. It can
+C      be hardcoded or 
+C     made more elegant using hashtables if its execution speed ever
+C      becomes a relevant
+C     factor. From two split order indices, it return the corresponding
+C      index in the squared 
+C     order canonical ordering.
+C     
+C     CONSTANTS
+C     
+
+      INTEGER    NSO, NSQUAREDSO, NAMPSO
+      PARAMETER (NSO=1, NSQUAREDSO=1, NAMPSO=1)
+C     
+C     ARGUMENTS
+C     
+      INTEGER ORDERINDEXA, ORDERINDEXB
+C     
+C     LOCAL VARIABLES
+C     
+      INTEGER I, SQORDERS(NSO)
+      INTEGER AMPSPLITORDERS(NAMPSO,NSO)
+      DATA (AMPSPLITORDERS(  1,I),I=  1,  1) /    0/
+      COMMON/AMPSPLITORDERS/AMPSPLITORDERS
+C     
+C     FUNCTION
+C     
+      INTEGER SOINDEX_FOR_SQUARED_ORDERS
+C     
+C     BEGIN CODE
+C     
+      DO I=1,NSO
+        SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)+AMPSPLITORDERS(ORDERI
+     $   NDEXB,I)
+      ENDDO
+      SQSOINDEX=SOINDEX_FOR_SQUARED_ORDERS(SQORDERS)
+      END
+
+      INTEGER FUNCTION SOINDEX_FOR_SQUARED_ORDERS(ORDERS)
+C     
+C     This functions returns the integer index identifying the squared
+C      split orders list passed in argument which corresponds to the
+C      values of the following list of couplings (and in this order).
+C     ['QCD']
+C     
+C     CONSTANTS
+C     
+      INTEGER    NSO, NSQSO, NAMPSO
+      PARAMETER (NSO=1, NSQSO=1, NAMPSO=1)
+C     
+C     ARGUMENTS
+C     
+      INTEGER ORDERS(NSO)
+C     
+C     LOCAL VARIABLES
+C     
+      INTEGER I,J
+      INTEGER SQSPLITORDERS(NSQSO,NSO)
+      DATA (SQSPLITORDERS(  1,I),I=  1,  1) /    0/
+      COMMON/SQPLITORDERS/SQPLITORDERS
+C     
+C     BEGIN CODE
+C     
+      DO I=1,NSQSO
+        DO J=1,NSO
+          IF (ORDERS(J).NE.SQSPLITORDERS(I,J)) GOTO 1009
+        ENDDO
+        SOINDEX_FOR_SQUARED_ORDERS = I
+        RETURN
+ 1009   CONTINUE
+      ENDDO
+
+      WRITE(*,*) 'ERROR:: Stopping in function'
+      WRITE(*,*) 'SOINDEX_FOR_SQUARED_ORDERS'
+      WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1,NSO)
+      STOP
+
+      END
+
+      SUBROUTINE GET_NSQSO_BORN(NSQSO)
+C     
+C     Simple subroutine returning the number of squared split order
+C     contributions returned when calling smatrix_split_orders 
+C     
+
+      INTEGER    NSQUAREDSO
+      PARAMETER  (NSQUAREDSO=1)
+
+      INTEGER NSQSO
+
+      NSQSO=NSQUAREDSO
+
+      END
+
+C     This is the inverse subroutine of SOINDEX_FOR_SQUARED_ORDERS.
+C      Not directly useful, but provided nonetheless.
+      SUBROUTINE GET_SQUARED_ORDERS_FOR_SOINDEX(SOINDEX,ORDERS)
+C     
+C     This functions returns the orders identified by the squared
+C      split order index in argument. Order values correspond to
+C      following list of couplings (and in this order):
+C     ['QCD']
+C     
+C     CONSTANTS
+C     
+      INTEGER    NSO, NSQSO
+      PARAMETER (NSO=1, NSQSO=1)
+C     
+C     ARGUMENTS
+C     
+      INTEGER SOINDEX, ORDERS(NSO)
+C     
+C     LOCAL VARIABLES
+C     
+      INTEGER I
+      INTEGER SQPLITORDERS(NSQSO,NSO)
+      COMMON/SQPLITORDERS/SQPLITORDERS
+C     
+C     BEGIN CODE
+C     
+      IF (SOINDEX.GT.0.AND.SOINDEX.LE.NSQSO) THEN
+        DO I=1,NSO
+          ORDERS(I) =  SQPLITORDERS(SOINDEX,I)
+        ENDDO
+        RETURN
+      ENDIF
+
+      WRITE(*,*) 'ERROR:: Stopping function GET_SQUARED_ORDERS_FOR_SOI'
+     $ //'NDEX'
+      WRITE(*,*) 'Could not find squared orders index ',SOINDEX
+      STOP
+
+      END SUBROUTINE
+
+C     This is the inverse subroutine of getting amplitude SO orders.
+C      Not directly useful, but provided nonetheless.
+      SUBROUTINE GET_ORDERS_FOR_AMPSOINDEX(SOINDEX,ORDERS)
+C     
+C     This functions returns the orders identified by the split order
+C      index in argument. Order values correspond to following list of
+C      couplings (and in this order):
+C     ['QCD']
+C     
+C     CONSTANTS
+C     
+      INTEGER    NSO, NAMPSO
+      PARAMETER (NSO=1, NAMPSO=1)
+C     
+C     ARGUMENTS
+C     
+      INTEGER SOINDEX, ORDERS(NSO)
+C     
+C     LOCAL VARIABLES
+C     
+      INTEGER I
+      INTEGER AMPSPLITORDERS(NAMPSO,NSO)
+      COMMON/AMPSPLITORDERS/AMPSPLITORDERS
+C     
+C     BEGIN CODE
+C     
+      IF (SOINDEX.GT.0.AND.SOINDEX.LE.NAMPSO) THEN
+        DO I=1,NSO
+          ORDERS(I) =  AMPSPLITORDERS(SOINDEX,I)
+        ENDDO
+        RETURN
+      ENDIF
+
+      WRITE(*,*) 'ERROR:: Stopping function GET_ORDERS_FOR_AMPSOINDEX'
+      WRITE(*,*) 'Could not find amplitude split orders index ',SOINDEX
+      STOP
+
+      END SUBROUTINE
+
+C     This function is not directly useful, but included for completene
+C     ss
+      INTEGER FUNCTION SOINDEX_FOR_AMPORDERS(ORDERS)
+C     
+C     This functions returns the integer index identifying the
+C      amplitude split orders passed in argument which correspond to
+C      the values of the following list of couplings (and in this
+C      order):
+C     ['QCD']
+C     
+C     CONSTANTS
+C     
+      INTEGER    NSO, NAMPSO
+      PARAMETER (NSO=1, NAMPSO=1)
+C     
+C     ARGUMENTS
+C     
+      INTEGER ORDERS(NSO)
+C     
+C     LOCAL VARIABLES
+C     
+      INTEGER I,J
+      INTEGER AMPSPLITORDERS(NAMPSO,NSO)
+      COMMON/AMPSPLITORDERS/AMPSPLITORDERS
+C     
+C     BEGIN CODE
+C     
+      DO I=1,NAMPSO
+        DO J=1,NSO
+          IF (ORDERS(J).NE.AMPSPLITORDERS(I,J)) GOTO 1009
+        ENDDO
+        SOINDEX_FOR_AMPORDERS = I
+        RETURN
+ 1009   CONTINUE
+      ENDDO
+
+      WRITE(*,*) 'ERROR:: Stopping function SOINDEX_FOR_AMPORDERS'
+      WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1,NSO)
+      STOP
+
+      END
+

=== added 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	1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa.f	2016-02-24 13:54:50 +0000
@@ -0,0 +1,756 @@
+      PROGRAM DRIVER
+C     *****************************************************************
+C     ********
+C     THIS IS THE DRIVER FOR CHECKING THE STANDALONE MATRIX ELEMENT.
+C     IT USES A SIMPLE PHASE SPACE GENERATOR
+C     *****************************************************************
+C     ********
+      IMPLICIT NONE
+C     
+C     CONSTANTS  
+C     
+      REAL*8 ZERO
+      PARAMETER (ZERO=0D0)
+
+      LOGICAL READPS
+      PARAMETER (READPS = .FALSE.)
+
+      INTEGER NPSPOINTS
+      PARAMETER (NPSPOINTS = 4)
+
+C     integer nexternal and number particles (incoming+outgoing) in
+C      the me 
+      INTEGER NEXTERNAL, NINCOMING
+      PARAMETER (NEXTERNAL=3,NINCOMING=2)
+
+      CHARACTER(512) MADLOOPRESOURCEPATH
+
+C     
+C     INCLUDE FILES
+C     
+C     the include file with the values of the parameters and masses   
+C        
+      INCLUDE 'coupl.inc'
+C     particle masses
+      REAL*8 PMASS(NEXTERNAL)
+C     integer    n_max_cg
+      INCLUDE 'ngraphs.inc'
+      INCLUDE 'nsqso_born.inc'
+      INCLUDE 'nsquaredSO.inc'
+
+C     
+C     LOCAL
+C     
+      INTEGER I,J,K
+C     four momenta. Energy is the zeroth component.
+      REAL*8 P(0:3,NEXTERNAL)
+      INTEGER MATELEM_ARRAY_DIM
+      REAL*8 , ALLOCATABLE :: MATELEM(:,:)
+      REAL*8 SQRTS,AO2PI,TOTMASS
+C     sqrt(s)= center of mass energy 
+      REAL*8 PIN(0:3), POUT(0:3)
+      CHARACTER*120 BUFF(NEXTERNAL)
+      INTEGER RETURNCODE, UNITS, TENS, HUNDREDS
+      INTEGER NSQUAREDSO_LOOP
+      REAL*8 , ALLOCATABLE :: PREC_FOUND(:)
+
+C     
+C     GLOBAL VARIABLES
+C     
+C     This is from ML code for the list of split orders selected by
+C     the process definition
+C     
+      INTEGER NLOOPCHOSEN
+      CHARACTER*20 CHOSEN_LOOP_SO_INDICES(NSQUAREDSO)
+      LOGICAL CHOSEN_LOOP_SO_CONFIGS(NSQUAREDSO)
+      COMMON/CHOSEN_LOOP_SQSO/CHOSEN_LOOP_SO_CONFIGS
+      INTEGER NBORNCHOSEN
+      CHARACTER*20 CHOSEN_BORN_SO_INDICES(NSQSO_BORN)
+      LOGICAL CHOSEN_BORN_SO_CONFIGS(NSQSO_BORN)
+      COMMON/CHOSEN_BORN_SQSO/CHOSEN_BORN_SO_CONFIGS
+
+C     
+C     SAVED VARIABLES
+C     
+      LOGICAL INIT
+      DATA INIT/.TRUE./
+      COMMON/INITCHECKSA/INIT
+C     
+C     EXTERNAL
+C     
+      REAL*8 DOT
+      EXTERNAL DOT
+
+C     
+C     BEGIN CODE
+C     
+C     
+
+      IF (INIT) THEN
+        INIT=.FALSE.
+        CALL GET_ANSWER_DIMENSION(MATELEM_ARRAY_DIM)
+        ALLOCATE(MATELEM(0:3,0:MATELEM_ARRAY_DIM))
+        CALL GET_NSQSO_LOOP(NSQUAREDSO_LOOP)
+        ALLOCATE(PREC_FOUND(0:NSQUAREDSO_LOOP))
+
+C       INITIALIZATION CALLS
+C       
+C       Call to initialize the values of the couplings, masses and
+C        widths 
+C       used in the evaluation of the matrix element. The primary
+C        parameters of the
+C       models are read from Cards/param_card.dat. The secondary
+C        parameters are calculated
+C       in Source/MODEL/couplings.f. The values are stored in common
+C        blocks that are listed
+C       in coupl.inc .
+C       first call to setup the paramaters
+        CALL SETPARA('param_card.dat')
+C       set up masses
+        INCLUDE 'pmass.inc'
+
+      ENDIF
+
+
+C     Start by initializing what is the squared split orders indices
+C      chosen
+      NLOOPCHOSEN=0
+      DO I=1,NSQUAREDSO
+        IF (CHOSEN_LOOP_SO_CONFIGS(I)) THEN
+          NLOOPCHOSEN=NLOOPCHOSEN+1
+          WRITE(CHOSEN_LOOP_SO_INDICES(NLOOPCHOSEN),'(I3,A2)') I,'L)'
+        ENDIF
+      ENDDO
+      NBORNCHOSEN=0
+      DO I=1,NSQSO_BORN
+        IF (CHOSEN_BORN_SO_CONFIGS(I)) THEN
+          NBORNCHOSEN=NBORNCHOSEN+1
+          WRITE(CHOSEN_BORN_SO_INDICES(NBORNCHOSEN),'(I3,A2)') I,'B)'
+        ENDIF
+      ENDDO
+
+      AO2PI=G**2/(8.D0*(3.14159265358979323846D0**2))
+
+      WRITE(*,*) 'AO2PI=',AO2PI
+C     Now use a simple multipurpose PS generator (RAMBO) just to get a 
+C     RANDOM set of four momenta of given masses pmass(i) to be used
+C      to evaluate 
+C     the madgraph matrix-element.       
+C     Alternatevely, here the user can call or set the four momenta at
+C      his will, see below.
+C     
+      IF(NINCOMING.EQ.1) THEN
+        SQRTS=PMASS(1)
+      ELSE
+        TOTMASS = 0.0D0
+        DO I=1,NEXTERNAL
+          TOTMASS = TOTMASS + PMASS(I)
+        ENDDO
+C       CMS energy in GEV
+        SQRTS=MAX(1000D0,2.0D0*TOTMASS)
+      ENDIF
+
+      CALL PRINTOUT()
+
+
+
+      DO K=1,NPSPOINTS
+
+        IF(READPS) THEN
+          OPEN(967, FILE='PS.input', ERR=976, STATUS='OLD', ACTION='RE'
+     $     //'AD')
+          DO I=1,NEXTERNAL
+            READ(967,*,END=978) P(0,I),P(1,I),P(2,I),P(3,I)
+          ENDDO
+          GOTO 978
+ 976      CONTINUE
+          STOP 'Could not read the PS.input phase-space point.'
+ 978      CONTINUE
+          CLOSE(967)
+        ELSE
+          IF ((NINCOMING.EQ.2).AND.((NEXTERNAL - NINCOMING .EQ.1))
+     $     ) THEN
+            IF (PMASS(3).EQ.0.0D0) THEN
+              STOP 'Cannot generate 2>1 kin. config. with m3=0.0d0'
+            ELSE
+C             deal with the case of only one particle in the final
+C              state
+              P(0,1) = PMASS(3)/2D0
+              P(1,1) = 0D0
+              P(2,1) = 0D0
+              P(3,1) = PMASS(3)/2D0
+              IF (PMASS(1).GT.0D0) THEN
+                P(3,1) = DSQRT(PMASS(3)**2/4D0 - PMASS(1)**2)
+              ENDIF
+              P(0,2) = PMASS(3)/2D0
+              P(1,2) = 0D0
+              P(2,2) = 0D0
+              P(3,2) = -PMASS(3)/2D0
+              IF (PMASS(2) > 0D0) THEN
+                P(3,2) = -DSQRT(PMASS(3)**2/4D0 - PMASS(1)**2)
+              ENDIF
+              P(0,3) = PMASS(3)
+              P(1,3) = 0D0
+              P(2,3) = 0D0
+              P(3,3) = 0D0
+            ENDIF
+          ELSE
+            CALL GET_MOMENTA(SQRTS,PMASS,P)
+          ENDIF
+        ENDIF
+
+        DO I=0,3
+          PIN(I)=0.0D0
+          DO J=1,NINCOMING
+            PIN(I)=PIN(I)+P(I,J)
+          ENDDO
+        ENDDO
+
+C       In standalone mode, always use sqrt_s as the renormalization
+C        scale.
+        SQRTS=DSQRT(DABS(DOT(PIN(0),PIN(0))))
+        MU_R=SQRTS
+
+C       Update the couplings with the new MU_R
+        CALL UPDATE_AS_PARAM()
+
+C       Optionally the user can set where to find the MadLoop5_resource
+C       s folder.
+C       Otherwise it will look for it automatically and find it if it
+C        has not
+C       been moved
+C       MadLoopResourcePath = '<MadLoop5_resources_path>'
+C       CALL SETMADLOOPPATH(MadLoopResourcePath)
+C       To force the stabiliy check to also be performed in the
+C        initialization phase
+C       CALL FORCE_STABILITY_CHECK(.TRUE.)
+C       To chose a particular tartget split order, SOTARGET is an
+C        integer labeling
+C       the possible squared order couplings contributions (only in
+C        optimized mode)
+C       CALL SET_COUPLINGORDERS_TARGET(SOTARGET)
+
+
+C       
+C       Now we can call the matrix element
+C       
+        CALL SLOOPMATRIX_THRES(P,MATELEM,-1.0D0,PREC_FOUND,RETURNCODE)
+
+C       
+C       write the information on the four momenta 
+C       
+        IF (K.EQ.NPSPOINTS) THEN
+          WRITE (*,*)
+          WRITE (*,*) ' Phase space point:'
+          WRITE (*,*)
+          WRITE (*,*) '---------------------------------'
+          WRITE (*,*)  'n  E  px py pz m'
+          DO I=1,NEXTERNAL
+            WRITE (*,'(i2,1x,5e15.7)') I, P(0,I),P(1,I),P(2,I),P(3,I)
+     $       ,DSQRT(DABS(DOT(P(0,I),P(0,I))))
+          ENDDO
+          WRITE (*,*) '---------------------------------'
+          WRITE (*,*) 'Detailed result for each coupling order'
+     $     //'s combination.'
+          WRITE(*,*) 'All Born contributions are of split order'
+     $     //'s (QCD=0)'
+          WRITE (*,*) '---------------------------------'
+          WRITE(*,*) 'All loop contributions are of split order'
+     $     //'s (QCD=2)'
+          WRITE (*,*) '---------------------------------'
+          UNITS=MOD(RETURNCODE,10)
+          TENS=(MOD(RETURNCODE,100)-UNITS)/10
+          HUNDREDS=(RETURNCODE-TENS*10-UNITS)/100
+          IF (HUNDREDS.EQ.1) THEN
+            IF (TENS.EQ.3.OR.TENS.EQ.4) THEN
+              WRITE(*,*) 'Unknown numerical stability because MadLoo'
+     $         //'p is in the initialization stage.'
+            ELSE
+              WRITE(*,*) 'Unknown numerical stability, check CTModeRu'
+     $         //'n value in MadLoopParams.dat.'
+            ENDIF
+          ELSEIF (HUNDREDS.EQ.2) THEN
+            WRITE(*,*) 'Stable kinematic configuration (SPS).'
+          ELSEIF (HUNDREDS.EQ.3) THEN
+            WRITE(*,*) 'Unstable kinematic configuration (UPS).'
+            WRITE(*,*) 'Quadruple precision rescue successful.'
+          ELSEIF (HUNDREDS.EQ.4) THEN
+            WRITE(*,*) 'Exceptional kinematic configuration (EPS).'
+            WRITE(*,*) 'Both double an quadruple precision computation'
+     $       //'s, are unstable.'
+          ENDIF
+          IF (TENS.EQ.2.OR.TENS.EQ.4) THEN
+            WRITE(*,*) 'Quadruple precision computation used.'
+          ENDIF
+          IF (HUNDREDS.NE.1) THEN
+            IF (PREC_FOUND(0).GT.0.0D0) THEN
+              WRITE(*,'(1x,a23,1x,1e10.2)') 'Relative accuracy     ='
+     $         ,PREC_FOUND(0)
+            ELSEIF (PREC_FOUND(0).EQ.0.0D0) THEN
+              WRITE(*,'(1x,a23,1x,1e10.2,1x,a30)') 'Relative accuracy'
+     $         //'     =',PREC_FOUND(0),'(i.e. beyond double precisio'
+     $         //'n)'
+            ELSE
+              WRITE(*,*) 'Estimated accuracy could not be computed fo'
+     $         //'r an unknown reason.'
+            ENDIF
+          ENDIF
+          WRITE (*,*) '---------------------------------'
+          IF (NBORNCHOSEN.EQ.0) THEN
+            WRITE (*,*) 'No Born contribution satisfied the square'
+     $       //'d order constraints.'
+          ELSE IF (NBORNCHOSEN.NE.NSQSO_BORN) THEN
+            WRITE (*,*) 'Selected squared coupling orders combinatio'
+     $       //'n for the Born summed result below:'
+            WRITE (*,*) (CHOSEN_BORN_SO_INDICES(I),I=1,NBORNCHOSEN)
+          ENDIF
+          IF (NLOOPCHOSEN.NE.NSQUAREDSO) THEN
+            WRITE (*,*) 'Selected squared coupling orders combinatio'
+     $       //'n for the loop summed result below:'
+            WRITE (*,*) (CHOSEN_LOOP_SO_INDICES(I),I=1,NLOOPCHOSEN)
+          ENDIF
+          WRITE (*,*) '---------------------------------'
+          WRITE (*,*) 'Matrix element born   = ', MATELEM(0,0)
+     $     , ' GeV^',-(2*NEXTERNAL-8)
+          WRITE (*,*) 'Matrix element finite = ', MATELEM(1,0)
+     $     , ' GeV^',-(2*NEXTERNAL-8)
+          WRITE (*,*) 'Matrix element 1eps   = ', MATELEM(2,0)
+     $     , ' GeV^',-(2*NEXTERNAL-8)
+          WRITE (*,*) 'Matrix element 2eps   = ', MATELEM(3,0)
+     $     , ' GeV^',-(2*NEXTERNAL-8)
+          WRITE (*,*) '---------------------------------'
+          IF (MATELEM(0,0).NE.0.0D0) THEN
+            WRITE (*,*) 'finite / (born*ao2pi) = ', MATELEM(1,0)
+     $       /MATELEM(0,0)/AO2PI
+            WRITE (*,*) '1eps   / (born*ao2pi) = ', MATELEM(2,0)
+     $       /MATELEM(0,0)/AO2PI
+            WRITE (*,*) '2eps   / (born*ao2pi) = ', MATELEM(3,0)
+     $       /MATELEM(0,0)/AO2PI
+          ELSE
+            WRITE (*,*) 'finite / ao2pi      = ', MATELEM(1,0)/AO2PI
+            WRITE (*,*) '1eps   / ao2pi      = ', MATELEM(2,0)/AO2PI
+            WRITE (*,*) '2eps   / ao2pi      = ', MATELEM(3,0)/AO2PI
+          ENDIF
+          WRITE (*,*) '---------------------------------'
+
+          OPEN(69, FILE='result.dat', ERR=976, ACTION='WRITE')
+          DO I=1,NEXTERNAL
+            WRITE (69,'(a2,1x,5e25.15)') 'PS',P(0,I),P(1,I),P(2,I),P(3
+     $       ,I)
+          ENDDO
+          WRITE (69,'(a3,1x,i2)') 'EXP',-(2*NEXTERNAL-8)
+          WRITE (69,'(a4,1x,1e25.15)') 'BORN',MATELEM(0,0)
+          IF (MATELEM(0,0).NE.0.0D0) THEN
+            WRITE (69,'(a3,1x,1e25.15)') 'FIN',MATELEM(1,0)/MATELEM(0
+     $       ,0)/AO2PI
+            WRITE (69,'(a4,1x,1e25.15)') '1EPS',MATELEM(2,0)/MATELEM(0
+     $       ,0)/AO2PI
+            WRITE (69,'(a4,1x,1e25.15)') '2EPS',MATELEM(3,0)/MATELEM(0
+     $       ,0)/AO2PI
+          ELSE
+            WRITE (69,'(a3,1x,1e25.15)') 'FIN',MATELEM(1,0)/AO2PI
+            WRITE (69,'(a4,1x,1e25.15)') '1EPS',MATELEM(2,0)/AO2PI
+            WRITE (69,'(a4,1x,1e25.15)') '2EPS',MATELEM(3,0)/AO2PI
+          ENDIF
+          WRITE (69,'(a6,1x,1e25.15)') 'ASO2PI',AO2PI
+          WRITE (69,*) 'Export_Format Default'
+          WRITE (69,'(a7,1x,i3)') 'RETCODE',RETURNCODE
+          WRITE (69,'(a3,1x,1e10.4)') 'ACC',PREC_FOUND(0)
+          WRITE (69,*) 'Born_kept',(CHOSEN_BORN_SO_CONFIGS(I),I=1
+     $     ,NSQSO_BORN)
+          WRITE (69,*) 'Loop_kept',(CHOSEN_LOOP_SO_CONFIGS(I),I=1
+     $     ,NSQUAREDSO)
+          WRITE (69,*) 'Born_SO_Results 0'
+          WRITE (69,*) 'SO_Born BORN ',MATELEM(0,1)
+          WRITE (69,*) 'Split_Orders_Names QCD'
+          WRITE (69,*) 'Loop_SO_Results 2'
+          WRITE (69,*) 'SO_Loop ACC  ',PREC_FOUND(1)
+          WRITE (69,*) 'SO_Loop FIN  ',MATELEM(1,1)
+          WRITE (69,*) 'SO_Loop 1EPS ',MATELEM(2,1)
+          WRITE (69,*) 'SO_Loop 2EPS ',MATELEM(3,1)
+          CLOSE(69)
+        ELSE
+          WRITE (*,*) 'PS Point #',K,' done.'
+        ENDIF
+      ENDDO
+
+C     C
+C     C      Copy down here (or read in) the four momenta as a string. 
+C     C      
+C     C
+C     buff(1)=" 1   0.5630480E+04  0.0000000E+00  0.0000000E+00 
+C      0.5630480E+04"
+C     buff(2)=" 2   0.5630480E+04  0.0000000E+00  0.0000000E+00
+C      -0.5630480E+04"
+C     buff(3)=" 3   0.5466073E+04  0.4443190E+03  0.2446331E+04
+C      -0.4864732E+04"
+C     buff(4)=" 4   0.8785819E+03 -0.2533886E+03  0.2741971E+03 
+C      0.7759741E+03"
+C     buff(5)=" 5   0.4916306E+04 -0.1909305E+03 -0.2720528E+04 
+C      0.4088757E+04"
+C     C
+C     C      Here the k,E,px,py,pz are read from the string into the
+C      momenta array.
+C     C      k=1,2          : incoming
+C     C      k=3,nexternal  : outgoing
+C     C
+C     do i=1,nexternal
+C     read (buff(i),*) k, P(0,i),P(1,i),P(2,i),P(3,i)
+C     enddo
+C     
+C     C print the momenta out
+C     
+C     do i=1,nexternal
+C     write (*,'(i2,1x,5e15.7)') i, P(0,i),P(1,i),P(2,i),P(3,i), 
+C     &dsqrt(dabs(DOT(p(0,i),p(0,i))))
+C     enddo
+C     
+C     CALL SLOOPMATRIX(P,MATELEM)
+C     
+C     write (*,*) "-------------------------------------------------"
+C     write (*,*) "Matrix element = ", MATELEM(1), " GeV^",-(2*nexterna
+C     l-8)      
+C     write (*,*) "-------------------------------------------------"
+
+      DEALLOCATE(MATELEM)
+      DEALLOCATE(PREC_FOUND)
+
+      END
+
+
+
+
+      DOUBLE PRECISION FUNCTION DOT(P1,P2)
+C     *************************************************************
+C     4-Vector Dot product
+C     *************************************************************
+      IMPLICIT NONE
+      DOUBLE PRECISION P1(0:3),P2(0:3)
+      DOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)
+      END
+
+
+      SUBROUTINE GET_MOMENTA(ENERGY,PMASS,P)
+C     auxiliary function to change convention between madgraph and
+C      rambo
+C     four momenta.         
+      IMPLICIT NONE
+      INTEGER NEXTERNAL, NINCOMING
+      PARAMETER (NEXTERNAL=3,NINCOMING=2)
+C     ARGUMENTS
+      REAL*8 ENERGY,PMASS(NEXTERNAL),P(0:3,NEXTERNAL),PRAMBO(4,10),WGT
+C     LOCAL
+      INTEGER I
+      REAL*8 ETOT2,MOM,M1,M2,E1,E2
+
+      ETOT2=ENERGY**2
+      M1=PMASS(1)
+      M2=PMASS(2)
+      MOM=(ETOT2**2 - 2*ETOT2*M1**2 + M1**4 - 2*ETOT2*M2**2 - 2*M1**2
+     $ *M2**2 + M2**4)/(4.*ETOT2)
+      MOM=DSQRT(MOM)
+      E1=DSQRT(MOM**2+M1**2)
+      E2=DSQRT(MOM**2+M2**2)
+C     write (*,*) e1+e2,mom
+
+      IF(NINCOMING.EQ.2) THEN
+
+        P(0,1)=E1
+        P(1,1)=0D0
+        P(2,1)=0D0
+        P(3,1)=MOM
+
+        P(0,2)=E2
+        P(1,2)=0D0
+        P(2,2)=0D0
+        P(3,2)=-MOM
+
+        CALL RAMBO(NEXTERNAL-2,ENERGY,PMASS(3),PRAMBO,WGT)
+        DO I=3, NEXTERNAL
+          P(0,I)=PRAMBO(4,I-2)
+          P(1,I)=PRAMBO(1,I-2)
+          P(2,I)=PRAMBO(2,I-2)
+          P(3,I)=PRAMBO(3,I-2)
+        ENDDO
+
+      ELSEIF(NINCOMING.EQ.1) THEN
+
+        P(0,1)=ENERGY
+        P(1,1)=0D0
+        P(2,1)=0D0
+        P(3,1)=0D0
+
+        CALL RAMBO(NEXTERNAL-1,ENERGY,PMASS(2),PRAMBO,WGT)
+        DO I=2, NEXTERNAL
+          P(0,I)=PRAMBO(4,I-1)
+          P(1,I)=PRAMBO(1,I-1)
+          P(2,I)=PRAMBO(2,I-1)
+          P(3,I)=PRAMBO(3,I-1)
+        ENDDO
+      ENDIF
+
+      RETURN
+      END
+
+
+      SUBROUTINE RAMBO(N,ET,XM,P,WT)
+C     *****************************************************************
+C     *****
+C     RAMBO                                         *
+C     RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)                 
+C      *
+C     *
+C     A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR               
+C      *
+C     AUTHORS:  S.D. ELLIS,  R. KLEISS,  W.J. STIRLING                
+C      *
+C     THIS IS VERSION 1.0 -  WRITTEN BY R. KLEISS                     
+C      *
+C     -- ADJUSTED BY HANS KUIJF, WEIGHTS ARE LOGARITHMIC (20-08-90)   
+C      *
+C     *
+C     N  = NUMBER OF PARTICLES                                        
+C      *
+C     ET = TOTAL CENTRE-OF-MASS ENERGY                                
+C      *
+C     XM = PARTICLE MASSES ( DIM=NEXTERNAL-nincoming )                
+C      *
+C     P  = PARTICLE MOMENTA ( DIM=(4,NEXTERNAL-nincoming) )           
+C      *
+C     WT = WEIGHT OF THE EVENT                                        
+C      *
+C     *****************************************************************
+C     *****
+      IMPLICIT REAL*8(A-H,O-Z)
+      INTEGER NEXTERNAL, NINCOMING
+      PARAMETER (NEXTERNAL=3,NINCOMING=2)
+      DIMENSION XM(NEXTERNAL-NINCOMING),P(4,NEXTERNAL-NINCOMING)
+      DIMENSION Q(4,NEXTERNAL-NINCOMING),Z(NEXTERNAL-NINCOMING),R(4)
+     $ ,B(3),P2(NEXTERNAL-NINCOMING),XM2(NEXTERNAL-NINCOMING)
+     $ ,E(NEXTERNAL-NINCOMING),V(NEXTERNAL-NINCOMING),IWARN(5)
+      SAVE ACC,ITMAX,IBEGIN,IWARN
+      DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/
+C     
+C     INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
+      IF(IBEGIN.NE.0) GOTO 103
+      IBEGIN=1
+      TWOPI=8.*DATAN(1.D0)
+      PO2LOG=LOG(TWOPI/4.)
+      Z(2)=PO2LOG
+      DO 101 K=3,(NEXTERNAL-NINCOMING)
+ 101  Z(K)=Z(K-1)+PO2LOG-2.*LOG(DFLOAT(K-2))
+      DO 102 K=3,(NEXTERNAL-NINCOMING)
+ 102  Z(K)=(Z(K)-LOG(DFLOAT(K-1)))
+C     
+C     CHECK ON THE NUMBER OF PARTICLES
+ 103  IF(N.GT.1.AND.N.LT.101) GOTO 104
+      PRINT 1001,N
+      STOP
+C     
+C     CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
+ 104  XMT=0.
+      NM=0
+      DO 105 I=1,N
+      IF(XM(I).NE.0.D0) NM=NM+1
+ 105  XMT=XMT+ABS(XM(I))
+      IF(XMT.LE.ET) GOTO 201
+      PRINT 1002,XMT,ET
+      STOP
+C     
+C     THE PARAMETER VALUES ARE NOW ACCEPTED
+C     
+C     GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
+ 201  DO 202 I=1,N
+      R1=RN(1)
+      C=2.*R1-1.
+      S=SQRT(1.-C*C)
+      F=TWOPI*RN(2)
+      R1=RN(3)
+      R2=RN(4)
+      Q(4,I)=-LOG(R1*R2)
+      Q(3,I)=Q(4,I)*C
+      Q(2,I)=Q(4,I)*S*COS(F)
+ 202  Q(1,I)=Q(4,I)*S*SIN(F)
+C     
+C     CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
+      DO 203 I=1,4
+ 203  R(I)=0.
+      DO 204 I=1,N
+      DO 204 K=1,4
+ 204  R(K)=R(K)+Q(K,I)
+      RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
+      DO 205 K=1,3
+ 205  B(K)=-R(K)/RMAS
+      G=R(4)/RMAS
+      A=1./(1.+G)
+      X=ET/RMAS
+C     
+C     TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
+      DO 207 I=1,N
+      BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
+      DO 206 K=1,3
+ 206  P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
+ 207  P(4,I)=X*(G*Q(4,I)+BQ)
+C     
+C     CALCULATE WEIGHT AND POSSIBLE WARNINGS
+      WT=PO2LOG
+      IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N)
+      IF(WT.GE.-180.D0) GOTO 208
+      IF(IWARN(1).LE.5) PRINT 1004,WT
+      IWARN(1)=IWARN(1)+1
+ 208  IF(WT.LE. 174.D0) GOTO 209
+      IF(IWARN(2).LE.5) PRINT 1005,WT
+      IWARN(2)=IWARN(2)+1
+C     
+C     RETURN FOR WEIGHTED MASSLESS MOMENTA
+ 209  IF(NM.NE.0) GOTO 210
+C     RETURN LOG OF WEIGHT
+      WT=WT
+      RETURN
+C     
+C     MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
+ 210  XMAX=SQRT(1.-(XMT/ET)**2)
+      DO 301 I=1,N
+      XM2(I)=XM(I)**2
+ 301  P2(I)=P(4,I)**2
+      ITER=0
+      X=XMAX
+      ACCU=ET*ACC
+ 302  F0=-ET
+      G0=0.
+      X2=X*X
+      DO 303 I=1,N
+      E(I)=SQRT(XM2(I)+X2*P2(I))
+      F0=F0+E(I)
+ 303  G0=G0+P2(I)/E(I)
+      IF(ABS(F0).LE.ACCU) GOTO 305
+      ITER=ITER+1
+      IF(ITER.LE.ITMAX) GOTO 304
+      PRINT 1006,ITMAX
+      GOTO 305
+ 304  X=X-F0/(X*G0)
+      GOTO 302
+ 305  DO 307 I=1,N
+      V(I)=X*P(4,I)
+      DO 306 K=1,3
+ 306  P(K,I)=X*P(K,I)
+ 307  P(4,I)=E(I)
+C     
+C     CALCULATE THE MASS-EFFECT WEIGHT FACTOR
+      WT2=1.
+      WT3=0.
+      DO 308 I=1,N
+      WT2=WT2*V(I)/E(I)
+ 308  WT3=WT3+V(I)**2/E(I)
+      WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET)
+C     
+C     RETURN FOR  WEIGHTED MASSIVE MOMENTA
+      WT=WT+WTM
+      IF(WT.GE.-180.D0) GOTO 309
+      IF(IWARN(3).LE.5) PRINT 1004,WT
+      IWARN(3)=IWARN(3)+1
+ 309  IF(WT.LE. 174.D0) GOTO 310
+      IF(IWARN(4).LE.5) PRINT 1005,WT
+      IWARN(4)=IWARN(4)+1
+C     RETURN LOG OF WEIGHT
+ 310  WT=WT
+      RETURN
+C     
+ 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
+ 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',' SMALLE'
+     $ //'R THAN TOTAL ENERGY =',D15.6)
+ 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
+ 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY  OVERFLOW')
+ 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE'
+     $ ,' DESIRED ACCURACY =',D15.6)
+      END
+
+      FUNCTION RN(IDUMMY)
+      REAL*8 RN,RAN
+      SAVE INIT
+      DATA INIT /1/
+      IF (INIT.EQ.1) THEN
+        INIT=0
+        CALL RMARIN(1802,9373)
+        END IF
+C       
+ 10     CALL RANMAR(RAN)
+        IF (RAN.LT.1D-16) GOTO 10
+        RN=RAN
+C       
+        END
+
+
+
+        SUBROUTINE RANMAR(RVEC)
+C       -----------------
+C       Universal random number generator proposed by Marsaglia and
+C        Zaman
+C       in report FSU-SCRI-87-50
+C       In this version RVEC is a double precision variable.
+        IMPLICIT REAL*8(A-H,O-Z)
+        COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
+        COMMON/ RASET2 / IRANMR,JRANMR
+        SAVE /RASET1/,/RASET2/
+        UNI = RANU(IRANMR) - RANU(JRANMR)
+        IF(UNI .LT. 0D0) UNI = UNI + 1D0
+        RANU(IRANMR) = UNI
+        IRANMR = IRANMR - 1
+        JRANMR = JRANMR - 1
+        IF(IRANMR .EQ. 0) IRANMR = 97
+        IF(JRANMR .EQ. 0) JRANMR = 97
+        RANC = RANC - RANCD
+        IF(RANC .LT. 0D0) RANC = RANC + RANCM
+        UNI = UNI - RANC
+        IF(UNI .LT. 0D0) UNI = UNI + 1D0
+        RVEC = UNI
+        END
+
+        SUBROUTINE RMARIN(IJ,KL)
+C       -----------------
+C       Initializing routine for RANMAR, must be called before
+C        generating
+C       any pseudorandom numbers with RANMAR. The input values should
+C        be in
+C       the ranges 0<=ij<=31328 ; 0<=kl<=30081
+        IMPLICIT REAL*8(A-H,O-Z)
+        COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
+        COMMON/ RASET2 / IRANMR,JRANMR
+        SAVE /RASET1/,/RASET2/
+C       This shows correspondence between the simplified input seeds
+C        IJ, KL
+C       and the original Marsaglia-Zaman seeds I,J,K,L.
+C       To get the standard values in the Marsaglia-Zaman paper
+C        (i=12,j=34
+C       k=56,l=78) put ij=1802, kl=9373
+        I = MOD( IJ/177 , 177 ) + 2
+        J = MOD( IJ     , 177 ) + 2
+        K = MOD( KL/169 , 178 ) + 1
+        L = MOD( KL     , 169 )
+        DO 300 II = 1 , 97
+        S =  0D0
+        T = .5D0
+        DO 200 JJ = 1 , 24
+        M = MOD( MOD(I*J,179)*K , 179 )
+        I = J
+        J = K
+        K = M
+        L = MOD( 53*L+1 , 169 )
+        IF(MOD(L*M,64) .GE. 32) S = S + T
+        T = .5D0*T
+ 200    CONTINUE
+        RANU(II) = S
+ 300    CONTINUE
+        RANC  =   362436D0 / 16777216D0
+        RANCD =  7654321D0 / 16777216D0
+        RANCM = 16777213D0 / 16777216D0
+        IRANMR = 97
+        JRANMR = 33
+        END
+
+
+
+
+
+
+

=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa_born_splitOrders.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa_born_splitOrders.f	1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa_born_splitOrders.f	2016-02-24 13:54:50 +0000
@@ -0,0 +1,524 @@
+      PROGRAM DRIVER
+C     *****************************************************************
+C     *********
+C     THIS IS THE DRIVER FOR CHECKING THE STANDALONE MATRIX ELEMENT.
+C     IT USES A SIMPLE PHASE SPACE GENERATOR
+C     Fabio Maltoni - 3rd Febraury 2007
+C     *****************************************************************
+C     *********
+      IMPLICIT NONE
+C     
+C     CONSTANTS  
+C     
+      REAL*8 ZERO
+      PARAMETER (ZERO=0D0)
+      INTEGER NSPLITORDERS
+      PARAMETER (NSPLITORDERS=1)
+C     
+C     INCLUDE FILES
+C     
+C     ---  the include file with the values of the parameters and
+C      masses	
+      INCLUDE 'coupl.inc'
+C     integer nexternal and number particles (incoming+outgoing) in
+C      the me 
+      INTEGER NEXTERNAL, NINCOMING
+      PARAMETER (NEXTERNAL=3,NINCOMING=2)
+C     ---  particle masses
+      REAL*8 PMASS(NEXTERNAL)
+      REAL*8 TOTALMASS
+C     ---  integer    n_max_cg
+      INCLUDE 'ngraphs.inc'  !how many diagrams (could be useful to know...)
+
+C     
+C     LOCAL
+C     
+      INTEGER I,J,K
+      REAL*8 P(0:3,NEXTERNAL)  ! four momenta. Energy is the zeroth component.
+      REAL*8 SQRTS  ! sqrt(s)= center of mass energy
+      REAL*8 MATELEM, MATELEMS(0:NSPLITORDERS)
+      REAL*8 PIN(0:3), POUT(0:3)
+      CHARACTER*120 BUFF(NEXTERNAL)
+
+      INTEGER NCHOSEN
+      CHARACTER*20 CHOSEN_SO_INDICES(NSPLITORDERS)
+
+C     
+C     EXTERNAL
+C     
+      REAL*8 DOT
+      EXTERNAL DOT
+
+      LOGICAL CHOSEN_SO_CONFIGS(NSPLITORDERS)
+      COMMON/CHOSEN_BORN_SQSO/CHOSEN_SO_CONFIGS
+
+C     -----
+C     BEGIN CODE
+C     -----
+C     
+C     Start by initializing what is the squared split orders indices
+C      chosen
+      NCHOSEN=0
+      DO I=1,NSPLITORDERS
+        IF (CHOSEN_SO_CONFIGS(I)) THEN
+          NCHOSEN=NCHOSEN+1
+          WRITE(CHOSEN_SO_INDICES(NCHOSEN),'(I3,A1)') I,')'
+        ENDIF
+      ENDDO
+
+C     ---  INITIALIZATION CALLS
+C     
+C     ---  Call to initialize the values of the couplings, masses and
+C      widths 
+C     used in the evaluation of the matrix element. The primary
+C      parameters of the
+C     models are read from Cards/param_card.dat. The secondary
+C      parameters are calculated
+C     in Source/MODEL/couplings.f. The values are stored in common
+C      blocks that are listed
+C     in coupl.inc .
+
+      CALL SETPARA('param_card.dat')  !first call to setup the paramaters
+      INCLUDE 'pmass.inc'  !set up masses
+
+      TOTALMASS = 0.0D0
+      DO I=1,NEXTERNAL
+        TOTALMASS = TOTALMASS + PMASS(I)
+      ENDDO
+
+C     ---  Now use a simple multipurpose PS generator (RAMBO) just to
+C      get a 
+C     RANDOM set of four momenta of given masses pmass(i) to be used
+C      to evaluate 
+C     the madgraph matrix-element.       
+C     Alternatevely, here the user can call or set the four momenta at
+C      his will, see below.
+C     
+      IF(NINCOMING.EQ.1) THEN
+        SQRTS=PMASS(1)
+      ELSE
+        SQRTS=1000D0  !CMS energy in GEV
+        IF (SQRTS.LE.2.0D0*TOTALMASS) THEN
+          SQRTS = 2.0D0*TOTALMASS
+        ENDIF
+      ENDIF
+
+      CALL PRINTOUT()
+
+      CALL GET_MOMENTA(SQRTS,PMASS,P)
+C     
+C     write the information on the four momenta 
+C     
+      WRITE (*,*)
+      WRITE (*,*) ' Phase space point:'
+      WRITE (*,*)
+      WRITE (*,*) '-----------------------------'
+      WRITE (*,*)  'n   E   px   py   pz   m '
+      DO I=1,NEXTERNAL
+        WRITE (*,'(i2,1x,5e15.7)') I, P(0,I),P(1,I),P(2,I),P(3,I)
+     $   ,DSQRT(DABS(DOT(P(0,I),P(0,I))))
+      ENDDO
+      WRITE (*,*) '-----------------------------'
+
+C     
+C     Now we can call the matrix element!
+C     
+      CALL SMATRIX_SPLITORDERS(P,MATELEMS)
+      MATELEM=MATELEMS(0)
+      WRITE(*,*) '1) Matrix element for (QCD=0) = ',MATELEMS(1)
+C     
+      IF (NCHOSEN.NE.NSPLITORDERS) THEN
+        WRITE (*,*) 'Selected squared coupling orders combination fo'
+     $   //'r the sum below:'
+        WRITE (*,*) (CHOSEN_SO_INDICES(I),I=1,NCHOSEN)
+      ENDIF
+      WRITE (*,*) 'Total Matrix element = ', MATELEM, ' GeV^',
+     $ -(2*NEXTERNAL-8)
+      WRITE (*,*) '-----------------------------'
+
+C     c
+C     c      Copy down here (or read in) the four momenta as a string. 
+C     c      
+C     c
+C     buff(1)=" 1   0.5630480E+04  0.0000000E+00  0.0000000E+00 
+C      0.5630480E+04"
+C     buff(2)=" 2   0.5630480E+04  0.0000000E+00  0.0000000E+00
+C      -0.5630480E+04"
+C     buff(3)=" 3   0.5466073E+04  0.4443190E+03  0.2446331E+04
+C      -0.4864732E+04"
+C     buff(4)=" 4   0.8785819E+03 -0.2533886E+03  0.2741971E+03 
+C      0.7759741E+03"
+C     buff(5)=" 5   0.4916306E+04 -0.1909305E+03 -0.2720528E+04 
+C      0.4088757E+04"
+C     c
+C     c      Here the k,E,px,py,pz are read from the string into the
+C      momenta array.
+C     c      k=1,2          : incoming
+C     c      k=3,nexternal  : outgoing
+C     c
+C     do i=1,nexternal
+C     read (buff(i),*) k, P(0,i),P(1,i),P(2,i),P(3,i)
+C     enddo
+C     
+C     c- print the momenta out
+C     
+C     do i=1,nexternal
+C     write (*,'(i2,1x,5e15.7)') i, P(0,i),P(1,i),P(2,i),P(3,i), 
+C     .dsqrt(dabs(DOT(p(0,i),p(0,i))))
+C     enddo
+C     
+C     CALL SMATRIX(P,MATELEM)
+C     
+C     write (*,*) "-------------------------------------------------"
+C     write (*,*) "Matrix element = ", MATELEM, " GeV^",-(2*nexternal-8
+C     )	
+C     write (*,*) "-------------------------------------------------"
+
+      END
+
+
+
+
+      DOUBLE PRECISION FUNCTION DOT(P1,P2)
+C     *****************************************************************
+C     ***********
+C     4-Vector Dot product
+C     *****************************************************************
+C     ***********
+      IMPLICIT NONE
+      DOUBLE PRECISION P1(0:3),P2(0:3)
+      DOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)
+      END
+
+
+      SUBROUTINE GET_MOMENTA(ENERGY,PMASS,P)
+C     ---- auxiliary function to change convention between madgraph
+C      and rambo
+C     ---- four momenta. 	  
+      IMPLICIT NONE
+C     integer nexternal and number particles (incoming+outgoing) in
+C      the me 
+      INTEGER NEXTERNAL, NINCOMING
+      PARAMETER (NEXTERNAL=3,NINCOMING=2)
+C     ARGUMENTS
+      REAL*8 ENERGY,PMASS(NEXTERNAL),P(0:3,NEXTERNAL),PRAMBO(4,10),WGT
+C     LOCAL
+      INTEGER I
+      REAL*8 ETOT2,MOM,M1,M2,E1,E2
+
+      ETOT2=ENERGY**2
+      M1=PMASS(1)
+      M2=PMASS(2)
+      MOM=(ETOT2**2 - 2*ETOT2*M1**2 + M1**4 - 2*ETOT2*M2**2 - 2*M1**2
+     $ *M2**2 + M2**4)/(4.*ETOT2)
+      MOM=DSQRT(MOM)
+      E1=DSQRT(MOM**2+M1**2)
+      E2=DSQRT(MOM**2+M2**2)
+      WRITE (*,*) E1+E2,MOM
+
+      IF(NINCOMING.EQ.2) THEN
+
+        P(0,1)=E1
+        P(1,1)=0D0
+        P(2,1)=0D0
+        P(3,1)=MOM
+
+        P(0,2)=E2
+        P(1,2)=0D0
+        P(2,2)=0D0
+        P(3,2)=-MOM
+
+        CALL RAMBO(NEXTERNAL-2,ENERGY,PMASS(3),PRAMBO,WGT)
+        DO I=3, NEXTERNAL
+          P(0,I)=PRAMBO(4,I-2)
+          P(1,I)=PRAMBO(1,I-2)
+          P(2,I)=PRAMBO(2,I-2)
+          P(3,I)=PRAMBO(3,I-2)
+        ENDDO
+
+      ELSEIF(NINCOMING.EQ.1) THEN
+
+        P(0,1)=ENERGY
+        P(1,1)=0D0
+        P(2,1)=0D0
+        P(3,1)=0D0
+
+        CALL RAMBO(NEXTERNAL-1,ENERGY,PMASS(2),PRAMBO,WGT)
+        DO I=2, NEXTERNAL
+          P(0,I)=PRAMBO(4,I-1)
+          P(1,I)=PRAMBO(1,I-1)
+          P(2,I)=PRAMBO(2,I-1)
+          P(3,I)=PRAMBO(3,I-1)
+        ENDDO
+      ENDIF
+
+      RETURN
+      END
+
+
+      SUBROUTINE RAMBO(N,ET,XM,P,WT)
+C     *****************************************************************
+C     ******
+C     *                       RAMBO                                   
+C           *
+C     *    RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)            
+C           *
+C     *                                                               
+C           *
+C     *    A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR          
+C           *
+C     *    AUTHORS:  S.D. ELLIS,  R. KLEISS,  W.J. STIRLING           
+C           *
+C     *    THIS IS VERSION 1.0 -  WRITTEN BY R. KLEISS                
+C           *
+C     *    -- ADJUSTED BY HANS KUIJF, WEIGHTS ARE LOGARITHMIC
+C      (20-08-90)    *
+C     *                                                               
+C           *
+C     *    N  = NUMBER OF PARTICLES                                   
+C           *
+C     *    ET = TOTAL CENTRE-OF-MASS ENERGY                           
+C           *
+C     *    XM = PARTICLE MASSES ( DIM=NEXTERNAL-nincoming )           
+C           *
+C     *    P  = PARTICLE MOMENTA ( DIM=(4,NEXTERNAL-nincoming) )      
+C           *
+C     *    WT = WEIGHT OF THE EVENT                                   
+C           *
+C     *****************************************************************
+C     ******
+      IMPLICIT REAL*8(A-H,O-Z)
+C     integer nexternal and number particles (incoming+outgoing) in
+C      the me 
+      INTEGER NEXTERNAL, NINCOMING
+      PARAMETER (NEXTERNAL=3,NINCOMING=2)
+      DIMENSION XM(NEXTERNAL-NINCOMING),P(4,NEXTERNAL-NINCOMING)
+      DIMENSION Q(4,NEXTERNAL-NINCOMING),Z(NEXTERNAL-NINCOMING),R(4)
+     $ , B(3),P2(NEXTERNAL-NINCOMING),XM2(NEXTERNAL-NINCOMING)
+     $ , E(NEXTERNAL-NINCOMING),V(NEXTERNAL-NINCOMING),IWARN(5)
+      SAVE ACC,ITMAX,IBEGIN,IWARN
+      DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/
+C     *
+C     * INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
+      IF(IBEGIN.NE.0) GOTO 103
+      IBEGIN=1
+      TWOPI=8.*DATAN(1.D0)
+      PO2LOG=LOG(TWOPI/4.)
+      Z(2)=PO2LOG
+      DO 101 K=3,NEXTERNAL-NINCOMING-1
+ 101  Z(K)=Z(K-1)+PO2LOG-2.*LOG(DFLOAT(K-2))
+      DO 102 K=3,NEXTERNAL-NINCOMING-1
+ 102  Z(K)=(Z(K)-LOG(DFLOAT(K-1)))
+C     *
+C     * CHECK ON THE NUMBER OF PARTICLES
+ 103  IF(N.GT.1.AND.N.LT.101) GOTO 104
+      PRINT 1001,N
+      STOP
+C     *
+C     * CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
+ 104  XMT=0.
+      NM=0
+      DO 105 I=1,N
+      IF(XM(I).NE.0.D0) NM=NM+1
+ 105  XMT=XMT+ABS(XM(I))
+      IF(XMT.LE.ET) GOTO 201
+      PRINT 1002,XMT,ET
+      STOP
+C     *
+C     * THE PARAMETER VALUES ARE NOW ACCEPTED
+C     *
+C     * GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
+ 201  DO 202 I=1,N
+      R1=RN(1)
+      C=2.*R1-1.
+      S=SQRT(1.-C*C)
+      F=TWOPI*RN(2)
+      R1=RN(3)
+      R2=RN(4)
+      Q(4,I)=-LOG(R1*R2)
+      Q(3,I)=Q(4,I)*C
+      Q(2,I)=Q(4,I)*S*COS(F)
+ 202  Q(1,I)=Q(4,I)*S*SIN(F)
+C     *
+C     * CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
+      DO 203 I=1,4
+ 203  R(I)=0.
+      DO 204 I=1,N
+      DO 204 K=1,4
+ 204  R(K)=R(K)+Q(K,I)
+      RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
+      DO 205 K=1,3
+ 205  B(K)=-R(K)/RMAS
+      G=R(4)/RMAS
+      A=1./(1.+G)
+      X=ET/RMAS
+C     *
+C     * TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
+      DO 207 I=1,N
+      BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
+      DO 206 K=1,3
+ 206  P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
+ 207  P(4,I)=X*(G*Q(4,I)+BQ)
+C     *
+C     * CALCULATE WEIGHT AND POSSIBLE WARNINGS
+      WT=PO2LOG
+      IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N)
+      IF(WT.GE.-180.D0) GOTO 208
+      IF(IWARN(1).LE.5) PRINT 1004,WT
+      IWARN(1)=IWARN(1)+1
+ 208  IF(WT.LE. 174.D0) GOTO 209
+      IF(IWARN(2).LE.5) PRINT 1005,WT
+      IWARN(2)=IWARN(2)+1
+C     *
+C     * RETURN FOR WEIGHTED MASSLESS MOMENTA
+ 209  IF(NM.NE.0) GOTO 210
+C     * RETURN LOG OF WEIGHT
+      WT=WT
+      RETURN
+C     *
+C     * MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
+ 210  XMAX=SQRT(1.-(XMT/ET)**2)
+      DO 301 I=1,N
+      XM2(I)=XM(I)**2
+ 301  P2(I)=P(4,I)**2
+      ITER=0
+      X=XMAX
+      ACCU=ET*ACC
+ 302  F0=-ET
+      G0=0.
+      X2=X*X
+      DO 303 I=1,N
+      E(I)=SQRT(XM2(I)+X2*P2(I))
+      F0=F0+E(I)
+ 303  G0=G0+P2(I)/E(I)
+      IF(ABS(F0).LE.ACCU) GOTO 305
+      ITER=ITER+1
+      IF(ITER.LE.ITMAX) GOTO 304
+      PRINT 1006,ITMAX
+      GOTO 305
+ 304  X=X-F0/(X*G0)
+      GOTO 302
+ 305  DO 307 I=1,N
+      V(I)=X*P(4,I)
+      DO 306 K=1,3
+ 306  P(K,I)=X*P(K,I)
+ 307  P(4,I)=E(I)
+C     *
+C     * CALCULATE THE MASS-EFFECT WEIGHT FACTOR
+      WT2=1.
+      WT3=0.
+      DO 308 I=1,N
+      WT2=WT2*V(I)/E(I)
+ 308  WT3=WT3+V(I)**2/E(I)
+      WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET)
+C     *
+C     * RETURN FOR  WEIGHTED MASSIVE MOMENTA
+      WT=WT+WTM
+      IF(WT.GE.-180.D0) GOTO 309
+      IF(IWARN(3).LE.5) PRINT 1004,WT
+      IWARN(3)=IWARN(3)+1
+ 309  IF(WT.LE. 174.D0) GOTO 310
+      IF(IWARN(4).LE.5) PRINT 1005,WT
+      IWARN(4)=IWARN(4)+1
+C     * RETURN LOG OF WEIGHT
+ 310  WT=WT
+      RETURN
+C     *
+ 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
+ 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT SMALLER THA'
+     $ //'N TOTAL ENERGY =',D15.6)
+ 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
+ 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY  OVERFLOW')
+ 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE DESIRE'
+     $ //'D ACCURACY =',D15.6)
+      END
+
+      FUNCTION RN(IDUMMY)
+      REAL*8 RN,RAN
+      SAVE INIT
+      DATA INIT /1/
+      IF (INIT.EQ.1) THEN
+        INIT=0
+        CALL RMARIN(1802,9373)
+        END IF
+C       *
+ 10     CALL RANMAR(RAN)
+        IF (RAN.LT.1D-16) GOTO 10
+        RN=RAN
+C       *
+        END
+
+
+
+        SUBROUTINE RANMAR(RVEC)
+C       *     -----------------
+C       * Universal random number generator proposed by Marsaglia and
+C        Zaman
+C       * in report FSU-SCRI-87-50
+C       * In this version RVEC is a double precision variable.
+        IMPLICIT REAL*8(A-H,O-Z)
+        COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
+        COMMON/ RASET2 / IRANMR,JRANMR
+        SAVE /RASET1/,/RASET2/
+        UNI = RANU(IRANMR) - RANU(JRANMR)
+        IF(UNI .LT. 0D0) UNI = UNI + 1D0
+        RANU(IRANMR) = UNI
+        IRANMR = IRANMR - 1
+        JRANMR = JRANMR - 1
+        IF(IRANMR .EQ. 0) IRANMR = 97
+        IF(JRANMR .EQ. 0) JRANMR = 97
+        RANC = RANC - RANCD
+        IF(RANC .LT. 0D0) RANC = RANC + RANCM
+        UNI = UNI - RANC
+        IF(UNI .LT. 0D0) UNI = UNI + 1D0
+        RVEC = UNI
+        END
+
+        SUBROUTINE RMARIN(IJ,KL)
+C       *     -----------------
+C       * Initializing routine for RANMAR, must be called before
+C        generating
+C       * any pseudorandom numbers with RANMAR. The input values
+C        should be in
+C       * the ranges 0<=ij<=31328 ; 0<=kl<=30081
+        IMPLICIT REAL*8(A-H,O-Z)
+        COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
+        COMMON/ RASET2 / IRANMR,JRANMR
+        SAVE /RASET1/,/RASET2/
+C       * This shows correspondence between the simplified input seeds
+C        IJ, KL
+C       * and the original Marsaglia-Zaman seeds I,J,K,L.
+C       * To get the standard values in the Marsaglia-Zaman paper
+C        (i=12,j=34
+C       * k=56,l=78) put ij=1802, kl=9373
+        I = MOD( IJ/177 , 177 ) + 2
+        J = MOD( IJ     , 177 ) + 2
+        K = MOD( KL/169 , 178 ) + 1
+        L = MOD( KL     , 169 )
+        DO 300 II = 1 , 97
+        S =  0D0
+        T = .5D0
+        DO 200 JJ = 1 , 24
+        M = MOD( MOD(I*J,179)*K , 179 )
+        I = J
+        J = K
+        K = M
+        L = MOD( 53*L+1 , 169 )
+        IF(MOD(L*M,64) .GE. 32) S = S + T
+        T = .5D0*T
+ 200    CONTINUE
+        RANU(II) = S
+ 300    CONTINUE
+        RANC  =   362436D0 / 16777216D0
+        RANCD =  7654321D0 / 16777216D0
+        RANCM = 16777213D0 / 16777216D0
+        IRANMR = 97
+        JRANMR = 33
+        END
+
+
+
+
+
+
+

=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%coef_specs.inc'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%coef_specs.inc	1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%coef_specs.inc	2016-02-24 13:54:50 +0000
@@ -0,0 +1,6 @@
+      INTEGER MAXLWFSIZE
+      PARAMETER (MAXLWFSIZE=4)
+      INTEGER LOOP_MAXCOEFS
+      PARAMETER (LOOP_MAXCOEFS=15)
+      INTEGER VERTEXMAXCOEFS
+      PARAMETER (VERTEXMAXCOEFS=5)

=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%improve_ps.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%improve_ps.f	1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%improve_ps.f	2016-02-24 13:54:50 +0000
@@ -0,0 +1,981 @@
+      SUBROUTINE IMPROVE_PS_POINT_PRECISION(P)
+      IMPLICIT NONE
+C     
+C     CONSTANTS
+C     
+      INTEGER    NEXTERNAL
+      PARAMETER (NEXTERNAL=3)
+C     
+C     ARGUMENTS 
+C     
+      DOUBLE PRECISION P(0:3,NEXTERNAL)
+      REAL*16 QP_P(0:3,NEXTERNAL)
+C     
+C     LOCAL VARIABLES 
+C     
+      INTEGER I,J
+
+C     ----------
+C     BEGIN CODE
+C     ----------
+
+      DO I=1,NEXTERNAL
+        DO J=0,3
+          QP_P(J,I)=P(J,I)
+        ENDDO
+      ENDDO
+
+      CALL MP_IMPROVE_PS_POINT_PRECISION(QP_P)
+
+      DO I=1,NEXTERNAL
+        DO J=0,3
+          P(J,I)=QP_P(J,I)
+        ENDDO
+      ENDDO
+
+      END
+
+
+      SUBROUTINE MP_IMPROVE_PS_POINT_PRECISION(P)
+      IMPLICIT NONE
+C     
+C     CONSTANTS
+C     
+      INTEGER    NEXTERNAL
+      PARAMETER (NEXTERNAL=3)
+C     
+C     ARGUMENTS 
+C     
+      REAL*16 P(0:3,NEXTERNAL)
+C     
+C     LOCAL VARIABLES 
+C     
+      INTEGER I,J
+      INTEGER ERRCODE,ERRCODETMP
+      REAL*16 NEWP(0:3,NEXTERNAL)
+C     
+C     FUNCTIONS
+C     
+      LOGICAL MP_IS_PHYSICAL
+C     
+C     SAVED VARIABLES
+C     
+      INCLUDE 'MadLoopParams.inc'
+C     
+C     SAVED VARIABLES
+C     
+      INTEGER WARNED
+      DATA WARNED/0/
+
+      LOGICAL TOLD_SUPPRESS
+      DATA TOLD_SUPPRESS/.FALSE./
+C     ----------
+C     BEGIN CODE
+C     ----------
+
+C     ERROR CODES CONVENTION
+C     
+C     1         ::  None physical PS point input
+C     100-1000  ::  Error in the origianl method for restoring
+C      precision
+C     1000-9999 ::  Error when restoring precision ala PSMC
+C     
+      ERRCODETMP=0
+      ERRCODE=0
+
+      DO J=1,NEXTERNAL
+        DO I=0,3
+          NEWP(I,J)=P(I,J)
+        ENDDO
+      ENDDO
+
+C     Check the sanity of the original PS point
+      IF (.NOT.MP_IS_PHYSICAL(NEWP,WARNED)) THEN
+        ERRCODE = 1
+        WRITE(*,*) 'ERROR:: The input PS point is not precise enough.'
+        GOTO 100
+      ENDIF
+
+C     Now restore the precision
+      IF (IMPROVEPSPOINT.EQ.1) THEN
+        CALL MP_PSMC_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODE,WARNED)
+      ELSEIF((IMPROVEPSPOINT.EQ.2).OR.(IMPROVEPSPOINT.LE.0)) THEN
+        CALL MP_ORIG_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODE,WARNED)
+      ENDIF
+      IF (ERRCODE.NE.0) THEN
+        IF (WARNED.LT.20) THEN
+          WRITE(*,*) 'INFO:: Attempting to rescue the precisio'
+     $     //'n improvement with an alternative method.'
+          WARNED=WARNED+1
+        ENDIF
+        IF (IMPROVEPSPOINT.EQ.1) THEN
+          CALL MP_ORIG_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODETMP
+     $     ,WARNED)
+        ELSEIF((IMPROVEPSPOINT.EQ.2).OR.(IMPROVEPSPOINT.LE.0)) THEN
+          CALL MP_PSMC_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODETMP
+     $     ,WARNED)
+        ENDIF
+        IF (ERRCODETMP.NE.0) GOTO 100
+      ENDIF
+
+C     Report to the user or update the PS point.
+
+      GOTO 101
+ 100  CONTINUE
+      IF (WARNED.LT.20) THEN
+        WRITE(*,*) 'WARNING:: This PS point could not be improved'
+     $   //'. Error code = ',ERRCODE,ERRCODETMP
+        CALL MP_WRITE_MOM(P)
+        WARNED = WARNED +1
+      ENDIF
+      GOTO 102
+ 101  CONTINUE
+      DO J=1,NEXTERNAL
+        DO I=0,3
+          P(I,J)=NEWP(I,J)
+        ENDDO
+      ENDDO
+ 102  CONTINUE
+
+      IF (WARNED.GE.20.AND..NOT.TOLD_SUPPRESS) THEN
+        WRITE(*,*) 'INFO:: Further warnings from the improve_p'
+     $   //'s routine will now be supressed.'
+        TOLD_SUPPRESS=.TRUE.
+      ENDIF
+
+      END
+
+
+      FUNCTION MP_IS_CLOSE(P,NEWP,WARNED)
+      IMPLICIT NONE
+C     
+C     CONSTANTS
+C     
+      INTEGER    NEXTERNAL
+      PARAMETER (NEXTERNAL=3)
+      REAL*16 ZERO
+      PARAMETER (ZERO=0.0E+00_16)
+      REAL*16 THRS_CLOSE
+      PARAMETER (THRS_CLOSE=1.0E-02_16)
+C     
+C     ARGUMENTS 
+C     
+      REAL*16 P(0:3,NEXTERNAL), NEWP(0:3,NEXTERNAL)
+      LOGICAL MP_IS_CLOSE
+      INTEGER WARNED
+C     
+C     LOCAL VARIABLES 
+C     
+      INTEGER I,J
+      REAL*16 REF,REF2
+      DOUBLE PRECISION BUFFDP
+