grey/BOSL2
1
0
Fork 0
mirror of https://github.com/BelfrySCAD/BOSL2.git synced 2024-12-29 16:29:40 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

changed order so linear_solve and submatrix are first.

Browse source
This commit is contained in:
Adrian Mariano 2020-02-29 22:59:39 -05:00
parent 5bd2cba0ff
commit 1e6cf426a9
1 changed files with 37 additions and 35 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

72
math.scad
View file

@ -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])