Some cleanup and marginal improvements in force calculation for US PP

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@144 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2003-04-04 14:36:00 +00:00
parent 0efed33427
commit 6e3853d809
11 changed files with 207 additions and 311 deletions

View File

@ -94,6 +94,7 @@ PWOBJS=../PW/pwcom.o \
../PW/multable.o \
../PW/n_plane_waves.o \
../PW/new_ns.o \
../PW/newd.o \
../PW/openfil.o \
../PW/potinit.o \
../PW/print_clock_pw.o \
@ -214,7 +215,6 @@ h_psi.o \
init_run.o \
interpolate.o \
mix_rho.o \
newd.o \
pw_gemm.o \
read_file.o \
rotate_wfc.o \

View File

@ -21,15 +21,8 @@ subroutine addusdens
! here the local variables
!
integer :: ig, na, nt, ih, jh, ijh, ir, is
! counter on G vectors
! counter on atoms
! the atom type
! counter on beta functions
! counter on beta functions
! composite index ih jh
! counter on mesh points
! counter on spin polarization
integer :: ig, na, nt, ih, jh, ijh, is
! counters
real(kind=DP), allocatable :: qmod (:), ylmk0 (:,:)
! the modulus of G
@ -41,14 +34,14 @@ subroutine addusdens
! work space for rho(G,nspin)
if (.not.okvan) return
call start_clock ('addusdens')
allocate (aux ( ngm, nspin))
allocate (qg( nrxx))
allocate (qmod( ngm))
allocate (ylmk0( ngm, lqx * lqx))
call setv (2 * ngm * nspin, 0.d0, aux, 1)
aux (:,:) = (0.d0, 0.d0)
call ylmr2 (lqx * lqx, ngm, g, gg, ylmk0)
do ig = 1, ngm
qmod (ig) = sqrt (gg (ig) )
@ -80,24 +73,22 @@ subroutine addusdens
endif
enddo
!
! convert aux to real space and add to the charge density
!
if (okvan) then
do is = 1, nspin
call setv (2 * nrxx, 0.d0, qg, 1)
do ig = 1, ngm
qg (nl (ig) ) = aux (ig, is)
qg (nlm(ig) ) = conjg(aux(ig, is))
enddo
call cft3 (qg, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
do ir = 1, nrxx
rho (ir, is) = rho (ir, is) + DREAL (qg (ir) )
enddo
enddo
endif
deallocate (ylmk0)
deallocate (qmod)
!
! convert aux to real space and add to the charge density
!
allocate (qg( nrxx))
do is = 1, nspin
qg(:) = (0.d0, 0.d0)
do ig = 1, ngm
qg (nl (ig) ) = aux (ig, is)
qg (nlm(ig) ) = conjg(aux(ig, is))
enddo
call cft3 (qg, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
rho (:, is) = rho (:, is) + DREAL (qg (:) )
enddo
deallocate (qg)
deallocate (aux)

View File

@ -107,7 +107,7 @@ subroutine force_us (forcenl)
!
call addusforce (forcenl)
#ifdef __PARA
! collect contributions across pool
! collect contributions across pools
call poolreduce (3 * nat, forcenl)
#endif
!

View File

@ -1,146 +0,0 @@
!
! Copyright (C) 2003 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!----------------------------------------------------------------------
subroutine newd
!----------------------------------------------------------------------
!
! This routine computes the integral of the effective potential with
! the Q function and adds it to the bare ionic D term which is used
! to compute the non-local term in the US scheme.
! LSDA case allowed: nspin = 2
!
#include "machine.h"
use pwcom
implicit none
integer :: ig, ir, nt, ih, jh, na, is
! counter on g vectors
! counter on mesh points
! the atom type
! counter on atomic beta functions
! counter on atomic beta functions
! counter on atoms
! counter on spin polarizations
complex(kind=DP), allocatable :: aux (:,:), vg (:)
! used to contain the potential
! used to contain the structure
real(kind=DP), allocatable :: ylmk0 (:,:), qmod (:)
! the spherical harmonics
! the modulus of G
real(kind=DP) :: DDOT
!
!
if (.not.okvan) then
! no ultrasoft potentials: use bare coefficients for projectors
do na = 1, nat
nt = ityp (na)
do is = 1, nspin
do ih = 1, nh (nt)
do jh = ih, nh (nt)
deeq (ih, jh, na, is) = dvan (ih, jh, nt)
deeq (jh, ih, na, is) = deeq (ih, jh, na, is)
enddo
enddo
enddo
end do
return
end if
call start_clock ('newd')
allocate (aux ( ngm , nspin))
allocate (vg ( nrxx))
allocate (qmod ( ngm))
allocate (ylmk0( ngm , lqx * lqx))
!
! set deeq to zero
!
call setv (nhm * nhm * nat * nspin, 0.d0, deeq, 1)
call ylmr2 (lqx * lqx, ngm, g, gg, ylmk0)
do ig = 1, ngm
qmod (ig) = sqrt (gg (ig) )
enddo
!
! fourier transform of the total effective potential
!
do is = 1, nspin
do ir = 1, nrxx
vg (ir) = vltot (ir) + vr (ir, is)
enddo
call cft3 (vg, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
do ig = 1, ngm
aux (ig, is) = vg (nl (ig) )
enddo
enddo
!
! here we compute the integral Q*V for each atom,
! I = sum_G exp(-iR.G) Q_nm v^*
!
do nt = 1, ntyp
if (tvanp (nt) ) then
do ih = 1, nh (nt)
do jh = ih, nh (nt)
call qvan2 (ngm, ih, jh, nt, qmod, qgm, ylmk0)
do na = 1, nat
if (ityp (na) .eq.nt) then
!
! The product of the potential and the structure factor
!
do is = 1, nspin
do ig = 1, ngm
vg (ig) = aux (ig, is) * &
conjg (eigts1 (ig1 (ig), na) &
* eigts2 (ig2 (ig), na) &
* eigts3 (ig3 (ig), na) )
enddo
!
! and the product with the Q functions
!
deeq (ih, jh, na, is) = omega * 2.d0 * &
DDOT (2 * ngm, vg, 1, qgm, 1)
if (gstart==2) deeq (ih, jh, na, is) = &
deeq (ih, jh, na, is)-omega*real(vg(1)*qgm(1))
enddo
endif
enddo
enddo
enddo
endif
enddo
#ifdef __PARA
call reduce (nhm * nhm * nat * nspin, deeq)
#endif
do na = 1, nat
nt = ityp (na)
do is = 1, nspin
! write(6,'( "dmatrix atom ",i4, " spin",i4)') na,is
! do ih = 1, nh(nt)
! write(6,'(8f9.4)') (deeq(ih,jh,na,is),jh=1,nh(nt))
! end do
do ih = 1, nh (nt)
do jh = ih, nh (nt)
deeq (ih, jh, na, is) = deeq (ih, jh, na, is) + dvan (ih, jh, nt)
deeq (jh, ih, na, is) = deeq (ih, jh, na, is)
enddo
enddo
enddo
! write(6,'( "dion pseudo ",i4)') nt
! do ih = 1, nh(nt)
! write(6,'(8f9.4)') (dvan(ih,jh,nt),jh=1,nh(nt))
! end do
enddo
deallocate (ylmk0)
deallocate (qmod)
deallocate (vg)
deallocate (aux)
call stop_clock ('newd')
return
end subroutine newd

View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
@ -20,15 +20,8 @@ subroutine addusdens
! here the local variables
!
integer :: ig, na, nt, ih, jh, ijh, ir, is
! counter on G vectors
! counter on atoms
! the atom type
! counter on beta functions
! counter on beta functions
! composite index ih jh
! counter on mesh points
! counter on spin polarization
integer :: ig, na, nt, ih, jh, ijh, is
! counters
real(kind=DP), allocatable :: qmod (:), ylmk0 (:,:)
! the modulus of G
@ -39,15 +32,15 @@ subroutine addusdens
! work space for FFT
! work space for rho(G,nspin)
if (.not.okvan) return
call start_clock ('addusdens')
allocate (aux ( ngm, nspin))
allocate (qg( nrxx))
allocate (qmod( ngm))
allocate (ylmk0( ngm, lqx * lqx))
call setv (2 * ngm * nspin, 0.d0, aux, 1)
aux (:,:) = (0.d0, 0.d0)
call ylmr2 (lqx * lqx, ngm, g, gg, ylmk0)
do ig = 1, ngm
qmod (ig) = sqrt (gg (ig) )
@ -79,24 +72,20 @@ subroutine addusdens
endif
enddo
!
! convert aux to real space and add to the charge density
!
if (okvan) then
do is = 1, nspin
call setv (2 * nrxx, 0.d0, qg, 1)
do ig = 1, ngm
qg (nl (ig) ) = aux (ig, is)
enddo
call cft3 (qg, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
do ir = 1, nrxx
rho (ir, is) = rho (ir, is) + DREAL (qg (ir) )
enddo
enddo
endif
deallocate (ylmk0)
deallocate (qmod)
!
! convert aux to real space and add to the charge density
!
allocate (qg( nrxx))
do is = 1, nspin
qg(:) = (0.d0, 0.d0)
do ig = 1, ngm
qg (nl (ig) ) = aux (ig, is)
enddo
call cft3 (qg, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
rho (:, is) = rho (:, is) + DREAL (qg (:) )
enddo
deallocate (qg)
deallocate (aux)

View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
@ -35,16 +35,6 @@ subroutine addusforce (forcenl)
fact = 1.d0
end if
allocate (aux(ngm,nspin))
allocate (aux1(ngm,3))
allocate (ddeeq( 3, (nhm*(nhm+1))/2,nat,nspin))
allocate (qmod( ngm))
allocate (ylmk0(ngm,lqx*lqx))
!
ddeeq(:,:,:,:) = 0.d0
!
call ylmr2 (lqx * lqx, ngm, g, gg, ylmk0)
!
qmod (:) = sqrt (gg (:) )
!
! fourier transform of the total effective potential
!
@ -56,6 +46,17 @@ subroutine addusforce (forcenl)
enddo
deallocate (vg)
!
allocate (aux1(ngm,3))
allocate (ddeeq( 3, (nhm*(nhm+1))/2,nat,nspin))
allocate (qmod( ngm))
allocate (ylmk0(ngm,lqx*lqx))
!
ddeeq(:,:,:,:) = 0.d0
!
call ylmr2 (lqx * lqx, ngm, g, gg, ylmk0)
!
qmod (:) = sqrt (gg (:) )
!
! here we compute the integral Q*V for each atom,
! I = sum_G i G_a exp(-iR.G) Q_nm v^*
!

View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
@ -24,24 +24,21 @@ subroutine force_corr (forcescc)
real(kind=DP) :: forcescc (3, nat)
!
real(kind=DP), allocatable :: rhocgnt (:), aux (:)
! work space
real(kind=DP) :: gx, arg, fact
integer :: ir, is, ig, igl0, nt, na, ipol
! counter on mesh points
! counter on spin polarizations
! counter on G vectors
! index of first shell with G != 0
! counter on atom types
! " " atoms
! " " cartesians components
!
!
call setv (2 * nrxx, 0.d0, psic, 1)
! temp factors
integer :: ir, isup, isdw, ig, igl0, nt, na, ipol
! counters
!
! vnew is V_out - V_in, psic is the temp space
!
do is = 1, nspin
call DAXPY (nrxx, 1.d0 / nspin, vnew (1, is), 1, psic, 2)
enddo
if (nspin == 1) then
psic(:) = vnew (:, 1)
else
isup = 1
isdw = 2
psic(:) = (vnew (:, isup) + vnew (:, isdw)) * 0.5d0
end if
!
allocate ( aux(ndm), rhocgnt(ngl) )
@ -70,9 +67,9 @@ subroutine force_corr (forcescc)
enddo
call simpson (msh (nt), aux, rab (1, nt), rhocgnt (ig) )
enddo
do ig = gstart, ngm
do na = 1, nat
if (nt.eq.ityp (na) ) then
do na = 1, nat
if (nt.eq.ityp (na) ) then
do ig = gstart, ngm
arg = (g (1, ig) * tau (1, na) + g (2, ig) * tau (2, na) &
+ g (3, ig) * tau (3, na) ) * tpi
do ipol = 1, 3
@ -80,8 +77,8 @@ subroutine force_corr (forcescc)
rhocgnt (igtongl(ig) ) * DCMPLX(sin(arg),cos(arg)) * &
g(ipol,ig) * tpiba * conjg(psic(nl(ig)))
enddo
endif
enddo
enddo
endif
enddo
enddo
#ifdef __PARA

View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
@ -23,36 +23,21 @@ subroutine force_us (forcenl)
real(kind=DP) :: forcenl (3, nat)
! output: the nonlocal contribution
complex(kind=DP), allocatable :: dbecp (:,:,:), work1 (:,:), aux1 (:)
complex(kind=DP), allocatable :: dbecp (:,:,:)
! auxiliary variable contains <dbeta|psi>
! auxiliary variable contains g*psi
! auxiliary variable contains beta(k+g)
complex(kind=DP) :: fact, ZDOTC
! a multiplicative factor
! the scalar product function
real(kind=DP) :: fac, ps
integer :: ik, ipol, ibnd, ig, ih, jh, na, nt, ikb, jkb, jkb2, &
ijkb0
! counter on k points
! counter on polarization
! counter on bands
! counter on G vectors
! counters on atomic beta functions
! counter on atoms
! counter on type of atoms
! counters on beta functions
complex(kind=DP), allocatable :: vkb1 (:,:)
! auxiliary variable contains g*|beta>
real(kind=DP) :: ps
integer :: ik, ipol, ibnd, ig, ih, jh, na, nt, ikb, jkb, ijkb0
! counters
!
forcenl(:,:) = 0.d0
allocate (dbecp( nkb, 3, nbnd))
allocate (work1( npwx, 3))
allocate (aux1( npwx))
allocate (dbecp( nkb, nbnd, 3))
allocate (vkb1( npwx, nkb))
if (nks.gt.1) rewind iunigk
!
! the forces are a sum over the K points and the bands
!
fact = DCMPLX (0.d0, tpiba)
do ik = 1, nks
if (lsda) current_spin = isk (ik)
if (nks.gt.1) then
@ -60,29 +45,19 @@ subroutine force_us (forcenl)
call davcio (evc, nwordwfc, iunwfc, ik, - 1)
call init_us_2 (npw, igk, xk (1, ik), vkb)
endif
!
call ccalbec (nkb, npwx, npw, nbnd, becp, vkb, evc)
jkb2 = 0
do nt = 1, ntyp
do na = 1, nat
if (ityp (na) .eq.nt) then
do ih = 1, nh (nt)
jkb2 = jkb2 + 1
do ig = 1, npw
do ipol = 1, 3
work1 (ig, ipol) = vkb(ig,jkb2)* g(ipol,igk(ig) )
enddo
enddo
do ibnd = 1, nbnd
do ipol = 1, 3
dbecp (jkb2, ipol, ibnd) = fact * &
ZDOTC (npw,work1(1, ipol),1,evc(1,ibnd),1)
enddo
enddo
enddo
endif
!
do ipol = 1, 3
do jkb = 1,nkb
do ig = 1, npw
vkb1 (ig, jkb) = vkb(ig,jkb) *(0.d0,-1.d0)*g(ipol,igk(ig) )
enddo
enddo
enddo
!
call ZGEMM ('C', 'N', nkb, nbnd, npw, (1.d0, 0.d0), vkb1, npwx, &
evc, npwx, (0.d0, 0.d0), dbecp(1,1,ipol), nkb)
end do
!
ijkb0 = 0
do nt = 1, ntyp
@ -93,11 +68,10 @@ subroutine force_us (forcenl)
do ibnd = 1, nbnd
ps = deeq (ih, ih, na, current_spin) - &
et (ibnd, ik) * qq (ih,ih, nt)
fac = wg (ibnd, ik)
do ipol = 1, 3
forcenl (ipol, na) = forcenl (ipol, na) - &
ps * fac * 2.d0 * &
DREAL (conjg (dbecp (ikb, ipol, ibnd) ) &
ps * wg (ibnd, ik) * 2.d0 * tpiba * &
DREAL (conjg (dbecp (ikb, ibnd, ipol) ) &
* becp (ikb, ibnd) )
enddo
enddo
@ -110,15 +84,14 @@ subroutine force_us (forcenl)
do jh = ih + 1, nh (nt)
jkb = ijkb0 + jh
do ibnd = 1, nbnd
fac = wg (ibnd, ik)
ps = deeq (ih, jh, na, current_spin) - &
et (ibnd, ik) * qq (ih, jh, nt)
do ipol = 1, 3
forcenl (ipol, na) = forcenl (ipol, na) - &
ps * fac * 2.d0 * &
DREAL (conjg (dbecp (ikb, ipol, ibnd) ) * &
ps * wg (ibnd, ik) * 2.d0 * tpiba * &
DREAL (conjg (dbecp (ikb, ibnd, ipol) ) * &
becp (jkb, ibnd) + &
dbecp (jkb, ipol, ibnd) * &
dbecp (jkb, ibnd, ipol) * &
conjg (becp (ikb, ibnd) ) )
enddo
enddo
@ -133,14 +106,17 @@ subroutine force_us (forcenl)
#ifdef __PARA
call reduce (3 * nat, forcenl)
#endif
deallocate (vkb1)
deallocate (dbecp)
!
! The total D matrix depends on the ionic position due to the augmentati
! The total D matrix depends on the ionic position via the augmentation
! part \int V_eff Q dr, the term deriving from the derivative of Q
! is added in the routine addusforce
!
call addusforce (forcenl)
!
#ifdef __PARA
! collect contributions across poo
! collect contributions across pools
call poolreduce (3 * nat, forcenl)
#endif
!
@ -161,11 +137,7 @@ subroutine force_us (forcenl)
!
do na = 1, nat
call trnvect (forcenl (1, na), at, bg, 1)
enddo
deallocate (aux1)
deallocate (work1)
deallocate (dbecp)
return
end subroutine force_us

View File

@ -48,7 +48,6 @@ subroutine forces
!
! The local contribution
!
call force_lc (nat, tau, ityp, alat, omega, ngm, ngl, igtongl, &
nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, g, rho, nl, nspin, &
gstart, gamma_only, vloc, forcelc)
@ -69,7 +68,6 @@ subroutine forces
! The SCF contribution
!
call force_corr (forcescc)
!
! here we sum all the contributions and compute the total force acting
! on the crstal

View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2003 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
@ -13,20 +13,18 @@ subroutine newd
! This routine computes the integral of the effective potential with
! the Q function and adds it to the bare ionic D term which is used
! to compute the non-local term in the US scheme.
! LSDA case allowed: nspin = 2
!
#include "machine.h"
use pwcom
implicit none
integer :: ig, ir, nt, ih, jh, na, is
! counters on g vectors, mesh points, atom type, beta functions x 2, atoms
! and spin polarizations
complex(kind=DP), allocatable :: aux (:,:), & ! used to contain the potential
vg (:) ! used to contain the structure
real(kind=DP), allocatable :: ylmk0 (:,:), & ! the spherical harmonics
qmod (:) ! the modulus of G
real(kind=DP) :: DDOT
integer :: ig, nt, ih, jh, na, is
! counters on g vectors, atom type, beta functions x 2, atoms, spin
complex(kind=DP), allocatable :: aux (:,:), vg (:)
! work space
real(kind=DP), allocatable :: ylmk0 (:,:), qmod (:)
! spherical harmonics, modulus of G
real(kind=DP) :: fact, DDOT
!
!
if (.not.okvan) then
@ -45,13 +43,16 @@ subroutine newd
return
end if
if (gamma_only) then
fact = 2.d0
else
fact = 1.d0
end if
call start_clock ('newd')
allocate ( aux(ngm,nspin), vg(nrxx), qmod(ngm), ylmk0(ngm, lqx*lqx) )
!
! set deeq to zero
!
deeq(:,:,:,:) = 0.d0
!
call ylmr2 (lqx * lqx, ngm, g, gg, ylmk0)
do ig = 1, ngm
qmod (ig) = sqrt (gg (ig) )
@ -60,9 +61,7 @@ subroutine newd
! fourier transform of the total effective potential
!
do is = 1, nspin
do ir = 1, nrxx
vg (ir) = vltot (ir) + vr (ir, is)
enddo
vg (:) = vltot (:) + vr (:, is)
call cft3 (vg, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
do ig = 1, ngm
aux (ig, is) = vg (nl (ig) )
@ -91,8 +90,11 @@ subroutine newd
!
! and the product with the Q functions
!
deeq (ih, jh, na, is) = omega * &
deeq (ih, jh, na, is) = fact * omega * &
DDOT (2 * ngm, vg, 1, qgm, 1)
if (gamma_only .and. gstart==2) &
deeq (ih, jh, na, is) = &
deeq (ih, jh, na, is) - omega*real(vg(1)*qgm(1))
enddo
endif
enddo

92
TODO Normal file
View File

@ -0,0 +1,92 @@
TODO LIST - 2 Apr 2003
INSTALLATION
- compile fftw by default on request for PW and CP as well
COMMON
- fix shift in eigenvalues between CP/FPMD and PW
The eigenvalues must be the same!!!
- fix shift in total energy between CP, FPMD and PW
for charged systems. Total energies for charged systems
must be the same!!!
- merge diagonalization routines, use parallel diag everywhere
- start merge of input routines
- merge spherical harmonics, bessel, integration routines, etc
- variable declarations : real(kind=DP) etc everywhere
- real vs dreal, cmplx vs dcmplx, dfloat etc
PW
- bfgs, md : atomic positions should be written in the same
format as they are read (but this should not spoil scripts
that extract coordinates from output file)
- remove analytical pps and related variables
- remove residual direct calls to MPI routines
use routines in mp.f90 instead
merge with existing wrapper if needed
- remove potential mixing, save and start from rho instead of V
- remove all calls to setv and to blas copy, scal
- save rho and/or V in G-space instead of R-space
- remove old punch, use new punch, add possibility to save
wavefunctions either 'distributed' or 'collected' in
parallel execution (requires identifier in the first line
of the file)
- add possibility to read atomic positions from file
- read a, b, c, cosab, cosbc, cosac instead of celldm
Use more standard choices for crystal axis
- fix problem with j'_l(x)Y_lm(x) for x -> 0 and l=1
- memory estimator tool, or memory report
POSTPROCESSING
- bands.x must either be finished or removed
- postprocessing in the parallel case and with new_punch
must be verified
- stm in non-scf calculations to be verified
- add more scripts that process output file
- dos, projected dos, etc: input data should be more uniform
PH
- Tone: ntyp in input needed for phonon GUI ?
- better algorithm for electron-phonon (Malgorzata?)
DOCUMENTATION
- expand, update, verify, etc
- add a list of FAQ, or AAQ (Already Answered Questions)
CPV:
- find reason for differences in empty states between CP and PW
with ultrasoft pps and box grids. If it is a bug, it is serious.
- check on input ipp and pp type. Even better: remove ipp
- remove eispack routines, replace with lapack or parallel diag
- remove all calls to zero and to blas copy, scal