diff --git a/.gitignore b/.gitignore
index c86a539e6..4d95c5bbe 100644
--- a/.gitignore
+++ b/.gitignore
@@ -19,6 +19,8 @@ sample/hubbard_2d
 sample/hubbard_2d-g
 sample/trg
 sample/trg-g
+sample/ctmrg
+sample/ctmrg-g
 this_dir.mk
 options.mk
 .exrc
diff --git a/sample/Makefile b/sample/Makefile
index b401df3bc..e4133f69c 100644
--- a/sample/Makefile
+++ b/sample/Makefile
@@ -19,11 +19,11 @@ LIBGFLAGS=-L'$(ITENSOR_LIBDIR)' $(ITENSOR_LIBGFLAGS)
 
 #Targets -----------------
 
-build: dmrg dmrg_table dmrgj1j2 exthubbard trg mixedspin hubbard_2d
+build: dmrg dmrg_table dmrgj1j2 exthubbard trg ctmrg mixedspin hubbard_2d
 
-debug: dmrg-g dmrg_table-g dmrgj1j2-g exthubbard-g trg-g mixedspin-g hubbard_2d-g
+debug: dmrg-g dmrg_table-g dmrgj1j2-g exthubbard-g trg-g ctmrg-g mixedspin-g hubbard_2d-g
 
-all: dmrg dmrg_table dmrgj1j2 exthubbard trg mixedspin hubbard_2d
+all: dmrg dmrg_table dmrgj1j2 exthubbard trg ctmrg mixedspin hubbard_2d
 
 dmrg: dmrg.o $(ITENSOR_LIBS) $(TENSOR_HEADERS)
 	$(CCCOM) $(CCFLAGS) dmrg.o -o dmrg $(LIBFLAGS)
@@ -61,6 +61,12 @@ trg: trg.o $(ITENSOR_LIBS) $(TENSOR_HEADERS)
 trg-g: mkdebugdir .debug_objs/trg.o $(ITENSOR_GLIBS) $(TENSOR_HEADERS)
 	$(CCCOM) $(CCGFLAGS) .debug_objs/trg.o -o trg-g $(LIBGFLAGS)
 
+ctmrg: ctmrg.o $(ITENSOR_LIBS) $(TENSOR_HEADERS)
+	$(CCCOM) $(CCFLAGS) ctmrg.o -o ctmrg $(LIBFLAGS)
+
+ctmrg-g: mkdebugdir .debug_objs/trg.o $(ITENSOR_GLIBS) $(TENSOR_HEADERS)
+	$(CCCOM) $(CCGFLAGS) .debug_objs/ctmrg.o -o ctmrg-g $(LIBGFLAGS)
+
 hubbard_2d: hubbard_2d.o $(ITENSOR_LIBS) $(TENSOR_HEADERS)
 	$(CCCOM) $(CCFLAGS) hubbard_2d.o -o hubbard_2d $(LIBFLAGS)
 
@@ -73,4 +79,4 @@ mkdebugdir:
 clean:
 	@rm -fr *.o .debug_objs dmrg dmrg-g \
 	dmrg_table dmrg_table-g dmrgj1j2 dmrgj1j2-g exthubbard exthubbard-g \
-    mixedspin mixedspin-g trg trg-g hubbard_2d hubbard_2d-g
+    mixedspin mixedspin-g trg trg-g ctmrg ctmrg-g hubbard_2d hubbard_2d-g
diff --git a/sample/README b/sample/README
index 088187738..ea115b3f2 100644
--- a/sample/README
+++ b/sample/README
@@ -43,3 +43,7 @@ mixedspin - demo code showing how to make and use
 trg - tensor renormalization group (TRG) algorithm
       for computing properties of large 2D classical
       stat mech systems
