Last active
August 3, 2018 17:47
-
-
Save jrus/0a5d827265516e4239b6 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
function [Q, R] = qr(A, econ) | |
%QR QR factorization of an array-valued CHEBFUN. | |
% [Q, R] = QR(A) or QR(A, 0), where A is a column CHEBFUN with n columns, | |
% produces a column CHEBFUN Q with n orthonormal columns and an n x n upper | |
% triangular matrix R such that A = Q*R. | |
% | |
% The algorithm used is described in L.N. Trefethen, "Householder | |
% triangularization of a quasimatrix", IMA J. Numer. Anal. (30), 887-897 | |
% (2010). | |
% | |
% See also SVD, MRDIVIDE, RANK. | |
% Copyright 2015 by The University of Oxford and The Chebfun Developers. | |
% See http://www.chebfun.org/ for Chebfun information. | |
% Check inputs | |
if ( (nargin == 2) && (econ ~= 0) ) | |
error('CHEBFUN:CHEBFUN:qr:twoargs',... | |
'Use qr(A) or qr(A, 0) for QR decomposition of an array-valued CHEBFUN.'); | |
end | |
if ( A(1).isTransposed ) | |
error('CHEBFUN:CHEBFUN:qr:transpose',... | |
'CHEBFUN QR works only for column CHEBFUN objects.') | |
end | |
if ( ~all(isfinite(A(1).domain)) ) | |
error('CHEBFUN:CHEBFUN:qr:infdomain', ... | |
'CHEBFUN QR does not support unbounded domains.'); | |
end | |
numCols = numColumns(A); | |
% If A has only one column we simply scale it. | |
if ( numCols == 1 ) | |
R = sqrt(innerProduct(A, A)); | |
Q = A./R; | |
return | |
end | |
% If A is a quasimatrix then try to convert it to an array-valued CHEBFUN and then | |
% call QR on that. Otherwise, use ABSTRACTQR with a Legendre-Vandermonde matrix | |
% as the second argument. | |
if ( numel(A) > 1 ) | |
[A, isSimple] = quasi2cheb(A); | |
if ( isSimple ) % Successfully converted to array-valued CHEBFUN. | |
[Q, R] = qr(A, 0); | |
return | |
end | |
% We can't convert the quasimatrix to an array-valued CHEBFUN. Tricky case. | |
% TODO: This case is not tested! | |
% DEVELOPER NOTE: | |
% Currently (5th Feb 2016) the only way this is reachable is in the | |
% case of a quasimatrix consisting of SINGFUN or DELTAFUN objects, | |
% neither of which return anything sensible when we attempt to compute | |
% a QR factorization. | |
% Legendre-Vandermonde matrix converted to also be a quasimatrix: | |
L = cheb2quasi(legpoly(0:numCols-1, domain(A), 'norm', 1)); | |
[Q, R] = abstractQR(A, L, @innerProduct, @normest); | |
return | |
end | |
% If A is an array-valued Chebfun, find the number of pieces. | |
numFuns = numel(A.funs); | |
% If there is only one piece (no breakpoints), call QR at the fun level. | |
if ( numFuns == 1 ) | |
[Q, R] = qr(A.funs{1}); | |
Q = chebfun({Q}); | |
return | |
end | |
% If there are multiple pieces, perform QR on each piece, stack the | |
% R matrices from the separate pieces and perform QR on that, yielding | |
% Qhat and Rhat, then break Qhat back down into sections and fold those | |
% back into the matrices Q. | |
% Perform QR on each piece: | |
R = cell(numFuns, 1); | |
Q = R; | |
for k = 1:numFuns | |
[Q{k}, R{k}] = qr(A.funs{k}); | |
end | |
% Compute [Qhat, Rhat] = qr(R): | |
[Qhat, Rhat] = qr(cell2mat(R)); | |
Rhat = Rhat(1:numCols,:); % Extract required entries. | |
Qhat = Qhat(:,1:numCols); | |
% Ensure that the diagonal of Rhat is non-negative: | |
neg = diag(Rhat) < 0; | |
Rhat(neg,:) = -1 * Rhat(neg,:); | |
Qhat(:,neg) = -1 * Qhat(:,neg); | |
% ... And fold Qhat back in to Q for each segment: | |
m = cellfun(@(v) size(v, 1), R); % Number of rows per segment. | |
Qhat = mat2cell(Qhat, m, numCols); | |
for k = 1:numFuns | |
Q{k} = Q{k}*Qhat{k,1}; | |
end | |
Q = chebfun(Q); | |
R = Rhat; | |
end |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
--- untitled 8 | |
+++ (clipboard) | |
@@ -13,14 +13,6 @@ | |
% Copyright 2015 by The University of Oxford and The Chebfun Developers. | |
% See http://www.chebfun.org/ for Chebfun information. | |
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
-% Developer note: | |
-% If A contains only a single FUN, then FUN/QR is used directly. If A has | |
-% multiple pieces but each of these are simple CHEBTECH objects, then | |
-% QRSIMPLE() is called. This violates OOP principles, but is _much_ more | |
-% efficient. | |
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
- | |
% Check inputs | |
if ( (nargin == 2) && (econ ~= 0) ) | |
error('CHEBFUN:CHEBFUN:qr:twoargs',... | |
@@ -37,74 +29,78 @@ | |
numCols = numColumns(A); | |
+% If A has only one column we simply scale it. | |
if ( numCols == 1 ) | |
- % If f has only one column we simply scale it. | |
R = sqrt(innerProduct(A, A)); | |
Q = A./R; | |
- | |
-elseif ( numel(A) > 1 ) | |
- % Quasimatrix case: | |
+ return | |
+end | |
+ | |
+% If A is a quasimatrix then try to convert it to an array-valued CHEBFUN and then | |
+% call QR on that. Otherwise, use ABSTRACTQR with a Legendre-Vandermonde matrix | |
+% as the second argument. | |
+if ( numel(A) > 1 ) | |
- % Try to convert to an array-valued CHEBFUN: | |
[A, isSimple] = quasi2cheb(A); | |
- | |
if ( isSimple ) % Successfully converted to array-valued CHEBFUN. | |
- | |
[Q, R] = qr(A, 0); | |
- | |
- else % We can't convert to array valued. Tricky case. | |
- | |
- % TODO: This case is not tested! | |
- % DEVELOPER NOTE: | |
- % Currently (5th Feb 2016) the only way this is reachable is in the | |
- % case of a quasimatrix consisting of SINGFUN or DELTAFUN objects, | |
- % neither of which return anything sensible when we attempt to compute | |
- % a QR factorization. | |
- | |
- % Legendre-Vandermonde matrix: | |
- L = legpoly(0:numCols-1, domain(A), 'norm', 1); | |
- % Convert so that L is also quasimatrix: | |
- L = cheb2quasi(L); | |
- % Call abstract QR: | |
- [Q, R] = abstractQR(A, L, @innerProduct, @normest); | |
- | |
+ return | |
end | |
- | |
-elseif ( numel(A.funs) == 1 ) | |
- % No breakpoints = easy case. | |
- | |
- % Call QR at the FUN level: | |
- [Q, R] = qr(A.funs{1}); | |
- Q = chebfun({Q}); | |
-else | |
+ % We can't convert the quasimatrix to an array-valued CHEBFUN. Tricky case. | |
+ % TODO: This case is not tested! | |
+ % DEVELOPER NOTE: | |
+ % Currently (5th Feb 2016) the only way this is reachable is in the | |
+ % case of a quasimatrix consisting of SINGFUN or DELTAFUN objects, | |
+ % neither of which return anything sensible when we attempt to compute | |
+ % a QR factorization. | |
+ | |
+ % Legendre-Vandermonde matrix converted to also be a quasimatrix: | |
+ L = cheb2quasi(legpoly(0:numCols-1, domain(A), 'norm', 1)); | |
+ [Q, R] = abstractQR(A, L, @innerProduct, @normest); | |
+ return | |
+end | |
- numFuns = numel(A.funs); | |
+% If A is an array-valued Chebfun, find the number of pieces. | |
+numFuns = numel(A.funs); | |
- % Perform QR on each piece: | |
- R = cell(numFuns, 1); | |
- Q = R; | |
- for k = 1:numFuns | |
- [Q{k}, R{k}] = qr(A.funs{k}); | |
- end | |
+% If there is only one piece (no breakpoints), call QR at the fun level. | |
+if ( numFuns == 1 ) | |
+ [Q, R] = qr(A.funs{1}); | |
+ Q = chebfun({Q}); | |
+ return | |
+end | |
- % Compute [Qhat, Rhat] = qr(Q): | |
- m = cellfun(@(v) size(v, 1), R); | |
- [Qhat, Rhat] = qr(cell2mat(R)); | |
- Rhat = Rhat(1:numCols,:); % Extract require entries. | |
- Qhat = mat2cell(Qhat(:,1:numCols), m, numCols); | |
- | |
- % Ensure diagonal has positive sign ... | |
- s = sign(diag(Rhat)); s(~s) = 1; | |
- S = spdiags(s, 0, numCols, numCols); | |
- R = S*Rhat; | |
- % ... and fold Qhat back in to Q: | |
- for k = 1:numFuns | |
- Q{k} = Q{k}*(Qhat{k,1}*S); | |
- end | |
- | |
- Q = chebfun(Q); | |
+% If there are multiple pieces, perform QR on each piece, stack the | |
+% R matrices from the separate pieces and perform QR on that, yielding | |
+% Qhat and Rhat, then break Qhat back down into sections and fold those | |
+% back into the matrices Q. | |
+ | |
+% Perform QR on each piece: | |
+R = cell(numFuns, 1); | |
+Q = R; | |
+for k = 1:numFuns | |
+ [Q{k}, R{k}] = qr(A.funs{k}); | |
+end | |
+% Compute [Qhat, Rhat] = qr(R): | |
+[Qhat, Rhat] = qr(cell2mat(R)); | |
+Rhat = Rhat(1:numCols,:); % Extract required entries. | |
+Qhat = Qhat(:,1:numCols); | |
+ | |
+% Ensure that the diagonal of Rhat is non-negative: | |
+neg = diag(Rhat) < 0; | |
+Rhat(neg,:) = -1 * Rhat(neg,:); | |
+Qhat(:,neg) = -1 * Qhat(:,neg); | |
+ | |
+% ... And fold Qhat back in to Q for each segment: | |
+m = cellfun(@(v) size(v, 1), R); % Number of rows per segment. | |
+Qhat = mat2cell(Qhat, m, numCols); | |
+for k = 1:numFuns | |
+ Q{k} = Q{k}*Qhat{k,1}; | |
end | |
+Q = chebfun(Q); | |
+R = Rhat; | |
+ | |
end |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Hello. I'm sorry. I'm weak in coding. Can you submit the code for the clenshaw method to my email?
Thank you very much
my email:
[email protected]