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
tntMpsExpecOutput.c
1 /*
2 Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch
3 $LastChangedDate: 2013-10-09 10:02:48 +0100 (Wed, 09 Oct 2013) $
4 (c) University of Oxford 2014
5 */
6 
12 /* Include the header for the TNT MPS library */
13 #include "tntMpsInternal.h"
14 
15 /* Internal function to avoid repeated code */
16 void _tnt_MpsExpecOutputOneOp(tntComplexArray *Exval, int printOutput, int saveOutput, char *savename, char *label, unsigned counter);
17 
18 
37 void tntMpsExpecOutput(tntNetwork mps,
38  tntMpsExOp *Op,
39  int orthCentre,
40  unsigned printOutput,
41  unsigned saveOutput,
42  char *savename,
43  unsigned counter)
46 {
47 
48  tntComplexArray Exval; /* array for holding expectation values */
49  double normsq; /* Norm squared of the wave function. Always calculated */
50  unsigned loop; /* Used to loop through the different expectation value types */
51  tntNodeArray op; /* Node array for holding operators for expectation value currently being calculated */
52  tntIntArray sitenum; /* Array for holding site number for expectation value currently being calculated */
53  int L; /* Length of the MPS */
54  tntNodeArray betap, gammap, *betas = &betap, *gammas = &gammap; /* Node arrays for precontracted parts of the network */
55 
56  if (!printOutput && !saveOutput) {
57  tntWarningPrint("No expectation values are being calculated as both output flags are set to off"); /* NO_COVERAGE */
58  return; /* NO_COVERAGE */
59  } /* NO_COVERAGE */
60 
62 
63  /* Find the length of the MPS */
64  L = (int) tntMpsLength(mps);
65 
66  /* Calculate the norm squared */
67  normsq = tntMpsSelfProduct(mps,orthCentre);
68 
69  if (printOutput) {
70  tntPrintf("\n~~~~~~~ Norm squared is % #6.6g ~~~~~~~~\n\n",normsq);
71  }
72  if (saveOutput) {
73  if (counter) {
74  tntDoubleParamsUpdate(savename, counter, normsq);
75  } else {
76  tntDoubleParamsSave(savename, normsq);
77  }
78 
79  }
80 
81  /* Calculate precontracted nodes since multiple types of the same contraction will be performed */
82  if (0 == orthCentre) {
83  tntMpsSelfInit(mps, betas, NULL);
84  gammas = NULL;
85  } else if (L-1 == orthCentre) {
86  tntMpsSelfInit(mps, NULL, gammas);
87  betas = NULL;
88  } else {
89  tntMpsSelfInit(mps, betas, gammas);
90  }
91 
92 
93  /* ----- Onsite expectation values ----- */
94 
95  if (Op->ExOpOs.sz > 0) {
96 
97  /* Allocate arrays of size one for onsite operators */
98  op = tntNodeArrayAlloc(1);
99  sitenum = tntIntArrayAlloc(1);
100 
101  if (printOutput) {
102  tntPrintf("\n~~~~~~~ Onsite expectation values ~~~~~~~~\n");
103  }
104 
105  }
106  /* Find the single site expectation values of the calculated ground state */
107  for (loop = 0; loop < Op->ExOpOs.sz; loop++) {
108 
109  /* Allocate memory for storing single-site expectation values */
110  Exval = tntComplexArrayAlloc(L);
111 
112  /* Copy over the required operator */
113  op.vals[0] = tntNodeCopy(Op->ExOpOs.vals[loop],0);
114 
115  for (sitenum.vals[0] = 0; sitenum.vals[0] < L; (sitenum.vals[0])++) {
116  /* Find the expectation value on each site using a product MPO. */
117  Exval.vals[sitenum.vals[0]] = tntMpsPmpoMpsProduct(mps, &op, &sitenum, orthCentre, betas, gammas);
118  Exval.vals[sitenum.vals[0]].re /= normsq;
119  Exval.vals[sitenum.vals[0]].im /= normsq;
120  }
121 
122  /* output and free values */
123  _tnt_MpsExpecOutputOneOp(&Exval, printOutput, saveOutput, savename, Op->LbOpOs.vals[loop], counter);
124 
125  /* Free the copy of the operator */
126  tntNodeFree(&(op.vals[0]));
127  }
128  if (Op->ExOpOs.sz > 0) {
129 
130  /* free the arrays */
131  tntIntArrayFree(&sitenum);
132  tntNodeArrayFree(&op);
133 
134  if (printOutput) {
135  tntPrintf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
136  }
137  }
138 
139 
140  /* ---------- Nearest-neighbour expectation values ----------- */
141 
142  if (Op->ExOp2nn.sz > 0) {
143 
144  /* Allocate arrays of size two for two-site operator site numbers */
145  sitenum = tntIntArrayAlloc(2);
146 
147  if (printOutput) {
148  tntPrintf("\n~~~~~~~ Nearest-neighbour expectation values ~~~~~~~~\n");
149  }
150 
151  }
152  for (loop = 0; loop < Op->ExOp2nn.sz; loop+=2) {
153 
154  /* Allocate memory for storing two-site expectation values */
155  Exval = tntComplexArrayAlloc(L-1);
156 
157  /* Allocate arrays of size two for two-site operators */
158  op = tntNodeArrayAlloc(2);
159 
160  /* Copy over the required operator */
161  op.vals[0] = tntNodeCopy(Op->ExOp2nn.vals[loop],0);
162  op.vals[1] = tntNodeCopy(Op->ExOp2nn.vals[loop+1],0);
163 
164  for (sitenum.vals[0] = 0; sitenum.vals[0] < L-1; (sitenum.vals[0])++) {
165  /* assign value for second site */
166  sitenum.vals[1] = sitenum.vals[0]+1;
167  /* Find the expectation value on each site using a product MPO. */
168  Exval.vals[sitenum.vals[0]] = tntMpsPmpoMpsProduct(mps, &op, &sitenum, orthCentre, betas, gammas);
169  Exval.vals[sitenum.vals[0]].re /= normsq;
170  Exval.vals[sitenum.vals[0]].im /= normsq;
171  }
172 
173  /* Free the operators */
174  tntNodeArrayFree(&op);
175 
176  /* output and free values */
177  _tnt_MpsExpecOutputOneOp(&Exval, printOutput, saveOutput, savename, Op->LbOp2nn.vals[loop/2], counter);
178 
179  }
180  if (Op->ExOp2nn.sz > 0) {
181 
182  /* free the arrays */
183  tntIntArrayFree(&sitenum);
184 
185  if (printOutput) {
186  tntPrintf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
187  }
188  }
189 
190  /* ---------- Centre-site expectation values ----------- */
191 
192  if (Op->ExOp2cs.sz > 0) {
193 
194  /* Allocate arrays of size two for two-site operator site numbers */
195  sitenum = tntIntArrayAlloc(2);
196 
197  if (printOutput) {
198  tntPrintf("\n~~~~~~~ Centre-site expectation values ~~~~~~~~\n");
199  }
200  }
201  for (loop = 0; loop < Op->ExOp2cs.sz; loop+=2) {
202 
203  /* Allocate memory for storing two-site expectation values */
204  Exval = tntComplexArrayAlloc(L);
205 
206  /* Allocate first site that expectation value will be found on */
207  sitenum.vals[0] = L/2;
208 
209  /* Allocate arrays of size two for two-site operators */
210  op = tntNodeArrayAlloc(2);
211 
212  /* Copy over the required operator */
213  op.vals[0] = tntNodeCopy(Op->ExOp2cs.vals[loop]);
214  op.vals[1] = tntNodeCopy(Op->ExOp2cs.vals[loop+1]);
215 
216  for (sitenum.vals[1] = 0; sitenum.vals[1] < L; (sitenum.vals[1])++) {
217  /* Find the expectation value on each site using a product MPO. */
218  Exval.vals[sitenum.vals[1]] = tntMpsPmpoMpsProduct(mps, &op, &sitenum, orthCentre, betas, gammas);
219  Exval.vals[sitenum.vals[1]].re /= normsq;
220  Exval.vals[sitenum.vals[1]].im /= normsq;
221  }
222 
223  /* Free the operators */
224  tntNodeArrayFree(&op);
225 
226  /* output and free values */
227  _tnt_MpsExpecOutputOneOp(&Exval, printOutput, saveOutput, savename, Op->LbOp2cs.vals[loop/2], counter);
228 
229  }
230  if (Op->ExOp2cs.sz > 0) {
231 
232  /* free the arrays */
233  tntIntArrayFree(&sitenum);
234 
235  if (printOutput) {
236  tntPrintf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
237  }
238  }
239 
240 
241  /* ---------- All pairs of two-site expectation values ----------- */
242 
243  if (Op->ExOp2ap.sz > 0) {
244 
245  /* Allocate arrays of size two for sitenumbers for two-site operators */
246  sitenum = tntIntArrayAlloc(2);
247 
248  if (printOutput) {
249  tntPrintf("\n~~~~~~~~~ All pairs of two-site expectation values ~~~~~~~~~\n");
250  }
251 
252  }
253  for (loop = 0; loop < Op->ExOp2ap.sz; loop+=2) {
254 
255  int ind; /* total matrix index */
256 
257  /* Allocate memory for storing two-site expectation values */
258  Exval = tntComplexArrayAlloc(L,L);
259 
260  /* Allocate arrays of size two for two-site operators */
261  op = tntNodeArrayAlloc(2);
262 
263  /* Copy over the required operator */
264  op.vals[0] = tntNodeCopy(Op->ExOp2ap.vals[loop]);
265  op.vals[1] = tntNodeCopy(Op->ExOp2ap.vals[loop+1]);
266 
267  /* Find the expectation value on each pair of sites using a product MPO. */
268  for (sitenum.vals[0] = 0; sitenum.vals[0] < L; (sitenum.vals[0])++) {
269  for (sitenum.vals[1] = 0; sitenum.vals[1] < L; (sitenum.vals[1])++) {
270  ind = sitenum.vals[0] + L*sitenum.vals[1];
271  Exval.vals[ind] = tntMpsPmpoMpsProduct(mps, &op, &sitenum, orthCentre, betas, gammas);
272  Exval.vals[ind].re /= normsq;
273  Exval.vals[ind].im /= normsq;
274  }
275  }
276 
277  /* Free the operators */
278  tntNodeArrayFree(&op);
279 
280  /* output and free values */
281  _tnt_MpsExpecOutputOneOp(&Exval, printOutput, saveOutput, savename, Op->LbOp2ap.vals[loop/2], counter);
282 
283  }
284  if (Op->ExOp2ap.sz > 0) {
285 
286  /* free the arrays */
287  tntIntArrayFree(&sitenum);
288 
289  if (printOutput) {
290  tntPrintf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
291  }
292 
293  }
294 
295  if (0 == orthCentre) {
296  tntNodeArrayFree(betas);
297  } else if (L-1 == orthCentre) {
298  tntNodeArrayFree(gammas);
299  } else {
300  tntNodeArrayFree(betas);
301  tntNodeArrayFree(gammas);
302  }
303 
305 }
306 
307 /* prints, saves and frees the complex array, converting to double if necessary, and taking account of flags and counters */
308 void _tnt_MpsExpecOutputOneOp(tntComplexArray *Exval, int printOutput, int saveOutput, char *savename, char *label, unsigned counter) {
309 
310  tntDoubleArray Exvalr; /* array for holding real part of expectation values */
311  tntComplexArray Exvalc; /* array for holding real part of expectation values */
312 
313  /* Now take different action depending on whether the expectation values are real or complex */
314  if (tntComplexArrayIsReal(Exval)) {
315  /* change complex array to real array */
316  Exvalr = tntComplexArrayToReal(Exval);
317  /* Print real output */
318  if (printOutput) {
319  tntNamedDoubleArrayPrint(label,&Exvalr);
320  }
321  /* Save this expectation value */
322  if (saveOutput) {
323  if (counter) tntDoubleArraysNamedUpdate(savename, counter,Exvalr, label);
324  else tntDoubleArraysNamedSave(savename, Exvalr, label);
325 
326  }
327 
328  /* free real array */
329  tntDoubleArrayFree(&Exvalr);
330  } else {
331  Exvalc = *Exval;
332  /* Print complex output */
333  if (printOutput) {
334  tntNamedComplexArrayPrint(label,Exval);
335  }
336  /* Save complex output */
337  if (saveOutput) {
338  if (counter) tntComplexArraysNamedUpdate(savename, counter, Exvalc, label);
339  else tntComplexArraysNamedSave(savename, Exvalc, label);
340  }
341  /* Free complex array */
342  tntComplexArrayFree(Exval);
343  }
344 }
345 
tntNode tntNodeCopy(tntNode A)
Definition: tntNodeUtil.c:304
#define tntMpsExOp
Definition: tntMps.h:75
void tntNodeArrayFree(tntNodeArray *arr)
Definition: tntArray.c:668
int tntComplexArrayIsReal(tntComplexArray *arr)
Definition: tntArray.c:570
void tntDoubleParamsSave(const char *fname, double varA, double varB, double varC, double varD, double varE)
Definition: tntSave.c:73
tntNodeArray tntNodeArrayAlloc(unsigned num_elems)
Definition: tntArray.c:428
tntDoubleArray tntComplexArrayToReal(tntComplexArray *arr)
Definition: tntArray.c:591
void tntDoubleArraysNamedUpdate(const char *fname, unsigned posn, tntDoubleArray arrA, const char *varnameA, tntDoubleArray arrB, const char *varnameB, tntDoubleArray arrC, const char *varnameC, tntDoubleArray arrD, const char *varnameD, tntDoubleArray arrE, const char *varnameE)
Definition: tntUpdate.c:302
void tntDoubleParamsUpdate(const char *fname, unsigned posn, double varA, double varB, double varC, double varD, double varE)
Definition: tntUpdate.c:88
double tntMpsSelfProduct(tntNetwork mps, int orth_centre)
Definition: tntMpsMps.c:544
void tntNodeFree(tntNode *A)
Definition: tntNodeUtil.c:275
tntComplexArray tntComplexArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntArray.c:392
unsigned tntMpsLength(tntNetwork wf)
Definition: tntMpsUtil.c:17
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray.c:620
void tntDoubleArrayFree(tntDoubleArray *arr)
Definition: tntArray.c:635
void tntMpsExpecOutput(tntNetwork mps, tntMpsExOp *Op, int orthCentre, unsigned printOutput, unsigned saveOutput, char *savename, unsigned counter)
tntIntArray tntIntArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntArray.c:29
tntComplex tntMpsPmpoMpsProduct(tntNetwork mps, tntNodeArray *op, tntIntArray *sitenum, int orth_centre, tntNodeArray *betas, tntNodeArray *gammas)
Definition: tntMpsPmpoMps.c:57
void tntComplexArraysNamedSave(const char *fname, tntComplexArray arrA, const char *varnameA, tntComplexArray arrB, const char *varnameB, tntComplexArray arrC, const char *varnameC, tntComplexArray arrD, const char *varnameD, tntComplexArray arrE, const char *varnameE)
Definition: tntSave.c:327
void tntComplexArrayFree(tntComplexArray *arr)
Definition: tntArray.c:650
void tntMpsSelfInit(tntNetwork mps, tntNodeArray *betas, tntNodeArray *gammas)
Definition: tntMpsMps.c:329
int tntPrintf(const char *format,...)
Definition: tntPrint.c:77
void tntSysQNClearWarnOn(void)
Definition: tntUtil.c:17
void tntSysQNClearWarnOff(void)
Definition: tntUtil.c:35
void tntDoubleArraysNamedSave(const char *fname, tntDoubleArray arrA, const char *varnameA, tntDoubleArray arrB, const char *varnameB, tntDoubleArray arrC, const char *varnameC, tntDoubleArray arrD, const char *varnameD, tntDoubleArray arrE, const char *varnameE)
Definition: tntSave.c:271
void tntComplexArraysNamedUpdate(const char *fname, unsigned posn, tntComplexArray arrA, const char *varnameA, tntComplexArray arrB, const char *varnameB, tntComplexArray arrC, const char *varnameC, tntComplexArray arrD, const char *varnameD, tntComplexArray arrE, const char *varnameE)
Definition: tntUpdate.c:371