+
+ctmrg - corner transfer matrix renormalization group (CTMRG) algorithm
+        for computing properties of large 2D classical
+        stat mech systems
diff --git a/sample/ctmrg.cc b/sample/ctmrg.cc
new file mode 100644
index 000000000..a3da45520
--- /dev/null
+++ b/sample/ctmrg.cc
@@ -0,0 +1,60 @@
+#include "src/ctmrg.h"
+#include "src/ising.h"
+#include "itensor/util/print_macro.h"
+
+int
+main()
+  {
+  Real betac = 0.5 * log(sqrt(2) + 1.0);
+  Real beta = 1.1 * betac;
+  int maxdim = 20;
+  int nsteps = 20;
+  
+  auto dim0 = 2;
+  
+  // Define an initial Index making up
+  // the Ising partition function
+  auto s = Index(dim0, "Site");
+  
+  // Define the indices of the scale-0
+  // Boltzmann weight tensor "A"
+  auto sh = addTags(s, "horiz");
+  auto sv = addTags(s, "vert");
+  
+  auto T = ising(sh, sv, beta);
+
+  auto l = Index(1, "Link");
+  auto lh = addTags(l, "horiz");
+  auto lv = addTags(l, "vert");
+  auto Clu0 = ITensor(lv, lh);
+  Clu0.set(1, 1, 1.0);
+  auto Al0 = ITensor(lv, prime(lv), sh);
+  Al0.set(lv = 1, prime(lv) = 1, sh = 1, 1.0);
+
+  auto [Clu, Al] = ctmrg(T, Clu0, Al0, maxdim, nsteps);
+
+  lv = commonIndex(Clu, Al);
+  lh = uniqueIndex(Clu, Al);
+
+  auto Au = replaceInds(Al, {lv, prime(lv), sh},
+                            {lh, prime(lh), sv});
+
+  auto ACl = Al * Clu * dag(prime(Clu));
+
+  auto ACTl = prime(ACl * dag(prime(Au)) * T * Au, -1);
+  auto kappa = elt(ACTl * dag(ACl));
+
+  auto Tsz = ising(sh, sv, beta, true);
+  auto ACTszl = prime(ACl * dag(prime(Au)) * Tsz * Au, -1);
+  auto m = elt(ACTszl * dag(ACl)) / kappa;
+
+  printfln("beta = %.12f", beta);
+  printfln("beta/betac = %.12f", beta / betac);
+  println("maxdim = ", maxdim);
+  println("niters = ", nsteps);
+  printfln("kappa = %.12f", kappa);
+  printfln("m = %.12f", m);
+
+  return 0;
+  }
+
diff --git a/sample/src/ctmrg.h b/sample/src/ctmrg.h
new file mode 100644
index 000000000..67a903e05
--- /dev/null
+++ b/sample/src/ctmrg.h
@@ -0,0 +1,60 @@
+#include "itensor/all.h"
+
+using namespace itensor;
+
+std::tuple<ITensor, ITensor>
+ctmrg(ITensor T,
+      ITensor Clu,
+      ITensor Al,
+      int maxdim,
+      int nsteps,
+      Real cutoff = 0.0)
+  {
+  auto sh = commonIndex(T, Al);
+  auto sv = uniqueIndex(T, {Al, prime(Al)}, "0");
+  auto lv = commonIndex(Clu, Al);
+  auto lh = uniqueIndex(Clu, Al);
+  auto Au = replaceInds(Al, {lv, prime(lv), sh},
+                            {lh, prime(lh), sv});
+  ITensor Uv;
+  for(auto i : range1(nsteps))
+    {
+    // Get the grown corner transfer matrix (CTM)
+    auto Clu_new = Al * Clu * Au * T;
+
+    Clu_new.noPrime();
+    Clu_new.replaceInds({lh, sh}, {prime(lv), prime(sv)});
+
+    // Diagonalize the grown CTM
+    std::tie(Uv, Clu) = diagPosSemiDef(Clu_new,
+                                       {"MaxDim = ", maxdim,
+                                        "Tags = ", tags(lv),
+                                        "Cutoff = ", cutoff});
+
+    lv = commonIndex(Clu, Uv);
+    lh = setTags(lv, tags(lh));
+
+    Clu = toDense(Clu);
+    Clu.replaceInds({prime(lv)}, {lh});
+
+    // The renormalized CTM is the diagonal matrix of eigenvalues
+    // Normalize the CTM
+    auto Cl = Clu * prime(dag(Clu), lh);
+    auto normC = pow(elt(Cl * dag(Cl)), 0.25);
+    Clu /= normC;
+
+    // Calculate the renormalized half row transfer matrix (HRTM)
+    Uv.noPrime();
+    Al = Al * Uv * T * dag(prime(Uv));
+    Al = replaceInds(Al, {prime(sh)}, {sh});
+
+    // Normalize the HRTM
+    auto ACl = Al * Clu * prime(dag(Clu));
+    auto normA = sqrt(elt(ACl * dag(ACl)));
+    Al /= normA;
+    Au = replaceInds(Al, {lv, prime(lv), sh},
+                         {lh, prime(lh), sv});
+    }
+  return std::tuple<ITensor, ITensor>({Clu, Al});
+  }
+
diff --git a/sample/src/ising.h b/sample/src/ising.h
index 710ec21bf..71c56b9a0 100644
--- a/sample/src/ising.h
+++ b/sample/src/ising.h
@@ -3,26 +3,68 @@
 using namespace itensor;
 
 ITensor
