diff --git a/lib/linalg/dlauu2.cpp b/lib/linalg/dlauu2.cpp
new file mode 100644
index 0000000000..d90a84798d
--- /dev/null
+++ b/lib/linalg/dlauu2.cpp
@@ -0,0 +1,77 @@
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "lmp_f2c.h"
+static doublereal c_b7 = 1.;
+static integer c__1 = 1;
+int dlauu2_(char *uplo, integer *n, doublereal *a, integer *lda, integer *info, ftnlen uplo_len)
+{
+    integer a_dim1, a_offset, i__1, i__2, i__3;
+    integer i__;
+    doublereal aii;
+    extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, integer *);
+    extern int dscal_(integer *, doublereal *, doublereal *, integer *);
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    extern int dgemv_(char *, integer *, integer *, doublereal *, doublereal *, integer *,
+                      doublereal *, integer *, doublereal *, doublereal *, integer *, ftnlen);
+    logical upper;
+    extern int xerbla_(char *, integer *, ftnlen);
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    *info = 0;
+    upper = lsame_(uplo, (char *)"U", (ftnlen)1, (ftnlen)1);
+    if (!upper && !lsame_(uplo, (char *)"L", (ftnlen)1, (ftnlen)1)) {
+        *info = -1;
+    } else if (*n < 0) {
+        *info = -2;
+    } else if (*lda < max(1, *n)) {
+        *info = -4;
+    }
+    if (*info != 0) {
+        i__1 = -(*info);
+        xerbla_((char *)"DLAUU2", &i__1, (ftnlen)6);
+        return 0;
+    }
+    if (*n == 0) {
+        return 0;
+    }
+    if (upper) {
+        i__1 = *n;
+        for (i__ = 1; i__ <= i__1; ++i__) {
+            aii = a[i__ + i__ * a_dim1];
+            if (i__ < *n) {
+                i__2 = *n - i__ + 1;
+                a[i__ + i__ * a_dim1] =
+                    ddot_(&i__2, &a[i__ + i__ * a_dim1], lda, &a[i__ + i__ * a_dim1], lda);
+                i__2 = i__ - 1;
+                i__3 = *n - i__;
+                dgemv_((char *)"No transpose", &i__2, &i__3, &c_b7, &a[(i__ + 1) * a_dim1 + 1], lda,
+                       &a[i__ + (i__ + 1) * a_dim1], lda, &aii, &a[i__ * a_dim1 + 1], &c__1,
+                       (ftnlen)12);
+            } else {
+                dscal_(&i__, &aii, &a[i__ * a_dim1 + 1], &c__1);
+            }
+        }
+    } else {
+        i__1 = *n;
+        for (i__ = 1; i__ <= i__1; ++i__) {
+            aii = a[i__ + i__ * a_dim1];
+            if (i__ < *n) {
+                i__2 = *n - i__ + 1;
+                a[i__ + i__ * a_dim1] =
+                    ddot_(&i__2, &a[i__ + i__ * a_dim1], &c__1, &a[i__ + i__ * a_dim1], &c__1);
+                i__2 = *n - i__;
+                i__3 = i__ - 1;
+                dgemv_((char *)"Transpose", &i__2, &i__3, &c_b7, &a[i__ + 1 + a_dim1], lda,
+                       &a[i__ + 1 + i__ * a_dim1], &c__1, &aii, &a[i__ + a_dim1], lda, (ftnlen)9);
+            } else {
+                dscal_(&i__, &aii, &a[i__ + a_dim1], lda);
+            }
+        }
+    }
+    return 0;
+}
+#ifdef __cplusplus
+}
+#endif
diff --git a/lib/linalg/dlauum.cpp b/lib/linalg/dlauum.cpp
new file mode 100644
index 0000000000..632bd4ba85
--- /dev/null
+++ b/lib/linalg/dlauum.cpp
@@ -0,0 +1,101 @@
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "lmp_f2c.h"
+static integer c__1 = 1;
+static integer c_n1 = -1;
+static doublereal c_b15 = 1.;
+int dlauum_(char *uplo, integer *n, doublereal *a, integer *lda, integer *info, ftnlen uplo_len)
+{
+    integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
+    integer i__, ib, nb;
+    extern int dgemm_(char *, char *, integer *, integer *, integer *, doublereal *, doublereal *,
+                      integer *, doublereal *, integer *, doublereal *, doublereal *, integer *,
+                      ftnlen, ftnlen);
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    extern int dtrmm_(char *, char *, char *, char *, integer *, integer *, doublereal *,
+                      doublereal *, integer *, doublereal *, integer *, ftnlen, ftnlen, ftnlen,
+                      ftnlen);
+    logical upper;
+    extern int dsyrk_(char *, char *, integer *, integer *, doublereal *, doublereal *, integer *,
+                      doublereal *, doublereal *, integer *, ftnlen, ftnlen),
+        dlauu2_(char *, integer *, doublereal *, integer *, integer *, ftnlen),
+        xerbla_(char *, integer *, ftnlen);
+    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *,
+                           ftnlen, ftnlen);
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    *info = 0;
+    upper = lsame_(uplo, (char *)"U", (ftnlen)1, (ftnlen)1);
+    if (!upper && !lsame_(uplo, (char *)"L", (ftnlen)1, (ftnlen)1)) {
+        *info = -1;
+    } else if (*n < 0) {
+        *info = -2;
+    } else if (*lda < max(1, *n)) {
+        *info = -4;
+    }
+    if (*info != 0) {
+        i__1 = -(*info);
+        xerbla_((char *)"DLAUUM", &i__1, (ftnlen)6);
+        return 0;
+    }
+    if (*n == 0) {
+        return 0;
+    }
+    nb = ilaenv_(&c__1, (char *)"DLAUUM", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
+    if (nb <= 1 || nb >= *n) {
+        dlauu2_(uplo, n, &a[a_offset], lda, info, (ftnlen)1);
+    } else {
+        if (upper) {
+            i__1 = *n;
+            i__2 = nb;
+            for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
+                i__3 = nb, i__4 = *n - i__ + 1;
+                ib = min(i__3, i__4);
+                i__3 = i__ - 1;
+                dtrmm_((char *)"Right", (char *)"Upper", (char *)"Transpose", (char *)"Non-unit", &i__3, &ib, &c_b15,
+                       &a[i__ + i__ * a_dim1], lda, &a[i__ * a_dim1 + 1], lda, (ftnlen)5, (ftnlen)5,
+                       (ftnlen)9, (ftnlen)8);
+                dlauu2_((char *)"Upper", &ib, &a[i__ + i__ * a_dim1], lda, info, (ftnlen)5);
+                if (i__ + ib <= *n) {
+                    i__3 = i__ - 1;
+                    i__4 = *n - i__ - ib + 1;
+                    dgemm_((char *)"No transpose", (char *)"Transpose", &i__3, &ib, &i__4, &c_b15,
+                           &a[(i__ + ib) * a_dim1 + 1], lda, &a[i__ + (i__ + ib) * a_dim1], lda,
+                           &c_b15, &a[i__ * a_dim1 + 1], lda, (ftnlen)12, (ftnlen)9);
+                    i__3 = *n - i__ - ib + 1;
+                    dsyrk_((char *)"Upper", (char *)"No transpose", &ib, &i__3, &c_b15,
+                           &a[i__ + (i__ + ib) * a_dim1], lda, &c_b15, &a[i__ + i__ * a_dim1], lda,
+                           (ftnlen)5, (ftnlen)12);
+                }
+            }
+        } else {
+            i__2 = *n;
+            i__1 = nb;
+            for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) {
+                i__3 = nb, i__4 = *n - i__ + 1;
+                ib = min(i__3, i__4);
+                i__3 = i__ - 1;
+                dtrmm_((char *)"Left", (char *)"Lower", (char *)"Transpose", (char *)"Non-unit", &ib, &i__3, &c_b15,
+                       &a[i__ + i__ * a_dim1], lda, &a[i__ + a_dim1], lda, (ftnlen)4, (ftnlen)5,
+                       (ftnlen)9, (ftnlen)8);
+                dlauu2_((char *)"Lower", &ib, &a[i__ + i__ * a_dim1], lda, info, (ftnlen)5);
+                if (i__ + ib <= *n) {
+                    i__3 = i__ - 1;
+                    i__4 = *n - i__ - ib + 1;
+                    dgemm_((char *)"Transpose", (char *)"No transpose", &ib, &i__3, &i__4, &c_b15,
+                           &a[i__ + ib + i__ * a_dim1], lda, &a[i__ + ib + a_dim1], lda, &c_b15,
+                           &a[i__ + a_dim1], lda, (ftnlen)9, (ftnlen)12);
+                    i__3 = *n - i__ - ib + 1;
+                    dsyrk_((char *)"Lower", (char *)"Transpose", &ib, &i__3, &c_b15, &a[i__ + ib + i__ * a_dim1],
+                           lda, &c_b15, &a[i__ + i__ * a_dim1], lda, (ftnlen)5, (ftnlen)9);
+                }
+            }
+        }
+    }
+    return 0;
+}
+#ifdef __cplusplus
+}
+#endif
diff --git a/lib/linalg/dpotri.cpp b/lib/linalg/dpotri.cpp
new file mode 100644
index 0000000000..9c0a609e1b
--- /dev/null
+++ b/lib/linalg/dpotri.cpp
@@ -0,0 +1,40 @@
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "lmp_f2c.h"
+int dpotri_(char *uplo, integer *n, doublereal *a, integer *lda, integer *info, ftnlen uplo_len)
+{
+    integer a_dim1, a_offset, i__1;
+    extern logical lsame_(char *, char *, ftnlen, ftnlen);
+    extern int xerbla_(char *, integer *, ftnlen),
+        dlauum_(char *, integer *, doublereal *, integer *, integer *, ftnlen),
+        dtrtri_(char *, char *, integer *, doublereal *, integer *, integer *, ftnlen, ftnlen);
+    a_dim1 = *lda;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+    *info = 0;
+    if (!lsame_(uplo, (char *)"U", (ftnlen)1, (ftnlen)1) && !lsame_(uplo, (char *)"L", (ftnlen)1, (ftnlen)1)) {
+        *info = -1;
+    } else if (*n < 0) {
+        *info = -2;
+    } else if (*lda < max(1, *n)) {
+        *info = -4;
+    }
+    if (*info != 0) {
+        i__1 = -(*info);
+        xerbla_((char *)"DPOTRI", &i__1, (ftnlen)6);
+        return 0;
+    }
+    if (*n == 0) {
+        return 0;
+    }
+    dtrtri_(uplo, (char *)"Non-unit", n, &a[a_offset], lda, info, (ftnlen)1, (ftnlen)8);
+    if (*info > 0) {
+        return 0;
+    }
+    dlauum_(uplo, n, &a[a_offset], lda, info, (ftnlen)1);
+    return 0;
+}
+#ifdef __cplusplus
+}
+#endif