sdmvar.m

这是matlab解2阶锥工具包

function [quiz,Xindex] = sdmvar(quiz,nbrow,nbcol,name)
% SDMPB/SDMVAR - declare a new matrix variable X
%
% [quiz,Xindex] = sdmvar(quiz,nbrow,nbcol,name);
%  
%   <quiz> is a variable describing a SDMPB object
%  
%   <nbrow> and <nbcol> are the dimensions of the matrix variable X.
%   When <nbcol> is ommitted, the variable is an <nbrow> square matrix.  
%  
%   Default the matrix variable has real entries. If the argument <nbrow>
%   is imaginary then the matrix variable is complex valued.
%   For example: sdmvar(quiz,2*i,3,name), where i stands for i=sqrt(-1),
%   declares a 2-by-3 complex valued matrix variable.
%  
%   Exceptions are if  
%   <nbcol='d'>  
%     X is square diagonal (and complex if <nbrow> is imaginary)
%   <nbcol='s'>  
%     X is square symmetric (and complex if <nbrow> is imaginary)
%   <nbcol='h'>  
%     X is square Hermitian (coherent only if <nbrow> is imaginary)
%   <nbcol='as'> 
%     X is square anti-symmetric (and complex if <nbrow> is imaginary)
%   <nbcol='ah'> 
%     X is square anti-Hermitian (coherent only if <nbrow> is imaginary)
%   <nbcol='st'> 
%     X has the same size as <nbrow> and depends on existing decision
%     variables or new decision variables s.t. (where i=sqrt(-1))
%       nbrow(k,j)=  0     if X(k,j)=0
%       nbrow(k,j)= +n     if X(k,j) real    = n-th decision variable  
%       nbrow(k,j)= -n     if X(k,j) real    = (-1)*n-th decision variable
%       nbrow(k,j)= +n+n*i if X(k,j) complex = n-th decision variable  
%       nbrow(k,j)= -n-n*i if X(k,j) complex = (-1)*n-th decision variable  
%       nbrow(k,j)= +n-n*i if X(k,j) complex = conjugate n-th decision variable  
%       nbrow(k,j)= -n-n*i if X(k,j) complex = (-1)*conjugate n-th decision variable  
%
%   <name> is an optional string describing the variable
%  
%   <Xindex> index of the new matrix variable 
%            =k if k-1 matrix variables have been declared so far.
%
% SEE ALSO sdmshow, sdmset, sdmlmi, sdmobj and sdmsol.