-ising(Index const& l, Index const& r,
-      Index const& u, Index const& d,
-      Real beta)
+ising(Index const& sh,
+      Index const& sv,
+      Real beta,
+      bool sz = false)
   {
-  int dim0 = 2;
-  auto A = ITensor(l, r, u, d);
-  auto T = 1/beta;
-
-  // Fill the A tensor with correct Boltzmann weights:
-  auto Sig = [](int s) { return 1.-2.*(s-1); };
-  for(auto sl : range1(dim0))
-  for(auto sd : range1(dim0))
-  for(auto sr : range1(dim0))
-  for(auto su : range1(dim0))
-    {
-    auto E = Sig(sl)*Sig(sd)+Sig(sd)*Sig(sr)
-            +Sig(sr)*Sig(su)+Sig(su)*Sig(sl);
-    auto P = exp(-E/T);
-    A.set(l(sl),r(sr),u(su),d(sd),P);
-    }
-  return A;
+  auto d = dim(sh);
+  auto shp = prime(sh);
+  auto svp = prime(sv);
+  auto T = ITensor(sh, shp, sv, svp);
+  for(auto i : range1(d)){ T.set(i, i, i, i, 1.0); }
+  if(sz) { T.set(1, 1, 1, 1, -1.0); }
+  auto th = sim(sh);
+  auto thp = sim(shp);
+  auto tv = sim(sv);
+  auto tvp = sim(svp);
+  auto Tp = T * delta(sh, th) *
+                delta(shp, thp) *
+                delta(sv, tv) *
+                delta(svp, tvp);
+  // Analytic square root of the bond matrix:
+  // [exp( beta) exp(-beta)
+  //  exp(-beta) exp( beta)]
+  auto lambda_p = sqrt(exp(beta) + exp(-beta));
+  auto lambda_m = sqrt(exp(beta) - exp(-beta));
+  auto x_p = 0.5 * (lambda_p + lambda_m);
+  auto x_m = 0.5 * (lambda_p - lambda_m);
+  auto Xh = ITensor(th, sh);
+  Xh.set(1, 1, x_p);
+  Xh.set(2, 1, x_m);
+  Xh.set(1, 2, x_m);
+  Xh.set(2, 2, x_p);
+  auto Xhp = replaceInds(Xh, {th, sh}, {thp, shp});
+  auto Xv  = replaceInds(Xh, {th, sh}, {tv,  sv });
+  auto Xvp = replaceInds(Xh, {th, sh}, {tvp, svp});
+  return Tp * Xhp * Xvp * Xh * Xv;
   }
 
