update: messing up with manual merge with devel_MPI.
UPDATE: CODE HAS BEEN MERGED WITH MPI-PARALLEL VERSION.
This new version of the code has MPI instructions
under precompilation flags _MPI.

The code has been tested for the Hubbard model 1-band
and BHZ model 2-bands with positive results.
The results agree within negligible error (10^-15) in parallel
mode. When executed serially it provide the same results at machine
precision.
We are still looking for the possible source of error (suspect is
MPI routines).

THE CODE IS STILL UNDER DEVELOPMENT.

Some minor issues have been solved:
- ED_HAMILTONIAN: sector setup routines have been completely removed.
- ED_HAMILTONIAN: Matrix vector routines have been moved to a new module ED_MATVEC.
- ARPACK_LANCZOS: the routines are now integrated in the module and not included.
- PLAIN_LANCZOS: the routines used to setup the abstract interface for mat-vec
multiplication have been successfully removed. mat-vec operation is passed as
last input to triadiagonalization and plain_lanczos routines.
- ED_GREENSFUNCTIONS: spin susceptibilty routine has been renamed.
NEW VERSION.

This version of the code is a substantial change with respect to previous version.

= parallel with MPI.
= sector symmetries have been implemented in the code to reduce memory usage and increase cpu speed.

The MPI version was implemented in last version but thoroughly tested here.
The small mismatch with MPI and Serial solution has been tracked down to sensitivity of P_ARPACK routines
with respect to tolerance, so that parallel execution possibly produce slightly different eigenvector (but
identical eigenvalues), leading to slightly different number on output. Error on Sigma function is 10^-15.
Error on observables is slightyl larger.

The sector symmetries have been implemented for the non-magnetic, normal channel. This symmetry
enforce the use of *twin* sectors with related quantum numbers (e.g. [nup,ndw]<==>[ndw,nup]).
The eigensolution for Sector A=[nup,ndw] is indentical to that of sector B=[ndw,nup] modulo a
rotation of the eigenbasis, corresponding to a given permutation of eigenvectors elements.

The *ed_twin=T,F variable select usage of symmetric sectors.
If True, only independet sectors are solved and eigenpairs stored. The state_list hosts "ghost"
elements, directly linked to their respective twins so to reduce memory storage.
Upon request the eigen-vector is returned by selectively pointing to the actual vector or its
reordering, depending on the sector contributing to the state_list.
This required some rewriting of the ED_EIGENSPACE routines so to avoid further changes in the
rest of the code.

Modules dependencies have been reshuffled for technical reasons related to implementation.

= New input variables have been included for
 + ground state threshold energy *gs_threshold
 + lanczos tolerance *lanc_tolerance

Minor bugs fixed here and there.

TODO:
1. implement the symmetry sectors in the Superconducting channel (sz==>-sz)
2. implement multi-orbital Superconducting channel, both irred and hybrid solution.