diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/gemma.cpp | 13 | ||||
| -rw-r--r-- | src/gemma_io.cpp | 2 | ||||
| -rw-r--r-- | src/lapack.cpp | 14 | ||||
| -rw-r--r-- | src/lmm.cpp | 9 |
4 files changed, 27 insertions, 11 deletions
diff --git a/src/gemma.cpp b/src/gemma.cpp index a183edc..78633e8 100644 --- a/src/gemma.cpp +++ b/src/gemma.cpp @@ -1624,7 +1624,7 @@ void GEMMA::BatchRun(PARAM &cPar) { clock_t time_begin, time_start; time_begin = clock(); - if (is_check_mode()) enable_segfpe(); // fast NaN checking + if (is_check_mode()) enable_segfpe(); // fast NaN checking by default // Read Files. cout << "Reading Files ... " << endl; @@ -1912,7 +1912,7 @@ void GEMMA::BatchRun(PARAM &cPar) { VARCOV cVarcov; cVarcov.CopyFromParam(cPar); - if (is_check_mode()) disable_segfpe(); // fast NaN checking + if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now if (!cPar.file_bfile.empty()) { cVarcov.AnalyzePlink(); @@ -2027,7 +2027,7 @@ void GEMMA::BatchRun(PARAM &cPar) { VARCOV cVarcov; cVarcov.CopyFromParam(cPar); - if (is_check_mode()) disable_segfpe(); // fast NaN checking + if (is_check_mode()) disable_segfpe(); // fast NaN checking for now if (!cPar.file_bfile.empty()) { cVarcov.AnalyzePlink(); @@ -2055,7 +2055,7 @@ void GEMMA::BatchRun(PARAM &cPar) { gsl_vector_view Y_col = gsl_matrix_column(Y, 0); - if (is_check_mode()) disable_segfpe(); // fast NaN checking + // if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now if (!cPar.file_gene.empty()) { cLm.AnalyzeGene(W, @@ -2763,7 +2763,7 @@ void GEMMA::BatchRun(PARAM &cPar) { LMM cLmm; cLmm.CopyFromParam(cPar); - if (is_check_mode()) disable_segfpe(); // fast NaN checking + if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now gsl_vector_view Y_col = gsl_matrix_column(Y, 0); gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0); @@ -2790,14 +2790,13 @@ void GEMMA::BatchRun(PARAM &cPar) { &Y_col.vector, env); } } - cLmm.WriteFiles(); cLmm.CopyToParam(cPar); } else { MVLMM cMvlmm; cMvlmm.CopyFromParam(cPar); - if (is_check_mode()) disable_segfpe(); // fast NaN checking + if (is_check_mode()) disable_segfpe(); // disable fast NaN checking if (!cPar.file_bfile.empty()) { if (cPar.file_gxe.empty()) { diff --git a/src/gemma_io.cpp b/src/gemma_io.cpp index 818c5e8..fad3e47 100644 --- a/src/gemma_io.cpp +++ b/src/gemma_io.cpp @@ -1483,12 +1483,14 @@ bool BimbamKin(const string file_geno, const set<string> ksnps, // and transpose // FIXME: the following is very slow + debug_msg("begin transpose"); for (size_t i = 0; i < ni_total; ++i) { for (size_t j = 0; j < i; ++j) { d = gsl_matrix_get(matrix_kin, j, i); gsl_matrix_set(matrix_kin, i, j, d); } } + debug_msg("end transpose"); // GSL is faster - and there are even faster methods // enforce_gsl(gsl_matrix_transpose(matrix_kin)); diff --git a/src/lapack.cpp b/src/lapack.cpp index 187daa8..b309a1b 100644 --- a/src/lapack.cpp +++ b/src/lapack.cpp @@ -193,6 +193,8 @@ void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, double WORK_temp[1]; int IWORK_temp[1]; + if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now + // DSYEVR - computes selected eigenvalues and, optionally, // eigenvectors of a real symmetric matrix // Here compute both (JOBZ is V), all eigenvalues (RANGE is A) @@ -200,7 +202,12 @@ void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK_temp, &LWORK, IWORK_temp, &LIWORK, &INFO); - enforce_msg(INFO != 0, "lapack_eigen_symmv failed"); + // If info = 0, the execution is successful. + // If info = -i, the i-th parameter had an illegal value. + // If info = i, an internal error has occurred. + + if (INFO != 0) cerr << "ERROR: value of INFO is " << INFO; + enforce_msg(INFO == 0, "lapack_eigen_symmv failed"); LWORK = (int)WORK_temp[0]; LIWORK = (int)IWORK_temp[0]; @@ -210,7 +217,10 @@ void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK, &LWORK, IWORK, &LIWORK, &INFO); - enforce_msg(INFO != 0, "lapack_eigen_symmv failed"); + if (INFO != 0) cerr << "ERROR: value of INFO is " << INFO; + enforce_msg(INFO == 0, "lapack_eigen_symmv failed"); + + if (is_check_mode()) disable_segfpe(); // reinstate fast NaN checking gsl_matrix_transpose(evec); diff --git a/src/lmm.cpp b/src/lmm.cpp index acd9667..7e4c3e1 100644 --- a/src/lmm.cpp +++ b/src/lmm.cpp @@ -97,8 +97,9 @@ void LMM::CopyToParam(PARAM &cPar) { } void LMM::WriteFiles() { - string file_str; + debug_msg("LMM::WriteFiles"); + file_str = path_out + "/" + file_out; file_str += ".assoc.txt"; @@ -1319,6 +1320,7 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp, auto batch_compute = [&](size_t l) { // using a C++ closure // Compute SNPs in batch, note the computations are independent per SNP + debug_msg("enter batch_compute"); gsl_matrix_view Xlarge_sub = gsl_matrix_submatrix(Xlarge, 0, 0, inds, l); gsl_matrix_view UtXlarge_sub = gsl_matrix_submatrix(UtXlarge, 0, 0, inds, l); @@ -1366,6 +1368,7 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp, p_wald, p_lrt, p_score, logl_H1}; sumStat.push_back(SNPs); } + debug_msg("exit batch_compute"); }; const auto num_snps = indicator_snp.size(); @@ -1441,9 +1444,11 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp, gsl_vector_safe_memcpy(&Xlarge_col.vector, x); c++; // count SNPs going in - if (c % msize == 0) + if (c % msize == 0) { batch_compute(msize); + } } + batch_compute(c % msize); ProgressBar("Reading SNPs", num_snps - 1, num_snps - 1); // cout << "Counted SNPs " << c << " sumStat " << sumStat.size() << endl; |