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.

update: added in the code the ed_verbose=-1 level for more verbosity than standard, to
be used when Ns>10 and the code is slow.

TODO:
1. implement the symmetry sectors in the Superconducting channel (sz==>-sz)
2. implement multi-orbital Superconducting channel, both irred and hybrid solution.
UPDATE:
- added Compressed Storage Row format sparse matrix in the MATRIX_SPARSE module.
the added routine allows to handle sparse matrices in CSR format.
The CSR matrix can be created from a LL matrix in order to include only the
non-zero element.
CSR matrices are expected to perform better in the matrix-vector multiplication.

MPI Mat-Vev multiplication for CSR is not yet included.
added printing of Green's function poles and weights.
a possible bug in ED_GREENS_FUNCTION related to ed_verbose has been fixed.

some updates to drivers/ed_bhz.f90 related to our work.
update of driver/ed_hm_bethe.f90 to test FFT against W2Dynamics.
update version following SciFor2.0 version.
- simplification in the Makefile with explicit linking libraries.
- small bug fixed in PLAIN_LANCZOS, including variables in a wrong way
from ED_VARS_GLOBAL and ED_INPUT_VARS.
internal energy calculation introduced.

The internal energy of the problem is evaluated using information from the
impurity problem only, so that user is not bothered by performing any additional
calculation.

Potential energy contribution is evaluated as:
<Hi>=Tr[ H_interacting ]
where H_interacting is the part of the model Hamiltonian describing full local interaction,
there including the Hartree terms (e.g. for example if H_interacting = U*(n-1/2)*(n-1/2)
rather than U*n*n).

Kinetic energy contribution is evaluated as:
<K> = Tr{Tr[H(k)*G(k)]} = Tr[Delta*G] + Tr[Eimp*N]
the internal trace being thermal one (sum_{iw,k}), the external one is  over spin-orbital degrees of freedom.
This is easy to check
Tr{ H(k)*G(k) } = -1 + Tr{ (iw + mu - Sigma(iw))*G(k)}
                = -1 + Tr{ (iw + mu - G0^-1 - Gimp)*G(k)}
		= -1 + Tr{ (Delta + Eimp - Gimp)*G(k)}
We separated the Eimp contribution (<E0>) from <K>, so that
<K> = Tr[Delta*Gimp]

The other terms are evaluated from thermal averages of suitable terms.

testing: a simple test has shown very good agreement with former results for the
1-band Bethe Hubbard model. More tests are in order.

The introduction of Energy calculations required to move the impurity Green's function
outside the ED_GREENS_FUNCTIONS module, so I opted to put them in VARS_GLOBAL so to be
able to access the as Self-energies.

With the occasion I updated the code here and there and renamed some functions.
ED_HAMILTONIAN now is unified for both double precision and double complex calculations.
Something similar can be done for the ED_DIAG.
Merge branch 'devel'

Conflicts:
	DMFT_ED.f90
	ED_DIAG.f90
	ED_GREENS_FUNCTIONS.f90
	ED_OBSERVABLES.f90
	Makefile
	ed_build_gf_normal.f90