anchor.diff:        sl(2)  = 1*nn + 1
         call dcopy(2*nn,0.0d0,0,lmbda,1)
   
  -      n1(1)=1
  -      n2(1)=ne(1)
  -      n1(2)=ne(1)+1
  -      n2(2)=ne(1)+ne(2)
         
   *::::::::::::::::::::::  Lagrangian multipliers  ::::::::::::::::::::::
         DO 640 ms=1,ispin
           IF(ne(ms).le.0) GO TO 640
   
   *       ***** compute the overlap matrices ****
  -        call D3dB_cc_Vector_dot(nfft3d,n,ne(ms),
  -     >                          psi2(1,n1(ms)),
  -     >                          psi2(1,n1(ms)),
  -     >                          tmp(s22))
  -        call D3dB_cc_Vector_dot(nfft3d,n,ne(ms),
  -     >                          psi2(1,n1(ms)),
  -     >                          psi1(1,n1(ms)),
  -     >                          tmp(s21))
  -        call D3dB_cc_Vector_dot(nfft3d,n,ne(ms),
  -     >                          psi1(1,n1(ms)),
  -     >                          psi2(1,n1(ms)),
  -     >                          tmp(s12))
  -        call D3dB_cc_Vector_dot(nfft3d,n,ne(ms),
  -     >                          psi1(1,n1(ms)),
  -     >                          psi1(1,n1(ms)),
  -     >                          tmp(s11))
  +        call Pack_cc_dot(1,psi2(1,ms),psi2(1,ms),s22)
  +        call Pack_cc_dot(1,psi2(1,ms),psi1(1,ms),s21)
  +        call Pack_cc_dot(1,psi1(1,ms),psi2(1,ms),s12)
  +        call Pack_cc_dot(1,psi1(1,ms),psi1(1,ms),s11)
   
   *       ***** scale the overlap matrices ****
  -        do i=1,ne(ms)
  -          index = (i-1) + (i-1)*n
  -      
  -          tmp(s22+index)=(1.0d0-tmp(s22+index))*0.5d0/dte
  -          tmp(s21+index)=(1.0d0-tmp(s21+index))*0.5d0
  -          tmp(s12+index)=(1.0d0-tmp(s12+index))*0.5d0
  -          tmp(s11+index)= (-tmp(s11+index))*0.5d0*dte
  -          
  -           do j=i+1,ne(ms)
  -             index  = (i-1) + (j-1)*n
  -             indext = (j-1) + (i-1)*n
  -   
  -             tmp(s22+index)= (-tmp(s22+index))*0.5d0/dte
  -             tmp(s21+index)= (-tmp(s21+index))*0.5d0
  -             tmp(s12+index)= (-tmp(s12+index))*0.5d0
  -             tmp(s11+index)= (-tmp(s11+index))*0.5d0*dte
  -
  -             tmp(s22+indext)=tmp(s22+index)
  -             tmp(s21+indext)=tmp(s12+index)
  -             tmp(s12+indext)=tmp(s21+index)
  -             tmp(s11+indext)=tmp(s11+index)
  -          end do
  -        end do
  +        s22=(1.0d0-s22)*0.5d0/dte
  +        s21=(1.0d0-s21)*0.5d0
  +        s12=(1.0d0-s12)*0.5d0
  +        s11= (-s11)*0.5d0*dte
   
           call dcopy(nn,tmp(s22),1,tmp(sa0),1)
  +        sa0 = s22
   
           do it=1,itrlmd
  -          CALL dcopy(nn,tmp(s22),1,tmp(sa1),1)
  -          CALL DMMUL(n,ne(MS), tmp(s21), tmp(sa0), tmp(st1))
  -          CALL DMMUL(n,ne(MS), tmp(sa0), tmp(s12), tmp(st2))
  -          CALL DMADD(n,ne(MS), tmp(st1), tmp(st2), tmp(st1))
  -          CALL DMADD(n,ne(MS), tmp(st1), tmp(sa1), tmp(sa1))
  -          CALL DMMUL(n,ne(MS), tmp(s11), tmp(sa0), tmp(st1))
  -          CALL DMMUL(n,ne(MS), tmp(sa0), tmp(st1), tmp(st2))
  -          CALL DMADD(n,ne(MS), tmp(st2), tmp(sa1), tmp(sa1))
  -          CALL DMSUB(n,ne(MS), tmp(sa1), tmp(sa0), tmp(st1))
  -          adiff=tmp(st1 - 1 + (idamax(n*ne(ms),tmp(st1),1)))
  -          if(adiff.lt.convg) GO TO 630
  -          call dcopy(n*ne(ms),tmp(sa1),1,tmp(sa0),1)
  +          sa1 = s22
  +          st1 = s21*sa0
  +          st2 = sa0*s12
  +
  +          st1 = st1+st2
  +          sa1 = st1+sa1
  +
  +          st1 = s11*sa0
  +          st2 = sa0*st1
  +
  +          sa1 = st2+sa1
  +          st1 = sa1-sa0
  +          adiff = st1 
  +          if (adiff.lt.convg) go to 630
  +          sa0 = sa1
  +
  +c          CALL dcopy(nn,tmp(s22),1,tmp(sa1),1)
  +c          CALL DMMUL(n,ne(MS), tmp(s21), tmp(sa0), tmp(st1))
  +c          CALL DMMUL(n,ne(MS), tmp(sa0), tmp(s12), tmp(st2))
  +c
  +c          CALL DMADD(n,ne(MS), tmp(st1), tmp(st2), tmp(st1))
  +c          CALL DMADD(n,ne(MS), tmp(st1), tmp(sa1), tmp(sa1))
  +c
  +c          CALL DMMUL(n,ne(MS), tmp(s11), tmp(sa0), tmp(st1))
  +c          CALL DMMUL(n,ne(MS), tmp(sa0), tmp(st1), tmp(st2))
  +c
  +c          CALL DMADD(n,ne(MS), tmp(st2), tmp(sa1), tmp(sa1))
  +c          CALL DMSUB(n,ne(MS), tmp(sa1), tmp(sa0), tmp(st1))
  +c          adiff=tmp(st1 - 1 + (idamax(n*ne(ms),tmp(st1),1)))
  +
  +c          if(adiff.lt.convg) GO TO 630
  +c          call dcopy(n*ne(ms),tmp(sa1),1,tmp(sa0),1)
           end do
           ierr=10
           WRITE(6,*) 'ierr=10',adiff
   C       return
     630   continue
  -        call dcopy(n*ne(ms),tmp(sa1),1,lmbda(sl(ms)),1)
  +c        call dcopy(n*ne(ms),tmp(sa1),1,lmbda(sl(ms)),1)
  +         lmbda(ms) = sa1
     640 continue
   
   *     **** correction due to the constraint ****
         do ms=1,ispin
  -
  -        do ii=n1(ms),n2(ms)
  -           i=ii-n1(ms)+1
  -
  -           do jj=n1(ms),n2(ms)
  -              j=jj-n1(ms)+1
  -
  -              index = (i-1) + (j-1)*n
  -              alpha = dte*lmbda(sl(ms) + index)
  -
  -c             do k=1,nfft3d
  -c                psi2(k,ii) = psi2(k,ii)  + alpha*psi1(k,jj)
  -c             end do
  -              call D3dB_cc_daxpy(alpha,psi1(1,jj),psi2(1,ii))
  -           end do
  -        end do
  +         alpha = dte*lmbda(ms)
  +         call Pack_cc_daxpy(1,alpha,psi1(1,ms),psi2(1,ms))
         end do
   
  -      call current_second(tim2)
  -      times(3) = times(3) + (tim2-tim1)
   
         return
         end
  
