anchor.diff:              !deallocate(tmp3)
               !deallocate(X)
               
  -            call zinit_sfun(ry(i), ry(i), ra(i), ry(i), ry(i), phinew(i)%p, phitmp)
  -           sz = ry(i)*ry(i)
  +            call zinit_sfun(ry(i), ry(i), ra(i), rnew, rnew, phinew(i)%p, phitmp)
  +            !sz = rnew * ry(i)
               !allocate(X(sz,sz))
               !allocate(tmp1(sz))
               !do j = 1,sz
  
positive.code:           else !i == 1 
              call zcopy(ry(i)*n(i)*ry(i+1),curcr,1,crnew(i)%p,1) 
            end if
            
         end if
  
  
  
negative.code:           
              !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,tau,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) 
              !allocate(tmp3(ry(i),ra(i),ry(i)))
  