u1casc-ordinary/obs_diagcorr.hpp

   1 #ifndef OBS_DIAGCORR_HPP
   2 #define OBS_DIAGCORR_HPP
   3
   4 #include "latlib/o815/o815.h"
   5
   6 #include "latlib/writeout.h"
   7
   8 #include "latlib/obstat.hpp"
   9
  10 #include <iostream>
  11 #include <complex>
  12
  13 #include <cmath>
  14
  15 #include <gsl/gsl_complex.h>
  16 #include <gsl/gsl_complex_math.h>
  17 #include <gsl/gsl_math.h>
  18 #include <gsl/gsl_eigen.h>
  19 #include <gsl/gsl_linalg.h>
  20 #include <gsl/gsl_cblas.h>
  21 #include <gsl/gsl_blas.h>
  22
  23 using namespace std;
  24
  25 class obs_diagcorr : public o815::obs {
  26
  27 public:
  28   obs_diagcorr(o815 *_O815);
  29
  30 private:
  31   void _start();
  32   void _meas(bool loadedobs, const int& nthmeas);
  33   void _finish();
  34
  35   void corrCompute();
  36   static complex<double> effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para);
  37   static void preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para);
  38
  39   sim *Sim;
  40
  41   int spatialV;
  42
  43   complex<double> *OM[16+1];
  44
  45   gsl_matrix_complex ***measurements;
  46
  47   gsl_vector_complex **measurements_disco;
  48
  49   static void cdiag (gsl_vector *v, gsl_matrix_complex *m);
  50
  51   static void invertcm (gsl_matrix_complex *invm, const gsl_matrix_complex *m);
  52 };
  53
  54 obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr",
  55                                                     _O815->paraQ->getParaNames() +
  56                                                     "tsep"
  57                                                     "...",
  58                                                     _O815, sizeof(complex<double>) * ( 16*_O815->comargs.lsize[1]/2 + 1*4 ) ) {
  59   for (int ivar = 0; ivar<16+1; ivar++)
  60     OM[ivar]  = (complex<double>*)( obsMem + ivar*sizeof(complex<double>)*( _O815->comargs.lsize[1]/2 ) );
  61
  62   Sim = (sim*)O815->Sim;
  63   spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0];
  64
  65   measurements = new gsl_matrix_complex**[O815->comargs.nmeas];
  66   for (int imeas = 0; imeas<O815->comargs.nmeas; imeas++) {
  67     measurements[imeas] = new gsl_matrix_complex*[_O815->comargs.lsize[1]/2];
  68     for (int itsep=0; itsep<_O815->comargs.lsize[1]/2; itsep++)
  69       measurements[imeas][itsep] = gsl_matrix_complex_alloc(4,4);
  70   }
  71
  72   measurements_disco = new gsl_vector_complex*[O815->comargs.nmeas];
  73   for (int imeas = 0; imeas<O815->comargs.nmeas; imeas++)
  74     measurements_disco[imeas] = gsl_vector_complex_alloc(4);
  75 }
  76
  77 void obs_diagcorr::_start() {
  78   //  *out << O815->comargs.nmeas << endl;
  79
  80   //*out << "OBS_test: start" << endl;
  81 };
  82
  83 void obs_diagcorr::_meas(bool loadedobs, const int& nthmeas) {
  84   if (!loadedobs)
  85     corrCompute();
  86
  87   for (int icorr=0; icorr<4; icorr++)
  88     for (int jcorr=0; jcorr<4; jcorr++)
  89       for (int itsep=0; itsep<O815->comargs.lsize[1]/2; itsep++) {
  90         gsl_matrix_complex_set(measurements[nthmeas][itsep], icorr, jcorr,
  91                                gsl_complex_rect( OM[icorr*4+jcorr][itsep].real(), OM[icorr*4+jcorr][itsep].imag() ));
  92       }
  93
  94   for (int icorr=0; icorr<4; icorr++) {
  95     gsl_vector_complex_set(measurements_disco[nthmeas], icorr,
  96                            gsl_complex_rect( OM[16][icorr].real(), OM[16][icorr].imag() ) );
  97   }
  98 };
  99
 100 void obs_diagcorr::invertcm (gsl_matrix_complex *invm, const gsl_matrix_complex *m)
 101 {
 102   int s;
 103   gsl_permutation * perm = gsl_permutation_alloc (m->size1);
 104
 105   gsl_matrix_complex *mtmp = gsl_matrix_complex_alloc(m->size1,m->size2);
 106   gsl_matrix_complex_memcpy(mtmp, m);
 107
 108   gsl_linalg_complex_LU_decomp (mtmp, perm, &s);
 109   gsl_linalg_complex_LU_invert (mtmp, perm, invm);
 110
 111   gsl_matrix_complex_free(mtmp);
 112 }
 113
 114 void obs_diagcorr::cdiag (gsl_vector *v, gsl_matrix_complex *m)
 115 {
 116   gsl_matrix_complex *evec = gsl_matrix_complex_alloc(4,4);
 117   gsl_eigen_hermv_workspace *wspace = gsl_eigen_hermv_alloc(100);
 118
 119   gsl_eigen_hermv(m, v, evec, wspace);
 120   gsl_eigen_hermv_sort(v, evec, GSL_EIGEN_SORT_VAL_ASC);
 121
 122   gsl_eigen_hermv_free(wspace);
 123   gsl_matrix_complex_free(evec);
 124 }
 125
 126 void obs_diagcorr::_finish() {
 127   gsl_matrix_complex *totalval = gsl_matrix_complex_alloc(4,4);
 128   gsl_matrix_complex *tmpmatrix = gsl_matrix_complex_alloc(4,4);
 129   gsl_matrix_complex *tmpmatrix2 = gsl_matrix_complex_alloc(4,4);
 130   gsl_vector *tmpvec = gsl_vector_alloc(4);
 131   gsl_vector_complex *tmpvecc = gsl_vector_complex_alloc(4);
 132   gsl_vector *tmpvec2 = gsl_vector_alloc(4);
 133   gsl_vector *jackres = gsl_vector_alloc(4);
 134   gsl_vector *jackerrornorm = gsl_vector_alloc(4);
 135   gsl_vector_complex *totaldisco = gsl_vector_complex_alloc(4);
 136   gsl_matrix_complex *invC0 = NULL;
 137
 138   gsl_vector_complex_set_zero(totaldisco);
 139   for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
 140     gsl_vector_complex_add( totaldisco, measurements_disco[imeas] );
 141
 142   for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
 143     gsl_matrix_complex_set_zero(totalval);
 144     gsl_vector_set_zero(jackerrornorm);
 145
 146     for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
 147       gsl_matrix_complex_add( totalval, measurements[imeas][itsep] );
 148
 149     gsl_matrix_complex_memcpy (tmpmatrix, totalval);
 150     gsl_matrix_complex_scale (tmpmatrix, gsl_complex_rect(1.0/O815->comargs.nmeas, 0.0) );
 151     for (int icorr=0; icorr<4; icorr++)
 152       for (int jcorr=0; jcorr<4; jcorr++) {
 153         gsl_complex discopart = gsl_complex_mul( gsl_vector_complex_get(totaldisco, icorr),
 154                                                  gsl_complex_conjugate( gsl_vector_complex_get(totaldisco, jcorr) )
 155                                                  );
 156         discopart = gsl_complex_mul_real(discopart, pow(1.0/O815->comargs.nmeas,2));
 157         gsl_matrix_complex_set( tmpmatrix, icorr, jcorr,
 158                                 gsl_complex_sub( gsl_matrix_complex_get( tmpmatrix, icorr, jcorr ),
 159                                                  discopart
 160                                                  )
 161                                 );
 162       }
 163
 164     if (invC0 == NULL) {
 165       invC0 = gsl_matrix_complex_alloc(4,4);
 166       invertcm(invC0, tmpmatrix);
 167     }
 168
 169     gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, gsl_complex_rect(1.0,0.0),
 170                     invC0, tmpmatrix, gsl_complex_rect(0.0,0.0), tmpmatrix2);
 171     cdiag(jackres, tmpmatrix2);
 172
 173     for (int imeas=0; imeas<O815->comargs.nmeas; imeas++) {
 174       gsl_matrix_complex_memcpy (tmpmatrix, totalval);
 175       gsl_matrix_complex_sub (tmpmatrix, measurements[imeas][itsep]);
 176       gsl_matrix_complex_scale ( tmpmatrix, gsl_complex_rect(1.0/(O815->comargs.nmeas-1), 0.0) );
 177       gsl_vector_complex_memcpy( tmpvecc, totaldisco );
 178       gsl_vector_complex_sub( tmpvecc, measurements_disco[imeas] );
 179       for (int icorr=0; icorr<4; icorr++)
 180         for (int jcorr=0; jcorr<4; jcorr++) {
 181           gsl_complex discopart = gsl_complex_mul( gsl_vector_complex_get(tmpvecc, icorr),
 182                                                    gsl_complex_conjugate( gsl_vector_complex_get(tmpvecc, jcorr) )
 183                                                    );
 184           discopart = gsl_complex_mul_real(discopart, pow(1.0/(O815->comargs.nmeas-1),2));
 185           gsl_matrix_complex_set( tmpmatrix, icorr, jcorr,
 186                                   gsl_complex_sub( gsl_matrix_complex_get( tmpmatrix, icorr, jcorr ),
 187                                                    discopart
 188                                                    )
 189                                   );
 190         }
 191
 192       gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, gsl_complex_rect(1.0,0.0),
 193                       invC0, tmpmatrix, gsl_complex_rect(0.0,0.0), tmpmatrix2);
 194
 195       cdiag(tmpvec, tmpmatrix2);
 196       gsl_vector_sub(tmpvec, jackres);
 197       gsl_vector_memcpy(tmpvec2, tmpvec);
 198       gsl_vector_mul(tmpvec, tmpvec2);
 199       gsl_vector_add(jackerrornorm, tmpvec);
 200     }
 201     gsl_vector_scale( jackerrornorm,
 202                       (double)(O815->comargs.nmeas-1)/O815->comargs.nmeas );
 203
 204     *out << O815->paraQ->getParaVals();
 205     *out << "\t" << itsep;
 206     for (int iev = 0; iev < 4; iev++) {
 207       *out << "\t" << gsl_vector_get(jackres,iev)
 208            << "\t" << sqrt(gsl_vector_get(jackerrornorm,iev));
 209     }
 210     *out << endl;
 211   }
 212
 213   gsl_matrix_complex_free(totalval);
 214   gsl_matrix_complex_free(tmpmatrix);
 215   gsl_vector_free(tmpvec);
 216   gsl_vector_complex_free(tmpvecc);
 217   gsl_vector_free(tmpvec2);
 218   gsl_vector_free(jackres);
 219   gsl_vector_free(jackerrornorm);
 220   gsl_vector_complex_free(totaldisco);
 221 };
 222
 223 void obs_diagcorr::corrCompute()
 224 {
 225   complex<double> phislice[4][O815->comargs.lsize[1]];
 226
 227   for (int icorr=0; icorr<4; icorr++)
 228     OM[16][icorr] = 0;
 229
 230   for (int it = 0; it < O815->comargs.lsize[1]; it++) {
 231     for (int icorr=0; icorr<4; icorr++)
 232       phislice[icorr][it] = 0;
 233
 234     for (int ix = 0; ix < spatialV; ix++) {
 235       phislice[0][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ];
 236       phislice[1][it] += Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ] * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
 237       phislice[2][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]);
 238       phislice[3][it] += conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
 239     }
 240
 241     for (int icorr=0; icorr<4; icorr++) {
 242       phislice[icorr][it] /= spatialV;
 243       OM[16][icorr] += phislice[icorr][it];
 244     }
 245   }
 246
 247   for (int icorr = 0; icorr < 4; icorr++) {
 248     for (int jcorr = 0; jcorr < 4; jcorr++) {
 249       for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
 250         OM[icorr*4+jcorr][itsep] = 0;
 251
 252         for (int it = 0; it < O815->comargs.lsize[1]; it++)
 253           OM[icorr*4+jcorr][itsep] += phislice[icorr][ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[jcorr][it] );
 254
 255         OM[icorr*4+jcorr][itsep] /= O815->comargs.lsize[1];
 256       }
 257     }
 258
 259     OM[16][icorr] /= O815->comargs.lsize[1];
 260   }
 261 }
 262
 263 #endif