about summary refs log tree commit diff
path: root/test
diff options
context:
space:
mode:
Diffstat (limited to 'test')
-rwxr-xr-xtest/dev_test_suite.sh48
-rw-r--r--test/src/unittests-math.cpp117
-rwxr-xr-xtest/test_suite.sh6
3 files changed, 165 insertions, 6 deletions
diff --git a/test/dev_test_suite.sh b/test/dev_test_suite.sh
index 0fc4423..284c9aa 100755
--- a/test/dev_test_suite.sh
+++ b/test/dev_test_suite.sh
@@ -11,7 +11,8 @@ testBXDStandardRelatednessMatrixKSingularError() {
            -c ../example/BXD_covariates.txt \
            -a ../example/BXD_snps.txt \
            -gk \
-           -debug -o $outn
+           -debug \
+           -o $outn
     assertEquals 22 $? # should show singular error
 }
 
@@ -44,7 +45,7 @@ testBXDLMMLikelihoodRatio() {
     assertEquals 0 $?
 
     outfn=output/$outn.assoc.txt
-    assertEquals "80498" `wc -w < $outfn`
+    assertEquals "73180" `wc -w < $outfn`
     assertEquals "3088458212.93" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
 }
 
@@ -83,6 +84,49 @@ testUnivariateLinearMixedModelLOCO1() {
     assertEquals "15465346.22" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
 }
 
+testPlinkCenteredRelatednessMatrixKLOCO1() {
+    return 0
+    outn=mouse_hs1940_Plink_LOCO1
+    rm -f output/$outn.*
+    $gemma -bfile ../example/mouse_hs1940 \
+           -a ../example/mouse_hs1940.anno.txt \
+           -snps ../example/mouse_hs1940_snps.txt \
+           -nind 400 \
+           -loco 1 \
+           -gk \
+           -debug \
+           -o $outn
+    assertEquals 0 $?
+    grep "total computation time" < output/$outn.log.txt
+    outfn=output/$outn.cXX.txt
+    assertEquals 0 $?
+    assertEquals "400" `wc -l < $outfn`
+    assertEquals "0.312" `head -c 5 $outfn`
+    assertEquals "71.03" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
+}
+
+
+testPlinkUnivariateLinearMixedModelLOCO1() {
+    return 0
+    outn=mouse_hs1940_CD8_Plink_LOCO1_lmm
+    rm -f output/$outn.*
+    $gemma -bfile ../example/mouse_hs1940 \
+	   -n 1 \
+	   -loco 1 \
+           -k ./output/mouse_hs1940_Plink_LOCO1.cXX.txt \
+           -a ../example/mouse_hs1940.anno.txt \
+	   -snps ../example/mouse_hs1940_snps.txt -lmm \
+	   -nind 400 \
+	   -debug \
+           -o $outn
+    assertEquals 0 $?
+    grep "total computation time" < output/$outn.log.txt
+    assertEquals 0 $?
+    outfn=output/$outn.assoc.txt
+    assertEquals "68" `wc -l < $outfn`
+    assertEquals "15465346.22" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
+}
+
 shunit2=`which shunit2`
 
 if [ -x "$shunit2" ]; then
diff --git a/test/src/unittests-math.cpp b/test/src/unittests-math.cpp
index ac4c180..757c2dc 100644
--- a/test/src/unittests-math.cpp
+++ b/test/src/unittests-math.cpp
@@ -1,14 +1,23 @@
 #include <catch.hpp>
 #include <iostream>
 #include "gsl/gsl_matrix.h"
-#include "mathfunc.h"
+#include <cblas.h>
+
 #include <algorithm>
 #include <limits>
 #include <numeric>
 