+//
+// Dual partition function
+//
+
+//ITensor
+//ising(Index const& sh, Index const& sv,
+//      Real beta)
+//  {
+//  int dim0 = 2;
+//  auto T = ITensor(sh, prime(sh), sv, prime(sv));
+//
+//  // Fill the T tensor with correct Boltzmann weights:
+//  auto Sig = [](int s) { return 1. - 2. * (s - 1); };
+//  for(auto ssh  : range1(dim0))
+//  for(auto ssvp : range1(dim0))
+//  for(auto sshp : range1(dim0))
+//  for(auto ssv  : range1(dim0))
+//    {
+//    auto E = Sig(ssh)  * Sig(ssvp) +
+//             Sig(ssvp) * Sig(sshp) +
+//             Sig(sshp) * Sig(ssv)  +
+//             Sig(ssv)  * Sig(ssh);
+//    auto P = exp(-beta * E);
+//    T.set(sh = ssh, prime(sh) = sshp, sv = ssv, prime(sv) = ssvp, P);
+//    }
+//  return T;
+//  }
+
diff --git a/sample/src/trg.h b/sample/src/trg.h
index 1589ab4a7..cc2c14478 100644
--- a/sample/src/trg.h
+++ b/sample/src/trg.h
@@ -4,73 +4,63 @@ using namespace itensor;
 
 std::tuple<ITensor, Real>
 trg(ITensor const& A0,
-    Index l, Index r,
-    Index u, Index d,
-    int maxdim, int topscale)
+    int maxdim, int topscale,
+    Real cutoff = 0.0)
   {
-  auto A = addTags(A0, "scale=0");
-  l.addTags("scale=0");
-  r.addTags("scale=0");
-  u.addTags("scale=0");
-  d.addTags("scale=0");
+  auto A = A0;
+  auto is = findInds(A, "0");
+  auto sh = is(1);
+  auto sv = is(2);
 
   // Keep track of partition function per site, z = Z^(1/N)
   Real z = 1.0;
 
   for(auto scale : range1(topscale))
-      {
-      //printfln("\n---------- Scale %d -> %d  ----------",scale-1,scale);
+    {
+    //printfln("\n---------- Scale %d -> %d  ----------",scale-1,scale);
 
-      // Get the upper-left and lower-right tensors
-      auto [Fl, Fr] = factor(A, {r, d}, {l, u},
-                             {"MaxDim = ", maxdim,
-                              "Tags = ", "left,scale=" + str(scale),
-                              "SVDMethod = ", "gesdd",
-                              "ShowEigs = ", false});
+    // Get the upper-left and lower-right tensors
+    auto [Fh, Fhp] = factor(A, {prime(sh), prime(sv)}, {sh, sv},
+                            {"MaxDim = ", maxdim,
+                             "Tags = ", "horiz",
+                             "SVDMethod = ", "gesdd",
+                             "Cutoff = ", cutoff,
+                             "ShowEigs = ", false});
 
-      // Grab the new left Index
-      auto l_new = commonIndex(Fl, Fr);
+    // Grab the new left Index
+    auto sh_new = commonIndex(Fh, Fhp);
+    Fhp *= delta(sh_new, prime(sh_new));
 
-      // Get the upper-right and lower-left tensors
-      auto [Fu, Fd] = factor(A, {l, d}, {u, r},
-                             {"MaxDim = ", maxdim,
-                              "Tags = ", "up,scale=" + str(scale),
-                              "SVDMethod = ", "gesdd",
-                              "ShowEigs = ", false});
+    // Get the upper-right and lower-left tensors
+    auto [Fv, Fvp] = factor(A, {sh, prime(sv)}, {prime(sh), sv},
+                            {"MaxDim = ", maxdim,
+                             "Tags = ", "vert",
+                             "SVDMethod = ", "gesdd",
+                             "Cutoff = ", cutoff,
+                             "ShowEigs = ", false});
 
-      // Grab the new up Index
-      auto u_new = commonIndex(Fu, Fd);
+    // Grab the new up Index
+    auto sv_new = commonIndex(Fv, Fvp);
+    Fvp *= delta(sv_new, prime(sv_new));
 
-      // Make the new index of Fl distinct
-      // from the new index of Fr by changing
-      // the tag from "left" to "right"
-      auto r_new = replaceTags(l_new, "left", "right");
-      Fr *= delta(l_new, r_new);
+    A = (Fh  * delta(prime(sh), sh)) *
+        (Fv  * delta(prime(sv), sv)) *
+        (Fhp * delta(sh, prime(sh))) *
+        (Fvp * delta(sv, prime(sv)));
+    
+    // Update the indices
+    sh = sh_new;
+    sv = sv_new;
 
-      // Make the new index of Fd distinct
-      // from the new index of Fu by changing the tag
-      // from "up" to "down"
-      auto d_new = replaceTags(u_new, "up", "down");
-      Fd *= delta(u_new, d_new);
+    // Normalize the current tensor and keep track of
+    // the total normalization
+    Real TrA = elt(A * delta(sh, prime(sh)) * delta(sv, prime(sv)));
+    A /= TrA;
+    z *= pow(TrA, 1.0 / pow(2, scale));
 
-      Fl *= delta(r, l);
-      Fu *= delta(d, u);
-      Fr *= delta(l, r);
-      Fd *= delta(u, d);
-      A = Fl * Fu * Fr * Fd;
-      
-      // Update the indices
-      l = l_new;
-      r = r_new;
-      u = u_new;
-      d = d_new;
-
-      // Normalize the current tensor and keep track of
-      // the total normalization
-      Real TrA = elt(A * delta(l, r) * delta(u, d));
-      A /= TrA;
-      z *= pow(TrA, 1.0 / pow(2, 1 + scale));
-      }
+    // If using the dual Ising partition function
+    //z *= pow(TrA, 1.0 / pow(2, 1 + scale));
+    }
 
   //printfln("log(Z)/N = %.12f",log(z));
 
diff --git a/sample/trg.cc b/sample/trg.cc
index 9b42f4a7e..cd4a9e428 100644
--- a/sample/trg.cc
+++ b/sample/trg.cc
@@ -18,13 +18,11 @@ main()
   
   // Define the indices of the scale-0
   // Boltzmann weight tensor "A"
-  auto l = addTags(s, "left");
-  auto r = addTags(s, "right");
-  auto u = addTags(s, "up");
-  auto d = addTags(s, "down");
+  auto sh = addTags(s, "horiz");
+  auto sv = addTags(s, "vert");
   
-  auto A0 = ising(l, r, u, d, beta);
-  auto [A, z] = trg(A0, l, r, u, d, maxdim, topscale);
+  auto A0 = ising(sh, sv, beta);
+  auto [A, z] = trg(A0, maxdim, topscale);
   (void)A;
 
   printfln("beta = %.12f", beta);