v2DM-DOCI  1.0
FourIndex.cpp
Go to the documentation of this file.
1 /*
2  CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
3  Copyright (C) 2013, 2014 Sebastian Wouters
4 
5  This program is free software; you can redistribute it and/or modify
6  it under the terms of the GNU General Public License as published by
7  the Free Software Foundation; either version 2 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License along
16  with this program; if not, write to the Free Software Foundation, Inc.,
17  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19 
20 #include <stdlib.h>
21 #include <assert.h>
22 #include <iostream>
23 #include <sstream>
24 #include <string>
25 #include <cstring>
26 
27 #include "FourIndex.h"
28 #include "Lapack.h"
29 #include "MyHDF5.h"
30 
31 using namespace std;
32 
33 CheMPS2::FourIndex::FourIndex(const int nGroup, const int * IrrepSizes){
34 
35  SymmInfo.setGroup(nGroup);
36 
37  Isizes = new int[SymmInfo.getNumberOfIrreps()];
38  storage = new long long****[SymmInfo.getNumberOfIrreps()];
39 
40  for (int Icenter=0; Icenter<SymmInfo.getNumberOfIrreps(); Icenter++){
41  Isizes[Icenter] = IrrepSizes[Icenter];
42  }
43 
44  arrayLength = calcNumberOfUniqueElements(true); //true means allocate the storage!
45  theElements = new double[arrayLength];
46 
47 }
48 
49 CheMPS2::FourIndex::FourIndex(const CheMPS2::FourIndex &orig)
50 {
51  SymmInfo = orig.SymmInfo;
52 
53  Isizes = new int[SymmInfo.getNumberOfIrreps()];
54  storage = new long long****[SymmInfo.getNumberOfIrreps()];
55 
56  memcpy(Isizes, orig.Isizes, sizeof(int)*SymmInfo.getNumberOfIrreps());
57 
58  arrayLength = calcNumberOfUniqueElements(true); //true means allocate the storage!
59  theElements = new double[arrayLength];
60 
61  memcpy(theElements, orig.theElements, sizeof(double)*arrayLength);
62 }
63 
64 long long CheMPS2::FourIndex::calcNumberOfUniqueElements(const bool allocate){
65 
66  //The object size: see text above storage in Fourindex.h
67  long long theTotalSize = 0;
68 
69  for (int Icenter=0; Icenter<SymmInfo.getNumberOfIrreps(); Icenter++){
70  if (allocate){ storage[Icenter] = new long long***[SymmInfo.getNumberOfIrreps()]; }
71  for (int I_i=0; I_i<SymmInfo.getNumberOfIrreps(); I_i++){
72  int I_j = Irreps::directProd(Icenter,I_i);
73  if ((Isizes[I_i]>0)&&(Isizes[I_j]>0)){
74  if (allocate){ storage[Icenter][I_i] = new long long**[SymmInfo.getNumberOfIrreps()]; }
75  for (int I_k=I_i; I_k<SymmInfo.getNumberOfIrreps(); I_k++){
76  int I_l = Irreps::directProd(Icenter,I_k);
77  if ((Isizes[I_k]>0)&&(Isizes[I_l]>0)){
78  if ((I_i <= I_j) && (I_j <= I_l)){
79  if (Icenter == 0){ // I_i = I_j and I_k = I_l
80  if (I_i == I_k){
81  if (allocate){ storage[Icenter][I_i][I_k] = new long long*[Isizes[I_i]*(Isizes[I_i]+1)/2]; }
82  for (int i=0; i<Isizes[I_i]; i++){
83  for (int k=i; k<Isizes[I_k]; k++){
84  if (allocate){
85  storage[Icenter][I_i][I_k][i + k*(k+1)/2] = new long long[Isizes[I_j]-i];
86  for (int j=i; j<Isizes[I_j]; j++){
87  storage[Icenter][I_i][I_k][i + k*(k+1)/2][j-i] = theTotalSize;
88  theTotalSize += Isizes[I_l] - ((i==j) ? k : j);
89  }
90  } else {
91  delete [] storage[Icenter][I_i][I_k][i + k*(k+1)/2];
92  }
93  }
94  }
95  if (!allocate){ delete [] storage[Icenter][I_i][I_k]; }
96  } else { // I_i < I_k
97  if (allocate){ storage[Icenter][I_i][I_k] = new long long*[Isizes[I_i]*Isizes[I_k]]; }
98  for (int i=0; i<Isizes[I_i]; i++){
99  for (int k=0; k<Isizes[I_k]; k++){
100  if (allocate){
101  storage[Icenter][I_i][I_k][i + k*Isizes[I_i]] = new long long[Isizes[I_j]-i];
102  for (int j=i; j<Isizes[I_j]; j++){
103  storage[Icenter][I_i][I_k][i+k*Isizes[I_i]][j-i] = theTotalSize;
104  theTotalSize += Isizes[I_l] - ((i==j) ? k : 0 );
105  }
106  } else {
107  delete [] storage[Icenter][I_i][I_k][i + k*Isizes[I_i]];
108  }
109  }
110  }
111  if (!allocate){ delete [] storage[Icenter][I_i][I_k]; }
112  }
113  } else { //Icenter !=0 ; I_i < I_j and I_k != I_l
114  if (I_i == I_k){
115  if (allocate){ storage[Icenter][I_i][I_k] = new long long*[Isizes[I_i]*(Isizes[I_i]+1)/2]; }
116  for (int i=0; i<Isizes[I_i]; i++){
117  for (int k=i; k<Isizes[I_k]; k++){
118  if (allocate){
119  storage[Icenter][I_i][I_k][i + k*(k+1)/2] = new long long[Isizes[I_j]];
120  for (int j=0; j<Isizes[I_j]; j++){
121  storage[Icenter][I_i][I_k][i + k*(k+1)/2][j] = theTotalSize;
122  theTotalSize += Isizes[I_l] - j;
123  }
124  } else {
125  delete [] storage[Icenter][I_i][I_k][i + k*(k+1)/2];
126  }
127  }
128  }
129  if (!allocate){ delete [] storage[Icenter][I_i][I_k]; }
130  } else { // I_i < I_k
131  if (allocate){ storage[Icenter][I_i][I_k] = new long long*[Isizes[I_i]*Isizes[I_k]]; }
132  for (int i=0; i<Isizes[I_i]; i++){
133  for (int k=0; k<Isizes[I_k]; k++){
134  if (allocate){
135  storage[Icenter][I_i][I_k][i + k*Isizes[I_i]] = new long long[Isizes[I_j]];
136  for (int j=0; j<Isizes[I_j]; j++){
137  storage[Icenter][I_i][I_k][i + k*Isizes[I_i]][j] = theTotalSize;
138  theTotalSize += Isizes[I_l];
139  }
140  } else {
141  delete [] storage[Icenter][I_i][I_k][i + k*Isizes[I_i]];
142  }
143  }
144  }
145  if (!allocate){ delete [] storage[Icenter][I_i][I_k]; }
146  }
147  }
148  }
149  }
150  }
151  if (!allocate){ delete [] storage[Icenter][I_i]; }
152  }
153  }
154  if (!allocate){ delete [] storage[Icenter]; }
155  }
156 
157  return theTotalSize;
158 
159 }
160 
161 CheMPS2::FourIndex::~FourIndex(){
162 
163  arrayLength = calcNumberOfUniqueElements(false); //false means delete the storage!
164  delete [] theElements;
165  delete [] Isizes;
166  delete [] storage;
167 
168 }
169 
170 void CheMPS2::FourIndex::set(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l, const double val){
171 
172  theElements[getPointer(irrep_i, irrep_j, irrep_k, irrep_l, i, j, k, l)] = val;
173 
174 }
175 
176 void CheMPS2::FourIndex::add(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l, const double val){
177 
178  theElements[getPointer(irrep_i, irrep_j, irrep_k, irrep_l, i, j, k, l)] += val;
179 
180 }
181 
182 double CheMPS2::FourIndex::get(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l) const{
183 
184  return theElements[getPointer(irrep_i, irrep_j, irrep_k, irrep_l, i, j, k, l)];
185 
186 }
187 
188 long long CheMPS2::FourIndex::getPointer(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l) const {
189 
190  if (Irreps::directProd(irrep_i,irrep_j)==Irreps::directProd(irrep_k,irrep_l)){
191 
192  if ((irrep_i <= irrep_j) && (irrep_i <= irrep_k) && (irrep_j <= irrep_l)){ // (ijkl irrep ordering)
193  return getPtrIrrepOrderOK(irrep_i,irrep_j,irrep_k,irrep_l,i,j,k,l);
194  }
195 
196  if ((irrep_j <= irrep_i) && (irrep_j <= irrep_l) && (irrep_i <= irrep_k)){ // (jilk irrep ordering)
197  return getPtrIrrepOrderOK(irrep_j,irrep_i,irrep_l,irrep_k,j,i,l,k);
198  }
199 
200  if ((irrep_k <= irrep_j) && (irrep_k <= irrep_i) && (irrep_j <= irrep_l)){ // (kjil irrep ordering)
201  return getPtrIrrepOrderOK(irrep_k,irrep_j,irrep_i,irrep_l,k,j,i,l);
202  }
203 
204  if ((irrep_j <= irrep_k) && (irrep_j <= irrep_l) && (irrep_k <= irrep_i)){ // (jkli irrep ordering)
205  return getPtrIrrepOrderOK(irrep_j,irrep_k,irrep_l,irrep_i,j,k,l,i);
206  }
207 
208  if ((irrep_i <= irrep_l) && (irrep_i <= irrep_k) && (irrep_l <= irrep_j)){ // (ilkj irrep ordering)
209  return getPtrIrrepOrderOK(irrep_i,irrep_l,irrep_k,irrep_j,i,l,k,j);
210  }
211 
212  if ((irrep_l <= irrep_i) && (irrep_l <= irrep_j) && (irrep_i <= irrep_k)){ // (lijk irrep ordering)
213  return getPtrIrrepOrderOK(irrep_l,irrep_i,irrep_j,irrep_k,l,i,j,k);
214  }
215 
216  if ((irrep_k <= irrep_l) && (irrep_k <= irrep_i) && (irrep_l <= irrep_j)){ // (klij irrep ordering)
217  return getPtrIrrepOrderOK(irrep_k,irrep_l,irrep_i,irrep_j,k,l,i,j);
218  }
219 
220  if ((irrep_l <= irrep_k) && (irrep_l <= irrep_j) && (irrep_k <= irrep_i)){ // (lkji irrep ordering)
221  return getPtrIrrepOrderOK(irrep_l,irrep_k,irrep_j,irrep_i,l,k,j,i);
222  }
223 
224  }
225 
226  return -1;
227 
228 }
229 
230 long long CheMPS2::FourIndex::getPtrIrrepOrderOK(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l) const {
231 
232  //I_i <= I_j <= I_l and I_i <= I_k
233  int Icenter = Irreps::directProd(irrep_i,irrep_j);
234 
235  if (Icenter>0){ // I_i < I_j and I_k != I_l
236 
237  if (irrep_i == irrep_k){ //and hence I_j == I_l
238 
239  if (k>=i){
240  if (l>=j) return getPtrAllOK(5, Icenter, irrep_i, irrep_k, i, j, k, l);
241  else return getPtrAllOK(5, Icenter, irrep_i, irrep_k, i, l, k, j);
242  } else {
243  if (l>=j) return getPtrAllOK(5, Icenter, irrep_k, irrep_i, k, j, i, l);
244  else return getPtrAllOK(5, Icenter, irrep_k, irrep_i, k, l, i, j);
245  }
246 
247  } else return getPtrAllOK(6, Icenter, irrep_i, irrep_k, i, j, k, l);
248 
249  } else {
250 
251  if (irrep_i == irrep_k){ // all irreps the same
252 
253  //i en j
254  if ((i < j) && (i <= k) && (j <= l)) return getPtrAllOK(2, Icenter, irrep_i, irrep_k, i, j, k, l); // (ijkl ordering)
255  if ((i == j) && (i <= k) && (j <= l)){
256  if (l>=k) return getPtrAllOK(1, Icenter, irrep_i, irrep_k, i, j, k, l); // (ijkl ordering)
257  else return getPtrAllOK(1, Icenter, irrep_j, irrep_l, j, i, l, k); // (jilk ordering)
258  }
259  if ((j < i) && (j <= l) && (i <= k)) return getPtrAllOK(2, Icenter, irrep_j, irrep_l, j, i, l, k); // (jilk ordering)
260 
261  //k en l
262  if ((k < l) && (k <= i) && (l <= j)) return getPtrAllOK(2, Icenter, irrep_k, irrep_i, k, l, i, j); // (klij ordering)
263  if ((k == l) && (k <= i) && (l <= j)){
264  if (j>=i) return getPtrAllOK(1, Icenter, irrep_k, irrep_i, k, l, i, j); // (klij ordering)
265  else return getPtrAllOK(1, Icenter, irrep_l, irrep_j, l, k, j, i); // (lkji ordering)
266  }
267  if ((l < k) && (l <= j) && (k <= i)) return getPtrAllOK(2, Icenter, irrep_l, irrep_j, l, k, j, i); // (lkji ordering)
268 
269  // k en j
270  if ((k < j) && (k <= i) && (j <= l)) return getPtrAllOK(2, Icenter, irrep_k, irrep_i, k, j, i, l); // (kjil ordering)
271  if ((k == j) && (k <= i) && (j <= l)){
272  if (l>=i) return getPtrAllOK(1, Icenter, irrep_k, irrep_i, k, j, i, l); // (kjil ordering)
273  else return getPtrAllOK(1, Icenter, irrep_j, irrep_l, j, k, l, i); // (jkli ordering)
274  }
275  if ((j < k) && (j <= l) && (k <= i)) return getPtrAllOK(2, Icenter, irrep_j, irrep_l, j, k, l, i); // (jkli ordering)
276 
277  // i en l
278  if ((i < l) && (i <= k) && (l <= j)) return getPtrAllOK(2, Icenter, irrep_i, irrep_k, i, l, k, j); // (ilkj ordering)
279  if ((i == l) && (i <= k) && (l <= j)){
280  if (j>=k) return getPtrAllOK(1, Icenter, irrep_i, irrep_k, i, l, k, j); // (ilkj ordering)
281  else return getPtrAllOK(1, Icenter, irrep_l, irrep_j, l, i, j, k); // (lijk ordering)
282  }
283  if ((l < i) && (l <= j) && (i <= k)) return getPtrAllOK(2, Icenter, irrep_l, irrep_j, l, i, j, k); // (lijk ordering)
284 
285  } else {
286 
287  if (j==i){
288  if (l>=k) return getPtrAllOK(3, Icenter, irrep_i, irrep_k, i, j, k, l);
289  else return getPtrAllOK(3, Icenter, irrep_j, irrep_l, j, i, l, k);
290  } else {
291  if (j>i) return getPtrAllOK(4, Icenter, irrep_i, irrep_k, i, j, k, l);
292  else return getPtrAllOK(4, Icenter, irrep_j, irrep_l, j, i, l, k);
293  }
294 
295  }
296 
297  }
298 
299  return -1;
300 
301 }
302 
303 long long CheMPS2::FourIndex::getPtrAllOK(const int number, const int Icent, const int irrep_i, const int irrep_k, const int i, const int j, const int k, const int l) const {
304 
305  switch (number){
306  case 1:
307  return storage[Icent][irrep_i][irrep_k][i + k*(k+1)/2][j-i] + l-k;
308  case 2:
309  return storage[Icent][irrep_i][irrep_k][i + k*(k+1)/2][j-i] + l-j;
310  case 3:
311  return storage[Icent][irrep_i][irrep_k][i + Isizes[irrep_i]*k][j-i] + l-k;
312  case 4:
313  return storage[Icent][irrep_i][irrep_k][i + Isizes[irrep_i]*k][j-i] + l;
314  case 5:
315  return storage[Icent][irrep_i][irrep_k][i + k*(k+1)/2][j] + l-j;
316  case 6:
317  return storage[Icent][irrep_i][irrep_k][i + Isizes[irrep_i]*k][j] + l;
318  }
319 
320  return -1;
321 
322 }
323 
324 void CheMPS2::FourIndex::save(const std::string name) const{
325 
326  //The hdf5 file
327  hid_t file_id = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
328 
329  //The metadata
330  hid_t group_id = H5Gcreate(file_id, "/MetaData", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
331 
332  //The IrrepSizes
333  hsize_t dimarray = SymmInfo.getNumberOfIrreps();
334  hid_t dataspace_id = H5Screate_simple(1, &dimarray, NULL);
335  hid_t dataset_id = H5Dcreate(group_id, "IrrepSizes", H5T_STD_I32LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
336  H5Dwrite(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, Isizes);
337 
338  //Attributes
339  hid_t attribute_space_id1 = H5Screate(H5S_SCALAR);
340  hid_t attribute_id1 = H5Acreate(dataset_id, "nGroup", H5T_STD_I32LE, attribute_space_id1, H5P_DEFAULT, H5P_DEFAULT);
341  int nGroup = SymmInfo.getGroupNumber();
342  H5Awrite(attribute_id1, H5T_NATIVE_INT, &nGroup);
343 
344  hid_t attribute_space_id2 = H5Screate(H5S_SCALAR);
345  hid_t attribute_id2 = H5Acreate(dataset_id, "nIrreps", H5T_STD_I32LE, attribute_space_id2, H5P_DEFAULT, H5P_DEFAULT);
346  int nIrreps = SymmInfo.getNumberOfIrreps();
347  H5Awrite(attribute_id2, H5T_NATIVE_INT, &nIrreps);
348 
349  hid_t attribute_space_id3 = H5Screate(H5S_SCALAR);
350  hid_t attribute_id3 = H5Acreate(dataset_id, "theTotalSize", H5T_STD_I64LE, attribute_space_id3, H5P_DEFAULT, H5P_DEFAULT);
351  H5Awrite(attribute_id3, H5T_NATIVE_LLONG, &arrayLength);
352 
353  H5Aclose(attribute_id1);
354  H5Aclose(attribute_id2);
355  H5Aclose(attribute_id3);
356  H5Sclose(attribute_space_id1);
357  H5Sclose(attribute_space_id2);
358  H5Sclose(attribute_space_id3);
359 
360  H5Dclose(dataset_id);
361  H5Sclose(dataspace_id);
362 
363  H5Gclose(group_id);
364 
365  //The object itself
366  hid_t group_id7 = H5Gcreate(file_id, "/FourIndexObject", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
367 
368  hsize_t dimarray7 = arrayLength; //hsize_t is defined by default as unsigned long long, so no problem
369  hid_t dataspace_id7 = H5Screate_simple(1, &dimarray7, NULL);
370  hid_t dataset_id7 = H5Dcreate(group_id7, "Matrix elements", H5T_IEEE_F64LE, dataspace_id7, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
371  H5Dwrite(dataset_id7, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, theElements);
372 
373  H5Dclose(dataset_id7);
374  H5Sclose(dataspace_id7);
375 
376  H5Gclose(group_id7);
377 
378  H5Fclose(file_id);
379 
380 }
381 
382 
383 void CheMPS2::FourIndex::save2(const std::string name) const
384 {
385  //The hdf5 file
386  hid_t file_id = H5Fopen(name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
387 
388  hid_t group_id2 = H5Gcreate(file_id, "FourIndex", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
389 
390  //The metadata
391  hid_t group_id = H5Gcreate(group_id2, "MetaData", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
392 
393  //The IrrepSizes
394  hsize_t dimarray = SymmInfo.getNumberOfIrreps();
395  hid_t dataspace_id = H5Screate_simple(1, &dimarray, NULL);
396  hid_t dataset_id = H5Dcreate(group_id, "IrrepSizes", H5T_STD_I32LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
397  H5Dwrite(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, Isizes);
398 
399  //Attributes
400  hid_t attribute_space_id1 = H5Screate(H5S_SCALAR);
401  hid_t attribute_id1 = H5Acreate(dataset_id, "nGroup", H5T_STD_I32LE, attribute_space_id1, H5P_DEFAULT, H5P_DEFAULT);
402  int nGroup = SymmInfo.getGroupNumber();
403  H5Awrite(attribute_id1, H5T_NATIVE_INT, &nGroup);
404 
405  hid_t attribute_space_id2 = H5Screate(H5S_SCALAR);
406  hid_t attribute_id2 = H5Acreate(dataset_id, "nIrreps", H5T_STD_I32LE, attribute_space_id2, H5P_DEFAULT, H5P_DEFAULT);
407  int nIrreps = SymmInfo.getNumberOfIrreps();
408  H5Awrite(attribute_id2, H5T_NATIVE_INT, &nIrreps);
409 
410  hid_t attribute_space_id3 = H5Screate(H5S_SCALAR);
411  hid_t attribute_id3 = H5Acreate(dataset_id, "theTotalSize", H5T_STD_I64LE, attribute_space_id3, H5P_DEFAULT, H5P_DEFAULT);
412  H5Awrite(attribute_id3, H5T_NATIVE_LLONG, &arrayLength);
413 
414  H5Aclose(attribute_id1);
415  H5Aclose(attribute_id2);
416  H5Aclose(attribute_id3);
417  H5Sclose(attribute_space_id1);
418  H5Sclose(attribute_space_id2);
419  H5Sclose(attribute_space_id3);
420 
421  H5Dclose(dataset_id);
422  H5Sclose(dataspace_id);
423 
424  H5Gclose(group_id);
425 
426  //The object itself
427  hid_t group_id7 = H5Gcreate(group_id2, "FourIndexObject", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
428 
429  hsize_t dimarray7 = arrayLength; //hsize_t is defined by default as unsigned long long, so no problem
430  hid_t dataspace_id7 = H5Screate_simple(1, &dimarray7, NULL);
431  hid_t dataset_id7 = H5Dcreate(group_id7, "Matrix elements", H5T_IEEE_F64LE, dataspace_id7, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
432  H5Dwrite(dataset_id7, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, theElements);
433 
434  H5Dclose(dataset_id7);
435  H5Sclose(dataspace_id7);
436 
437  H5Gclose(group_id7);
438 
439  H5Gclose(group_id2);
440 
441  H5Fclose(file_id);
442 
443 }
444 
445 void CheMPS2::FourIndex::read(const std::string name){
446 
447  //The hdf5 file
448  hid_t file_id = H5Fopen(name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
449 
450  //The metadata
451  hid_t group_id = H5Gopen(file_id, "/MetaData",H5P_DEFAULT);
452 
453  //The IrrepSizes
454  hid_t dataset_id = H5Dopen(group_id, "IrrepSizes", H5P_DEFAULT);
455 
456  //Attributes
457  hid_t attribute_id1 = H5Aopen_by_name(group_id,"IrrepSizes", "nGroup", H5P_DEFAULT, H5P_DEFAULT);
458  int nGroup;
459  H5Aread(attribute_id1, H5T_NATIVE_INT, &nGroup);
460  assert( nGroup==SymmInfo.getGroupNumber() );
461 
462  hid_t attribute_id2 = H5Aopen_by_name(group_id,"IrrepSizes", "nIrreps", H5P_DEFAULT, H5P_DEFAULT);
463  int nIrreps;
464  H5Aread(attribute_id2, H5T_NATIVE_INT, &nIrreps);
465  assert( nIrreps==SymmInfo.getNumberOfIrreps() );
466 
467  hid_t attribute_id3 = H5Aopen_by_name(group_id,"IrrepSizes", "theTotalSize", H5P_DEFAULT, H5P_DEFAULT);
468  long long theTotalSize;
469  H5Aread(attribute_id3, H5T_NATIVE_LLONG, &theTotalSize);
470  assert( theTotalSize==arrayLength );
471 
472  H5Aclose(attribute_id1);
473  H5Aclose(attribute_id2);
474  H5Aclose(attribute_id3);
475 
476  int * IsizesAgain = new int[SymmInfo.getNumberOfIrreps()];
477  H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, IsizesAgain);
478  for (int cnt=0; cnt<SymmInfo.getNumberOfIrreps(); cnt++){
479  assert( IsizesAgain[cnt]==Isizes[cnt] );
480  }
481  delete [] IsizesAgain;
482  H5Dclose(dataset_id);
483 
484  H5Gclose(group_id);
485 
486  std::cout << "FourIndex::read : loading " << arrayLength << " doubles." << std::endl;
487 
488  //The object itself.
489  hid_t group_id7 = H5Gopen(file_id, "/FourIndexObject", H5P_DEFAULT);
490 
491  hid_t dataset_id7 = H5Dopen(group_id7, "Matrix elements", H5P_DEFAULT);
492  H5Dread(dataset_id7, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, theElements);
493  H5Dclose(dataset_id7);
494 
495  H5Gclose(group_id7);
496 
497  H5Fclose(file_id);
498 
499  std::cout << "FourIndex::read : everything loaded!" << std::endl;
500 
501 }
502 
503 
504 void CheMPS2::FourIndex::read2(const std::string name)
505 {
506  //The hdf5 file
507  hid_t file_id = H5Fopen(name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
508 
509  //The metadata
510  hid_t group_id = H5Gopen(file_id, "/FourIndex/MetaData",H5P_DEFAULT);
511 
512  //The IrrepSizes
513  hid_t dataset_id = H5Dopen(group_id, "IrrepSizes", H5P_DEFAULT);
514 
515  //Attributes
516  hid_t attribute_id1 = H5Aopen_by_name(group_id,"IrrepSizes", "nGroup", H5P_DEFAULT, H5P_DEFAULT);
517  int nGroup;
518  H5Aread(attribute_id1, H5T_NATIVE_INT, &nGroup);
519  assert( nGroup==SymmInfo.getGroupNumber() );
520 
521  hid_t attribute_id2 = H5Aopen_by_name(group_id,"IrrepSizes", "nIrreps", H5P_DEFAULT, H5P_DEFAULT);
522  int nIrreps;
523  H5Aread(attribute_id2, H5T_NATIVE_INT, &nIrreps);
524  assert( nIrreps==SymmInfo.getNumberOfIrreps() );
525 
526  hid_t attribute_id3 = H5Aopen_by_name(group_id,"IrrepSizes", "theTotalSize", H5P_DEFAULT, H5P_DEFAULT);
527  long long theTotalSize;
528  H5Aread(attribute_id3, H5T_NATIVE_LLONG, &theTotalSize);
529  assert( theTotalSize==arrayLength );
530 
531  H5Aclose(attribute_id1);
532  H5Aclose(attribute_id2);
533  H5Aclose(attribute_id3);
534 
535  int * IsizesAgain = new int[SymmInfo.getNumberOfIrreps()];
536  H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, IsizesAgain);
537  for (int cnt=0; cnt<SymmInfo.getNumberOfIrreps(); cnt++){
538  assert( IsizesAgain[cnt]==Isizes[cnt] );
539  }
540  delete [] IsizesAgain;
541  H5Dclose(dataset_id);
542 
543  H5Gclose(group_id);
544 
545  std::cout << "FourIndex::read : loading " << arrayLength << " doubles." << std::endl;
546 
547  //The object itself.
548  hid_t group_id7 = H5Gopen(file_id, "/FourIndex/FourIndexObject", H5P_DEFAULT);
549 
550  hid_t dataset_id7 = H5Dopen(group_id7, "Matrix elements", H5P_DEFAULT);
551  H5Dread(dataset_id7, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, theElements);
552  H5Dclose(dataset_id7);
553 
554  H5Gclose(group_id7);
555 
556  H5Fclose(file_id);
557 
558  std::cout << "FourIndex::read : everything loaded!" << std::endl;
559 }
560 
561 void CheMPS2::FourIndex::reset()
562 {
563  memset(theElements, 0, sizeof(double)*arrayLength);
564 }
CheMPS2::FourIndex::save
void save(const std::string name) const
Save the FourIndex object.
Definition: FourIndex.cpp:324
CheMPS2::FourIndex::getPtrIrrepOrderOK
long long getPtrIrrepOrderOK(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l) const
Definition: FourIndex.cpp:230
CheMPS2::FourIndex::save2
void save2(const std::string name) const
Definition: FourIndex.cpp:383
CheMPS2::FourIndex::read2
void read2(const std::string name)
Definition: FourIndex.cpp:504
CheMPS2::FourIndex::reset
void reset()
set everything to zero
Definition: FourIndex.cpp:561
std
STL namespace.
CheMPS2::FourIndex::get
double get(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l) const
Get an element.
Definition: FourIndex.cpp:182
CheMPS2::FourIndex::calcNumberOfUniqueElements
long long calcNumberOfUniqueElements(const bool allocateStorage)
Definition: FourIndex.cpp:64
CheMPS2::FourIndex::getPtrAllOK
long long getPtrAllOK(const int number, const int Icent, const int irrep_i, const int irrep_k, const int i, const int j, const int k, const int l) const
Definition: FourIndex.cpp:303
CheMPS2::FourIndex::~FourIndex
virtual ~FourIndex()
Destructor.
Definition: FourIndex.cpp:161
CheMPS2::FourIndex::set
void set(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l, const double val)
Set an element.
Definition: FourIndex.cpp:170
CheMPS2::FourIndex::theElements
double * theElements
Definition: FourIndex.h:136
CheMPS2::FourIndex::getPointer
long long getPointer(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l) const
Definition: FourIndex.cpp:188
CheMPS2::FourIndex::FourIndex
FourIndex(const int nGroup, const int *IrrepSizes)
Constructor.
Definition: FourIndex.cpp:33
CheMPS2::FourIndex::add
void add(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l, const double val)
Add a double to an element.
Definition: FourIndex.cpp:176
CheMPS2::FourIndex::read
void read(const std::string name)
Load the FourIndex object.
Definition: FourIndex.cpp:445