source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/third_party/opt/yalmip/@sdpvar/plus.m @ 37

Last change on this file since 37 was 37, checked in by (none), 14 years ago

Added original make3d

File size: 5.3 KB
Line 
1function y = plus(X,Y)
2%PLUS (overloaded)
3
4% Author Johan Löfberg
5% $Id: plus.m,v 1.20 2006/07/26 20:17:58 joloef Exp $
6
7X_is_spdvar = isa(X,'sdpvar');
8Y_is_spdvar = isa(Y,'sdpvar');
9
10% Convert block objects
11if ~X_is_spdvar
12    if isa(X,'blkvar')
13        X = sdpvar(X);
14        X_is_spdvar = isa(X,'sdpvar');
15    elseif ~isa(X,'double')
16        error(['Cannot add SDPVAR object and ' upper(class(X)) ' object']);
17    end
18end
19
20if ~Y_is_spdvar
21    if isa(Y,'blkvar')
22        Y = sdpvar(Y);
23        Y_is_spdvar = isa(Y,'sdpvar');;
24    elseif ~isa(Y,'double')
25         error(['Cannot add SDPVAR object and ' upper(class(Y)) ' object']);
26    end
27end
28
29
30switch 2*X_is_spdvar+Y_is_spdvar
31    case 1
32        if isempty(X)
33            try
34                y = full(X - reshape(Y.basis(:,1),Y.dim(1),Y.dim(2)));
35            catch
36                error(lasterr);
37            end
38            return
39        end
40
41        y = Y;
42        n_Y = Y.dim(1);
43        m_Y = Y.dim(2);
44        [n_X,m_X] = size(X);
45        x_isscalar = (n_X*m_X==1);
46        y_isscalar = (n_Y*m_Y==1);
47        any_scalar = x_isscalar | y_isscalar;
48       
49        if x_isscalar & y_isscalar           
50             y.basis(1) = y.basis(1)+X;
51             % Reset info about conic terms
52             y.conicinfo = [0 0];
53             return
54         end
55         
56        if any_scalar | ([n_Y m_Y]==[n_X m_X])
57            if y_isscalar
58                y.basis = repmat(y.basis,n_X*m_X,1);
59                y.dim(1) = n_X;
60                y.dim(2) = m_X;
61            end
62            y.basis(:,1) = y.basis(:,1)+X(:);
63        else
64            error('Matrix dimensions must agree.');
65        end
66        % Reset info about conic terms
67        y.conicinfo = [0 0];
68
69    case 2
70
71        if isempty(Y)
72            try
73                y = full(reshape(X.basis(:,1),X.dim(1),X.dim(2))-Y);
74            catch
75                error(lasterr);
76            end
77            return
78        end
79
80        y = X;
81        n_X = X.dim(1);
82        m_X = X.dim(2);
83        [n_Y,m_Y] = size(Y);
84        x_isscalar = (n_X*m_X==1);
85        y_isscalar = (n_Y*m_Y==1);
86        any_scalar = x_isscalar | y_isscalar;
87       
88         % Special special case...
89         if x_isscalar & y_isscalar
90             y.basis(1) = y.basis(1)+Y;
91             % Reset info about conic terms
92             y.conicinfo = [0 0];
93             return
94         end
95         
96        if any_scalar | ([n_Y m_Y]==[n_X m_X])
97            if x_isscalar
98                y.basis = repmat(y.basis,n_Y*m_Y,1);
99                y.dim(1) = n_Y;
100                y.dim(2) = m_Y;
101            end
102            y.basis(:,1) = y.basis(:,1)+Y(:);
103        else
104            error('Matrix dimensions must agree.');
105        end
106        % Reset info about conic terms
107        y.conicinfo = [0 0];
108
109    case 3
110
111        n_X = X.dim(1);
112        m_X = X.dim(2);
113        n_Y = Y.dim(1);
114        m_Y = Y.dim(2);
115        x_isscalar = (n_X*m_X==1);
116        y_isscalar = (n_Y*m_Y==1);
117        any_scalar = x_isscalar | y_isscalar;
118
119        if (~((n_X==n_Y) & (m_X==m_Y))) & ~any_scalar
120            error('Matrix dimensions must agree.')
121        end
122
123        all_lmi_variables = uniquestripped([X.lmi_variables Y.lmi_variables]);
124        y = X;
125        X.basis = [];
126        y.lmi_variables = all_lmi_variables;
127
128        % ismembc faster (buggy?)
129        in_X = find(ismembc(all_lmi_variables,X.lmi_variables));
130        in_Y = find(ismembc(all_lmi_variables,Y.lmi_variables));
131        % in_X = find(ismember(all_lmi_variables,X.lmi_variables));
132        % in_Y = find(ismember(all_lmi_variables,Y.lmi_variables));
133
134        if isequal(X.lmi_variables,Y.lmi_variables) & n_Y==n_X & m_Y==m_X
135            y.basis = y.basis + Y.basis;
136             if length(X.lmi_variables)==1
137                 if all(y.basis(:,2)==0)
138                     y = full(y.basis(1));
139                 else
140                     % Reset info about conic terms
141                     y.conicinfo = [0 0];
142                 end
143                return
144            end
145        else
146            if 1
147                [ix,jx,sx] = find(y.basis);y.basis = [];
148                [iy,jy,sy] = find(Y.basis);Y.basis = [];
149                mapX = [1 1+in_X];
150                mapY = [1 1+in_Y];
151                basis_X = sparse(ix,mapX(jx),sx,n_X*m_X,1+length(all_lmi_variables));ix=[];jx=[];sx=[];
152                basis_Y = sparse(iy,mapY(jy),sy,n_Y*m_Y,1+length(all_lmi_variables));iy=[];jy=[];sy=[];
153            else
154                % MATLAB sparse fails on this for huge problems at a certain size
155                basis_X = spalloc(n_X*m_X,1+length(all_lmi_variables),nnz(X.basis));
156                basis_Y = spalloc(n_Y*m_Y,1+length(all_lmi_variables),nnz(Y.basis));
157                basis_X(:,[1 1+in_X])=y.basis;y.basis = [];
158                basis_Y(:,[1 1+in_Y])=Y.basis;Y.basis = [];
159            end
160
161            % Fix addition of matrix+scalar
162            if n_X*m_X<n_Y*m_Y
163                y.dim(1) = n_Y;
164                y.dim(2) = m_Y;
165                basis_X = repmat(basis_X,n_Y*m_Y,1);
166            end
167            if n_Y*m_Y<n_X*m_X
168                y.dim(1) = n_X;
169                y.dim(2) = m_X;
170                basis_Y = repmat(basis_Y,n_X*m_X,1);
171            end
172            % OK, solution is...
173            y.basis = basis_X;basis_X = [];
174            y.basis = y.basis+basis_Y;basis_Y = [];
175        end
176        % Reset info about conic terms
177        y.conicinfo = [0 0];
178        y = clean(y);
179
180    otherwise
181end
182
183
Note: See TracBrowser for help on using the repository browser.