Commit 5417b27a authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

Fixed indentation of matlab files.

parent 3d91cf97
...@@ -26,107 +26,107 @@ function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,ka ...@@ -26,107 +26,107 @@ function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,ka
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
k = size(DT,3); % number of structural parameters k = size(DT,3); % number of structural parameters
smpl = size(Y,2); % Sample size. smpl = size(Y,2); % Sample size.
pp = size(Y,1); % Maximum number of observed variables. pp = size(Y,1); % Maximum number of observed variables.
mm = size(T,2); % Number of state variables. mm = size(T,2); % Number of state variables.
a = zeros(mm,1); % State vector. a = zeros(mm,1); % State vector.
Om = R*Q*transpose(R); % Variance of R times the vector of structural innovations. Om = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
t = 0; % Initialization of the time index. t = 0; % Initialization of the time index.
oldK = 0; oldK = 0;
notsteady = 1; % Steady state flag. notsteady = 1; % Steady state flag.
F_singular = 1; F_singular = 1;
lik = zeros(smpl,1); % Initialization of the vector gathering the densities. lik = zeros(smpl,1); % Initialization of the vector gathering the densities.
LIK = Inf; % Default value of the log likelihood. LIK = Inf; % Default value of the log likelihood.
if nargout > 1 if nargout > 1
DLIK = zeros(k,1); % Initialization of the score. DLIK = zeros(k,1); % Initialization of the score.
end end
AHess = zeros(k,k); % Initialization of the Hessian AHess = zeros(k,k); % Initialization of the Hessian
Da = zeros(mm,k); % State vector. Da = zeros(mm,k); % State vector.
Dv = zeros(length(mf),k); Dv = zeros(length(mf),k);
% for ii = 1:k % for ii = 1:k
% DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii)); % DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
% end % end
while notsteady && t<smpl while notsteady && t<smpl
t = t+1; t = t+1;
v = Y(:,t)-a(mf); v = Y(:,t)-a(mf);
F = P(mf,mf) + H; F = P(mf,mf) + H;
if rcond(F) < kalman_tol if rcond(F) < kalman_tol
if ~all(abs(F(:))<kalman_tol) if ~all(abs(F(:))<kalman_tol)
return return
else
a = T*a;
P = T*P*transpose(T)+Om;
end
else else
F_singular = 0; a = T*a;
iF = inv(F); P = T*P*transpose(T)+Om;
K = P(:,mf)*iF;
lik(t) = log(det(F))+transpose(v)*iF*v;
[DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
for ii = 1:k
Dv(:,ii) = -Da(mf,ii) - DYss(mf,ii);
Da(:,ii) = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
if t>=start && nargout > 1
DLIK(ii,1) = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
end
end
vecDPmf = reshape(DP(mf,mf,:),[],k);
% iPmf = inv(P(mf,mf));
if t>=start
AHess = AHess + Dv'*iF*Dv + .5*(vecDPmf' * kron(iF,iF) * vecDPmf);
end
a = T*(a+K*v);
P = T*(P-K*P(mf,:))*transpose(T)+Om;
DP = DP1;
end end
notsteady = max(max(abs(K-oldK))) > riccati_tol; else
oldK = K; F_singular = 0;
end iF = inv(F);
K = P(:,mf)*iF;
lik(t) = log(det(F))+transpose(v)*iF*v;
if F_singular [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
error('The variance of the forecast error remains singular until the end of the sample')
end
for ii = 1:k
if t < smpl Dv(:,ii) = -Da(mf,ii) - DYss(mf,ii);
t0 = t+1; Da(:,ii) = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
while t < smpl if t>=start && nargout > 1
t = t+1; DLIK(ii,1) = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
v = Y(:,t)-a(mf); end
for ii = 1:k end
Dv(:,ii) = -Da(mf,ii)-DYss(mf,ii); vecDPmf = reshape(DP(mf,mf,:),[],k);
Da(:,ii) = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii)); % iPmf = inv(P(mf,mf));
if t>=start && nargout >1 if t>=start
DLIK(ii,1) = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v; AHess = AHess + Dv'*iF*Dv + .5*(vecDPmf' * kron(iF,iF) * vecDPmf);
end
end
if t>=start
AHess = AHess + Dv'*iF*Dv;
end
a = T*(a+K*v);
lik(t) = transpose(v)*iF*v;
end end
AHess = AHess + .5*(smpl-t0+1)*(vecDPmf' * kron(iF,iF) * vecDPmf); a = T*(a+K*v);
if nargout > 1 P = T*(P-K*P(mf,:))*transpose(T)+Om;
DP = DP1;
end
notsteady = max(max(abs(K-oldK))) > riccati_tol;
oldK = K;
end
if F_singular
error('The variance of the forecast error remains singular until the end of the sample')
end
if t < smpl
t0 = t+1;
while t < smpl
t = t+1;
v = Y(:,t)-a(mf);
for ii = 1:k for ii = 1:k
% DLIK(ii,1) = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) ); Dv(:,ii) = -Da(mf,ii)-DYss(mf,ii);
Da(:,ii) = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
if t>=start && nargout >1
DLIK(ii,1) = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
end
end
if t>=start
AHess = AHess + Dv'*iF*Dv;
end end
a = T*(a+K*v);
lik(t) = transpose(v)*iF*v;
end
AHess = AHess + .5*(smpl-t0+1)*(vecDPmf' * kron(iF,iF) * vecDPmf);
if nargout > 1
for ii = 1:k
% DLIK(ii,1) = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) );
end end
lik(t0:smpl) = lik(t0:smpl) + log(det(F)); end
% for ii = 1:k; lik(t0:smpl) = lik(t0:smpl) + log(det(F));
% for jj = 1:ii % for ii = 1:k;
% H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)')); % for jj = 1:ii
% end % H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
% end % end
end % end
end
AHess = -AHess;
AHess = -AHess;
if nargout > 1 if nargout > 1
DLIK = DLIK/2; DLIK = DLIK/2;
end end
...@@ -134,15 +134,15 @@ end ...@@ -134,15 +134,15 @@ end
lik = (lik + pp*log(2*pi))/2; lik = (lik + pp*log(2*pi))/2;
LIK = sum(lik(start:end)); % Minus the log-likelihood. LIK = sum(lik(start:end)); % Minus the log-likelihood.
% end of main function % end of main function
function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K) function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
k = size(DT,3); k = size(DT,3);
tmp = P-K*P(mf,:); tmp = P-K*P(mf,:);
for ii = 1:k for ii = 1:k
DF(:,:,ii) = DP(mf,mf,ii) + DH(:,:,ii); DF(:,:,ii) = DP(mf,mf,ii) + DH(:,:,ii);
DiF(:,:,ii) = -iF*DF(:,:,ii)*iF; DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
DK(:,:,ii) = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii); DK(:,:,ii) = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
Dtmp = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii); Dtmp = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
...@@ -150,6 +150,3 @@ for ii = 1:k ...@@ -150,6 +150,3 @@ for ii = 1:k
end end
% end of computeDKalman % end of computeDKalman
\ No newline at end of file
...@@ -29,7 +29,7 @@ function e = SPAimerr(c) ...@@ -29,7 +29,7 @@ function e = SPAimerr(c)
% Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3, % Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3,
% pages 472-489 % pages 472-489
if(c==1) e='Aim: unique solution.'; if(c==1) e='Aim: unique solution.';
elseif(c==2) e='Aim: roots not correctly computed by real_schur.'; elseif(c==2) e='Aim: roots not correctly computed by real_schur.';
elseif(c==3) e='Aim: too many big roots.'; elseif(c==3) e='Aim: too many big roots.';
elseif(c==35) e='Aim: too many big roots, and q(:,right) is singular.'; elseif(c==35) e='Aim: too many big roots, and q(:,right) is singular.';
......
function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ... function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
SPAmalg(h,neq,nlag,nlead,condn,uprbnd) SPAmalg(h,neq,nlag,nlead,condn,uprbnd)
% [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ... % [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
% SPAmalg(h,neq,nlag,nlead,condn,uprbnd) % SPAmalg(h,neq,nlag,nlead,condn,uprbnd)
% %
...@@ -8,9 +8,9 @@ function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ... ...@@ -8,9 +8,9 @@ function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
% roots. This procedure will fail if the companion matrix is % roots. This procedure will fail if the companion matrix is
% defective and does not have a linearly independent set of % defective and does not have a linearly independent set of
% eigenvectors associated with the big roots. % eigenvectors associated with the big roots.
% %
% Input arguments: % Input arguments:
% %
% h Structural coefficient matrix (neq,neq*(nlag+1+nlead)). % h Structural coefficient matrix (neq,neq*(nlag+1+nlead)).
% neq Number of equations. % neq Number of equations.
% nlag Number of lags. % nlag Number of lags.
...@@ -19,9 +19,9 @@ function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ... ...@@ -19,9 +19,9 @@ function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
% by numeric_shift and reduced_form. % by numeric_shift and reduced_form.
% uprbnd Inclusive upper bound for the modulus of roots % uprbnd Inclusive upper bound for the modulus of roots
% allowed in the reduced form. % allowed in the reduced form.
% %
% Output arguments: % Output arguments:
% %
% b Reduced form coefficient matrix (neq,neq*nlag). % b Reduced form coefficient matrix (neq,neq*nlag).
% rts Roots returned by eig. % rts Roots returned by eig.
% ia Dimension of companion matrix (number of non-trivial % ia Dimension of companion matrix (number of non-trivial
...@@ -57,7 +57,7 @@ function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ... ...@@ -57,7 +57,7 @@ function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
% pages 472-489 % pages 472-489
b=[];rts=[];ia=[];nexact=[];nnumeric=[];lgroots=[];aimcode=[]; b=[];rts=[];ia=[];nexact=[];nnumeric=[];lgroots=[];aimcode=[];
if(nlag<1 || nlead<1) if(nlag<1 || nlead<1)
error('Aim_eig: model must have at least one lag and one lead'); error('Aim_eig: model must have at least one lag and one lead');
end end
% Initialization. % Initialization.
...@@ -75,17 +75,17 @@ if (iq>qrows) ...@@ -75,17 +75,17 @@ if (iq>qrows)
end end
[a,ia,js] = SPBuild_a(h,qcols,neq); [a,ia,js] = SPBuild_a(h,qcols,neq);
if (ia ~= 0) if (ia ~= 0)
if any(any(isnan(a))) || any(any(isinf(a))) if any(any(isnan(a))) || any(any(isinf(a)))
disp('A is NAN or INF') disp('A is NAN or INF')
aimcode=63; aimcode=63;
return return
end end
[w,rts,lgroots,flag_trouble]=SPEigensystem(a,uprbnd,min(length(js),qrows-iq+1)); [w,rts,lgroots,flag_trouble]=SPEigensystem(a,uprbnd,min(length(js),qrows-iq+1));
if flag_trouble==1 if flag_trouble==1
disp('Problem in SPEIG'); disp('Problem in SPEIG');
aimcode=64; aimcode=64;
return return
end end
q = SPCopy_w(q,w,js,iq,qrows); q = SPCopy_w(q,w,js,iq,qrows);
end end
test=nexact+nnumeric+lgroots; test=nexact+nnumeric+lgroots;
......
...@@ -41,9 +41,9 @@ hs(:,left) = -hs(:,right)\hs(:,left); ...@@ -41,9 +41,9 @@ hs(:,left) = -hs(:,right)\hs(:,left);
a = zeros(qcols,qcols); a = zeros(qcols,qcols);
if(qcols > neq) if(qcols > neq)
eyerows = 1:qcols-neq; eyerows = 1:qcols-neq;
eyecols = neq+1:qcols; eyecols = neq+1:qcols;
a(eyerows,eyecols) = eye(qcols-neq); a(eyerows,eyecols) = eye(qcols-neq);
end end
hrows = qcols-neq+1:qcols; hrows = qcols-neq+1:qcols;
a(hrows,:) = hs(:,left); a(hrows,:) = hs(:,left);
...@@ -51,14 +51,14 @@ a(hrows,:) = hs(:,left); ...@@ -51,14 +51,14 @@ a(hrows,:) = hs(:,left);
% Delete inessential lags and build index array js. js indexes the % Delete inessential lags and build index array js. js indexes the
% columns in the big transition matrix that correspond to the % columns in the big transition matrix that correspond to the
% essential lags in the model. They are the columns of q that will % essential lags in the model. They are the columns of q that will
% get the unstable left eigenvectors. % get the unstable left eigenvectors.
js = 1:qcols; js = 1:qcols;
zerocols = sum(abs(a)) == 0; zerocols = sum(abs(a)) == 0;
while( any(zerocols) ) while( any(zerocols) )
a(:,zerocols) = []; a(:,zerocols) = [];
a(zerocols,:) = []; a(zerocols,:) = [];
js(zerocols) = []; js(zerocols) = [];
zerocols = sum(abs(a)) == 0; zerocols = sum(abs(a)) == 0;
end end
ia = length(js); ia = length(js);
...@@ -2,7 +2,7 @@ function q = SPCopy_w(q,w,js,iq,qrows) ...@@ -2,7 +2,7 @@ function q = SPCopy_w(q,w,js,iq,qrows)
% q = SPCopy_w(q,w,js,iq,qrows) % q = SPCopy_w(q,w,js,iq,qrows)
% %
% Copy the eigenvectors corresponding to the largest roots into the % Copy the eigenvectors corresponding to the largest roots into the
% remaining empty rows and columns js of q % remaining empty rows and columns js of q
% Author: Gary Anderson % Author: Gary Anderson
% Original file downloaded from: % Original file downloaded from:
...@@ -30,7 +30,7 @@ function q = SPCopy_w(q,w,js,iq,qrows) ...@@ -30,7 +30,7 @@ function q = SPCopy_w(q,w,js,iq,qrows)
if(iq < qrows) if(iq < qrows)
lastrows = iq+1:qrows; lastrows = iq+1:qrows;
wrows = 1:length(lastrows); wrows = 1:length(lastrows);
q(lastrows,js) = w(:,wrows)'; q(lastrows,js) = w(:,wrows)';
end end
...@@ -30,7 +30,7 @@ function [w,rts,lgroots,flag_trouble] = SPEigensystem(a,uprbnd,rowsLeft) ...@@ -30,7 +30,7 @@ function [w,rts,lgroots,flag_trouble] = SPEigensystem(a,uprbnd,rowsLeft)
% Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3, % Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3,
% pages 472-489 % pages 472-489
opts.disp=0; opts.disp=0;
% next block is commented out because eigs() intermitently returns different rts % next block is commented out because eigs() intermitently returns different rts
%try %try
% [w,d] = eigs(a',rowsLeft,'LM',opts); % [w,d] = eigs(a',rowsLeft,'LM',opts);
...@@ -39,24 +39,24 @@ opts.disp=0; ...@@ -39,24 +39,24 @@ opts.disp=0;
% [mag,k] = sort(-mag); % [mag,k] = sort(-mag);
% rts = rts(k); % rts = rts(k);
%catch %catch
%disp('Catch in SPE'); %disp('Catch in SPE');
%pause(0.5); %pause(0.5);
%aStr=datestr(clock); %aStr=datestr(clock);
%eval(['save ' regexprep(aStr,' ','') ' a']); %eval(['save ' regexprep(aStr,' ','') ' a']);
try try
[w,d]=eig(a'); [w,d]=eig(a');
catch catch
lasterr lasterr
w=[];rts=[];lgroots=[]; w=[];rts=[];lgroots=[];
flag_trouble=1; flag_trouble=1;
return return
end end
rts = diag(d); rts = diag(d);
mag = abs(rts); mag = abs(rts);
[mag,k] = sort(-mag); [mag,k] = sort(-mag);
rts = rts(k); rts = rts(k);
%end %end
flag_trouble=0; flag_trouble=0;
%ws=SPSparse(w); %ws=SPSparse(w);
ws=sparse(w); ws=sparse(w);
...@@ -65,7 +65,7 @@ ws = ws(:,k); ...@@ -65,7 +65,7 @@ ws = ws(:,k);
% Given a complex conjugate pair of vectors W = [w1,w2], there is a % Given a complex conjugate pair of vectors W = [w1,w2], there is a
% nonsingular matrix D such that W*D = real(W) + imag(W). That is to % nonsingular matrix D such that W*D = real(W) + imag(W). That is to
% say, W and real(W)+imag(W) span the same subspace, which is all % say, W and real(W)+imag(W) span the same subspace, which is all
% that aim cares about. % that aim cares about.
ws = real(ws) + imag(ws); ws = real(ws) + imag(ws);
......
...@@ -36,12 +36,11 @@ right = qcols+1:qcols+neq; ...@@ -36,12 +36,11 @@ right = qcols+1:qcols+neq;
zerorows = find( sum(abs( hs(:,right)' ))==0 ); zerorows = find( sum(abs( hs(:,right)' ))==0 );
while( any(zerorows) && iq <= qrows ) while( any(zerorows) && iq <= qrows )
nz = length(zerorows); nz = length(zerorows);
q(iq+1:iq+nz,:) = hs(zerorows,left); q(iq+1:iq+nz,:) = hs(zerorows,left);
hs(zerorows,:) = SPShiftright(hs(zerorows,:),neq); hs(zerorows,:) = SPShiftright(hs(zerorows,:),neq);
iq = iq + nz; iq = iq + nz;
nexact = nexact + nz; nexact = nexact + nz;
zerorows = find( sum(abs( hs(:,right)' ))==0 ); zerorows = find( sum(abs( hs(:,right)' ))==0 );
end end
h=full(hs); h=full(hs);
...@@ -37,14 +37,14 @@ right = qcols+1:qcols+neq; ...@@ -37,14 +37,14 @@ right = qcols+1:qcols+neq;
zerorows = find( abs(diag(R)) <= condn ); zerorows = find( abs(diag(R)) <= condn );
while( any(zerorows) && iq <= qrows ) while( any(zerorows) && iq <= qrows )
h=sparse(h); h=sparse(h);
Q=sparse(Q); Q=sparse(Q);
h = Q'*h; h = Q'*h;
nz = length(zerorows); nz = length(zerorows);
q(iq+1:iq+nz,:) = h(zerorows,left); q(iq+1:iq+nz,:) = h(zerorows,left);
h(zerorows,:) = SPShiftright( h(zerorows,:), neq ); h(zerorows,:) = SPShiftright( h(zerorows,:), neq );
iq = iq + nz; iq = iq + nz;
nnumeric = nnumeric + nz; nnumeric = nnumeric + nz;
[Q,R,E] = qr( full(h(:,right)) ); [Q,R,E] = qr( full(h(:,right)) );
zerorows = find( abs(diag(R)) <= condn ); zerorows = find( abs(diag(R)) <= condn );
end end
...@@ -2,7 +2,7 @@ function scof = SPObstruct(cof,cofb,neq,nlag,nlead) ...@@ -2,7 +2,7 @@ function scof = SPObstruct(cof,cofb,neq,nlag,nlead)
% scof = SPObstruct(cof,cofb,neq,nlag,nlead) % scof = SPObstruct(cof,cofb,neq,nlag,nlead)
% %
% Construct the coefficients in the observable structure. % Construct the coefficients in the observable structure.
% %
% Input arguments: % Input arguments:
% %
% cof structural coefficients % cof structural coefficients
...@@ -51,17 +51,17 @@ qs=sparse(q); ...@@ -51,17 +51,17 @@ qs=sparse(q);
qs(1:rc,1:cc) = sparse(cofb); qs(1:rc,1:cc) = sparse(cofb);
qcols = neq*(nlag+nlead); qcols = neq*(nlag+nlead);
if( nlead > 1 ) if( nlead > 1 )
for i = 1:nlead-1 for i = 1:nlead-1
rows = i*neq + (1:neq); rows = i*neq + (1:neq);
qs(rows,:) = SPShiftright( qs((rows-neq),:), neq ); qs(rows,:) = SPShiftright( qs((rows-neq),:), neq );
end end
end end
l = (1: neq*nlag); l = (1: neq*nlag);
r = (neq*nlag+1: neq*(nlag+nlead)); r = (neq*nlag+1: neq*(nlag+nlead));
qs(:,l) = - qs(:,r) \ qs(:,l); qs(:,l) = - qs(:,r) \ qs(:,l);
minus = 1: neq*(nlag+1); minus = 1: neq*(nlag+1);
plus = neq*(nlag+1)+1: neq*(nlag+1+nlead); plus = neq*(nlag+1)+1: neq*(nlag+1+nlead);
......
...@@ -38,7 +38,7 @@ if(nonsing) ...@@ -38,7 +38,7 @@ if(nonsing)
b = qs(1:neq,1:bcols); b = qs(1:neq,1:bcols);
b = full(b); b = full(b);
else %rescale by dividing row by maximal qr element else %rescale by dividing row by maximal qr element
%'inverse condition number small, rescaling ' %'inverse condition number small, rescaling '
themax=max(abs(qs(:,right)),[],2); themax=max(abs(qs(:,right)),[],2);
oneover = diag(1 ./ themax); oneover = diag(1 ./ themax);
nonsing = rcond(full(oneover *qs(:,right))) > condn; nonsing = rcond(full(oneover *qs(:,right))) > condn;
...@@ -48,4 +48,3 @@ else %rescale by dividing row by maximal qr element ...@@ -48,4 +48,3 @@ else %rescale by dividing row by maximal qr element
b = full(b); b = full(b);
end