diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9061aec..5e154ad 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,3 +1,5 @@
+cmake_minimum_required(VERSION 3.20)
+
 set(LinAlgebra_Sources
     Linalg.h
     LUDecomp.cpp
@@ -6,18 +8,29 @@ set(LinAlgebra_Sources
     Mtx.h
     Range.h
     Vect.cpp
-    Vect.h)
+    Vect.h
+)
 
-#set(LinAlgebra_Required_Libs
-#   debug MSVCRTD.LIB
-#   debug MSVCPRTD.LIB
-#   optimized MSVCRT.LIB
-#   optimized MSVCPRT.LIB)
+set(LinAlgebra_Required_Libs
+)
 
+if (APPLE)
+    find_library(ACCELERATE_FRAMEWORK Accelerate)
+    if (ACCELERATE_FRAMEWORK)
+        message(STATUS "Accelerate framework found")
+        add_compile_definitions(LINALG_USE_ACCELERATE ACCELERATE_LAPACK_ILP64=1)
+        add_compile_options(-DACCELERATE_NEW_LAPACK)
+    else()
+        message(WARNING "Accelerate framework not found")
+    endif()
+endif()
 
-add_library(LinAlgebra ${LinAlgebra_Sources})
-# libraries that needs to be linked to the present library
-#target_link_libraries(LinAlgebra ${LinAlgebra_Required_Libs})
+find_package(BLAS REQUIRED)
+include_directories(${BLAS_INCLUDE_DIRS})
+set(LinAlgebra_Required_Libs ${LinAlgebra_Required_Libs} ${BLAS_LIBRARIES})
 
+add_library(LinAlgebra ${LinAlgebra_Sources})
 
+# libraries that needs to be linked to the present library
+target_link_libraries(LinAlgebra ${LinAlgebra_Required_Libs})
 
diff --git a/Linalg.h b/Linalg.h
index b3bfdeb..aa0d2aa 100644
--- a/Linalg.h
+++ b/Linalg.h
@@ -74,6 +74,11 @@ inline void CLin_assert(bool a) { _ASSERT(a); }
 inline void CLin_assert(bool ) { }
 #endif
 
+#if defined(LINALG_USE_ACCELERATE)
+#   include <Accelerate/Accelerate.h>
+#else
+#   include <cblas.h>
+#endif
 
 #endif // LINALG_H
 
diff --git a/Mtx.h b/Mtx.h
index b228c87..43caf12 100644
--- a/Mtx.h
+++ b/Mtx.h
@@ -88,6 +88,9 @@ class CLin_Matrix
 
 	inline CLin_subscript size() const { return m_iMN; }
 
+    inline double *array() { return m_pV; }
+
+    inline const double *array() const { return m_pV; }
 	// returns # of rows or cols according to parameter d (1=rows, 2=cols, other=0)
 	inline CLin_subscript dim(CLin_subscript d)
 	{
@@ -365,90 +368,94 @@ inline CLin_Vector matmult(const CLin_Matrix  &A, const CLin_Vector &x)
 #endif
 
 	CLin_subscript M, N;
-	CLin_subscript i, j;
 
 	M = A.num_rows();
 	N = A.num_cols();
 
 	CLin_Vector tmp(M);
-	double sum;
-
-	for (i=0; i<M; i++)
-	{
-		sum = 0;
-		const double* rowi = A[i];
-		for (j=0; j<N; j++)
-			sum = sum +  rowi[j] * x[j];
-
-		tmp[i] = sum;
-	}
+    
+    // Y <- alpha*A*X + beta*Y
+    cblas_dgemv(CblasRowMajor,
+                CblasNoTrans,
+                M,
+                N,
+                1.0, // Scaling factor for the product of matrix A and vector X
+                A.array(),
+                M,
+                x.array(),
+                1,
+                0.0, // Scaling factor for vector Y
+                tmp.array(),
+                1
+                );
 
 	return tmp;
 }
 
