minverse   xtlang

(bind-func minverse
  "matrix inverse - matrix must be square (= nrows ncols)"
  (lambda (m1:double* nrows:i64 result:double*)
    (let ((m2:double* (salloc (* nrows nrows)))
          (fac:double* (salloc (* nrows nrows)))
          (d:double 0.0)
          (i 0) (j 0) (m 0) (n 0) (q:i64 0) (p:i64 0))
      (dotimes (q nrows)
        (dotimes (p nrows)
          (set! m 0)
          (set! n 0)
          (dotimes (i nrows)
            (dotimes (j nrows)
              (pset! m2 (+ (* i nrows) j) 0.0)
              (if (and (<> i q) (<> j p))
                  (begin
                    (pset! m2 (+ (* m (- nrows 1)) n)
                           (pref m1 (+ (* i nrows) j)))
                    (if (< n (- nrows 2))
                        (set! n (+ n 1))
                        (begin
                          (set! n 0)
                          (set! m (+ m 1)))))
                  (begin 1))))
          (pset! fac (+ (* q nrows) p)
                 (* (pow -1.0 (i64tod (+ q p)))
                    (mdeterminant m2 (- nrows 1))))))
      (dotimes (i nrows)
        (dotimes (j nrows)
          (pset! m2 (+ (* i nrows) j)
                 (pref fac (+ (* j nrows) i)))))
      (set! d (mdeterminant m1 nrows))
      (dotimes (i nrows)
        (dotimes (j nrows)
          (pset! result (+ (* i nrows) j)
                 (/ (pref m2 (+ (* i nrows) j)) d))))
      void)))