gsl_vector_complex **measurements_disco;
static void cdiag (gsl_vector *v, gsl_matrix_complex *m);
-
- static void invertcm (gsl_matrix_complex *invm, const gsl_matrix_complex *m);
};
obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr",
}
};
-void obs_diagcorr::invertcm (gsl_matrix_complex *invm, const gsl_matrix_complex *m)
-{
- int s;
- gsl_permutation * perm = gsl_permutation_alloc (m->size1);
-
- gsl_matrix_complex *mtmp = gsl_matrix_complex_alloc(m->size1,m->size2);
- gsl_matrix_complex_memcpy(mtmp, m);
-
- gsl_linalg_complex_LU_decomp (mtmp, perm, &s);
- gsl_linalg_complex_LU_invert (mtmp, perm, invm);
-
- gsl_matrix_complex_free(mtmp);
-}
-
void obs_diagcorr::cdiag (gsl_vector *v, gsl_matrix_complex *m)
{
gsl_matrix_complex *evec = gsl_matrix_complex_alloc(4,4);
void obs_diagcorr::_finish() {
gsl_matrix_complex *totalval = gsl_matrix_complex_alloc(4,4);
gsl_matrix_complex *tmpmatrix = gsl_matrix_complex_alloc(4,4);
- gsl_matrix_complex *tmpmatrix2 = gsl_matrix_complex_alloc(4,4);
gsl_vector *tmpvec = gsl_vector_alloc(4);
gsl_vector_complex *tmpvecc = gsl_vector_complex_alloc(4);
gsl_vector *tmpvec2 = gsl_vector_alloc(4);
gsl_vector *jackres = gsl_vector_alloc(4);
gsl_vector *jackerrornorm = gsl_vector_alloc(4);
gsl_vector_complex *totaldisco = gsl_vector_complex_alloc(4);
- gsl_matrix_complex *invC0 = NULL;
gsl_vector_complex_set_zero(totaldisco);
for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
)
);
}
-
- if (invC0 == NULL) {
- invC0 = gsl_matrix_complex_alloc(4,4);
- invertcm(invC0, tmpmatrix);
- }
-
- gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, gsl_complex_rect(1.0,0.0),
- invC0, tmpmatrix, gsl_complex_rect(0.0,0.0), tmpmatrix2);
- cdiag(jackres, tmpmatrix2);
+ cdiag(jackres, tmpmatrix);
for (int imeas=0; imeas<O815->comargs.nmeas; imeas++) {
gsl_matrix_complex_memcpy (tmpmatrix, totalval);
discopart
)
);
- }
-
- gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, gsl_complex_rect(1.0,0.0),
- invC0, tmpmatrix, gsl_complex_rect(0.0,0.0), tmpmatrix2);
-
- cdiag(tmpvec, tmpmatrix2);
+ }
+ cdiag(tmpvec, tmpmatrix);
gsl_vector_sub(tmpvec, jackres);
gsl_vector_memcpy(tmpvec2, tmpvec);
gsl_vector_mul(tmpvec, tmpvec2);