www.gusucode.com > rctobsolete 工具箱 matlab源码程序 > rctobsolete/mutools/subs/ruqvsolb.m
% function [q,gopt] = ruqvsolb(r0,r1,uu,vv) % % Copyright 1991-2004 MUSYN Inc. and The MathWorks, Inc. function [maxlam] = ruqvsolb(r0,r1,uu,vv) [nrr ncr] = size(r0); [nru ncu] = size(uu); [nrv ncv] = size(vv); if nrr ~= nru | ncr ~= ncv | nru < ncu | ncv < nrv error('Dimension problems') return end [Uu,Su,Vu] = svd(uu); utol = max(size(uu))*max(max(Su))*eps; rku = sum(diag(Su)>utol); Vu1 = Vu(:,1:rku); Su1 = Su(1:rku,1:rku); Ueq = Uu(:,1:rku); if rku<nru uperp = Uu(:,rku+1:nru); else uperp = []; end [Uv,Sv,Vv] = svd(vv); vtol = max(size(vv))*max(max(Sv))*eps; rkv = sum(diag(Sv)>vtol); Uv1 = Uv(:,1:rkv); Sv1 = Sv(1:rkv,1:rkv); Veq = Vv(:,1:rkv)'; if rkv<ncv vperp = Vv(:,rkv+1:ncv)'; else vperp = []; end if ~isempty(uperp) upr0 = uperp'*r0; upr1 = uperp'*r1; b = [zeros(nru-ncu,nru-ncu) upr0;upr0' zeros(ncr,ncr)]; a = -[eye(nru-ncu) upr1;upr1' eye(ncr)]; lam = eig(a,b); keep = find(abs(imag(lam))<1e-8); keep = 1:length(lam); vals = real(lam(keep)); vals = sort(vals); keep = find(~isinf(vals) & ~isnan(vals)); vals = vals(keep); lowu = []; highu = []; cnt = 0; tmpu = zeros(length(vals),1); for i=1:length(vals) tmpu(i) = norm(uperp'*(vals(i)*r0+r1)); if abs(norm(uperp'*(vals(i)*r0+r1))-1)<1e-6 if isempty(lowu) lowu = vals(i); cnt = cnt + 1; else highu = vals(i); cnt = cnt + 1; end end end if cnt>2 % disp('Warning in RUQVSOLB: Too many U-Points found'); % tmpu end else highu = inf; end if ~isempty(vperp) r0vp = r0*vperp'; r1vp = r1*vperp'; b = [zeros(nrr,nrr) r0vp;r0vp' zeros(ncv-nrv,ncv-nrv)]; a = -[eye(nrr) r1vp;r1vp' eye(ncv-nrv)]; lam = eig(a,b); keep = find(abs(imag(lam))<1e-8); keep = 1:length(lam); vals = real(lam(keep)); vals = sort(vals); keep = find(~isinf(vals) & ~isnan(vals)); vals = vals(keep); lowv = []; highv = []; cnt = 0; tmpv = zeros(length(vals),1); for i=1:length(vals) tmpv(i) = norm((vals(i)*r0+r1)*vperp'); if abs(norm((vals(i)*r0+r1)*vperp')-1)<1e-6 if isempty(lowv) lowv = vals(i); cnt = cnt + 1; else highv = vals(i); cnt = cnt + 1; end end end if cnt>2 % disp('Warning in RUQVSOLB: Too many V-Points found'); % tmpv end else highv = inf; end if ~isempty(highu) & ~isempty(highv) if isinf(highu) maxlam = highv; elseif isinf(highv) maxlam = highu; else if highu<=highv & highu>=lowv maxlam = highu; elseif highv<=highu & highv>=lowu maxlam = highv; else maxlam = []; end end end