diff --git a/Tests/MathLib/TestGlobalMatrixInterface.cpp b/Tests/MathLib/TestGlobalMatrixInterface.cpp
index 71c17d97c1e995bed2d3f7794f537c7f4c959ccf..f157984c4b2c2d506027e89ed014018bb4e34668 100644
--- a/Tests/MathLib/TestGlobalMatrixInterface.cpp
+++ b/Tests/MathLib/TestGlobalMatrixInterface.cpp
@@ -15,6 +15,8 @@
 
 #include <gtest/gtest.h>
 
+#include "MathLib/LinAlg/BLAS.h"
+
 #if defined(USE_LIS)
 #include "MathLib/LinAlg/Lis/LisMatrix.h"
 #elif defined(USE_PETSC)
@@ -30,6 +32,8 @@
 
 #include "NumLib/NumericsConfig.h"
 
+using namespace MathLib::BLAS;
+
 namespace
 {
 
@@ -94,12 +98,25 @@ void checkGlobalMatrixInterfaceMPI(T_MATRIX &m, T_VECTOR &v)
 
     MathLib::finalizeMatrixAssembly(m);
 
+    // Multiply by a vector
+    // v = 1.;
+    set(v, 1.);
+    const bool deep_copy = false;
+    T_VECTOR y(v, deep_copy);
+    matMult(m, v, y);
+
+    ASSERT_EQ(sqrt(3*(3*3 + 7*7)), norm2(y));
+
     // set a value
     m.set(2 * mrank, 2 * mrank, 5.0);
     MathLib::finalizeMatrixAssembly(m);
     // add a value
     m.add(2 * mrank+1, 2 * mrank+1, 5.0);
     MathLib::finalizeMatrixAssembly(m);
+
+    matMult(m, v, y);
+
+    ASSERT_EQ(sqrt((3*7*7 + 3*12*12)), norm2(y));
 }
 
 // Rectanglular matrix
@@ -142,6 +159,12 @@ void checkGlobalRectangularMatrixInterfaceMPI(T_MATRIX &m, T_VECTOR &v)
 
     MathLib::finalizeMatrixAssembly(m);
 
+    // Multiply by a vector
+    set(v, 1);
+    T_VECTOR y(m.getNRows());
+    matMult(m, v, y);
+
+    ASSERT_NEAR(6.*sqrt(6.), norm2(y), 1.e-10);
 }
 
 #endif // end of: ifdef USE_PETSC // or MPI
diff --git a/Tests/MathLib/TestGlobalVectorInterface.cpp b/Tests/MathLib/TestGlobalVectorInterface.cpp
index 0a63772956d4057a6debb9490eadb966de4509a0..e51b04eef2d8b0e0261f172f0da5bb6297feb318 100644
--- a/Tests/MathLib/TestGlobalVectorInterface.cpp
+++ b/Tests/MathLib/TestGlobalVectorInterface.cpp
@@ -16,6 +16,8 @@
 #include <gtest/gtest.h>
 #include "../TestTools.h"
 
+#include "MathLib/LinAlg/BLAS.h"
+
 #if defined(USE_LIS)
 #include "MathLib/LinAlg/Lis/LisVector.h"
 #elif defined(USE_PETSC)
@@ -29,6 +31,8 @@
 #include "MathLib/LinAlg/FinalizeVectorAssembly.h"
 #include "NumLib/NumericsConfig.h"
 
+using namespace MathLib::BLAS;
+
 namespace
 {
 
@@ -50,15 +54,27 @@ void checkGlobalVectorInterface()
     x.add(0, 1.0);
     ASSERT_EQ(2.0, x.get(0));
 
-    std::vector<double> local_vec(2, 1.0);
-    std::vector<GlobalIndexType> vec_pos(2);
-    vec_pos[0] = 0;
-    vec_pos[1] = 3;
-
     T_VECTOR y(x);
     ASSERT_EQ(2.0, y.get(0));
     ASSERT_EQ(0.0, y.get(1));
+    // y += x;
+    axpy(y, 1., x);
 
+    ASSERT_EQ(4.0, y.get(0));
+    //y -= x;
+    axpy(y, -1., x);
+    ASSERT_EQ(2.0, y.get(0));
+    //y = 1.0;
+    set(y, 1.0);
+    ASSERT_EQ(1.0, y.get(0));
+    // y = x;
+    copy(x, y);
+    ASSERT_EQ(2.0, y.get(0));
+
+    std::vector<double> local_vec(2, 1.0);
+    std::vector<GlobalIndexType> vec_pos(2);
+    vec_pos[0] = 0;
+    vec_pos[1] = 3;
     y.add(vec_pos, local_vec);
     ASSERT_EQ(3.0, y.get(0));
     ASSERT_EQ(0.0, y.get(1));
@@ -85,11 +101,27 @@ void checkGlobalVectorInterfacePETSc()
     //x.get(0) is expensive, only get local value. Use it for test purpose
     ASSERT_EQ(.0, x.get(r0));
 
+    set(x, 10.);
+
+
     // Value of x is not copied to y
     const bool deep_copy = false;
     T_VECTOR y(x, deep_copy);
     ASSERT_EQ(0, y.get(r0));
 
+    set(y, 10.0);
+    ASSERT_EQ(10, y.get(r0));
+
+    // y += x
+    axpy(y, 1., x);
+    ASSERT_EQ(20, y.get(r0));
+    ASSERT_EQ(80., norm2(y));
+
+    // y -= x
+    axpy(y, -1., x);
+    ASSERT_EQ(10, y.get(r0));
+    ASSERT_EQ(40., norm2(y));
+
     std::vector<double> local_vec(2, 10.0);
     std::vector<GlobalIndexType> vec_pos(2);
 
@@ -98,6 +130,10 @@ void checkGlobalVectorInterfacePETSc()
 
     y.add(vec_pos, local_vec);
 
+    double normy = std::sqrt(6.0*400+10.0*100);
+
+    ASSERT_NEAR(0.0, normy-norm2(y), 1.e-10);
+
     double x0[16];
     double z[] =
     {
@@ -163,6 +199,7 @@ void checkGlobalVectorInterfacePETSc()
     // Deep copy
     MathLib::finalizeVectorAssembly(x_fixed_p);
     T_VECTOR x_deep_copied(x_fixed_p);
+    ASSERT_NEAR(sqrt(3.0*5), norm2(x_deep_copied), 1.e-10);
 
     // -----------------------------------------------------------------
     // Vector with ghost entries