positive.code: *     ***************************
  *     *				*
  *     *		psi_lmbda	*
  *     *				*
  *     ***************************
  
  *     Routine now uses D3dB_Vector routines to reduce latency
  *
  *
        subroutine psi_lmbda(ispin,ne,nemax,
       >                     psi1,psi2,
       >                     dte,
       >                     lmbda,tmp,ierr)
  
        implicit none
        include 'aimd.h'
        integer nfft3d
        parameter (nfft3d=(nfft1/2+1)*nfft2*nslabs) 
        integer ispin,ne(2),nemax
        complex*16 psi1(nfft3d,nemax)
        complex*16 psi2(nfft3d,nemax)
        real*8     dte
        real*8     lmbda(*)
        real*8     tmp(*)
        integer	 ierr
  
  
  *     ***** timing variables *****
  *     1 - fft_time
  *     2 - dot_time
  *     3 - ABC_time
  *     4 - vxc_time
  *     5 - vlocal_time
  *     6 - vnl_time
  *     7 - vh_time
  *     8 - strfac_time
  *     9 - maskC_time
        real*8 tim1,tim2
        real*8 times(50)
        common / timing / times
  
  
  
  *     **** local variables ****
        integer n1(2),n2(2)
        integer i,j,ii,jj,ms,it,k
        integer index,indext,n,nn
        integer s11,s12,s21,s22,st1,st2,sa1,sa0
  
negative.code: 
        call current_second(tim1)
  
        n    = 1
        nn   = n**2
        
        s11  = 0*nn + 1
  