+#include "debug.h"
+#include "mathfunc.h"
+#include "fastblas.h"
+#include "fastopenblas.h"
+
 using namespace std;
 
 TEST_CASE( "Math functions", "[math]" ) {
+  debug_set_debug_mode(true);
+  debug_set_no_check_mode(false);
+  debug_set_strict_mode(true);
   double data[] = { 2,-1, 0,
                    -1, 2,-1,
                     0,-1, 2};
@@ -51,3 +60,109 @@ TEST_CASE( "Math functions", "[math]" ) {
   REQUIRE (std::isnan(v3[2]));
   REQUIRE(has_nan(v3));
 }
+
+TEST_CASE("cblas_dgemm", "[math]") {
+   double *A, *B, *C;
+   int m, n, k, i, j;
+   double alpha, beta;
+
+   printf ("\n This example computes real matrix C=alpha*A*B+beta*C using \n"
+           " Intel(R) MKL function dgemm, where A, B, and  C are matrices and \n"
+           " alpha and beta are double precision scalars\n\n");
+
+   m = 2000, k = 200, n = 1000;
+   printf (" Initializing data for matrix multiplication C=A*B for matrix \n"
+           " A(%ix%i) and matrix B(%ix%i)\n\n", m, k, k, n);
+   alpha = 1.0; beta = 0.0;
+
+   printf (" Allocating memory for matrices aligned on 64-byte boundary for better \n"
+           " performance \n\n");
+   A = (double *)malloc( m*k*sizeof( double ));
+   B = (double *)malloc( k*n*sizeof( double ));
+   C = (double *)malloc( m*n*sizeof( double ));
+
+   printf (" Intializing matrix data \n\n");
+   for (i = 0; i < (m*k); i++) {
+     A[i] = (double)(i+1);
+   }
+
+   for (i = 0; i < (k*n); i++) {
+     B[i] = (double)(-i-1);
+   }
+
+   for (i = 0; i < (m*n); i++) {
+     C[i] = 0.0;
+   }
+
+   printf (" Computing matrix product using Intel(R) MKL dgemm function via CBLAS interface \n\n");
+   assert(m==2000);
+   assert(k==200);
+   assert(n==1000);
+   //cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+   //            m, n, k, alpha, A, k, B, n, beta, C, n);
+   fast_cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+                    m, n, k, alpha, A, k, B, n, beta, C, n);
+
+   REQUIRE(trunc(C[0]) == -2666620100.0 );
+   REQUIRE(trunc(C[1]) == -2666640200.0 );
+   REQUIRE(trunc(C[2003]) == -10627000400.0 );
+
+}
+
+TEST_CASE("fast_dgemm", "[math]") {
+   double *A, *B, *C;
+   int m, n, k, i, j;
+   double alpha, beta;
+
+   printf ("\n This example computes real matrix C=alpha*A*B+beta*C using \n"
+           " Intel(R) MKL function dgemm, where A, B, and  C are matrices and \n"
+           " alpha and beta are double precision scalars\n\n");
+
+   m = 2000, k = 200, n = 1000;
+   printf (" Initializing data for matrix multiplication C=A*B for matrix \n"
+           " A(%ix%i) and matrix B(%ix%i)\n\n", m, k, k, n);
+   alpha = 1.0; beta = 0.0;
+
+   printf (" Allocating memory for matrices aligned on 64-byte boundary for better \n"
+           " performance \n\n");
+   A = (double *)malloc( m*k*sizeof( double ));
+   B = (double *)malloc( k*n*sizeof( double ));
+   C = (double *)malloc( m*n*sizeof( double ));
+
+   printf (" Intializing matrix data \n\n");
+   for (i = 0; i < (m*k); i++) {
+     A[i] = (double)(i+1);
+   }
+
+   for (i = 0; i < (k*n); i++) {
+     B[i] = (double)(-i-1);
+   }
+
+   for (i = 0; i < (m*n); i++) {
+     C[i] = 0.0;
+   }
+
+   printf (" Computing matrix product using Intel(R) MKL dgemm function via CBLAS interface \n\n");
+   // cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+   //            m, n, k, alpha, A, k, B, n, beta, C, n);
+   // eigenlib_dgemm(const char *TransA, const char *TransB, const double alpha,
+   //                const gsl_matrix *A, const gsl_matrix *B, const double beta,
+   //                gsl_matrix *C) {
+   gsl_matrix *AM = gsl_matrix_safe_alloc(m,k); // rows x cols
+   gsl_matrix *BM = gsl_matrix_safe_alloc(k,n);
+   gsl_matrix *CM = gsl_matrix_calloc(m,n);
+
+   fast_copy(AM,A);
+   fast_copy(BM,B);
+   fast_copy(CM,C);
+   fast_dgemm("N","N",alpha,AM,BM,beta,CM);
+   printf ("\n Computations completed.\n\n");
+   A = AM->data;
+   B = BM->data;
+   C = CM->data;
+
+   REQUIRE(trunc(C[0]) == -2666620100.0 );
+   REQUIRE(trunc(C[1]) == -2666640200.0 );
+   REQUIRE(trunc(C[2003]) == -10627000400.0 );
+
+}
diff --git a/test/test_suite.sh b/test/test_suite.sh
index 350fc27..dc6053a 100755
--- a/test/test_suite.sh
+++ b/test/test_suite.sh
@@ -62,7 +62,7 @@ testUnivariateLinearMixedModel() {
     assertEquals "4038540440.86" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn`
 }
 
-testMultivariateLinearMixedModel() {
+testLinearMixedModelPhenotypes() {
     $gemma -g ../example/mouse_hs1940.geno.txt.gz \
            -p ../example/mouse_hs1940.pheno.txt \
            -n 1 6 \
@@ -92,8 +92,8 @@ testPlinkStandardRelatednessMatrixK() {
 
 # Test for https://github.com/genetics-statistics/GEMMA/issues/58
 # fixed GSLv2 NaN's that appeared with covariates.
-testPlinkMultivariateLinearMixedModel() {
-    testname=testPlinkMultivariateLinearMixedModel
+testPlinkLinearMixedModelCovariates() {
+    testname=testPlinkLinearMixedModelCovariates
     datadir=../example
     $gemma -bfile $datadir/HLC \
            -k output/testPlinkStandardRelatednessMatrixK.sXX.txt \