001 function [x0,pmx]=mindist3(g,b,u)
002
003
004
005
006
007
008
009 g=g(:);
010 b=b(:);
011 d0=length(g);
012 if length(b)~=d0
013 error('g and b should have the same length')
014 end
015
016 if u==0
017 x0=zeros(d0,1);
018 pmx=zeros(d0,1);
019 return
020 end
021
022 nulldim=(g==0)&(b==0);
023
024 if sum(nulldim)==d0
025 if u==0
026 x0=zeros(d0,1);
027 pmx=zeros(d0,1);
028 return
029 else
030 error('u is out of range')
031 end
032 end
033
034 g(nulldim)=[];
035 b(nulldim)=[];
036 d=length(b);
037 if (all(g>0)|all(g<0))&u==-.25*sum(b.^2./g)
038 x00=-.5*b./g;
039 pmx=zeros(d,1);
040 x0=zeros(d0,1);
041 x0(~nulldim)=x00;
042 end
043 if all(g>0)&u<-.25*sum(b.^2./g)
044 error('u out of range')
045 end
046 if all(g<0)&u>-.25*sum(b.^2./g)
047 error('u out of range')
048 end
049 pmx0=zeros(d,1);
050 x00=[inf;zeros(d-1,1)];
051 I0=b==0;
052 I1=b~=0;
053 gI0=g(I0);
054 gI0unique=unique(gI0);
055 for k=1:length(gI0unique)
056 if all(g(I1)-gI0unique(k))
057 xp=zeros(d,1);
058 xp(I1)=.5*b(I1)./(gI0unique(k)-g(I1));
059 up=u-sum(b.*xp+g.*xp.^2,1);
060 if sign(up)*sign(gI0unique(k))>=0
061 xp(min(find(g==gI0unique(k))))=sqrt(up/gI0unique(k));
062 if sum(xp.^2)<sum(x00.^2)
063 x00=xp;
064 pmx0=zeros(d,1);
065 pmx0(g==gI0unique(k))=1;
066 end
067 end
068 end
069 end
070
071 bet=b(I1);
072 gam=g(I1);
073 if sum(I1)>0
074 R=[bet(1)^2*[0 2 -gam(1)];1 -2*gam(1) gam(1)^2];
075 for k=2:sum(I1);
076 R=rfadd(R,[bet(k)^2*[0 2 -gam(k)];1 -2*gam(k) gam(k)^2]);
077 end
078 lam=roots(R(1,:)-4*u*R(2,:));
079 lam=real(lam(abs(real(lam))>(1e6)*abs(imag(lam))));
080 x1=zeros(d,1);
081 for k=1:length(lam)
082 x1(I1)=.5*bet./(lam(k)-gam);
083 if sum(x1.^2)<sum(x00.^2);
084 x00=x1;
085 pmx0=zeros(size(pmx0));
086 end
087 end
088 end
089
090
091 x0=zeros(d0,1);
092 x0(~nulldim)=x00;
093 pmx=zeros(d0,1);
094 pmx(~nulldim)=pmx0;
095