Tensor Network Theory Library  Beta release 1.2.1
A library of routines for performing TNT-based operations
 All Data Structures Functions Variables Groups Pages
tntMpsPmpo.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 void tntMpsPmpoProduct(tntNetwork mps,
36  tntNodeArray *op,
39  tntIntArray *sitenum)
42 {
43  /* Networks and nodes during the contraction */
44  tntNode A; /* The current MPS node */
45  tntNode opn; /* A copy of the operator to use for contracting */
46 
47  /* Counters and site information used during the contraction */
48  int sitecount; /* the current site count */
49  int L; /* The number of MPS nodes in the network */
50  unsigned opnum=0; /* Counter for the operators */
51  tntIntArray siterank; /* Rank for each site-number, to make then in ascending order */
52  tntIntArray sitenum_cp; /* Copy of the site numbers which will be used for sorting */
53  unsigned i, j; /* Used for looping */
54  int sitenummax, sitemax; /* Position of current maximum site, current maximum site */
55 
56  /* Determine the length of the MPS. */
57  L = (int) tntMpsLength(mps);
58 
59  /* Check that operators if they have been provided */
60  if ((NULL == op) || (NULL == sitenum) || (0 == op->sz) || (0 == sitenum->sz)) {
61  tntWarningPrint("Leaving function tntMpsPmpoProduct() with no action taken as no operators or sitenumbers have been provided"); /* NO_COVERAGE */
62  return; /* NO_COVERAGE */
63  } /* NO_COVERAGE */
64 
65  /* Allocate arrays for copies of the site numbers, and rank of site numbers */
66  sitenum_cp = tntIntArrayAlloc(op->sz);
67  siterank = tntIntArrayAlloc(op->sz);
68 
69  /* Check that the sites and orthogonality centre given are sensible, and at the same time copy over site numbers */
70  for (opnum = 0; opnum < op->sz; opnum++) {
71  if (sitenum->vals[opnum] >= L) tntErrorPrint("The site for the operator %d is %d, which is greater than the length of the MPS (%d)", opnum, sitenum->vals[opnum], L);
72 
73  if (sitenum->vals[opnum] < 0) tntErrorPrint("The site for the operator %d is %d - site numbers should be zero or more", opnum, sitenum->vals[opnum]); /* NO_COVERAGE */
74 
75  sitenum_cp.vals[opnum] = sitenum->vals[opnum];
76  }
77 
78  /* First give each element of sitenum a rank, where siterank i gives the element of ith biggest site */
79  /* Loop through all legs from end to beginning. For each index bl, find the blth largest leg dimension */
80  for (i = opnum; i-->0; ) {
81  /* Reset marker for largest site number */
82  sitemax = -1;
83  for (j = opnum; j-- > 0; ) {
84  if (sitenum_cp.vals[j] > sitemax) {
85  sitemax = sitenum_cp.vals[j];
86  sitenummax = j;
87  }
88  }
89  /* 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 */
90  sitenum_cp.vals[sitenummax] = -1;
91  siterank.vals[i] = sitenummax;
92  }
93 
94  /* Copy of site numbers no longer required */
95  tntIntArrayFree(&sitenum_cp);
96 
97  /* Reset counter for number of operators */
98  opnum = 0;
99 
100  /* 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 */
101 
102  /* identify the first node */
103  A = tntNodeFindFirst(mps);
104 
105  /* Go through remaining sites, contracting them with alpha one by one */
106  for (sitecount = 0; sitecount < L; sitecount++) {
107 
108  /* If there are operators on this site, join to the physical leg and contract with A */
109  while ((opnum < op->sz) && (sitenum->vals[siterank.vals[opnum]] == sitecount)) {
110  /* Make a copy of the operator */
111  opn = tntNodeCopy(op->vals[siterank.vals[opnum]]);
112 
113  /* join the operator to the physical leg */
114  tntNodeJoin(opn, "U", A, "D");
115 
116  /* contract the first node with the operator on the physical leg */
117  A = tntNodeContract(A,opn);
118 
119  /* increment the counter for the operator number */
120  opnum++;
121  }
122 
123  /* Break if all the operators have been applied */
124  if (opnum == op->sz) {
125  break;
126  }
127 
128  /* Find the next node along in the MPS */
129  A = tntNodeFindConn(A, "R");
130  }
131 
132  /* Free the site rank */
133  tntIntArrayFree(&siterank);
134 
135  return;
136 
137 }
138 
139 
140 
146 void tntMpsNodeProduct(tntNetwork mps,
147  tntNode op,
148  unsigned sitenum)
149 {
150  tntIntArray sitenum_arr = tntIntArrayAlloc(1); /* Current site number */
151  tntNodeArray op_arr = tntNodeArrayAlloc(1); /* Current operator */
152 
153  sitenum_arr.vals[0] = sitenum;
154  op_arr.vals[0] = tntNodeCopy(op);
155 
156  /* Create a site-wide MPO to preserve QN information if the mps is quantum number conserving */
158  tntNetwork pmpo = tntMpsCreateProductMpo(tntMpsLength(mps), &op_arr, &sitenum_arr);
159  double err = tntMpsMpoProduct(mps, pmpo, -1);
160  err = (double) err;
161  } else {
162  tntMpsPmpoProduct(mps, &op_arr, &sitenum_arr);
163  }
164 
165  tntIntArrayFree(&sitenum_arr);
166  tntNodeArrayFree(&op_arr);
167 
168  return;
169 }
170 
tntNode tntNodeCopy(tntNode A)
Definition: tntNodeUtil.c:304
void tntNodeArrayFree(tntNodeArray *arr)
Definition: tntArray.c:668
void tntNodeJoin(tntNode A, tntLegLabel legA, tntNode B, tntLegLabel legB)
Definition: tntNodeConn.c:52
tntNodeArray tntNodeArrayAlloc(unsigned num_elems)
Definition: tntArray.c:428
unsigned tntMpsLength(tntNetwork wf)
Definition: tntMpsUtil.c:17
unsigned tntNodeIsCovariant(tntNode A)
Definition: tntNodeInfo.c:535
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray.c:620
double tntMpsMpoProduct(tntNetwork mps, tntNetwork mpo, int chi)
Definition: tntMpsMpo.c:225
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
void tntMpsPmpoProduct(tntNetwork mps, tntNodeArray *op, tntIntArray *sitenum)
Definition: tntMpsPmpo.c:35
tntNetwork tntMpsCreateProductMpo(unsigned L, tntNodeArray *op, tntIntArray *sitenum)
void tntMpsNodeProduct(tntNetwork mps, tntNode op, unsigned sitenum)
Definition: tntMpsPmpo.c:146