diff --git a/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/vectorsCodegen.kt b/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/vectorsCodegen.kt
index 018d84773..e227dc721 100644
--- a/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/vectorsCodegen.kt
+++ b/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/vectorsCodegen.kt
@@ -18,8 +18,8 @@ private fun KtPsiFactory.createVectorClass(
     val structName = sn("gsl_vectorR", cTypeName)
 
     @Language("kotlin") val text =
-        """internal class $className(override val nativeHandle: CPointer<$structName>, scope: DeferScope) 
-            : GslVector<$kotlinTypeName, $structName>(scope) {
+        """internal class $className(override val nativeHandle: CPointer<$structName>, scope: DeferScope) : 
+    GslVector<$kotlinTypeName, $structName>(scope) {
     override val size: Int
         get() = nativeHandle.pointed.size.toInt()
 
diff --git a/kmath-gsl/build.gradle.kts b/kmath-gsl/build.gradle.kts
index df1e635be..d7d53e797 100644
--- a/kmath-gsl/build.gradle.kts
+++ b/kmath-gsl/build.gradle.kts
@@ -12,7 +12,7 @@ kotlin {
 
     val nativeTarget = when (System.getProperty("os.name")) {
         "Mac OS X" -> macosX64()
-        "Linux" -> linuxX64()
+        "Linux" -> linuxX64("native")
 
         else -> {
             logger.warn("Current OS cannot build any of kmath-gsl targets.")
@@ -29,7 +29,7 @@ kotlin {
     val test by nativeTarget.compilations.getting
 
     sourceSets {
-        val nativeMain by creating {
+        val nativeMain by getting {
             val codegen by tasks.creating {
                 matricesCodegen(kotlin.srcDirs.first().absolutePath + "/kscience/kmath/gsl/_Matrices.kt")
                 vectorsCodegen(kotlin.srcDirs.first().absolutePath + "/kscience/kmath/gsl/_Vectors.kt")
@@ -42,11 +42,11 @@ kotlin {
             }
         }
 
-        val nativeTest by creating {
+        val nativeTest by getting {
             dependsOn(nativeMain)
         }
 
-        main.defaultSourceSet.dependsOn(nativeMain)
-        test.defaultSourceSet.dependsOn(nativeTest)
+//        main.defaultSourceSet.dependsOn(nativeMain)
+//        test.defaultSourceSet.dependsOn(nativeTest)
     }
 }
diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslException.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslException.kt
new file mode 100644
index 000000000..65e466015
--- /dev/null
+++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslException.kt
@@ -0,0 +1,70 @@
+package kscience.kmath.gsl
+
+import kotlinx.cinterop.staticCFunction
+import kotlinx.cinterop.toKString
+import org.gnu.gsl.gsl_set_error_handler
+import org.gnu.gsl.gsl_set_error_handler_off
+import kotlin.native.concurrent.AtomicInt
+
+private object Container {
+    val isKmathHandlerRegistered = AtomicInt(0)
+}
+
+internal enum class GslErrnoValue(val code: Int, val text: String) {
+    GSL_SUCCESS(org.gnu.gsl.GSL_SUCCESS, ""),
+    GSL_FAILURE(org.gnu.gsl.GSL_FAILURE, ""),
+    GSL_CONTINUE(org.gnu.gsl.GSL_CONTINUE, "iteration has not converged"),
+    GSL_EDOM(org.gnu.gsl.GSL_EDOM, "input domain error, e.g sqrt(-1)"),
+    GSL_ERANGE(org.gnu.gsl.GSL_ERANGE, "output range error, e.g. exp(1e100)"),
+    GSL_EFAULT(org.gnu.gsl.GSL_EFAULT, "invalid pointer"),
+    GSL_EINVAL(org.gnu.gsl.GSL_EINVAL, "invalid argument supplied by user"),
+    GSL_EFAILED(org.gnu.gsl.GSL_EFAILED, "generic failure"),
+    GSL_EFACTOR(org.gnu.gsl.GSL_EFACTOR, "factorization failed"),
+    GSL_ESANITY(org.gnu.gsl.GSL_ESANITY, "sanity check failed - shouldn't happen"),
+    GSL_ENOMEM(org.gnu.gsl.GSL_ENOMEM, "malloc failed"),
+    GSL_EBADFUNC(org.gnu.gsl.GSL_EBADFUNC, "problem with user-supplied function"),
+    GSL_ERUNAWAY(org.gnu.gsl.GSL_ERUNAWAY, "iterative process is out of control"),
+    GSL_EMAXITER(org.gnu.gsl.GSL_EMAXITER, "exceeded max number of iterations"),
+    GSL_EZERODIV(org.gnu.gsl.GSL_EZERODIV, "tried to divide by zero"),
+    GSL_EBADTOL(org.gnu.gsl.GSL_EBADTOL, "user specified an invalid tolerance"),
+    GSL_ETOL(org.gnu.gsl.GSL_ETOL, "failed to reach the specified tolerance"),
+    GSL_EUNDRFLW(org.gnu.gsl.GSL_EUNDRFLW, "underflow"),
+    GSL_EOVRFLW(org.gnu.gsl.GSL_EOVRFLW, "overflow "),
+    GSL_ELOSS(org.gnu.gsl.GSL_ELOSS, "loss of accuracy"),
+    GSL_EROUND(org.gnu.gsl.GSL_EROUND, "failed because of roundoff error"),
+    GSL_EBADLEN(org.gnu.gsl.GSL_EBADLEN, "matrix, vector lengths are not conformant"),
+    GSL_ENOTSQR(org.gnu.gsl.GSL_ENOTSQR, "matrix not square"),
+    GSL_ESING(org.gnu.gsl.GSL_ESING, "apparent singularity detected"),
+    GSL_EDIVERGE(org.gnu.gsl.GSL_EDIVERGE, "integral or series is divergent"),
+    GSL_EUNSUP(org.gnu.gsl.GSL_EUNSUP, "requested feature is not supported by the hardware"),
+    GSL_EUNIMPL(org.gnu.gsl.GSL_EUNIMPL, "requested feature not (yet) implemented"),
+    GSL_ECACHE(org.gnu.gsl.GSL_ECACHE, "cache limit exceeded"),
+    GSL_ETABLE(org.gnu.gsl.GSL_ETABLE, "table limit exceeded"),
+    GSL_ENOPROG(org.gnu.gsl.GSL_ENOPROG, "iteration is not making progress towards solution"),
+    GSL_ENOPROGJ(org.gnu.gsl.GSL_ENOPROGJ, "jacobian evaluations are not improving the solution"),
+    GSL_ETOLF(org.gnu.gsl.GSL_ETOLF, "cannot reach the specified tolerance in F"),
+    GSL_ETOLX(org.gnu.gsl.GSL_ETOLX, "cannot reach the specified tolerance in X"),
+    GSL_ETOLG(org.gnu.gsl.GSL_ETOLG, "cannot reach the specified tolerance in gradient"),
+    GSL_EOF(org.gnu.gsl.GSL_EOF, "end of file");
+
+    override fun toString(): String = "${name}('$text')"
+
+    companion object {
+        fun valueOf(code: Int): GslErrnoValue? = values().find { it.code == code }
+    }
+}
+
+internal class GslException internal constructor(file: String, line: Int, reason: String, errno: Int) :
+    RuntimeException("$file:$line: $reason. errno - $errno, ${GslErrnoValue.valueOf(errno)}") {
+}
+
+internal fun ensureHasGslErrorHandler() {
+    if (Container.isKmathHandlerRegistered.value == 1) return
+    gsl_set_error_handler_off()
+
+    gsl_set_error_handler(staticCFunction { reason, file, line, errno ->
+        throw GslException(checkNotNull(file).toKString(), line, checkNotNull(reason).toKString(), errno)
+    })
+
+    Container.isKmathHandlerRegistered.value = 1
+}
diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContext.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContext.kt
index 8b0f12883..762cf1252 100644
--- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContext.kt
+++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContext.kt
@@ -1,7 +1,175 @@
 package kscience.kmath.gsl
 
 import kotlinx.cinterop.CStructVar
+import kotlinx.cinterop.DeferScope
+import kotlinx.cinterop.pointed
 import kscience.kmath.linear.MatrixContext
 import kscience.kmath.linear.Point
+import kscience.kmath.operations.Complex
+import kscience.kmath.operations.ComplexField
+import kscience.kmath.operations.toComplex
 import kscience.kmath.structures.Matrix
+import org.gnu.gsl.*
 
+internal inline fun <T : Any, H : CStructVar> GslMatrix<T, H>.fill(initializer: (Int, Int) -> T): GslMatrix<T, H> =
+    apply {
+        (0 until rowNum).forEach { row -> (0 until colNum).forEach { col -> this[row, col] = initializer(row, col) } }
+    }
+
+internal inline fun <T : Any, H : CStructVar> GslVector<T, H>.fill(initializer: (Int) -> T): GslVector<T, H> =
+    apply { (0 until size).forEach { index -> this[index] = initializer(index) } }
+
+public abstract class GslMatrixContext<
+        T : Any,
+        H1 : CStructVar,
+        H2 : CStructVar> internal constructor(internal val scope: DeferScope) : MatrixContext<T> {
+    init {
+        ensureHasGslErrorHandler()
+    }
+
+    @Suppress("UNCHECKED_CAST")
+    public fun Matrix<T>.toGsl(): GslMatrix<T, H1> = (if (this is GslMatrix<*, *>)
+        this as GslMatrix<T, H1>
+    else
+        produce(rowNum, colNum) { i, j -> this[i, j] }).copy()
+
+    @Suppress("UNCHECKED_CAST")
+    public fun Point<T>.toGsl(): GslVector<T, H2> =
+        (if (this is GslVector<*, *>) this as GslVector<T, H2> else produceDirtyVector(size).fill { this[it] }).copy()
+
+    internal abstract fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix<T, H1>
+    internal abstract fun produceDirtyVector(size: Int): GslVector<T, H2>
+
+    public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): GslMatrix<T, H1> =
+        produceDirtyMatrix(rows, columns).fill(initializer)
+}
+
+public class GslRealMatrixContext(scope: DeferScope) : GslMatrixContext<Double, gsl_matrix, gsl_vector>(scope) {
+    override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix<Double, gsl_matrix> =
+        GslRealMatrix(nativeHandle = requireNotNull(gsl_matrix_alloc(rows.toULong(), columns.toULong())), scope = scope)
+
+    override fun produceDirtyVector(size: Int): GslVector<Double, gsl_vector> =
+        GslRealVector(nativeHandle = requireNotNull(gsl_vector_alloc(size.toULong())), scope = scope)
+
+    public override fun Matrix<Double>.dot(other: Matrix<Double>): GslMatrix<Double, gsl_matrix> {
+        val x = toGsl().nativeHandle
+        val a = other.toGsl().nativeHandle
+        val result = requireNotNull(gsl_matrix_calloc(a.pointed.size1, a.pointed.size2))
+        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, x, a, 1.0, result)
+        return GslRealMatrix(result, scope = scope)
+    }
+
+    public override fun Matrix<Double>.dot(vector: Point<Double>): GslVector<Double, gsl_vector> {
+        val x = toGsl().nativeHandle
+        val a = vector.toGsl().nativeHandle
+        val result = requireNotNull(gsl_vector_calloc(a.pointed.size))
+        gsl_blas_dgemv(CblasNoTrans, 1.0, x, a, 1.0, result)
+        return GslRealVector(result, scope = scope)
+    }
+
+    public override fun Matrix<Double>.times(value: Double): GslMatrix<Double, gsl_matrix> {
+        val g1 = toGsl()
+        gsl_matrix_scale(g1.nativeHandle, value)
+        return g1
+    }
+
+    public override fun add(a: Matrix<Double>, b: Matrix<Double>): GslMatrix<Double, gsl_matrix> {
+        val g1 = a.toGsl()
+        gsl_matrix_add(g1.nativeHandle, b.toGsl().nativeHandle)
+        return g1
+    }
+
+    public override fun multiply(a: Matrix<Double>, k: Number): GslMatrix<Double, gsl_matrix> {
+        val g1 = a.toGsl()
+        gsl_matrix_scale(g1.nativeHandle, k.toDouble())
+        return g1
+    }
+}
+
+public class GslFloatMatrixContext(scope: DeferScope) :
+    GslMatrixContext<Float, gsl_matrix_float, gsl_vector_float>(scope) {
+    override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix<Float, gsl_matrix_float> =
+        GslFloatMatrix(requireNotNull(gsl_matrix_float_alloc(rows.toULong(), columns.toULong())), scope = scope)
+
+    override fun produceDirtyVector(size: Int): GslVector<Float, gsl_vector_float> =
+        GslFloatVector(requireNotNull(gsl_vector_float_alloc(size.toULong())), scope)
+
+    public override fun Matrix<Float>.dot(other: Matrix<Float>): GslMatrix<Float, gsl_matrix_float> {
+        val x = toGsl().nativeHandle
+        val a = other.toGsl().nativeHandle
+        val result = requireNotNull(gsl_matrix_float_calloc(a.pointed.size1, a.pointed.size2))
+        gsl_blas_sgemm(CblasNoTrans, CblasNoTrans, 1f, x, a, 1f, result)
+        return GslFloatMatrix(nativeHandle = result, scope = scope)
+    }
+
+    public override fun Matrix<Float>.dot(vector: Point<Float>): GslVector<Float, gsl_vector_float> {
+        val x = toGsl().nativeHandle
+        val a = vector.toGsl().nativeHandle
+        val result = requireNotNull(gsl_vector_float_calloc(a.pointed.size))
+        gsl_blas_sgemv(CblasNoTrans, 1f, x, a, 1f, result)
+        return GslFloatVector(nativeHandle = result, scope = scope)
+    }
+
+    public override fun Matrix<Float>.times(value: Float): GslMatrix<Float, gsl_matrix_float> {
+        val g1 = toGsl()
+        gsl_matrix_float_scale(g1.nativeHandle, value.toDouble())
+        return g1
+    }
+
+    public override fun add(a: Matrix<Float>, b: Matrix<Float>): GslMatrix<Float, gsl_matrix_float> {
+        val g1 = a.toGsl()
+        gsl_matrix_float_add(g1.nativeHandle, b.toGsl().nativeHandle)
+        return g1
+    }
+
+    public override fun multiply(a: Matrix<Float>, k: Number): GslMatrix<Float, gsl_matrix_float> {
+        val g1 = a.toGsl()
+        gsl_matrix_float_scale(g1.nativeHandle, k.toDouble())
+        return g1
+    }
+}
+
+public class GslComplexMatrixContext(scope: DeferScope) :
+    GslMatrixContext<Complex, gsl_matrix_complex, gsl_vector_complex>(scope) {
+    override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix<Complex, gsl_matrix_complex> = GslComplexMatrix(
+        nativeHandle = requireNotNull(gsl_matrix_complex_alloc(rows.toULong(), columns.toULong())),
+        scope = scope
+    )
+
+    override fun produceDirtyVector(size: Int): GslVector<Complex, gsl_vector_complex> =
+        GslComplexVector(requireNotNull(gsl_vector_complex_alloc(size.toULong())), scope)
+
+    public override fun Matrix<Complex>.dot(other: Matrix<Complex>): GslMatrix<Complex, gsl_matrix_complex> {
+        val x = toGsl().nativeHandle
+        val a = other.toGsl().nativeHandle
+        val result = requireNotNull(gsl_matrix_complex_calloc(a.pointed.size1, a.pointed.size2))
+        gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, ComplexField.one.toGsl(), x, a, ComplexField.one.toGsl(), result)
+        return GslComplexMatrix(nativeHandle = result, scope = scope)
+    }
+
+    public override fun Matrix<Complex>.dot(vector: Point<Complex>): GslVector<Complex, gsl_vector_complex> {
+        val x = toGsl().nativeHandle
+        val a = vector.toGsl().nativeHandle
+        val result = requireNotNull(gsl_vector_complex_calloc(a.pointed.size))
+        gsl_blas_zgemv(CblasNoTrans, ComplexField.one.toGsl(), x, a, ComplexField.one.toGsl(), result)
+        return GslComplexVector(result, scope)
+    }
+
+    public override fun Matrix<Complex>.times(value: Complex): GslMatrix<Complex, gsl_matrix_complex> {
+        val g1 = toGsl()
+        gsl_matrix_complex_scale(g1.nativeHandle, value.toGsl())
+        return g1
+    }
+
+    public override fun add(a: Matrix<Complex>, b: Matrix<Complex>): GslMatrix<Complex, gsl_matrix_complex> {
+        val g1 = a.toGsl()
+        gsl_matrix_complex_add(g1.nativeHandle, b.toGsl().nativeHandle)
+        return g1
+    }
+
+    public override fun multiply(a: Matrix<Complex>, k: Number): GslMatrix<Complex, gsl_matrix_complex> {
+        val g1 = a.toGsl()
+        gsl_matrix_complex_scale(g1.nativeHandle, k.toComplex().toGsl())
+        return g1
+    }
+}
diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContexts.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContexts.kt
deleted file mode 100644
index 3ac6582eb..000000000
--- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMatrixContexts.kt
+++ /dev/null
@@ -1,170 +0,0 @@
-package kscience.kmath.gsl
-
-import kotlinx.cinterop.CStructVar
-import kotlinx.cinterop.DeferScope
-import kotlinx.cinterop.pointed
-import kscience.kmath.linear.MatrixContext
-import kscience.kmath.linear.Point
-import kscience.kmath.operations.Complex
-import kscience.kmath.operations.ComplexField
-import kscience.kmath.operations.toComplex
-import kscience.kmath.structures.Matrix
-import org.gnu.gsl.*
-
-internal inline fun <T : Any, H : CStructVar> GslMatrix<T, H>.fill(initializer: (Int, Int) -> T): GslMatrix<T, H> =
-    apply {
-        (0 until rowNum).forEach { row -> (0 until colNum).forEach { col -> this[row, col] = initializer(row, col) } }
-    }
-
-internal inline fun <T : Any, H : CStructVar> GslVector<T, H>.fill(initializer: (Int) -> T): GslVector<T, H> =
-    apply { (0 until size).forEach { index -> this[index] = initializer(index) } }
-
-public abstract class GslMatrixContext<T : Any, H1 : CStructVar, H2 : CStructVar> internal constructor(
-    internal val scope: DeferScope
-) : MatrixContext<T> {
-    @Suppress("UNCHECKED_CAST")
-    public fun Matrix<T>.toGsl(): GslMatrix<T, H1> = (if (this is GslMatrix<*, *>)
-        this as GslMatrix<T, H1>
-    else
-        produce(rowNum, colNum) { i, j -> this[i, j] }).copy()
-
-    @Suppress("UNCHECKED_CAST")
-    public fun Point<T>.toGsl(): GslVector<T, H2> =
-        (if (this is GslVector<*, *>) this as GslVector<T, H2> else produceDirtyVector(size).fill { this[it] }).copy()
-
-    internal abstract fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix<T, H1>
-    internal abstract fun produceDirtyVector(size: Int): GslVector<T, H2>
-
-    public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): GslMatrix<T, H1> =
-        produceDirtyMatrix(rows, columns).fill(initializer)
-}
-
-public class GslRealMatrixContext(scope: DeferScope) : GslMatrixContext<Double, gsl_matrix, gsl_vector>(scope) {
-    override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix<Double, gsl_matrix> =
-        GslRealMatrix(nativeHandle = requireNotNull(gsl_matrix_alloc(rows.toULong(), columns.toULong())), scope = scope)
-
-    override fun produceDirtyVector(size: Int): GslVector<Double, gsl_vector> =
-        GslRealVector(nativeHandle = requireNotNull(gsl_vector_alloc(size.toULong())), scope = scope)
-
-    public override fun Matrix<Double>.dot(other: Matrix<Double>): GslMatrix<Double, gsl_matrix> {
-        val x = toGsl().nativeHandle
-        val a = other.toGsl().nativeHandle
-        val result = requireNotNull(gsl_matrix_calloc(a.pointed.size1, a.pointed.size2))
-        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, x, a, 1.0, result)
-        return GslRealMatrix(result, scope = scope)
-    }
-
-    public override fun Matrix<Double>.dot(vector: Point<Double>): GslVector<Double, gsl_vector> {
-        val x = toGsl().nativeHandle
-        val a = vector.toGsl().nativeHandle
-        val result = requireNotNull(gsl_vector_calloc(a.pointed.size))
-        gsl_blas_dgemv(CblasNoTrans, 1.0, x, a, 1.0, result)
-        return GslRealVector(result, scope = scope)
-    }
-
-    public override fun Matrix<Double>.times(value: Double): GslMatrix<Double, gsl_matrix> {
-        val g1 = toGsl()
-        gsl_matrix_scale(g1.nativeHandle, value)
-        return g1
-    }
-
-    public override fun add(a: Matrix<Double>, b: Matrix<Double>): GslMatrix<Double, gsl_matrix> {
-        val g1 = a.toGsl()
-        gsl_matrix_add(g1.nativeHandle, b.toGsl().nativeHandle)
-        return g1
-    }
-
-    public override fun multiply(a: Matrix<Double>, k: Number): GslMatrix<Double, gsl_matrix> {
-        val g1 = a.toGsl()
-        gsl_matrix_scale(g1.nativeHandle, k.toDouble())
-        return g1
-    }
-}
-
-public class GslFloatMatrixContext(scope: DeferScope) :
-    GslMatrixContext<Float, gsl_matrix_float, gsl_vector_float>(scope) {
-    override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix<Float, gsl_matrix_float> =
-        GslFloatMatrix(requireNotNull(gsl_matrix_float_alloc(rows.toULong(), columns.toULong())), scope = scope)
-
-    override fun produceDirtyVector(size: Int): GslVector<Float, gsl_vector_float> =
-        GslFloatVector(requireNotNull(gsl_vector_float_alloc(size.toULong())), scope)
-
-    public override fun Matrix<Float>.dot(other: Matrix<Float>): GslMatrix<Float, gsl_matrix_float> {
-        val x = toGsl().nativeHandle
-        val a = other.toGsl().nativeHandle
-        val result = requireNotNull(gsl_matrix_float_calloc(a.pointed.size1, a.pointed.size2))
-        gsl_blas_sgemm(CblasNoTrans, CblasNoTrans, 1f, x, a, 1f, result)
-        return GslFloatMatrix(nativeHandle = result, scope = scope)
-    }
-
-    public override fun Matrix<Float>.dot(vector: Point<Float>): GslVector<Float, gsl_vector_float> {
-        val x = toGsl().nativeHandle
-        val a = vector.toGsl().nativeHandle
-        val result = requireNotNull(gsl_vector_float_calloc(a.pointed.size))
-        gsl_blas_sgemv(CblasNoTrans, 1f, x, a, 1f, result)
-        return GslFloatVector(nativeHandle = result, scope = scope)
-    }
-
-    public override fun Matrix<Float>.times(value: Float): GslMatrix<Float, gsl_matrix_float> {
-        val g1 = toGsl()
-        gsl_matrix_float_scale(g1.nativeHandle, value.toDouble())
-        return g1
-    }
-
-    public override fun add(a: Matrix<Float>, b: Matrix<Float>): GslMatrix<Float, gsl_matrix_float> {
-        val g1 = a.toGsl()
-        gsl_matrix_float_add(g1.nativeHandle, b.toGsl().nativeHandle)
-        return g1
-    }
-
-    public override fun multiply(a: Matrix<Float>, k: Number): GslMatrix<Float, gsl_matrix_float> {
-        val g1 = a.toGsl()
-        gsl_matrix_float_scale(g1.nativeHandle, k.toDouble())
-        return g1
-    }
-}
-
-public class GslComplexMatrixContext(scope: DeferScope) :
-    GslMatrixContext<Complex, gsl_matrix_complex, gsl_vector_complex>(scope) {
-    override fun produceDirtyMatrix(rows: Int, columns: Int): GslMatrix<Complex, gsl_matrix_complex> = GslComplexMatrix(
-        nativeHandle = requireNotNull(gsl_matrix_complex_alloc(rows.toULong(), columns.toULong())),
-        scope = scope
-    )
-
-    override fun produceDirtyVector(size: Int): GslVector<Complex, gsl_vector_complex> =
-        GslComplexVector(requireNotNull(gsl_vector_complex_alloc(size.toULong())), scope)
-
-    public override fun Matrix<Complex>.dot(other: Matrix<Complex>): GslMatrix<Complex, gsl_matrix_complex> {
-        val x = toGsl().nativeHandle
-        val a = other.toGsl().nativeHandle
-        val result = requireNotNull(gsl_matrix_complex_calloc(a.pointed.size1, a.pointed.size2))
-        gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, ComplexField.one.toGsl(), x, a, ComplexField.one.toGsl(), result)
-        return GslComplexMatrix(nativeHandle = result, scope = scope)
-    }
-
-    public override fun Matrix<Complex>.dot(vector: Point<Complex>): GslVector<Complex, gsl_vector_complex> {
-        val x = toGsl().nativeHandle
-        val a = vector.toGsl().nativeHandle
-        val result = requireNotNull(gsl_vector_complex_calloc(a.pointed.size))
-        gsl_blas_zgemv(CblasNoTrans, ComplexField.one.toGsl(), x, a, ComplexField.one.toGsl(), result)
-        return GslComplexVector(result, scope)
-    }
-
-    public override fun Matrix<Complex>.times(value: Complex): GslMatrix<Complex, gsl_matrix_complex> {
-        val g1 = toGsl()
-        gsl_matrix_complex_scale(g1.nativeHandle, value.toGsl())
-        return g1
-    }
-
-    public override fun add(a: Matrix<Complex>, b: Matrix<Complex>): GslMatrix<Complex, gsl_matrix_complex> {
-        val g1 = a.toGsl()
-        gsl_matrix_complex_add(g1.nativeHandle, b.toGsl().nativeHandle)
-        return g1
-    }
-
-    public override fun multiply(a: Matrix<Complex>, k: Number): GslMatrix<Complex, gsl_matrix_complex> {
-        val g1 = a.toGsl()
-        gsl_matrix_complex_scale(g1.nativeHandle, k.toComplex().toGsl())
-        return g1
-    }
-}
diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMemoryHolder.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMemoryHolder.kt
index d3544874d..d52b78ccc 100644
--- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMemoryHolder.kt
+++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/GslMemoryHolder.kt
@@ -8,6 +8,7 @@ public abstract class GslMemoryHolder<H : CStructVar> internal constructor(inter
     internal abstract val nativeHandle: CPointer<H>
 
     init {
+        ensureHasGslErrorHandler()
         scope.defer(::close)
     }
 
diff --git a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Vectors.kt b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Vectors.kt
index 8990416b3..7ef8c24e7 100644
--- a/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Vectors.kt
+++ b/kmath-gsl/src/nativeMain/kotlin/kscience/kmath/gsl/_Vectors.kt
@@ -1,10 +1,12 @@
 package kscience.kmath.gsl
 
-import kotlinx.cinterop.*
+import kotlinx.cinterop.CPointer
+import kotlinx.cinterop.DeferScope
+import kotlinx.cinterop.pointed
 import org.gnu.gsl.*
 
-internal class GslRealVector(override val nativeHandle: CPointer<gsl_vector>, scope: DeferScope) 
-            : GslVector<Double, gsl_vector>(scope) {
+internal class GslRealVector(override val nativeHandle: CPointer<gsl_vector>, scope: DeferScope) : 
+    GslVector<Double, gsl_vector>(scope) {
     override val size: Int
         get() = nativeHandle.pointed.size.toInt()
 
@@ -25,8 +27,8 @@ internal class GslRealVector(override val nativeHandle: CPointer<gsl_vector>, sc
     override fun close(): Unit = gsl_vector_free(nativeHandle)
 }
 
-internal class GslFloatVector(override val nativeHandle: CPointer<gsl_vector_float>, scope: DeferScope) 
-            : GslVector<Float, gsl_vector_float>(scope) {
+internal class GslFloatVector(override val nativeHandle: CPointer<gsl_vector_float>, scope: DeferScope) : 
+    GslVector<Float, gsl_vector_float>(scope) {
     override val size: Int
         get() = nativeHandle.pointed.size.toInt()
 
@@ -47,8 +49,8 @@ internal class GslFloatVector(override val nativeHandle: CPointer<gsl_vector_flo
     override fun close(): Unit = gsl_vector_float_free(nativeHandle)
 }
 
-internal class GslShortVector(override val nativeHandle: CPointer<gsl_vector_short>, scope: DeferScope) 
-            : GslVector<Short, gsl_vector_short>(scope) {
+internal class GslShortVector(override val nativeHandle: CPointer<gsl_vector_short>, scope: DeferScope) : 
+    GslVector<Short, gsl_vector_short>(scope) {
     override val size: Int
         get() = nativeHandle.pointed.size.toInt()
 
@@ -69,8 +71,8 @@ internal class GslShortVector(override val nativeHandle: CPointer<gsl_vector_sho
     override fun close(): Unit = gsl_vector_short_free(nativeHandle)
 }
 
-internal class GslUShortVector(override val nativeHandle: CPointer<gsl_vector_ushort>, scope: DeferScope) 
-            : GslVector<UShort, gsl_vector_ushort>(scope) {
+internal class GslUShortVector(override val nativeHandle: CPointer<gsl_vector_ushort>, scope: DeferScope) : 
+    GslVector<UShort, gsl_vector_ushort>(scope) {
     override val size: Int
         get() = nativeHandle.pointed.size.toInt()
 
@@ -91,8 +93,8 @@ internal class GslUShortVector(override val nativeHandle: CPointer<gsl_vector_us
     override fun close(): Unit = gsl_vector_ushort_free(nativeHandle)
 }
 
-internal class GslLongVector(override val nativeHandle: CPointer<gsl_vector_long>, scope: DeferScope) 
-            : GslVector<Long, gsl_vector_long>(scope) {
+internal class GslLongVector(override val nativeHandle: CPointer<gsl_vector_long>, scope: DeferScope) : 
+    GslVector<Long, gsl_vector_long>(scope) {
     override val size: Int
         get() = nativeHandle.pointed.size.toInt()
 
@@ -113,8 +115,8 @@ internal class GslLongVector(override val nativeHandle: CPointer<gsl_vector_long
     override fun close(): Unit = gsl_vector_long_free(nativeHandle)
 }
 
-internal class GslULongVector(override val nativeHandle: CPointer<gsl_vector_ulong>, scope: DeferScope) 
-            : GslVector<ULong, gsl_vector_ulong>(scope) {
+internal class GslULongVector(override val nativeHandle: CPointer<gsl_vector_ulong>, scope: DeferScope) : 
+    GslVector<ULong, gsl_vector_ulong>(scope) {
     override val size: Int
         get() = nativeHandle.pointed.size.toInt()
 
@@ -135,8 +137,8 @@ internal class GslULongVector(override val nativeHandle: CPointer<gsl_vector_ulo
     override fun close(): Unit = gsl_vector_ulong_free(nativeHandle)
 }
 
-internal class GslIntVector(override val nativeHandle: CPointer<gsl_vector_int>, scope: DeferScope) 
-            : GslVector<Int, gsl_vector_int>(scope) {
+internal class GslIntVector(override val nativeHandle: CPointer<gsl_vector_int>, scope: DeferScope) : 
+    GslVector<Int, gsl_vector_int>(scope) {
     override val size: Int
         get() = nativeHandle.pointed.size.toInt()
 
@@ -157,8 +159,8 @@ internal class GslIntVector(override val nativeHandle: CPointer<gsl_vector_int>,
     override fun close(): Unit = gsl_vector_int_free(nativeHandle)
 }
 
-internal class GslUIntVector(override val nativeHandle: CPointer<gsl_vector_uint>, scope: DeferScope) 
-            : GslVector<UInt, gsl_vector_uint>(scope) {
+internal class GslUIntVector(override val nativeHandle: CPointer<gsl_vector_uint>, scope: DeferScope) : 
+    GslVector<UInt, gsl_vector_uint>(scope) {
     override val size: Int
         get() = nativeHandle.pointed.size.toInt()
 
diff --git a/kmath-gsl/src/nativeTest/kotlin/ErrorHandler.kt b/kmath-gsl/src/nativeTest/kotlin/ErrorHandler.kt
new file mode 100644
index 000000000..1259f58b9
--- /dev/null
+++ b/kmath-gsl/src/nativeTest/kotlin/ErrorHandler.kt
@@ -0,0 +1,23 @@
+package kscience.kmath.gsl
+
+import kotlinx.cinterop.memScoped
+import org.gnu.gsl.gsl_block_calloc
+import kotlin.test.Test
+import kotlin.test.assertFailsWith
+
+internal class ErrorHandler {
+    @Test
+    fun blockAllocation() {
+        assertFailsWith<GslException> {
+            ensureHasGslErrorHandler()
+            gsl_block_calloc(ULong.MAX_VALUE)
+        }
+    }
+
+    @Test
+    fun matrixAllocation() {
+        assertFailsWith<GslException> {
+            memScoped { GslRealMatrixContext(this).produce(Int.MAX_VALUE, Int.MAX_VALUE) { _, _ -> 0.0 } }
+        }
+    }
+}