anchor.diff:                 !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)*rnew, 4, zsfun_matvec, zsfun_matvec_transp)
  -               ! anorm_t = znorm_true(ry(i)*rnew, zsfun_matvec, zsfun_matvec_transp)
  -               ! print *, 'anorm, anorm_t = ', anorm, anorm_t
                  call zexp_mv(ry(i)*rnew, order, -tau0, R, Stmp, eps, anorm, zsfun_matvec)
                  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))
  
positive.code:              ra(i+1), phinew(i)%p, crA(pa(i)), crnew(i)%p, crnew(i)%p, phitmp)
               call dinit_sfun(ry(i+1), ry(i+1), ra(i+1), ry(i+1), ry(i+1), phitmp, phinew(i+1)%p)
               anorm = normest(ry(i+1)*ry(i+1),4,dsfun_matvec,dsfun_matvec_transp)
               anorm_t = norm_true(ry(i+1)*ry(i+1), dsfun_matvec,dsfun_matvec_transp)
               print *, 'anorm, anorm_t = ', anorm, anorm_t 
               call dexp_mv(ry(i+1)*ry(i+1), order, -tau/2, R, Stmp, eps, anorm, dsfun_matvec)
               call dgemm('n','n',ry(i)*n(i),ry(i+1),ry(i+1),1d0,crnew(i)%p,ry(i)*n(i),Stmp,ry(i+1),0d0,curcr,ry(i)*n(i))
            else
  
negative.code:              ra(i), ra(i+1), phinew(i)%p, crA(pa(i)), curcr, curcr, phitmp)
               call dinit_sfun(ry(i+1), ry(i+1), ra(i+1), ry(i+1), ry(i+1), phitmp, phinew(i+1)%p)
               anorm = normest(ry(i+1)*ry(i+1), 4, dsfun_matvec, dsfun_matvec_transp)
               call dexp_mv(ry(i+1)*ry(i+1), order, -tau/2, R, Stmp, eps, anorm, dsfun_matvec)
               call dgemm('n','n',ry(i)*n(i),ry(i+1),ry(i+1),1d0,curcr, ry(i)*n(i), Stmp, ry(i+1), 0d0, crnew(i)%p, ry(i)*n(i))
               call dcopy(ry(i)*n(i)*ry(i+1),crnew(i)%p,1,curcr,1)
  