From 1e6cf426a9067f03c60436af86ff2caddb6663ca Mon Sep 17 00:00:00 2001
From: Adrian Mariano <avm4@cornell.edu>
Date: Sat, 29 Feb 2020 22:59:39 -0500
Subject: [PATCH] changed order so linear_solve and submatrix are first.
---
math.scad | 72 ++++++++++++++++++++++++++++---------------------------
1 file changed, 37 insertions(+), 35 deletions(-)
diff --git a/math.scad b/math.scad
index d270831..010c307 100644
--- a/math.scad
+++ b/math.scad
@@ -532,6 +532,43 @@ function mean(v) = sum(v)/len(v);
// Section: Matrix math
+// Function: linear_solve()
+// Usage: linear_solve(A,b)
+// Description:
+// Solves the linear system Ax=b. If A is square and non-singular the unique solution is returned. If A is overdetermined
+// the least squares solution is returned. If A is underdetermined, the minimal norm solution is returned.
+// If A is rank deficient or singular then linear_solve returns `undef`.
+function linear_solve(A,b) =
+ let(
+ dim = array_dim(A),
+ m=dim[0], n=dim[1]
+ )
+ assert(len(b)==m,str("Incompatible matrix and vector",dim,len(b)))
+ let (
+ qr = m<n ? qr_factor(transpose(A)) : qr_factor(A),
+ maxdim = max(n,m),
+ mindim = min(n,m),
+ Q = submatrix(qr[0],[0:maxdim-1], [0:mindim-1]),
+ R = submatrix(qr[1],[0:mindim-1], [0:mindim-1]),
+ zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i]
+ )
+ zeros != [] ? undef :
+ m<n ? Q*back_substitute(R,b,transpose=true) :
+ back_substitute(R, transpose(Q)*b);
+
+
+// Function: submatrix()
+// Usage: submatrix(M, ind1, ind2)
+// Description:
+// Returns a submatrix with the specified index ranges or index sets.
+function submatrix(M,ind1,ind2) =
+ [for(i=ind1)
+ [for(j=ind2)
+ M[i][j]
+ ]
+ ];
+
+
// Function: qr_factor()
// Usage: qr = qr_factor(A)
// Description:
@@ -563,41 +600,6 @@ function _qr_factor(A,Q, column, m, n) =
_qr_factor(Qf*A, Q*Qf, column+1, m, n);
-// Function: submatrix()
-// Usage: submatrix(M, ind1, ind2)
-// Description:
-// Returns a submatrix with the specified index ranges or index sets.
-function submatrix(M,ind1,ind2) =
- [for(i=ind1)
- [for(j=ind2)
- M[i][j]
- ]
- ];
-
-
-// Function: linear_solve()
-// Usage: linear_solve(A,b)
-// Description:
-// Solves the linear system Ax=b. If A is square and non-singular the unique solution is returned. If A is overdetermined
-// the least squares solution is returned. If A is underdetermined, the minimal norm solution is returned.
-// If A is rank deficient or singular then linear_solve returns `undef`.
-function linear_solve(A,b) =
- let(
- dim = array_dim(A),
- m=dim[0], n=dim[1]
- )
- assert(len(b)==m,str("Incompatible matrix and vector",dim,len(b)))
- let (
- qr = m<n ? qr_factor(transpose(A)) : qr_factor(A),
- maxdim = max(n,m),
- mindim = min(n,m),
- Q = submatrix(qr[0],[0:maxdim-1], [0:mindim-1]),
- R = submatrix(qr[1],[0:mindim-1], [0:mindim-1]),
- zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i]
- )
- zeros != [] ? undef :
- m<n ? Q*back_substitute(R,b,transpose=true) :
- back_substitute(R, transpose(Q)*b);
// Function: back_substitute()
// Usage: back_substitute(R, b, [transpose])