From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001
From: Joseph Hunkeler <jhunkeler@gmail.com>
Date: Wed, 8 Jul 2015 20:46:52 -0400
Subject: Initial commit

---
 sys/mwcs/mwinvertd.x | 40 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)
 create mode 100644 sys/mwcs/mwinvertd.x

(limited to 'sys/mwcs/mwinvertd.x')

diff --git a/sys/mwcs/mwinvertd.x b/sys/mwcs/mwinvertd.x
new file mode 100644
index 00000000..e2744821
--- /dev/null
+++ b/sys/mwcs/mwinvertd.x
@@ -0,0 +1,40 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+# MW_INVERTD -- Invert a square matrix, double precision version.  The matrix
+# need not be symmetric.  The input and output matrices cannot be the same.
+
+procedure mw_invertd (o_ltm, n_ltm, ndim)
+
+double	o_ltm[ndim,ndim]	#I input matrix
+double	n_ltm[ndim,ndim]	#O output (inverted) matrix
+int	ndim			#I dimensionality of system
+
+pointer	sp, ix, ltm
+int	nelem, i, j
+
+begin
+	call smark (sp)
+
+	nelem = ndim * ndim
+	call salloc (ix, ndim, TY_INT)
+	call salloc (ltm, nelem, TY_DOUBLE)
+
+	# Make scratch copy (to be modified) of input matrix.
+	call amovd (o_ltm, Memd[ltm], nelem)
+
+	# Set up identity matrix.
+	do i = 1, ndim {
+	    do j = 1, ndim
+		n_ltm[i,j] = 0.0
+	    n_ltm[i,i] = 1.0
+	}
+
+	# Perform the LU decomposition.
+	call mw_ludecompose (Memd[ltm], Memi[ix], ndim)
+
+	# Compute the inverse matrix by backsubstitution.
+	do j = 1, ndim
+	    call mw_lubacksub (Memd[ltm], Memi[ix], n_ltm[1,j], ndim)
+
+	call sfree (sp)
+end
-- 
cgit