diff --git a/examples/Su2ExactDiagonalization/FT_Su2CoarseGrain2To1Site.npy b/examples/Su2ExactDiagonalization/FT_Su2CoarseGrain2To1Site.npy
new file mode 100644
index 000000000..cc16d159d
Binary files /dev/null and b/examples/Su2ExactDiagonalization/FT_Su2CoarseGrain2To1Site.npy differ
diff --git a/examples/Su2ExactDiagonalization/Su2Engine.py b/examples/Su2ExactDiagonalization/Su2Engine.py
new file mode 100644
index 000000000..021481ec3
--- /dev/null
+++ b/examples/Su2ExactDiagonalization/Su2Engine.py
@@ -0,0 +1,263 @@
+# Copyright 2020 The TensorNetwork Authors
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Basic SU(2) representation theory data"""
+
+from math import floor
+from math import factorial
+from math import sqrt
+from collections import defaultdict
+import numpy as np
+from tensornetwork.ncon_interface import ncon
+from Su2Tensor import Su2Tensor
+from Su2Tensor import FusionTree
+
+def myisinteger(
+        num: int) -> bool:
+    """
+    Checks if num is an integer
+    """
+
+    val = 1 if num == floor(num) else 0
+    return val
+
+#######################################################################################
+def compatible(
+        j1: float, j2: float, j12: float) -> bool:
+    """
+    Checks if total spins j1,j2,j12 are compatible with SU(2) fusion rules
+    """
+
+    val = True if j12 in np.arange(abs(j1-j2),j1+j2+1,1).tolist() else False
+    return val
+
+#######################################################################################
+def wigner3j(
+        j1: float, m1: float, j2: float, m2: float, j12: float, m12: float) -> float:
+    """
+    Computes the Wigner 3j symbol. j's at the total spin values and m's are the
+    corresponding spin projection values
+    """
+    # compute Wigner 3j symbol
+    if not myisinteger(2*j1) or not myisinteger(2*j2) or not myisinteger(2*j12) \
+            or not myisinteger(2*m1) or not myisinteger(2*m2) or not myisinteger(2*m12):
+        raise ValueError('Input values must be (semi-)integers!')
+    if not myisinteger(2*(j1-m1)):
+        raise ValueError('j1-m1 must be an integer!')
+    if not myisinteger(2*(j2-m2)):
+        raise ValueError('j2-m2 must be an integer!')
+    if not myisinteger(2*(j12-m12)):
+        raise ValueError('j12-m12 must be an integer!')
+    if not compatible(j1,j2,j12):
+        raise ValueError('The spins j1,j2,j12 are not compatible with the fusion rules!')
+    if abs(m1) > j1:
+        raise ValueError('m1 is out of bounds.')
+    if abs(m2) > j2:
+        raise ValueError('m2 is out of bounds.')
+    if abs(m12) > j12:
+        raise ValueError('m12 is out of bounds.')
+    if m1 + m2 + m12 != 0:
+        return 0    # m's are not compatible. So coefficient must be 0.
+    t1 = j2 - m1 - j12
+    t2 = j1 + m2 - j12
+    t3 = j1 + j2 - j12
+    t4 = j1 - m1
+    t5 = j2 + m2
+
+    tmin = max([0, t1, t2])
+    tmax = min([t3, t4, t5])
+
+    wigner = 0
+    for t in np.arange(tmin, tmax+1, 1).tolist():
+        wigner = wigner + (-1)**t / (factorial(t) * factorial(t - t1) * factorial(t - t2)\
+                  * factorial(t3 - t) * factorial(t4 - t) * factorial(t5 - t))
+
+    wigner = wigner * (-1)**(j1-j2-m12) * sqrt(factorial(j1+j2-j12) * factorial(j1-j2+j12)\
+              * factorial(-j1+j2+j12) / factorial(j1+j2+j12+1) * factorial(j1+m1)\
+              * factorial(j1-m1) * factorial(j2+m2) * factorial(j2-m2) * factorial(j12+m12)\
+              * factorial(j12-m12))
+
+    return wigner
+
+#######################################################################################
+def clebschgordan(
+        j1: float, m1: float, j2: float, m2: float, j12: float, m12: float) -> float:
+    """
+    Computes a Clebsch-Gordan coefficient. j's at the total spin values and m's are the corresponding
+    spin projection values
+    """
+    return (-1)**(j1-j2+m12) * sqrt(2*j12+1) * wigner3j(j1,m1,j2,m2,j12,-m12) # Clebsch-Gordan coefficient from Wigner 3j
+
+#######################################################################################
+def createCGtensor(
+        j1: float, j2: float, j12: float) -> np.array:
+    """
+    Computes a Clebsch-Gordan tensor for the fusion j1xj2 -> j12
+    """
+
+    C = np.zeros((floor(2*j1+1), floor(2*j2+1), floor(2*j12+1)))
+    if not compatible(j1,j2,j12):
+        return C
+    for m1 in np.arange(-j1, j1+1, 1).tolist():
+        for m2 in np.arange(-j2, j2+1, 1).tolist():
+            for m12 in np.arange(-j12, j12+1, 1).tolist():
+                C[floor(j1+m1)][floor(j2+m2)][floor(j12+m12)] = clebschgordan(j1,m1,j2,m2,j12,m12) if m1+m2 == m12 else 0
+    return C
+
+#######################################################################################
+def rsymbol(
+        j1: float, j2: float, j12: float) -> float:
+    """
+    Returns a R-symbol of SU(2)
+    """
+
+    C1 = createCGtensor(j1,j2,j12)
+    C2 = createCGtensor(j2,j1,j12)
+    C = ncon([C1, C2],[[1, 2, -1],[2, 1, -2]])
+    return C[0][0]
+
+#######################################################################################
+def racahDenominator(
+        t: float, j1: float, j2: float, j3: float, J1: float, J2: float, J3: float) -> float:
+    """
+    Helper function used in wigner6j().
+    """
+
+    return factorial(t-j1-j2-j3) * factorial(t-j1-J2-J3) * factorial(t-J1-j2-J3) * factorial(t-J1-J2-j3) \
+    * factorial(j1+j2+J1+J2-t) * factorial(j2+j3+J2+J3-t) * factorial(j3+j1+J3+J1-t)
+
+#######################################################################################
+def triangleFactor(
+        a: int, b: int, c: int) -> float:
+    """
+    Helper function used in wigner6j().
+    """
+
+    return factorial(a+b-c) * factorial(a-b+c) * factorial(-a+b+c)/factorial(a+b+c+1)
+
+#######################################################################################
+def wigner6j(
+        j1: float, j2: float, j3: float, J1: float, J2: float, J3: float) -> float:
+    """
+    Returns a Wigner 6j-symbol of SU(2).
+    """
+
+    tri1 = triangleFactor(j1, j2, j3)
+    tri2 = triangleFactor(j1, J2, J3)
+    tri3 = triangleFactor(J1, j2, J3)
+    tri4 = triangleFactor(J1, J2, j3)
+
+    # Calculate the summation range in Racah formula
+    tmin = max([j1+j2+j3, j1+J2+J3, J1+j2+J3, J1+J2+j3])
+    tmax = min([j1+j2+J1+J2, j2+j3+J2+J3, j3+j1+J3+J1])
+
+    Wigner = 0
+
+    for t in np.arange(tmin, tmax+1, 1).tolist():
+        Wigner = Wigner + sqrt(tri1*tri2*tri3*tri4) * ((-1)**t) * factorial(t+1)/racahDenominator(t,j1,j2,j3,J1,J2,J3)
+
+    return Wigner
+
+#######################################################################################
+def fsymbol(
+        j1: float, j2: float, j3: float, j12: float, j23: float, j: float) -> float:
+    """
+    Returns he recouping coefficient of SU(2), related to the Wigner 6j symbol by an overall factor.
+    """
+
+    if not compatible(j1,j2,j12) or not compatible(j2,j3,j23) or not compatible(j1,j23,j) or not compatible(j12,j3,j):
+        return 0
+
+    return wigner6j(j1,j2,j12,j3,j,j23) * sqrt(2*j12+1) * sqrt(2*j23+1) * (-1)**(j1+j2+j3+j)
+
+#######################################################################################
+def isvalidindex(
+        index: dict) -> None:
+    """
+    Checks for a valid SU(2) index, and raises an exception if its not.
+    index: a Dictionary of the type {total spin value: degeneracy dimension}.
+    """
+
+    if not isinstance(index, dict):
+        raise ValueError('SU(2) index should be a dictionary of the type: {total spin, degeneracy dimension}')
+
+    for key in index:
+        if not myisinteger(2*key):
+            raise ValueError('Spin values in an index should be (semi-)integers.')
+        if not myisinteger(index[key]):
+            raise ValueError('Degeneracy values in an index should be integers.')
+
+#########################################################################################
+def fuseindices(
+        ia: dict, ib: dict) -> dict:
+    """
+    Fuses two SU(2) indices ia, ib according to the fusion rules.
+    ia, ib: Dictionaries of the type {total spin value: degeneracy dimension}.
+    """
+
+    if len(ia) == 0 or len(ib) == 0:
+        raise ValueError('Input indices cannot be empty.')
+    isvalidindex(ia)
+    isvalidindex(ib)
+
+    iab = defaultdict(int) # initialize fused index to empty default dictionary
+    for ja in ia:
+        da = ia[ja]
+        for jb in ib:
+            db = ib[jb]
+            for jjab in np.arange(abs(ja-jb), ja+jb+1, 1).tolist():
+                if len(iab) == 0:
+                    iab[jjab] = da*db
+                else:
+                    iab[jjab] += da*db
+
+    return iab
+
+########################################################################################
+def fuse(
+        ia: dict, ib: dict, totalSectors: list=[]) -> Su2Tensor:
+    """
+    Creates a 3-index tensor X that two SU(2) indices according to the fusion rules.
+    ia, ib: Dictionaries of the type {total spin value: degeneracy dimension}.
+    totalSectors (optional): list of total sectors to keep. [total spin values]
+    """
+    iab = fuseindices(ia, ib)
+
+    # remove some total sectors if targetSectors are specified
+    if totalSectors:
+        iabCopy = dict(iab) # shallow copy should suffice
+        for jab in iab:
+            if not jab in totalSectors:
+                iabCopy.pop(jab) # remove unwanted total sectors
+        iab = iabCopy # update iab
+
+    dabGone = defaultdict(int) # dictionary index by jab to keep track to degeneracies as they are filled in X
+
+    X = Su2Tensor([iab, ia, ib])
+    for ja in ia:
+        da = ia[ja]
+        for jb in ib:
+            db = ib[jb]
+            for jab in iab:
+                if not compatible(ja,jb,jab):
+                    continue
+                dab = iab[jab]
+                temp = np.zeros((da*db,dab))
+                temp[0:db*da, dabGone[jab]:db*da+dabGone[jab]] = np.identity(db*da)
+                temp = temp.reshape((db, da, dab), order='F')  # Specifying order=F means a column-major reshape, similar to MATLAB. Default in python is row-major
+                temp = temp.transpose(2,1,0)
+                X.addblock(FusionTree((jab, ja, jb)),temp)
+                dabGone[jab] += da*db
+
+    return X
\ No newline at end of file
diff --git a/examples/Su2ExactDiagonalization/Su2ExactDiagonalization.py b/examples/Su2ExactDiagonalization/Su2ExactDiagonalization.py
new file mode 100644
index 000000000..fbcf9f9e6
--- /dev/null
+++ b/examples/Su2ExactDiagonalization/Su2ExactDiagonalization.py
@@ -0,0 +1,277 @@
+# Copyright 2020 The TensorNetwork Authors
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+""" Performs exact diagonalization of SU(2)-symmetric Hamiltonians by exploiting the symmetry.
+    - Diagonalization can be restricted to specified total spin sectors.
+    - Uses pre-computation of contraction factors resulting from contracting Clebsch-Gordan tensors;
+      factors are computed on first run of the code and stored in a file.
+    - Provides a preset of popular SU(2)-symmetric Hamiltonian densities
+    - Assumptions: 1) Hamiltonian density must be a 2-site operator. 2) Presently only implements open
+                   boundary condition. 3) Only uses ncon() facility of the TensorNetwork library."""
+
+import numpy as np
+from numpy import linalg as LA
+from math import floor
+from scipy.sparse.linalg import eigs
+
+from tensornetwork.ncon_interface import ncon
+
+from Su2Tensor import Su2Tensor
+from Su2Tensor import FusionTree
+import Su2Engine
+
+# globals
+total, left, right = (0, 1, 2)  # some constants (index values)
+
+########################################################################################################################
+def Su2ExactDiagonalization(
+        N: 'int', d: 'dict', h: Su2Tensor, targetSectors: list=[],
+        tableFileName: 'string' = 'FT_Su2CoarseGrain2To1Site.npy',
+        doCreateTable: 'bool' =False) -> tuple[dict, dict, np.array]:
+    """
+    Performs exact diagonalization in specified targetSectors. Returns the low-energy states and their energies in
+    these sectors.
+    N: number of lattice sites.
+    d: SU(2) index describing each lattice site. Must be a dictionary of the type {total spin: degeneracy dimension}
+    h: 2-site Hamiltonian density. Must be an object of SU2Tensor class
+    targetSectors: Specifies which total spin sectors (of the entire lattice) to diagonalize and how many states
+                   required in each sector. Must be dictionary of the type: {total spin, number of states}
+                   If None diagonalizes all the sectors
+    doCreateTable: set to True if the factor tables required by this function have not been created yet.
+    tableFileName: file to store the table when table creation is demanded.
+    """
+
+    # check validity of inputs
+    if not Su2Engine.myisinteger(N):
+        raise ValueError('Lattice size must be a natural number!')
+    if N == 0 or N == 1:
+        raise ValueError('Lattice size must be greater than equal to 2!')
+    Su2Engine.isvalidindex(d)
+    if not isinstance(h,Su2Tensor):
+        raise ValueError('Hamiltonian density must be a SU(2) tensor.')
+
+    # create blocking isometries from bottom (physical index) to top(total index)
+    w = [] # list of blocking isometries
+    bond = d
+    for k in range(N-2):
+        X = Su2Engine.fuse(bond, d)
+        w.append(X)
+        bond = X.getIndices()[0] # set bond to new total index
+
+    # Create the final blocking isometry specially
+    if targetSectors: # the user has asked to truncate final index
+        X = Su2Engine.fuse(bond, d, list(targetSectors.keys())) # truncate total spin sectors
+    else:
+        X = Su2Engine.fuse(bond, d) # otherwise the full blocking tensor
+    w.append(X)
+    bond = X.getIndices()[0] # bond is now the total spin basis on the entire lattice
+
+    if doCreateTable:
+        factorTable = createFactorTable_Su2CoarseGrain_2To1Site([[w[N-4],w[N-3]], [w[N-3],w[N-2]]],w[0], tableFileName)
+    else:
+        factorTable = np.load(tableFileName, allow_pickle='TRUE').item()
+
+    su2H = Su2Tensor([bond, bond]) # Total Hamiltonian.
+    wd = w[0]
+    for term in range(N-1):
+        # first coarse-grain local term to a 1-site operator. Fist term has to
+        # be treated differently from the others.
+
+        if term == 0:
+            hh = h # already blocked to a one-site operator
+        else:
+            hh = Su2CoarseGrain_2To1site(w[term-1], w[term], wd, h, factorTable)
+
+        # now coarse-grain this one-site operator through the remaining linear fusion tree
+        for m in range(term+1,N-1):
+            hh = Su2CoarseGrain_1To1site(w[m],hh)
+
+        # hh is now the one-site version of the local hamiltonian term. Add to current sum
+        hhBlocks = hh.getAllBlocks()
+        for fusionTree in hhBlocks:
+            su2H.addblock(fusionTree,hhBlocks[fusionTree])
+
+    # Now let's diagonalize
+    spectrum = dict() # dictionary of the type: {total spin: spectra}
+    eigenvectors = dict() # dictionary of the type: {total spin: eigenvectors}
+    fullspectrum = None # full spectrum (including degneracies) as a NumPy array
+    su2HBlocks = su2H.getAllBlocks()
+    for fusionTree in su2HBlocks:
+        degdim = 2*fusionTree[0]+1
+        if targetSectors: # diagonalize only the sectors that the user wants
+            temp, ev = eigs(su2HBlocks[fusionTree], k=targetSectors[fusionTree[0]], which='SR')
+        else: # do full digonalization of all the blocks
+            temp = LA.eigvals(su2HBlocks[fusionTree]) # full diagonalization
+            ev = None
+        spectrum[(fusionTree[0])] = temp
+        eigenvectors[(fusionTree[0])] = ev
+        temp = np.kron(temp,np.ones((1,floor(degdim))))
+        if fullspectrum is None:
+            fullspectrum = temp
+        else:
+            fullspectrum = np.append(fullspectrum, temp)
+    fullspectrum.sort()
+
+    return spectrum, eigenvectors, fullspectrum
+
+########################################################################################################################
+def Su2CoarseGrain_2To1site(
+        w1: Su2Tensor, w2: Su2Tensor, wd: Su2Tensor, h: Su2Tensor, factorTable: dict) -> Su2Tensor:
+    """
+    Coarse-grains 2-site Ham h to a 1-site term through isometries w1,w2 inside a linear fusion tree
+    w2 is the blocking isometry from the higher level. h,hL,w1,w2 are SU(2)-symmetric tensors.
+    """
+
+    hh = Su2Tensor(h.getIndices())  # the coarse-grained hamiltonian. Even though h is a 4-index tensor it is
+                                    # stored as a 2-index matrix. So indices of h are the blocked 2-site indices
+    w2Blocks = w2.getAllBlocks()
+    w1Blocks = w1.getAllBlocks()
+    wdBlocks = wd.getAllBlocks()
+    hBlocks = h.getAllBlocks()
+    for ft in factorTable: # iterate through all fusion trees in table
+        c1 = FusionTree((ft[0], ft[1], ft[4]))
+        c2 = FusionTree((ft[1], ft[2], ft[3]))
+        c3 = FusionTree((ft[5], ft[3], ft[4]))
+        c4 = FusionTree((ft[5], ft[6], ft[7]))
+        c5 = FusionTree((ft[8], ft[2], ft[6]))
+        c6 = FusionTree((ft[0], ft[8], ft[7]))
+        ch = FusionTree((ft[5], ft[5]))
+
+        if not c1 in w2Blocks or not c2 in w1Blocks or not c3 in wdBlocks or not c4 in wdBlocks or not c5 in w1Blocks \
+                or not c6 in w2Blocks or not ch in hBlocks:
+            continue
+        temp = ncon([w2Blocks[c1],w1Blocks[c2],wdBlocks[c3],hBlocks[ch],wdBlocks[c4],w1Blocks[c5],w2Blocks[c6]],\
+                    [[-1,2,4],[2,8,3],[1,3,4],[1,9],[9,6,7],[5,8,6],[-2,5,7]])
+        hh.addblock(FusionTree((ft[0],ft[0])),factorTable[ft]*temp)
+
+    return hh
+
+########################################################################################################################
+def Su2CoarseGrain_1To1site(
+    w: Su2Tensor, h: Su2Tensor) -> Su2Tensor:
+    """
+    Coarse-grains (blocks) 1-site operator h through isometry w.
+    w: the blocking isometry
+    h: the 1-site SU(2)-symmetric operator
+    """
+
+    hh = Su2Tensor(h.getIndices()) # the coarse-grained hamiltonian. Even though h is a 4-index tensor it is
+                                   # stored as a 2-index matrix. So indices of h are the blocked 2-site indices
+    wBlocks = w.getAllBlocks() # a dictionary {FusionTree:block}
+    hBlocks = h.getAllBlocks() # a dictionary {FusionTree:block}
+    for ft in wBlocks:
+        ch = FusionTree((ft[1],ft[1]))
+        if not ch in hBlocks:
+            continue
+        temp = ncon([wBlocks[ft],hBlocks[ch],wBlocks[ft]],[[-1, 1, 3],[1, 2],[-2, 2, 3]])
+        hh.addblock(FusionTree((ft[0],ft[0])),temp)
+
+    return hh
+
+########################################################################################################################
+def createFactorTable_Su2CoarseGrain_2To1Site(
+        W: 'list of pairs of Su2Tensors', wd: Su2Tensor, tableFileName: str) -> dict:
+    """
+    (Precomputation.) Creates the factors and compatible chargesectors for the function Su2CoarseGrain_2To1Site()
+    W: a list of pairs of SU(2) tensors
+    wd: a SU(2) tensor (the isometry that blocks 2 lattice sites)
+    tableFileName: file name to store the created factor table
+    """
+
+    FactorTable = dict() # Dictionary from fusion trees to factors
+    for k in range(len(W)):
+        w1 = W[k][0]
+        w2 = W[k][1]
+        for c1 in w2.getAllBlocks():
+            for c2 in w1.getAllBlocks():
+                if c2[total] != c1[left]:
+                    continue
+                for c3 in wd.getAllBlocks():
+                    if c2[right] != c3[left] or c1[right] != c3[right] or not Su2Engine.compatible(c2[left], c3[total],\
+                                                                                                             c1[total]):
+                        continue
+                    for c4 in wd.getAllBlocks():
+                        if c4[total] != c3[total]:
+                            continue
+                        for c5 in w1.getAllBlocks():
+                            if c5[right] != c4[left] or c5[left] != c2[left] or not Su2Engine.compatible(c5[total],\
+                                                                                                c4[right],c1[total]):
+                                continue
+                            for c6 in w2.getAllBlocks():
+                                if c6[left] != c5[total] or c6[right] != c4[right] or c6[total] != c1[total]:
+                                    continue
+
+                                C1 = Su2Engine.createCGtensor(c1[left], c1[right], c1[total])
+                                C2 = Su2Engine.createCGtensor(c2[left], c2[right], c2[total])
+                                C3 = Su2Engine.createCGtensor(c3[left], c3[right], c3[total])
+                                C4 = Su2Engine.createCGtensor(c4[left], c4[right], c4[total])
+                                C5 = Su2Engine.createCGtensor(c5[left], c5[right], c5[total])
+                                C6 = Su2Engine.createCGtensor(c6[left], c6[right], c6[total])
+                                Cmat = ncon([C1,C2,C3,C4,C5,C6],[[1,3,-1],[7,2,1],[2,3,8],[5,6,8],[7,5,4],[4,6,-2]])
+                                FactorTable[FusionTree((c1[total],c1[left],c2[left],c2[right],c1[right],c3[total],\
+                                                                           c4[left],c4[right],c5[total]))] = Cmat[0][0]
+
+    np.save(tableFileName, FactorTable)
+    return FactorTable
+
+########################################################################################################################
+def Su2Create2SiteHamiltonianDensity(
+        whichHam: str) -> Su2Tensor:
+    """
+    Creates 2-site Hamiltonian density as a SU(2) tensor. d is the SU(2) index for each lattice site.
+    whichHam is a string to select from a preset of SU(2)-symmetry Hamiltonian densities:
+    'afmHeisenberg', 'spin1Heisenberg', 'aklt'"""
+    if whichHam == 'afmHeisenberg':
+        d2 = {0:1,1:1}
+        h = Su2Tensor([d2, d2])
+        block0 = np.ones((1,1))
+        block0[0][0] = -3
+        block1 = np.ones((1,1))
+        h.addblock(FusionTree((0, 0)), block0)
+        h.addblock(FusionTree((1, 1)), block1)
+    elif whichHam == 'identity': # for testing purposes
+        d2 = {0: 1, 1: 1}
+        h = Su2Tensor([d2, d2])
+        block0 = np.ones((1, 1))
+        block1 = np.ones((1, 1))
+        h.addblock(FusionTree((0, 0)), block0)
+        h.addblock(FusionTree((1, 1)), block1)
+    elif whichHam == 'spin1Heisenberg':
+        d2 = {0: 1, 1: 1, 2: 1}
+        h = Su2Tensor([d2, d2])
+        block0 = np.ones((1, 1))
+        block0[0][0] = -2
+        block1 = np.ones((1, 1))
+        block1[0][0] = -1
+        block2 = np.ones((1, 1))
+        block2[0][0] = 1
+        h.addblock(FusionTree((0, 0)), block0)
+        h.addblock(FusionTree((1, 1)), block1)
+        h.addblock(FusionTree((2, 2)), block2)
+    elif whichHam == 'aklt':
+        d2 = {0: 1, 1: 1, 2: 1}
+        h = Su2Tensor([d2, d2])
+        block0 = np.ones((1, 1))
+        block0[0][0] = -2/3
+        block1 = np.ones((1, 1))
+        block1[0][0] = -2/3
+        block2 = np.ones((1, 1))
+        block2[0][0] = 4/3
+        h.addblock(FusionTree((0, 0)), block0)
+        h.addblock(FusionTree((1, 1)), block1)
+        h.addblock(FusionTree((2, 2)), block2)
+    elif whichHam == 'blbqHeisenberg':
+        pass
+
+    return h
diff --git a/examples/Su2ExactDiagonalization/Su2Tensor.py b/examples/Su2ExactDiagonalization/Su2Tensor.py
new file mode 100644
index 000000000..f38ad62bb
--- /dev/null
+++ b/examples/Su2ExactDiagonalization/Su2Tensor.py
@@ -0,0 +1,90 @@
+# Copyright 2020 The TensorNetwork Authors
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Definition of classes FusionTree and Su2Tensor. (Very basic at present, only handles 2- and 3- index tensors)"""
+import numpy as np
+
+class FusionTree:
+    """
+    Class to encapsulate a FusionTree. Currently used only for 2-index and 3-index Su2Tensors. In these cases FusionTree
+    is essentially a wrapper on a list of charges (2 or 3 charges)
+    """
+    def __init__(self, charges):
+        self.charges = charges # charges is simply a tuple for now. (Should be replaced by a decorated tree in future)
+    def __hash__(self):
+        """
+        To use objects of FusionTree as a key in a dictionary
+        """
+        return hash(self.charges)
+
+    def __eq__(self, other):
+        """
+        To use objects of FusionTree as a key in a dictionary.
+        """
+        return self.charges == other.charges
+
+    def __ne__(self, other):
+        """
+        To use objects of FusionTree as a key in a dictionary. Not strictly necessary, but to avoid having
+        both x==y and x!=y True at the same time.
+        """
+        return not(self == other)
+
+    def __getitem__(self, key):
+        """
+        Override the built-in [] operator, so we can access elements of charges directly using the object
+        Currently, the class FusionTree is a wrapper to a list charges.
+        """
+        return self.charges[key]
+
+class Su2Tensor:
+    """
+    Only for 2- and 3-index SU(2) tensors for now! Currently essentially a wrapper for a dict: {FusionTree,block}
+    """
+    def __init__(
+            self, indices: list):
+        self.blocks = dict() # blocks is a list indexed by keys which are FusionTrees
+        self.indices = indices
+
+    def addblock(
+            self, fusiontree: FusionTree, block: np.array):
+        """
+        Adds a block indexed by the given fusiontree.
+        """
+        if fusiontree in self.blocks:
+            self.blocks[fusiontree] += block
+        else:
+            self.blocks[fusiontree] = block
+
+    def getblock(
+            self, fusiontree: FusionTree):
+        """
+        Retrieves a block indexed by the given fusiontree.
+        """
+        if fusiontree in self.blocks:
+            return self.blocks[fusiontree]
+        else:
+            raise KeyError('Fusion tree not found in the tensor.')
+
+    def getAllBlocks(self):
+        """
+        Retrieves the whole block dictionary.
+        """
+        return self.blocks
+
+    def getIndices(self):
+        """
+        Retrieves all the indices of the tensor.
+        """
+        return self.indices
diff --git a/examples/Su2ExactDiagonalization/main.py b/examples/Su2ExactDiagonalization/main.py
new file mode 100644
index 000000000..b3bd7bfb9
--- /dev/null
+++ b/examples/Su2ExactDiagonalization/main.py
@@ -0,0 +1,113 @@
+"""
+Test cases for SU(2) Exact Diagonalization
+"""
+import numpy as np
+from tensornetwork.ncon_interface import ncon
+import Su2Engine
+from Su2Tensor import FusionTree
+import Su2ExactDiagonalization as SD
+from scipy.sparse.linalg import eigs
+
+########################################################################################################################
+if __name__ == '__main__':
+
+    ###### Test1 : diagonalize spin-1 Heisenberg model on a lattice of N=4 sites #######################################
+    print('Test1: spin-1 Heisenberg model on a lattice of N=4 sites')
+    N = 4 # number of lattice sites
+    d = {1:1} # SU(2) charge on each lattice site = spin 1
+    h = SD.Su2Create2SiteHamiltonianDensity('spin1Heisenberg')
+
+    # full diagonalization. When running Su2ExactDiagonalization() for the first time for given N,d set doCreateTable
+    # to True in order to generate the factor tables for use in subsequent calls (even for different ham, see Test2).
+    # The table for N also work for any small M < N.
+    spec = SD.Su2ExactDiagonalization(N,d,h,tableFileName = 'spinOneFactorTable.npy',doCreateTable = True)
+    SD.printSpectrum(spec[0]) # prints sector-wise spectrum, but degeneracies are suppressed
+    print('===========================================================')
+    ###### End of Test1 : diagonalize spin-1 Heisenberg model on a lattice of N=4 sites ################################
+
+    ###### Test2 : diagonalize spin-1 AKLT model on a lattice of N=4 sites #############################################
+    print('Test2: spin-1 AKLT model on a lattice of N=4 sites')
+    N = 4 # number of lattice sites
+    d = {1:1} # SU(2) charge on each lattice site = spin 1
+    h = SD.Su2Create2SiteHamiltonianDensity('aklt')
+
+    # full diagonalization. Reuse factor table from Test1
+    spec = SD.Su2ExactDiagonalization(N,d,h,tableFileName = 'spinOneFactorTable.npy',doCreateTable = False)
+    SD.printSpectrum(spec[0])  # prints sector-wise spectrum, but degeneracies are suppressed
+
+    # diagonalize to find 2 lowest states in total spin = 0, 3 states in total spin = 1, 2 states in total spin 2
+    spec = SD.Su2ExactDiagonalization(N, d, h, targetSectors = {0:2, 1:3, 2:2}, tableFileName='spinOneFactorTable.npy',
+                                                        doCreateTable=False) # reuse factor tables from previous runs
+    SD.printSpectrum(spec[0])  # prints sector-wise spectrum, but degeneracies are suppressed
+    print('===========================================================')
+    ###### End of Test2 : diagonalize spin-1 AKLT model on a lattice of N=4 sites ######################################
+
+    ###### Test3 : diagonalize spin-1/2 AFM Heisenberg model on a lattice of N=8 sites #################################
+    print('Test2: spin-1/2 AFM Heisenberg model on a lattice of N=8 sites')
+    N = 8  # number of lattice sites
+    d = {0.5: 1}  # SU(2) charge on each lattice site = spin 1/2
+    h = SD.Su2Create2SiteHamiltonianDensity('afmHeisenberg')
+
+    # full diagonalization. Create factor tables for spin 1/2 chain with N=8 sites and store in specified file.
+    spec = SD.Su2ExactDiagonalization(N, d, h, tableFileName='spinHalfFactorTable.npy', doCreateTable=True)
+    SD.printSpectrum(spec[0])  # prints sector-wise spectrum, but degeneracies are suppressed
+    print('===========================================================')
+    ###### End of Test3 : diagonalize spin-1 AKLT model on a lattice of N=4 sites ######################################
+
+    ###### Test4 : Calculate Wigner and Clebsch-Gordan coefficients of SU(2) ###########################################
+    print('Demonstration of functions to calculate Wigner and Clebsch-Gordan coefficients of SU(2)')
+    w3j = Su2Engine.wigner3j(0.5, 0.5, 0.5, -0.5, 0, 0)
+    cg = Su2Engine.clebschgordan(0.5, 0.5, 0.5, -0.5, 0, 0)
+    cg1 = Su2Engine.clebschgordan(1, -1, 1, 0, 1, -1)
+    cg2 = Su2Engine.clebschgordan(1, -1, 1, 0, 2, -2)
+    C = Su2Engine.createCGtensor(0.5, 0.5, 1)
+    factor = Su2Engine.fsymbol(1, 1, 1, 1, 1, 1)
+    factor1 = Su2Engine.fsymbol(1, 2, 0.5, 3, 1.5, 2.5)
+    print('===========================================================')
+    ###### End of Test4 : Calculate Wigner and Clebsch-Gordan coefficients of SU(2) ####################################
+
+    ###### Test5 : Check fusion rules of SU(2) #########################################################################
+    print('Demonstration of functions to check fusion rules of SU(2)')
+    res = Su2Engine.compatible(j1=0.5, j2=0.5, j12=0.5)
+    res = Su2Engine.compatible(j1=1, j2=1, j12=1)
+    res = Su2Engine.compatible(j1=1, j2=0.5, j12=1)
+    print('===========================================================')
+    ###### End of Test5 : Check fusion rules of SU(2) ##################################################################
+
+    ###### Test5 : Fusing two SU(2) indices ############################################################################
+    print('Demonstration of fusing two SU(2) indices')
+    ia = {0.5:1}
+    ib = {0.5:1}
+    ic = {0.5:3, 2.5:1.3}
+
+    try:
+        Su2Engine.isvalidindex(ia)
+    except:
+        print('Invalid index')
+
+    try:
+        Su2Engine.isvalidindex(ic)
+    except:
+        print('Invalid index')
+
+    iab = Su2Engine.fuseindices(ia,ib) # fuses indices ia and ib. index iab has the fused spin and degeneracy values
+    iiab = Su2Engine.fuseindices(iab,iab) # another example
+
+    ia = {0:1, 1:1}
+    ib = {0:1, 1:1}
+    X = Su2Engine.fuse(ia,ib,[0,1]) # X is the fusion tensor. Its a Su2Tensor that can be contracted with another
+                                    # Su2Tensor inorder to fuse two indices and reshape that tensor
+    block000 = X.getblock(FusionTree((0, 0, 0))) # should be some assignments of 0's and 1's to give an injective map
+                                                 # from input indices to fused index
+    block011 = X.getblock(FusionTree((0, 1, 1)))
+    block101 = X.getblock(FusionTree((1, 0, 1)))
+    block110 = X.getblock(FusionTree((1, 1, 0)))
+    block111 = X.getblock(FusionTree((1, 1, 1)))
+
+    ia = {0.5:1, 1:2, 3:4}
+    ib = {0.5:1, 5:5, 2:1}
+    X = Su2Engine.fuse(ia,ib)
+    block0 = X.getblock(FusionTree((0, 0.5, 0.5)))
+    block1 = X.getblock(FusionTree((1, 1, 2)))
+    print('===========================================================')
+    ###### End of Test5 : Fusing two SU(2) indices #####################################################################
diff --git a/examples/Su2ExactDiagonalization/spinOneFactorTable.npy b/examples/Su2ExactDiagonalization/spinOneFactorTable.npy
new file mode 100644
index 000000000..cc16d159d
Binary files /dev/null and b/examples/Su2ExactDiagonalization/spinOneFactorTable.npy differ