diff --git a/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/codegenUtilities.kt b/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/codegenUtilities.kt
new file mode 100644
index 000000000..0df7141d7
--- /dev/null
+++ b/buildSrc/src/main/kotlin/kscience/kmath/gsl/codegen/codegenUtilities.kt
@@ -0,0 +1,61 @@
+package kscience.kmath.gsl.codegen
+
+import org.jetbrains.kotlin.cli.common.CLIConfigurationKeys
+import org.jetbrains.kotlin.cli.common.messages.MessageCollector
+import org.jetbrains.kotlin.cli.jvm.compiler.EnvironmentConfigFiles
+import org.jetbrains.kotlin.cli.jvm.compiler.KotlinCoreEnvironment
+import org.jetbrains.kotlin.com.intellij.mock.MockProject
+import org.jetbrains.kotlin.com.intellij.openapi.extensions.ExtensionPoint
+import org.jetbrains.kotlin.com.intellij.openapi.extensions.Extensions
+import org.jetbrains.kotlin.com.intellij.openapi.util.UserDataHolderBase
+import org.jetbrains.kotlin.com.intellij.pom.PomModel
+import org.jetbrains.kotlin.com.intellij.pom.PomModelAspect
+import org.jetbrains.kotlin.com.intellij.pom.PomTransaction
+import org.jetbrains.kotlin.com.intellij.pom.impl.PomTransactionBase
+import org.jetbrains.kotlin.com.intellij.pom.tree.TreeAspect
+import org.jetbrains.kotlin.com.intellij.psi.PsiElement
+import org.jetbrains.kotlin.com.intellij.psi.impl.source.tree.TreeCopyHandler
+import org.jetbrains.kotlin.config.CompilerConfiguration
+import org.jetbrains.kotlin.psi.KtFile
+import sun.reflect.ReflectionFactory
+
+internal fun createProject(): MockProject {
+    val project = KotlinCoreEnvironment.createForProduction(
+        {},
+        CompilerConfiguration().apply { put(CLIConfigurationKeys.MESSAGE_COLLECTOR_KEY, MessageCollector.NONE) },
+        EnvironmentConfigFiles.JVM_CONFIG_FILES
+    ).project as MockProject
+
+    val extensionPoint = "org.jetbrains.kotlin.com.intellij.treeCopyHandler"
+
+    arrayOf(project.extensionArea, Extensions.getRootArea())
+        .asSequence()
+        .filterNot { it.hasExtensionPoint(extensionPoint) }
+        .forEach {
+            it.registerExtensionPoint(extensionPoint, TreeCopyHandler::class.java.name, ExtensionPoint.Kind.INTERFACE)
+        }
+
+    project.registerService(PomModel::class.java, object : UserDataHolderBase(), PomModel {
+        override fun runTransaction(transaction: PomTransaction) = (transaction as PomTransactionBase).run()
+
+        @Suppress("UNCHECKED_CAST")
+        override fun <T : PomModelAspect> getModelAspect(aspect: Class<T>): T? {
+            if (aspect == TreeAspect::class.java) {
+                val constructor = ReflectionFactory.getReflectionFactory().newConstructorForSerialization(
+                    aspect,
+                    Any::class.java.getDeclaredConstructor(*arrayOfNulls(0))
+                )
+
+                return constructor.newInstance() as T
+            }
+
+            return null
+        }
+    })
+
+    return project
+}
+
+internal operator fun PsiElement.plusAssign(e: PsiElement) {
+    add(e)
+}
diff --git a/kmath-gsl/src/nativeMain/kotlin/generated/Matrices.kt b/kmath-gsl/src/nativeMain/kotlin/generated/Matrices.kt
new file mode 100644
index 000000000..47cb145c8
--- /dev/null
+++ b/kmath-gsl/src/nativeMain/kotlin/generated/Matrices.kt
@@ -0,0 +1,279 @@
+package kscience.kmath.gsl
+
+import kotlinx.cinterop.*
+import kscience.kmath.linear.*
+import kscience.kmath.operations.*
+import org.gnu.gsl.*
+
+internal class GslRealMatrix(
+    override val nativeHandle: CPointer<gsl_matrix>,
+    features: Set<MatrixFeature> = emptySet()
+) : GslMatrix<Double, gsl_matrix>() {
+    override val rowNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size1.toInt() }
+
+    override val colNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size2.toInt() }
+
+    override val features: Set<MatrixFeature> = features
+
+    override fun suggestFeature(vararg features: MatrixFeature): GslRealMatrix =
+        GslRealMatrix(nativeHandle, this.features + features)
+
+    override operator fun get(i: Int, j: Int): Double = gsl_matrix_get(nativeHandle, i.toULong(), j.toULong())
+
+    override operator fun set(i: Int, j: Int, value: Double): Unit =
+        gsl_matrix_set(nativeHandle, i.toULong(), j.toULong(), value)
+
+    override fun copy(): GslRealMatrix = memScoped {
+        val new = requireNotNull(gsl_matrix_alloc(rowNum.toULong(), colNum.toULong()))
+        gsl_matrix_memcpy(new, nativeHandle)
+        GslRealMatrix(new, features)
+    }
+
+    override fun close(): Unit = gsl_matrix_free(nativeHandle)
+
+    override fun equals(other: Any?): Boolean {
+        if (other is GslRealMatrix) gsl_matrix_equal(nativeHandle, other.nativeHandle)
+        return super.equals(other)
+    }
+}
+
+internal class GslFloatMatrix(
+    override val nativeHandle: CPointer<gsl_matrix_float>,
+    features: Set<MatrixFeature> = emptySet()
+) : GslMatrix<Float, gsl_matrix_float>() {
+    override val rowNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size1.toInt() }
+
+    override val colNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size2.toInt() }
+
+    override val features: Set<MatrixFeature> = features
+
+    override fun suggestFeature(vararg features: MatrixFeature): GslFloatMatrix =
+        GslFloatMatrix(nativeHandle, this.features + features)
+
+    override operator fun get(i: Int, j: Int): Float = gsl_matrix_float_get(nativeHandle, i.toULong(), j.toULong())
+
+    override operator fun set(i: Int, j: Int, value: Float): Unit =
+        gsl_matrix_float_set(nativeHandle, i.toULong(), j.toULong(), value)
+
+    override fun copy(): GslFloatMatrix = memScoped {
+        val new = requireNotNull(gsl_matrix_float_alloc(rowNum.toULong(), colNum.toULong()))
+        gsl_matrix_float_memcpy(new, nativeHandle)
+        GslFloatMatrix(new, features)
+    }
+
+    override fun close(): Unit = gsl_matrix_float_free(nativeHandle)
+
+    override fun equals(other: Any?): Boolean {
+        if (other is GslFloatMatrix) gsl_matrix_float_equal(nativeHandle, other.nativeHandle)
+        return super.equals(other)
+    }
+}
+
+internal class GslShortMatrix(
+    override val nativeHandle: CPointer<gsl_matrix_short>,
+    features: Set<MatrixFeature> = emptySet()
+) : GslMatrix<Short, gsl_matrix_short>() {
+    override val rowNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size1.toInt() }
+
+    override val colNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size2.toInt() }
+
+    override val features: Set<MatrixFeature> = features
+
+    override fun suggestFeature(vararg features: MatrixFeature): GslShortMatrix =
+        GslShortMatrix(nativeHandle, this.features + features)
+
+    override operator fun get(i: Int, j: Int): Short = gsl_matrix_short_get(nativeHandle, i.toULong(), j.toULong())
+
+    override operator fun set(i: Int, j: Int, value: Short): Unit =
+        gsl_matrix_short_set(nativeHandle, i.toULong(), j.toULong(), value)
+
+    override fun copy(): GslShortMatrix = memScoped {
+        val new = requireNotNull(gsl_matrix_short_alloc(rowNum.toULong(), colNum.toULong()))
+        gsl_matrix_short_memcpy(new, nativeHandle)
+        GslShortMatrix(new, features)
+    }
+
+    override fun close(): Unit = gsl_matrix_short_free(nativeHandle)
+
+    override fun equals(other: Any?): Boolean {
+        if (other is GslShortMatrix) gsl_matrix_short_equal(nativeHandle, other.nativeHandle)
+        return super.equals(other)
+    }
+}
+
+internal class GslUShortMatrix(
+    override val nativeHandle: CPointer<gsl_matrix_ushort>,
+    features: Set<MatrixFeature> = emptySet()
+) : GslMatrix<UShort, gsl_matrix_ushort>() {
+    override val rowNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size1.toInt() }
+
+    override val colNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size2.toInt() }
+
+    override val features: Set<MatrixFeature> = features
+
+    override fun suggestFeature(vararg features: MatrixFeature): GslUShortMatrix =
+        GslUShortMatrix(nativeHandle, this.features + features)
+
+    override operator fun get(i: Int, j: Int): UShort = gsl_matrix_ushort_get(nativeHandle, i.toULong(), j.toULong())
+
+    override operator fun set(i: Int, j: Int, value: UShort): Unit =
+        gsl_matrix_ushort_set(nativeHandle, i.toULong(), j.toULong(), value)
+
+    override fun copy(): GslUShortMatrix = memScoped {
+        val new = requireNotNull(gsl_matrix_ushort_alloc(rowNum.toULong(), colNum.toULong()))
+        gsl_matrix_ushort_memcpy(new, nativeHandle)
+        GslUShortMatrix(new, features)
+    }
+
+    override fun close(): Unit = gsl_matrix_ushort_free(nativeHandle)
+
+    override fun equals(other: Any?): Boolean {
+        if (other is GslUShortMatrix) gsl_matrix_ushort_equal(nativeHandle, other.nativeHandle)
+        return super.equals(other)
+    }
+}
+
+internal class GslLongMatrix(
+    override val nativeHandle: CPointer<gsl_matrix_long>,
+    features: Set<MatrixFeature> = emptySet()
+) : GslMatrix<Long, gsl_matrix_long>() {
+    override val rowNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size1.toInt() }
+
+    override val colNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size2.toInt() }
+
+    override val features: Set<MatrixFeature> = features
+
+    override fun suggestFeature(vararg features: MatrixFeature): GslLongMatrix =
+        GslLongMatrix(nativeHandle, this.features + features)
+
+    override operator fun get(i: Int, j: Int): Long = gsl_matrix_long_get(nativeHandle, i.toULong(), j.toULong())
+
+    override operator fun set(i: Int, j: Int, value: Long): Unit =
+        gsl_matrix_long_set(nativeHandle, i.toULong(), j.toULong(), value)
+
+    override fun copy(): GslLongMatrix = memScoped {
+        val new = requireNotNull(gsl_matrix_long_alloc(rowNum.toULong(), colNum.toULong()))
+        gsl_matrix_long_memcpy(new, nativeHandle)
+        GslLongMatrix(new, features)
+    }
+
+    override fun close(): Unit = gsl_matrix_long_free(nativeHandle)
+
+    override fun equals(other: Any?): Boolean {
+        if (other is GslLongMatrix) gsl_matrix_long_equal(nativeHandle, other.nativeHandle)
+        return super.equals(other)
+    }
+}
+
+internal class GslULongMatrix(
+    override val nativeHandle: CPointer<gsl_matrix_ulong>,
+    features: Set<MatrixFeature> = emptySet()
+) : GslMatrix<ULong, gsl_matrix_ulong>() {
+    override val rowNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size1.toInt() }
+
+    override val colNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size2.toInt() }
+
+    override val features: Set<MatrixFeature> = features
+
+    override fun suggestFeature(vararg features: MatrixFeature): GslULongMatrix =
+        GslULongMatrix(nativeHandle, this.features + features)
+
+    override operator fun get(i: Int, j: Int): ULong = gsl_matrix_ulong_get(nativeHandle, i.toULong(), j.toULong())
+
+    override operator fun set(i: Int, j: Int, value: ULong): Unit =
+        gsl_matrix_ulong_set(nativeHandle, i.toULong(), j.toULong(), value)
+
+    override fun copy(): GslULongMatrix = memScoped {
+        val new = requireNotNull(gsl_matrix_ulong_alloc(rowNum.toULong(), colNum.toULong()))
+        gsl_matrix_ulong_memcpy(new, nativeHandle)
+        GslULongMatrix(new, features)
+    }
+
+    override fun close(): Unit = gsl_matrix_ulong_free(nativeHandle)
+
+    override fun equals(other: Any?): Boolean {
+        if (other is GslULongMatrix) gsl_matrix_ulong_equal(nativeHandle, other.nativeHandle)
+        return super.equals(other)
+    }
+}
+
+internal class GslIntMatrix(
+    override val nativeHandle: CPointer<gsl_matrix_int>,
+    features: Set<MatrixFeature> = emptySet()
+) : GslMatrix<Int, gsl_matrix_int>() {
+    override val rowNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size1.toInt() }
+
+    override val colNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size2.toInt() }
+
+    override val features: Set<MatrixFeature> = features
+
+    override fun suggestFeature(vararg features: MatrixFeature): GslIntMatrix =
+        GslIntMatrix(nativeHandle, this.features + features)
+
+    override operator fun get(i: Int, j: Int): Int = gsl_matrix_int_get(nativeHandle, i.toULong(), j.toULong())
+
+    override operator fun set(i: Int, j: Int, value: Int): Unit =
+        gsl_matrix_int_set(nativeHandle, i.toULong(), j.toULong(), value)
+
+    override fun copy(): GslIntMatrix = memScoped {
+        val new = requireNotNull(gsl_matrix_int_alloc(rowNum.toULong(), colNum.toULong()))
+        gsl_matrix_int_memcpy(new, nativeHandle)
+        GslIntMatrix(new, features)
+    }
+
+    override fun close(): Unit = gsl_matrix_int_free(nativeHandle)
+
+    override fun equals(other: Any?): Boolean {
+        if (other is GslIntMatrix) gsl_matrix_int_equal(nativeHandle, other.nativeHandle)
+        return super.equals(other)
+    }
+}
+
+internal class GslUIntMatrix(
+    override val nativeHandle: CPointer<gsl_matrix_uint>,
+    features: Set<MatrixFeature> = emptySet()
+) : GslMatrix<UInt, gsl_matrix_uint>() {
+    override val rowNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size1.toInt() }
+
+    override val colNum: Int
+        get() = memScoped { nativeHandle.getPointer(this).pointed.size2.toInt() }
+
+    override val features: Set<MatrixFeature> = features
+
+    override fun suggestFeature(vararg features: MatrixFeature): GslUIntMatrix =
+        GslUIntMatrix(nativeHandle, this.features + features)
+
+    override operator fun get(i: Int, j: Int): UInt = gsl_matrix_uint_get(nativeHandle, i.toULong(), j.toULong())
+
+    override operator fun set(i: Int, j: Int, value: UInt): Unit =
+        gsl_matrix_uint_set(nativeHandle, i.toULong(), j.toULong(), value)
+
+    override fun copy(): GslUIntMatrix = memScoped {
+        val new = requireNotNull(gsl_matrix_uint_alloc(rowNum.toULong(), colNum.toULong()))
+        gsl_matrix_uint_memcpy(new, nativeHandle)
+        GslUIntMatrix(new, features)
+    }
+
+    override fun close(): Unit = gsl_matrix_uint_free(nativeHandle)
+
+    override fun equals(other: Any?): Boolean {
+        if (other is GslUIntMatrix) gsl_matrix_uint_equal(nativeHandle, other.nativeHandle)
+        return super.equals(other)
+    }
+}
+