%   This file is part of SeDuMi Interface 1.04 (JUL2002)
%   Last Update 6th September 2002
%   Copyright (C) 2002 Dimitri Peaucelle & Krysten Taitz
%   LAAS-CNRS, Toulouse, France
% 
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
% 
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
% 
%   You should have received a copy of the GNU General Public License
%   along with this program; if not, write to the Free Software
%   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 
  
  %%% check inputs
  if nargin<2
    error('not enough input arguments');
  elseif nargin>4
    error('too many input arguments');
  elseif ~isa(quiz,'sdmpb')
    error('1st input argument must be a SDMPB object');
  elseif nargin<3
    nbcol=abs(nbrow);
    name='NO NAME';
  elseif nargin<4
    name='NO NAME';
  end;
  if isempty(nbcol)
    nbcol=abs(nbrow);
  end;
  %%% clear the previous solution if necessary
  quiz = sdmclear(quiz);
  %%% create new variable
  quiz.var.nb=quiz.var.nb+1;
  Xindex=quiz.var.nb;
  if isempty(quiz.var.M)
    previousM=0;
  else
    previousM=quiz.var.M(Xindex-1);
  end;
  %%% Index of the last decision variable
  startdecnb=get(quiz,'vardecnb');
  %%% initialize ycomplex defining the complex decision variables
  ycomplex=[];
  %%% dimensions and structure of the variable
  nl=abs(nbrow);
  % if diagonal
  if isequal(nbcol,'d') 
    nc=nl;
    struc=sparse(nl*nc,nl);
    strucconj=sparse(nl*nc,nl);
    struc=sparse((1:nc)'*(nc+1)-nc,(1:nc)',ones(size((1:nc)')));
    struc=[sparse(size(struc,1),size(quiz.var.struc,2)),struc];
    strucconj=[sparse(size(struc,1),size(quiz.var.struc,2)),strucconj];
    if ~isreal(nbrow)
      ycomplex=1:nl;
    end;
    quiz.K.ycomplex=[quiz.K.ycomplex,startdecnb+ycomplex];
  % if symmetric
  elseif isequal(nbcol,'s') 
    nc=nl;
    struc=sparse(nl*nc,nl*(nl+1)/2);
   
    strucconj=sparse(nl*nc,nl*(nl+1)/2);
    Rindex=0;
    for jj=1:nc
      Rindex=Rindex+1;
      struc((nc+1)*jj-nc,Rindex)=1;
      for ii=(jj+1):nl
	Rindex=Rindex+1;
	struc(nl*(jj-1)+ii,Rindex)=1;
	struc(nl*(ii-1)+jj,Rindex)=1;
      end;
    end;
    struc=[sparse(size(struc,1),size(quiz.var.struc,2)),struc];
    strucconj=[sparse(size(struc,1),size(quiz.var.struc,2)),strucconj];
    if ~isreal(nbrow)
      ycomplex=1:(nl*(nl+1)/2);
    end;
    quiz.K.ycomplex=[quiz.K.ycomplex,startdecnb+ycomplex];
  % if Hermitian
  elseif isequal(nbcol,'h') 
    if isreal(nbrow)
      error(['2nd input argument must be imaginary: only complex matrices' ...
	     ' are Hermitian']);
    end;
    nc=nl;
    struc=sparse(nl*nc,nl*(nl+1)/2);
    strucconj=sparse(nl*nc,nl*(nl+1)/2);
    Rindex=0;
    for jj=1:nc
      Rindex=Rindex+1;
      struc((nc+1)*jj-nc,Rindex)=1;
      for ii=(jj+1):nl
	Rindex=Rindex+1;
	struc(nl*(jj-1)+ii,Rindex)=1;
	strucconj(nl*(ii-1)+jj,Rindex)=1;
	ycomplex=[ycomplex,Rindex];
      end;
    end;
    struc=[sparse(size(struc,1),size(quiz.var.struc,2)),struc];
    strucconj=[sparse(size(struc,1),size(quiz.var.struc,2)),strucconj];
    quiz.K.ycomplex=[quiz.K.ycomplex,startdecnb+ycomplex];
  % if anti-symmetric
  elseif isequal(nbcol,'as') 
    nc=nl;
    struc=sparse(nl*nc,nl*(nl-1)/2);
    strucconj=sparse(nl*nc,nl*(nl-1)/2);
    Rindex=0;
    for jj=1:nc
      for ii=(jj+1):nl
	Rindex=Rindex+1;
	struc(nl*(jj-1)+ii,Rindex)=1;
	struc(nl*(ii-1)+jj,Rindex)=-1;
      end;
    end;
    struc=[sparse(size(struc,1),size(quiz.var.struc,2)),struc];
    strucconj=[sparse(size(struc,1),size(quiz.var.struc,2)),strucconj];
    if ~isreal(nbrow)
      ycomplex=1:(nl*(nl-1)/2);
    end;
    quiz.K.ycomplex=[quiz.K.ycomplex,startdecnb+ycomplex];
  % if anti-Hermitian
  elseif isequal(nbcol,'ah') 
    if isreal(nbrow)
      error(['2nd input argument must be imaginary: only complex matrices' ...
	     ' are anti-Hermitian']);
    end;
    nc=nl;
    struc=sparse(nl*nc,nl*(nl-1)/2);
    strucconj=sparse(nl*nc,nl*(nl-1)/2);
    Rindex=0;
    for jj=1:nc
      for ii=(jj+1):nl
	Rindex=Rindex+1;
	struc(nl*(jj-1)+ii,Rindex)=1;
	strucconj(nl*(ii-1)+jj,Rindex)=-1;
      end;
    end;
    struc=[sparse(size(struc,1),size(quiz.var.struc,2)),struc];
    strucconj=[sparse(size(struc,1),size(quiz.var.struc,2)),strucconj];
    ycomplex=1:(nl*(nl-1)/2);
    quiz.K.ycomplex=[quiz.K.ycomplex,startdecnb+ycomplex];
  % if structured
  elseif isequal(nbcol,'st')
    Xdec=nbrow;
    ndec=max(startdecnb,max(max(abs(real(Xdec)))));
    [nl,nc]=size(Xdec);
    struci=[];
    strucj=[];
    strucs=[];
    strucconji=[];
    strucconjj=[];
    strucconjs=[];
    maxdecnb=startdecnb;
    for jj=1:nc
      for ii=1:nl
	if Xdec(ii,jj)~=0
	  rXdecij=real(Xdec(ii,jj));
	  iXdecij=imag(Xdec(ii,jj));
	  if ((mod(rXdecij,2)~=1 & mod(rXdecij,2)~=0) | (mod(iXdecij,2)~=1 ...
							   & mod(iXdecij,2)~=0))
	    errormsg=sprintf(['the (%i,%i) element of the 2nd input' ...
			       ' argument must be a real or complex integer'],ii,jj);
	    error(errormsg);
	  elseif (iXdecij & abs(rXdecij)~=abs(iXdecij))
	    errormsg=sprintf(['inconsistent (%i,%i) element of the 2nd input argument'],ii,jj);
	    error(errormsg);
	  elseif abs(rXdecij)<=maxdecnb
	    if isempty(quiz.K.ycomplex)
	      Kycomplex=0;
	    else
	      Kycomplex=quiz.K.ycomplex;
	    end;
	    if ( iXdecij & all(abs(rXdecij)~=Kycomplex) )
	      errormsg=sprintf(['the (%i,%i) element of the 2nd input' ...
				 ' argument must be real'],ii,jj);
	      error(errormsg);	    
	    elseif ( ~iXdecij & any(abs(rXdecij)==Kycomplex) )
	      errormsg=sprintf(['the (%i,%i) element of the 2nd input' ...
				 ' argument must be complex'],ii,jj);
	      error(errormsg);	    	    
	    end;
	  elseif (iXdecij & abs(rXdecij)>maxdecnb)
	    quiz.K.ycomplex=[quiz.K.ycomplex,abs(iXdecij)];	  
	    maxdecnb=abs(iXdecij);
	  end;
	  if rXdecij==-iXdecij
	    strucconji=[strucconji;nl*(jj-1)+ii];
	    strucconjj=[strucconjj;abs(rXdecij)];
	    strucconjs=[strucconjs;sign(rXdecij)];	  
	  else
	    struci=[struci;nl*(jj-1)+ii];
	    strucj=[strucj;abs(rXdecij)];
	    strucs=[strucs;sign(rXdecij)];
	  end;
	end;
	struc=sparse(struci,strucj,strucs,nl*nc,ndec);
	strucconj=sparse(strucconji,strucconjj,strucconjs,nl*nc,ndec);
      end;        
    end;
  % if full (rectangular) bloc
  else
    nl=abs(nbrow);
    nc=abs(nbcol);
    if ~isreal(nbrow)
      ycomplex=1:(nl*nc);
    end;
    quiz.K.ycomplex=[quiz.K.ycomplex,startdecnb+ycomplex];
    struc=speye(nl*nc);
    struc=[sparse(size(struc,1),size(quiz.var.struc,2)),struc];
    strucconj=sparse(nl*nc,nl*nc);
    strucconj=[sparse(size(struc,1),size(quiz.var.struc,2)),strucconj];
  end;  
  quiz.var.struc=[quiz.var.struc,...
		 sparse(size(quiz.var.struc,1),size(struc,2)-size(quiz.var.struc,2));...
		 struc];  
  quiz.var.strucconj=[quiz.var.strucconj,...
		     sparse(size(quiz.var.strucconj,1),size(strucconj,2)-size(quiz.var.strucconj,2));...
		     strucconj];  
  %%% define some fields used in other functions
  quiz.var.name{Xindex}=name;
  quiz.var.value{Xindex}=[];
  quiz.var.m=[quiz.var.m,previousM+1];
  quiz.var.M=[quiz.var.M,previousM+nl*nc];  
  quiz.var.row=[quiz.var.row,nl];
  quiz.var.col=[quiz.var.col,nc];
 
 
  %%% augment the matrices At and bt accordingly to the number of new
  %decision variables 
  newdecnb=get(quiz,'vardecnb');
  quiz.obj.bt=[quiz.obj.bt,sparse(1,nl*nc)];
  quiz.obj.btconj=[quiz.obj.btconj,sparse(1,nl*nc)];
  quiz.ineq.At=[quiz.ineq.At,sparse(size(quiz.ineq.At,1),nl*nc)];
  quiz.ineq.Atconj=[quiz.ineq.Atconj,sparse(size(quiz.ineq.Atconj,1),nl*nc)];
  quiz.eq.At=[quiz.eq.At,sparse(size(quiz.eq.At,1),nl*nc)];
  quiz.eq.Atconj=[quiz.eq.Atconj,sparse(size(quiz.eq.Atconj,1),nl*nc)];