Tensor Network Theory Library  Beta release 1.2.1 A library of routines for performing TNT-based operations
tntMpoPmpo.c
1 /*
2  Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch
3  Date: $LastChangedDate$
4  (c) University of Oxford 2014
5  */
6
11 /* Include the internal header for the TNT library */
12 #include "tntMpsInternal.h"
13
14
35 tntComplex tntMpoPmpoTrace(tntNetwork mpo,
36  tntNodeArray *op,
39  tntIntArray *sitenum)
41 {
42
43  tntNetwork mpoc; /* The full MPO-PMPO network */
44  tntComplex trace; /* The value of the trace to return */
45
46  /* the full network will be altered then deleted during the network contraction, so use a copy so that the original mpo is untouched. */
47  mpoc = tntNetworkCopy(mpo);
48
49  /* Take the product of the MPS with the PMPO */
50  tntMpoPmpoProduct(mpoc,op,sitenum);
51
52  /* Find the trace */
53  trace = tntMpoTrace(mpoc);
54
55  /* Free the network representing the product */
56  tntNetworkFree(&mpoc);
57
58  /* Return the trace */
59  return trace;
60 }
61
62
84 void tntMpoPmpoProduct(tntNetwork mpo,
85  tntNodeArray *op,
86  tntIntArray *sitenum)
87 {
88  /* Networks and nodes during the contraction */
89  tntNode M; /* The current MPO node */
90  tntNode opn; /* A copy of the operator to use for contracting */
91
92  /* Counters and site information used during the contraction */
93  int sitecount; /* the current site count */
94  int L; /* The number of MPS nodes in the network */
95  unsigned opnum=0; /* Counter for the operators */
96  tntIntArray siterank; /* Rank for each site-number, to make then in ascending order */
97  tntIntArray sitenum_cp; /* Copy of the site numbers which will be used for sorting */
98  unsigned i, j; /* Used for looping */
99  int sitenummax, sitemax; /* Position of current maximum site, current maximum site */
100
101  /* Determine the length of the MPS. */
102  L = (int) tntMpoLength(mpo);
103
104  /* Check that operators if they have been provided */
105  if ((NULL == op) || (NULL == sitenum) || (0 == op->sz) || (0 == sitenum->sz)) {
106  tntWarningPrint("Leaving function tntMpoPmpoProduct() with no action taken as no operators or sitenumbers have been provided"); /* NO_COVERAGE */
107  return; /* NO_COVERAGE */
108  } /* NO_COVERAGE */
109
110  /* Allocate arrays for copies of the site numbers, and rank of site numbers */
111  sitenum_cp = tntIntArrayAlloc(op->sz);
112  siterank = tntIntArrayAlloc(op->sz);
113
114  /* Check that the sites and orthogonality centre given are sensible, and at the same time copy over site numbers */
115  for (opnum = 0; opnum < op->sz; opnum++) {
116  if (sitenum->vals[opnum] >= L) {
117  tntErrorPrint("The site for the operator %d is %d, which is greater than the length of the MPS (%d).", opnum, sitenum->vals[opnum], L); /* NO_COVERAGE */
118  } /* NO_COVERAGE */
119
120  if (sitenum->vals[opnum] < 0) {
121  tntErrorPrint("The site for the operator %d is %d, which is not allowed - site numbers should be zero or more", opnum, sitenum->vals[opnum]); /* NO_COVERAGE */
122  } /* NO_COVERAGE */
123
124  sitenum_cp.vals[opnum] = sitenum->vals[opnum];
125  }
126
127  /* First give each element of sitenum a rank, where siterank i gives the element of ith biggest site */
128  /* Loop through all legs from end to beginning. For each index bl, find the blth largest leg dimension */
129  for (i = opnum; i-->0; ) {
130  /* Reset marker for largest site number */
131  sitemax = -1;
132  for (j = opnum; j-- > 0; ) {
133  if (sitenum_cp.vals[j] > sitemax) {
134  sitemax = sitenum_cp.vals[j];
135  sitenummax = j;
136  }
137  }
138  /* Marks current largest site as being in jth position in array, and set its dimension to 0 to indicate that it no longer requires sorting */
139  sitenum_cp.vals[sitenummax] = -1;
140  siterank.vals[i] = sitenummax;
141  }
142
143  /* Copy of site numbers no longer required */
144  tntIntArrayFree(&sitenum_cp);
145
146  /* Reset counter for number of operators */
147  opnum = 0;
148
149  /* Move through the network from left to right. On each site, check whether there is an operator to apply, and if contract it with the physical legs */
150
151  /* identify the first node */
152  M = tntNodeFindFirst(mpo);
153
154  /* Go through remaining sites, contracting them with alpha one by one */
155  for (sitecount = 0; sitecount < L; sitecount++) {
156
157  /* If there are operators on this site, join to the physical leg and contract with A */
158  while ((opnum < op->sz) && (sitenum->vals[siterank.vals[opnum]] == sitecount)) {
159  /* Make a copy of the operator */
160  opn = tntNodeCopy(op->vals[siterank.vals[opnum]]);
161
162  /* join the operator to the physical leg */
163  tntNodeJoin(opn, "U", M, "D");
164
165  /* contract the first node with the operator on the physical leg */
166  M = tntNodeContract(M,opn);
167
168  /* increment the counter for the operator number */
169  opnum++;
170  }
171
172  /* Break if all the operators have been applied */
173  if (opnum == op->sz) {
174  break;
175  }
176
177  /* Find the next node along in the MPS */
178  M = tntNodeFindConn(M, "R");
179  }
180
181  /* Free the site rank */
182  tntIntArrayFree(&siterank);
183
184  return;
185
186 }
187
void tntMpoPmpoProduct(tntNetwork mpo, tntNodeArray *op, tntIntArray *sitenum)
Definition: tntMpoPmpo.c:84
tntNode tntNodeCopy(tntNode A)
Definition: tntNodeUtil.c:304
void tntNetworkFree(tntNetwork *nwp)
Definition: tntNetwork.c:118
void tntNodeJoin(tntNode A, tntLegLabel legA, tntNode B, tntLegLabel legB)
Definition: tntNodeConn.c:52
tntComplex tntMpoTrace(tntNetwork mpo)
Definition: tntMpoTrace.c:20
Definition: tnt.h:65
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray.c:620
tntIntArray tntIntArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntArray.c:29
tntNode tntNodeFindConn(tntNode A, tntLegLabel legA)
Definition: tntNodeInfo.c:29
tntNode tntNodeContract(tntNode A, tntNode B, tntLegLabel legMapAC, tntLegLabel legMapBC)
tntNode tntNodeFindFirst(tntNetwork nw)
Definition: tntNodeInfo.c:50
tntNetwork tntNetworkCopy(tntNetwork nw)
Definition: tntNetwork.c:56
tntComplex tntMpoPmpoTrace(tntNetwork mpo, tntNodeArray *op, tntIntArray *sitenum)
Definition: tntMpoPmpo.c:35