-inline CLin_Matrix matmult(const CLin_Matrix &A, const CLin_Matrix &B)
-{
-
-#ifdef LINALG_BOUNDS_CHECK
-	CLin_assert(A.num_cols() == B.num_rows());
-#endif
-
+inline CLin_Matrix matmult(const CLin_Matrix &A, const CLin_Matrix &B)
+{
+
+#ifdef LINALG_BOUNDS_CHECK
+	CLin_assert(A.num_cols() == B.num_rows());
+#endif
+
 	CLin_subscript M, N, K;
-	CLin_subscript i, j, k;
-	double sum;
 
 	M = A.num_rows();
-	N = A.num_cols();
-	K = B.num_cols();
-
-	CLin_Matrix tmp(M,K);
-
-	for (i=0; i<M; i++)
-		for (k=0; k<K; k++)
-		{
-			sum = 0;
-			for (j=0; j<N; j++)
-				sum = sum +  A[i][j] * B[j][k];
-
-			tmp[i][k] = sum;
-		}
-
-	return tmp;
+    N = B.num_cols();
+	K = A.num_cols();
+
+	CLin_Matrix tmp(M, N);
+
+    // C <- alpha*A*B + beta*C
+    cblas_dgemm(CblasRowMajor,
+                CblasNoTrans,
+                CblasNoTrans,
+                M,
+                N,
+                K,
+                1.0,      // ALPHA, scaling factor for the product of matrices A and B.
+                A.array(),
+                M,
+                B.array(),
+                N,
+                0.0,      // BETA, Scaling factor for matrix C
+                tmp.array(),
+                M);
+
+	return tmp;
 }
-
-inline int matmult(CLin_Matrix &C, const CLin_Matrix &A, const CLin_Matrix &B)
-{
+
+inline int matmult(CLin_Matrix &C, const CLin_Matrix &A, const CLin_Matrix &B)
+{
 	CLin_subscript M, N, K;
-	CLin_subscript i, j, k;
-	double sum;
-	const double* row_i;
-	const double* col_k;
 
 	CLin_assert(A.num_cols() == B.num_rows());
 
 	M = A.num_rows();
-	N = A.num_cols();
-	K = B.num_cols();
-
-	C.newsize(M,K);
-
-	for (i=0; i<M; i++)
-		for (k=0; k<K; k++)
-		{
-			row_i  = &(A[i][0]);
-			col_k  = &(B[0][k]);
-			sum = 0;
-			for (j=0; j<N; j++)
-			{
-				sum  += *row_i * *col_k;
-				row_i++;
-				col_k += K;
-			}
-			C[i][k] = sum;
-		}
-
-	return 0;
-}
+    N = B.num_cols();
+	K = A.num_cols();
+
+	C.newsize(M, N);
+
+    // C <- alpha*A*B + beta*C
+    cblas_dgemm(CblasRowMajor,
+                CblasNoTrans,
+                CblasNoTrans,
+                M,
+                N,
+                K,
+                1.0,      // ALPHA, scaling factor for the product of matrices A and B.
+                A.array(),
+                M,
+                B.array(),
+                N,
+                0.0,      // BETA, Scaling factor for matrix C
+                C.array(),
+                M);
+
+	return 0;
+}
 
 inline CLin_Vector operator*(const CLin_Matrix	&A, const CLin_Vector &x)
 {
diff --git a/Vect.h b/Vect.h
index 37b9e4d..69ae8a5 100644
--- a/Vect.h
+++ b/Vect.h
@@ -181,6 +181,7 @@ class CLin_Vector
     inline CLin_subscript dim() const { return  m_iN; }
     inline CLin_subscript size() const { return  m_iN; }
     inline double* array() { return m_dV; }
+    inline const double* array() const { return m_dV; }
 
 	//
 	// methods
@@ -417,13 +418,15 @@ inline double dot_prod(const CLin_Vector &A, const CLin_Vector &B)
     CLin_subscript N = A.dim();
     CLin_assert(N == B.dim());
 
-    CLin_subscript i;
-    double sum = 0;
-
-    for (i=0; i<N; i++)
-		sum += A[i] * B[i];
-
-    return sum;
+    double result = cblas_ddot(
+        N,        // The number of elements in the vectors
+        A.array(), // Vector X
+        1,        // Stride within X. For example, if incX is 7, every 7th element is used
+        B.array(), // Vector Y
+        1         // Stride within Y. For example, if incY is 7, every 7th element is used
+    );
+    
+    return result;
 }
 
 inline double mod(const CLin_Vector &A)