sdmeq.m

这是matlab解2阶锥工具包

function quiz = sdmeq(quiz,lmedata,Xindex,inL,inR,inK,kside);  
% SDMPB/SDMEQ - add term in some equality constraint
%  
% quiz = sdmeq(quiz,[lmeindex blrow blcol],Xindex,L,R);
%
%   if <R> is non empty adds a term : L*X*R to the (blrow,blcol) block
%
%   else if <R=[]> (or ommited) adds the term (L*X*L') to the (blrow,blcol) block
%
%   INPUTS: 
%     lmeindex = +n   : left-hand side of the n-th LME (left = right)
%              = -n   : right-hand side of the n-th LME (left = right)
%     blrow           : row index of the block 
%     blcol           : column index of the block 
%     Xindex   =  0   : constant term <L*R>
%              = +m   : term on the m-th variable <X>
%              = -m   : term on the transpose of the m-th variable <X.'>
%              = +i*m : term on the conjugate of the m-th variable <conj(X)>
%              = -i*m : term on the conjugate transpose of the m-th variable <X'>
%     L : matrix of appropriate dimensions, scalars accepted (default L=1 ).
%     R : matrix of appropriate dimensions, scalars accepted (default R=L').
%
%
% quiz = sdmeq(quiz,lmeindex,Xindex,L,R);
%
%   Same use, but the terms are added regardless to the block decomposition.
%   This can be used for LMEs composed of a single block, but also for any other LME. 
%   The terms (L*X*R) or (L*X*L') are added to the whole LME.
%   <L> and <R> must have respectively their number of rows and columns
%   equal to the LMI dimensions (sum of the blocks rows and columns)  
%     
%
% quiz = sdmeq(quiz,[lmeindex blrow blcol],Xindex,L,R,K);
%
%   same use, but adds terms such that :
%    (L*(K kron X)*R)  or  (L*(K kron X)*L')  
%   where (K kron X) is the Kronecker product of K and X  
%
% quiz = sdmeq(quiz,[lmeindex blrow blcol],Xindex,L,R,K,-1);
%
%   same use, but adds terms such that :
%    (L*(X kron K)*R)  or  (L*(X kron K)*L')  
%
%  
% Note that complex valued L, R and K matrices are accepted. 
%  
% SEE ALSO sdmvar, sdmineq, sdmlme 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<3
    error('not enough input arguments');
  elseif nargin>7
    error('too many input arguments');
  elseif ~isa(quiz,'sdmpb')
    error('1st input argument must be a SDMPB object');
  end;      
  %%% check the input <lmedata>
  if isempty(lmedata)
    error('2nd input argument (LME and block indexes) must contain at least the index of the LME equality')
  elseif min(size(lmedata))~=1
    error('2nd input argument (LME and blocks indexes) must be a scalar or a vector');
  elseif any(round(lmedata)-lmedata)
    error('2nd input argument (LME and blocks indexes) must be a vector of integers');
  end;
  lmeindex=lmedata(1);

  %%% check the input <lmeindex>
  if lmeindex<0
    side=-1;
  else
    side=1;
  end;  
  lmeindex=abs(lmeindex);
  if lmeindex>quiz.eq.nb
    error('2nd input argument exceeds the number of equality constraints');
  end;
  lmem=quiz.eq.m(lmeindex);
  lmeM=quiz.eq.M(lmeindex);
  nbrow=quiz.eq.row(lmeindex);
  nbcol=quiz.eq.col(lmeindex);
  
  %%% check the inputs <blrow> and <blcol> if any
  if length(lmedata)<2
    % first case : no block partitionning
    lmedata(2)=0.5;
    lmedata(3)=0.5;
    % these two values are only added to remind that if there is no block 
    % then it means that all terms write "on the diagonal"
    mrow = 1;
    Mrow = nbrow;
    mcol = 1;
    Mcol = nbcol;
    % these four values are set to indicate that the block has the full
    % size of the LME constraint
  else
    % second case : block partitionning
    if length(lmedata)>3
      warning('2nd input argument (LME and blocks indexes) is too long! Only the first 3 will be used');
    elseif length(lmedata)==2
      error('2nd input argument (LMI and blocks indexes) must be either a scalar or a vector of 3 elements');
    end;
    if any(lmedata(2:3) <= 0)
      error('2nd and 3rd parameters of the 2nd input argument (line and column indexes of the block) must be positive');
    end;
    % Check if blocks indexes are defined
    if lmedata(2)>size(quiz.eq.blockMrow{lmeindex})
      error('2nd parameter of the 2nd input argument (blocks'' row index) refers to a non-existing block.');
    elseif lmedata(3)>size(quiz.eq.blockMcol{lmeindex})
      error('3rd parameter of the 2nd input argument (bolcks'' column index) refers to a non-existing block.');
    end;
    % Definition of the position of the block in the LME rows and columns
    mrow = quiz.eq.blockmrow{lmeindex};    mrow = mrow(lmedata(2));
    Mrow = quiz.eq.blockMrow{lmeindex};    Mrow = Mrow(lmedata(2));
    mcol = quiz.eq.blockmcol{lmeindex};    mcol = mcol(lmedata(3));
    Mcol = quiz.eq.blockMcol{lmeindex};    Mcol = Mcol(lmedata(3));
  end;
 
  %%% check the input <Xindex>
  if ( ~isnumeric(Xindex) | (mod((abs(Xindex)),2)~=1 & mod((abs(Xindex)),2)~=0) | any(size(Xindex)-[1,1]))
    error('3rd input argument (index of the variable) does not fit the convensions');
  end;
  % check if the variable is transposed and/or conjugated
  if isreal(Xindex)
    conjugue=0;
    if Xindex>0, trans=0; else, trans=1; end;
  else
    conjugue=1;
    if imag(Xindex)>0, trans=0; else, trans=1; end;
  end
  Xindex=abs(Xindex);
  % check if the variable exists
  if Xindex>quiz.var.nb
    error('3rd input argument exceeds the number of variables');
  elseif Xindex>0
    varm=quiz.var.m(Xindex);
    varM=quiz.var.M(Xindex);
    varrow=quiz.var.row(Xindex);
    varcol=quiz.var.col(Xindex);    
  end;
  
  %%% check the input <L>
  if nargin<4, inL=1;  end;
  if isempty(inL), inL=1;  end;
  L=inL;
  if ~isnumeric(L)
    error('4th input argument (left multiplier) must be a matrix or a scalar');
  % if scalar convert into a matrix with row block dimension
  elseif max(size(L))==1
    L=L*speye(Mrow-mrow+1);
  % check the dimensions fit with the specified block
  elseif ~isequal(size(L,1),(Mrow-mrow+1))
    error('4th input argument (left multiplier) must have the same number of rows as the dimension of the specified block');
  end;
  % convert to a left multiplier regarless of blocks
  [ii,jj,s]=find(L);
  ii=ii+mrow-1;
  L=sparse(ii,jj,s,nbrow,size(L,2));

  %%% check the input <R>
  if nargin<5, inR=[]; end;
  % Default, if R is not defined
  if isempty(inR), 
    %%% in case <L> was not a scalar
    if max(size(inL))>1
      %%% check that the specified block is square
      if (Mrow-mrow)~=(Mcol-mcol)
	error('The specified block is not square : inconsistent with this syntax');
      end;
      %%% check that the variable is square
      if Xindex
	decvar=get(quiz,'vardec',Xindex);
	if size(decvar,1)~=size(decvar,2)
	  error('the variable referenced by the 3rd input argument is not square : inconsistent with this syntax');	
	end;
      end;    
    end;
    inR=inL';
  end;
  R=inR;
  if ~isnumeric(R)
    error('5th input argument (right multiplier) must be a matrix or a scalar');
  % convert into a matrix with column block dimension
  elseif max(size(R))==1
    R=R*speye(Mcol-mcol+1);  
  % check the dimensions fit with the specified block
  elseif ~isequal(size(R,2),(Mcol-mcol+1))
    error('5th input argument (right multiplier) must have the same number of columns as the dimension of the specified block');
  end;
  % convert to a right multiplier regardless of blocks
  [ii,jj,s]=find(R);
  jj=jj+mcol-1;
  R=sparse(ii,jj,s,size(R,1),nbcol);
 
  %%% check the inputs <inK> and <kside> (only for variable terms)
  if Xindex    
    if nargin<6, inK=[]; end;
    if isempty(inK), inK=1; end;
    K=inK;
    if ~isnumeric(K)
      error('6th input argument (matrix of the Kronecker product) must be a matrix');
    end;
    if nargin<7, kside=1; end;
    if abs(kside)~=1
      error('7th input argument (to invert Kronecker product) must be 1 or -1');
    end;
    % in case of a scalar variable convert the terms into
    % x*k*I=((k*I) kron x) to make it convenient with respect to the matrix dimensions
    if varrow==1 & varcol==1 & max(size(K))==1
      K=K*speye(size(L,2));
    end;
    [krow,kcol]=size(K);      
    %%% if L(X kron K)R , permute in L(K kron X)R
    if kside<0      
      L=L*sdmupq(varrow,krow);
      R=sdmupq(kcol,varcol)*R;
    end;
  end;
  
  %%% check that all multiplications and Kronecker products fit in size
  % for variable terms
  if Xindex
    % first in the case the variable is transposed
    if trans
      if size(L,2)~=varcol*krow
	if max(size(inK))==1
	  error('Term cannot be added: inconsistent number of columns of the variable');
	else
	  error('Term cannot be added: inconsistent number of rows in the resulting Kronecker product');
	end;
      elseif size(R,1)~=varrow*kcol
	if max(size(inK))==1
	  error('Term cannot be added: inconsistent number of rows of the variable');
	else
	  error('Term cannot be added: inconsistent number of columns in the resulting Kronecker product');
	end;
      end;
    % second in the case the variable is not transposed
    else
      if size(L,2)~=varrow*krow
	if max(size(inK))==1
	  error('Term cannot be added: inconsistent number of rows of the variable');
	else
	  error('Term cannot be added: inconsistent number of rows in the resulting Kronecker product');
	end;
      elseif size(R,1)~=varcol*kcol
	if max(size(inK))==1
	  error('Term cannot be added: inconsistent number of columns of the variable');
	else
	  error('Term cannot be added: inconsistent number of columns in the resulting Kronecker product');
	end;
      end;
    end;
  % for constant terms  
  else
    if size(L,2)~=size(R,1)
      error(['Number of columns of 4th input argument (left multiplier)' ...
	     ' does not fit the number of rows of 5th input argument (right multiplier)']);
      end;
  end;
  
  %%% All the tests are accomplished 
  %%% the remaining is generic
        
  %%% clear the previous solution if necessary
  quiz = sdmclear(quiz);
  
  %%% if constant term
  if Xindex==0
    sec=length(quiz.eq.c);
    [ii,jj,s]=find(side*sdmvec(L*R));
    ii=ii+lmem-1;
    add=sparse(ii,jj,s,sec,1);
    quiz.eq.c = quiz.eq.c-add;
  %%% if variable term
  else
    %%% define the term for each value of K(ii,jj)
    add=sparse(kcol*varcol*krow*varrow,varcol*varrow);
    for ii=1:krow
      ER=kron(K(ii,:),speye(varcol));
      EL=kron(sparse(ii,1,1,krow,1),speye(varrow));
      add=add+kron(ER.',EL);      
    end;
    add=side*kron(R.',L)*add;
    %%% if X.' , permute decision variables of X
    if trans
      add=add*sdmupq(varrow,varcol);
    end;
    %%% add the term to the lme 
    if conjugue
      seAtc=size(quiz.eq.Atconj);
      [ii,jj,s]=find(add);
      ii=ii+lmem-1;
      jj=jj+varm-1;
      add=sparse(ii,jj,s,seAtc(1),seAtc(2));
      quiz.eq.Atconj=quiz.eq.Atconj+add;
    else
      seAt=size(quiz.eq.At);
      [ii,jj,s]=find(add);
      ii=ii+lmem-1;
      jj=jj+varm-1;
      add=sparse(ii,jj,s,seAt(1),seAt(2));
      quiz.eq.At=quiz.eq.At+add;
    end;
  end;