anchor.diff:      integer :: kickrank0, nswp0, verb0
       complex(8), intent(in) :: crA(*), crY0(*)
       real(8), intent(in) :: tau
  +    real(8) :: tau0
       type(zpointd) :: crnew(d+1)
       type(zpointd) :: phinew(d+1)
       complex(8),allocatable, target :: curcr(:)
  
positive.code:     call init_seed()
      ermax = 0d0
      swp = 1
  
      do while (swp .eq. 1)
         !True iteration when started from the left:
         !move (US), move S, next core
  
negative.code:        if ( dir < 0 ) then
             call init_bfun_sizes(ry(i),n(i),ry(i+1),ry(i),n(i),ry(i+1),ra(i),ra(i+1),ry(i)*n(i)*ry(i+1),ry(i)*n(i)*ry(i+1))
             call zinit_bfun_main(phinew(i)%p,crA(pa(i):pa(i+1)-1),phinew(i+1)%p)
  
             !anorm = znormest(ry(i)*n(i)*ry(i+1),4, zmatvec, zmatvec_transp)
             anorm = 1d0
             !test matvec
             !allocate(tmp1(ry(i)*n(i)*ry(i+1)))
             !call zmatvec(crnew(i)%p,tmp1)
             !print *,tmp1(1:4)/crnew(i)%p(1:4)
             !deallocate(tmp1)
             !read(*,*)
             call zexp_mv(ry(i) * n(i) * ry(i+1), 30, tau0, crnew(i)%p, curcr, eps, anorm, zmatvec)
             !print *,'And:',curcr(1:4)/crnew(i)%p(1:4) 
             !print *,'step in core',i
             !read(*,*)
             !call zcopy(ry(i)*n(i)*ry(i+1),crnew(i)%p,1,curcr,1)
             if ( i > 1 ) then
                 call ztransp(ry(i), n(i)*ry(i+1), curcr)
                 rnew = min(n(i)*ry(i+1), ry(i))
                 call zqr(n(i)*ry(i+1), ry(i), curcr, R) 
                 !curcr is n(i)*ry(i+1)xrnew
                 !Test code
                 !call zgemm('n','n',n(i)*ry(i+1),ry(i),ry(i),ONE,curcr,n(i)*ry(i+1),R,ry(i),ZERO,crnew(i)%p,n(i)*ry(i+1))
                 !call zcopy(ry(i)*n(i)*ry(i+1),curcr,1,crnew(i)%p,1)
                 !call ztransp(n(i)*ry(i+1),ry(i),crnew(i)%p)
                 !goto 104
                 !End test code
                 call ztransp(n(i) * ry(i+1), rnew, curcr)
                 call zcopy(rnew * n(i) * ry(i+1), curcr, 1, crnew(i)%p, 1)
                 call ztransp(rnew, ry(i), R)
  
                 call zphi_right(rnew, n(i), ry(i+1), rnew, n(i), ry(i+1), &
                 ra(i), ra(i+1), phinew(i+1)%p, crA(pa(i)), curcr, curcr, phitmp)
                 !phitmp is now ry(i) x ra(i) x ry(i) 
  
                 call zinit_sfun(ry(i), ry(i), ra(i), rnew, rnew, phinew(i)%p, phitmp)
                 !anorm = znormest(ry(i+1)*ry(i+1), 4, zsfun_matvec, zsfun_matvec_transp)
                 anorm = 1d0
                 call zexp_mv(ry(i)*rnew, 30, -tau0, R, Stmp, eps, anorm, zsfun_matvec)
                 !print *,Stmp(1:4)/R(1:4)
                 !read(*,*)
                 !print *,'step in R-matrix',i
                 !call zcopy(ry(i)*rnew,R,1,Stmp,1)
                 call zgemm('n','n',ry(i-1)*n(i-1),rnew,ry(i),ONE,crnew(i-1)%p,ry(i-1)*n(i-1),Stmp,ry(i),ZERO,curcr,ry(i-1)*n(i-1))
                 ry(i) = rnew
                 call zcopy(ry(i-1)*n(i-1)*ry(i), curcr, 1, crnew(i-1)%p, 1)
                 if ( size(phinew(i)%p) < ry(i)*ry(i)*ra(i) ) then
                     deallocate(phinew(i)%p)
                     allocate(phinew(i)%p(ry(i)*ry(i)*ra(i)*2))
                 end if
                 call zcopy(ry(i)*ra(i)*ry(i),phitmp,1,phinew(i)%p,1) !Update phi  
             else !i == 1 
                 call zcopy(ry(i)*n(i)*ry(i+1),curcr,1,crnew(i)%p,1) 
             end if
         else   ! dir > 0
            call init_bfun_sizes(ry(i),n(i),ry(i+1),ry(i),n(i),ry(i+1),ra(i),ra(i+1),ry(i)*n(i)*ry(i+1),ry(i)*n(i)*ry(i+1))
            call zinit_bfun_main(phinew(i)%p,crA(pa(i):pa(i+1)-1),phinew(i+1)%p)
            
            !anorm = znormest(ry(i)*n(i)*ry(i+1),4, zmatvec, zmatvec_transp)
            anorm = 1d0
            call zexp_mv(ry(i) * n(i) * ry(i+1), 30, tau0, crnew(i)%p, curcr, eps, anorm, zmatvec)
            if ( i < d ) then ! Update S
                call zqr(ry(i) * n(i), ry(i+1), curcr, R) 
                rnew = min(ry(i) * n(i), ry(i+1))
                call zcopy(ry(i) * n(i) * rnew, curcr, 1, crnew(i)%p, 1)
  
                call zphi_left(ry(i), n(i), rnew, ry(i), n(i), rnew, &
                ra(i), ra(i+1), phinew(i)%p, crA(pa(i)), curcr, curcr, phitmp)
                call zinit_sfun(rnew, rnew, ra(i+1), ry(i+1), ry(i+1), phitmp, phinew(i+1)%p)
                !anorm = znormest(ry(i+1)*ry(i+1), 4, zsfun_matvec, zsfun_matvec_transp)
                anorm = 1d0
                call zexp_mv(rnew * ry(i+1), 30, -tau0, R, Stmp, eps, anorm, zsfun_matvec)
                call zgemm('n', 'n', rnew, n(i+1) * ry(i+2), ry(i+1), ONE, &
                            Stmp, rnew, crnew(i+1)%p, ry(i+1), &
                            ZERO, curcr, rnew)
                call zcopy(rnew * n(i+1) * ry(i+2), curcr, 1, crnew(i+1) % p, 1) 
  
                ry(i+1) = rnew
                call zcopy(ry(i+1) * ra(i+1) * ry(i+1), phitmp, 1, phinew(i+1) % p, 1)
            else  ! Last core
                call zcopy(ry(i) * n(i) * ry(i+1), curcr, 1, crnew(i) % p, 1 ) 
            endif
  
        end if
         
  
  
  
         if ((dir>0) .and. (i==d )) then
             dir = -1
             i = d
             swp = swp + 1
             ermax = 0d0
         else if ((dir < 0) .and. (i == 1 )) then
             dir = 1
             i = 1
             if ( typ .eq. 1 ) then
                 goto 105
             endif
         else
             i = i + dir
         end if
      end do
  
  105 continue
      print *, 'Goint to exit'
      nn = sum(ry(2:d+1)*ry(1:d)*n(1:d))
      if ( allocated(zresult_core)) then
         if ( size(zresult_core) < nn ) then
  