1function SJplot(matrix,tol,option) 2% SJplot(matrix,tol,option) creates a plot of the calculated singular 3% value spectrum or partial singular value spectrum for a matrix from 4% the SJsingular database. 5% 6% Input: 7% matrix -- either: 8% a problem structure downloaded from the SJsingular database 9% a number: 1 to the # of matrices in the collection 10% a string: the name of a matrix in the collection 11% In the last two cases SJget must be installed first from 12% http://www.math.sjsu.edu/singular/matrices/SJget.html 13% tol (optional) -- draw a horizontal line at tol 14% (the default value is max(size(A))*eps(norm of A)) 15% option (optional) -- if option = 1 create high resolution image, 16% if option = 2 create a low resolution image, suitable for 17% thumbnail picture. The default is option = 1. 18% Note: if blue error bars appear in the plot then a singular 19% value of A is guaranteed to be in the interval pictured by the 20% blue bars around each of the calculated singular values. 21% Examples: 22% %if matrix information has already been downloaded 23% SJplot(Problem), shg 24% % or (assuming that SJget is installed) 25% SJplot('Meszaros/model6'), shg 26% % or 27% tol = 1.e-5; 28% SJplot(403, tol), shg 29% 30% See also SJrank. 31 32% Written by Nikolay Botev and Leslie Foster, 9/5/2008, Version 1.0 33% Copyright 2008, Leslie Foster 34 35 36 % check that the input is ok 37 if isstruct(matrix) 38 if isfield(matrix,'svals') 39 svals = matrix.svals; 40 if isfield(matrix,'svals_err') 41 svals_err = matrix.svals_err; 42 end 43 else 44 error(['Invalid matrix. The structure ', ... 45 inputname(1),' is missing a svals field.']) 46 end 47 else 48 if isempty(which('SJget')) 49 error(['SJget is not installed. It is available at ',... 50 'http://www.math.sjsu.edu/singular/matrices/SJget.html .']) 51 end 52 matrix = SJget(matrix); 53 svals = matrix.svals; 54 if isfield(matrix,'svals_err') 55 svals_err = matrix.svals_err; 56 end 57 end 58 [m,n] = size(matrix.A); 59 default_tol = 0; 60 if nargin < 2 || isempty(tol) 61 tol = max(m,n)*eps(svals(1)); 62 default_tol = 1; 63 end 64 if nargin >= 2 && tol < 0 65 error('tolerance is negative') 66 end 67 if nargin <= 2 || isempty(option) 68 option = 1; 69 end 70 if option ~= 1 && option ~= 2 71 error('option is not 1 or 2') 72 end 73 74 fullname = matrix.name; 75 if ~isempty(find(svals < 0, 1)) 76 full_svd = 0; 77 else 78 full_svd = 1; 79 end 80 81 %------------------------------------------------------------ 82 % Nikolay and LF: generate svd plots (thumbnail or full size) 83 %------------------------------------------------------------ 84 85 clf; % LF -- to reset graph settings to default 86 if full_svd % full svd plots 87 88 % Create plots for full singular value spectrum 89 % Using: 90 % svals: a (column) vector with partial information about 91 % singular values where any entry with a -1 indicates a 92 % sing. value that has not been calculated 93 % svals_err: a (column) vector with partial information about the 94 % accuracy of the calculated singular values where any 95 % entry with a -1 indicates a sing. value that has not been 96 % calculated. Currently svals_err is not used but code to 97 % use svals_err is included but not accessed 98 % result: plots of singular value spectrum 99 100 warning('off','MATLAB:Axes:NegativeDataInLogAxis') 101 if option == 2 102 % create svd image as thumbnail 103 semilogy(svals,'--rs','LineWidth',10,'MarkerEdgeColor','g', ... 104 'MarkerFaceColor','g','MarkerSize',15); 105 line([0 length(svals)], [1 1]*tol, 'Color', [0 0 .5]); 106 % marker for tolerance 107 set(gca,'XTick',[],'YTick',[]); 108 xlim([0 length(svals)]); 109 else 110 % prepare full-res svd image 111 %semilogy(svals,'--rs','LineWidth',3,'MarkerEdgeColor','g', ... 112 % 'MarkerFaceColor','b','MarkerSize',5); 113 semilogy(1:length(svals),svals,'--rs','LineWidth',3, ... 114 'MarkerEdgeColor','g', ... 115 'MarkerFaceColor','g','MarkerSize',5) 116 xlim([0 length(svals)]); 117 colorm = [0 .5 0]; 118 line(xlim, [1 1]*tol, 'Color', colorm,'Linewidth',1); 119 % marker for rank with MATLAB tolerance 120 121 if ( 1 == 1) 122 % to make the resolution sharper create dots at each svd 123 % location 124 hold on 125 semilogy(1:length(svals),svals,'.','MarkerSize',1); 126 hold off 127 end 128 129 % fix missing marker in Okunbor/aft01 plot 130 axisv1 = axis; 131 axisv1(4) = max( max(svals), axisv1(4) ); 132 axis(axisv1); 133 134 set(gca, 'FontSize', 17); 135 136 % Note in some cases (eg Pajek/SmaGri)the xlabel results 137 % are written to the pdf file or the matlab plot (unless the 138 % plot is full screen). This may occur when the y scale goes 139 % to the underflow limit. Is this a Matlab bug? 140 nsvals_zero = sum( svals == 0 ) ; 141 if ( nsvals_zero == 1 ) 142 xlabel([' singular value number ',... 143 '\newline', ... 144 ' Note: one calculated singular value is 0 and ',... 145 'is not shown'],'FontSize', 14); 146 elseif (nsvals_zero > 1 ) 147 xlabel([' singular value ',... 148 'number \newline ',... 149 ' Note: ', int2str(nsvals_zero), ... 150 ' calculated singular values are 0 ', ... 151 'and are not shown'],'FontSize', 14); 152 else 153 xlabel('singular value number', 'FontSize', 14); 154 end 155 156 % finish full-res svd image 157 grid on 158 line(xlim, [1 1]*tol, 'Color', colorm,'Linewidth',1); 159 % marker for rank with MATLAB tolerance 160 if default_tol 161 text(1, tol, ' max(size(A)) * eps(max(svals))', 'Color',... 162 colorm,'FontSize', 14, 'VerticalAlignment', 'bottom'); 163 else 164 text(1, tol,['tol = ',num2str(tol,'%0.5g')], 'Color',... 165 colorm,'FontSize', 14, 'VerticalAlignment', 'bottom'); 166 end 167 title(['Singular Value Spectrum for ' fullname],... 168 'FontSize',14, 'FontWeight', 'b', 'Interpreter', 'none'); 169 axisv = axis; 170 171 if ( exist('svals_err','var') && 1 == 0) 172 % skip this - the error bounds make the graph very 173 % cluttered near smaller singular values and the errors 174 % are small relative to the default tol 175 % include errbounds in plot 176 isvp=find(svals >= 0 ); 177 isvp = isvp(:)'; 178 bound_lw = 1; % width of line in error bounds 179 for it = isvp 180 if ( svals_err(it) > 0 ) 181 lowbound = max( axisv(3),svals(it)-svals_err(it)); 182 upbound = min(axisv(4),svals(it)+svals_err(it)); 183 line([it,it],[lowbound,upbound],'Color','b',... 184 'Linewidth',bound_lw); 185 boundlength = .15; 186 if ( lowbound ~= axisv(3) ) 187 line([it - boundlength,it + boundlength],... 188 [lowbound,lowbound],... 189 'Color','b','Linewidth',bound_lw); 190 end 191 if (upbound ~= axisv(4) ) 192 line([it - boundlength, it+boundlength],... 193 [upbound, upbound],... 194 'Color','b','Linewidth',bound_lw); 195 end 196 end 197 end 198 end 199 ylabel('calculated singular value', 'FontSize', 14); 200 end 201 warning('on','MATLAB:Axes:NegativeDataInLogAxis') 202 203 else % partial svd plots 204 205 % Create: plot of partial singular value spectrum (lvf + nb) 206 % Using: 207 % svals: a (column) vector with partial information about 208 % singular values where any entry with a -1 indicates 209 % a sing. value that has not been calculated 210 % svals_err: a (column) vector with partial information about the 211 % accuracy of the calculated singular values where any 212 % entry with a -1 indicates a sing. value that has not been 213 % calculated 214 % k_m1: the number of -1's to place in each string of 215 % uncalculated sing. values (set to 3 below) 216 % Produce: 217 % svalp: a vector with each string of -1's in svals replaced by 218 % a string of k_m1 -1's 219 % svalp_err: a vector with each string of -1's in svals_err 220 % replaced by a string of k_m1 -1's 221 % isvalp: a vector whose length is the number of strings of 222 % consecutive calculated singular values in svalp 223 % (and svals). isvalp(i) contains the starting index in 224 % svalp of the ith string of calculated singular values in 225 % svalp 226 % nsvalp: a vector whose length is the number of strings of 227 % consecutive calculated singular values in svalp 228 % (and svals). nsvalp(i) contains the number of calculated 229 % singular values in the ith string of consecutive 230 % calculated singular values 231 % isvals: a vector whose length is the number of strings of 232 % consecutive calculated singular values in svals 233 % (and svalp). isvals(i) contains the starting index in 234 % svals of the ith string of consecutive calculated 235 % singular values. Therefore the calculated singular 236 % values in svals (i.e. those not set equal to -1) have 237 % indices isvals(i): isvals(i) + nsvals(i) - 1 238 % for i = 1: length(isvals) 239 % AND plot of partial singular value spectrum 240 241 % have A, svals and svals_err from earlier code 242 %[m,n]=size(A); 243 isv = find( svals >= 0 ); 244 disv = diff(isv); 245 il1 = find( disv > 1) ; 246 k_m1=3; 247 248 %create svalp and svalp_err 249 clear isvals isvalp nsvalp 250 svalp = []; 251 svalp_err = []; 252 for i = 1:length(il1) 253 if ( i == 1 ) 254 isvalsv = isv(1):isv(il1(i)); 255 isvals(i) = isvalsv(1); 256 if ( isv(1) == 1 ) 257 svalp = svals(isvalsv); 258 svalp_err = svals_err(isvalsv); 259 isvalp(1)=1; 260 nsvalp(1) = isv(il1(i)); 261 else 262 svalp = [ - ones(k_m1,1) ;svals(isvalsv) ]; 263 svalp_err = [ - ones(k_m1,1) ;svals_err(isvalsv) ]; 264 isvalp(1) = k_m1+1; 265 nsvalp(1) = length(isvalsv); 266 end 267 else 268 isvalp(i) = length(svalp)+1; 269 isvalsv = isv(il1(i-1) + 1): isv(il1(i)); 270 isvals(i) = isvalsv(1); 271 nsvalp(i) = length( svals( isv(il1(i-1) + 1): isv(il1(i))) ); 272 svalp = [ svalp; svals( isv(il1(i-1) + 1): isv(il1(i)))]; 273 svalp_err = [ svalp_err; svals_err( isv(il1(i-1) + 1): ... 274 isv(il1(i)))]; 275 end 276 svalp = [svalp; - ones(k_m1,1) ]; 277 svalp_err = [svalp_err; - ones(k_m1,1) ]; 278 end 279 % create isvalp, nsvalp, isvals and finish svalp 280 if ( length( il1 ) == 0 ) 281 ib =1; 282 else 283 ib = isv( il1( length( il1 )) + 1 ); 284 end 285 isvalp(length(il1)+1) = length(svalp) + 1; 286 nsvalp(length(il1)+1) = length(ib:isv(end)); 287 isvalsv = ib :isv(end) ; 288 isvals(length(il1)+1) = isvalsv(1); 289 svalp = [ svalp ; svals( ib :isv(end)) ]; 290 svalp_err = [ svalp_err ; svals_err( ib :isv(end)) ]; 291 if(svals(end) < 0 ) 292 svalp = [svalp; - ones(k_m1,1)]; 293 svalp_err = [svalp_err; - ones(k_m1,1)]; 294 end 295 296 if ( 1 == 0 ) 297 % display results (for debugging) 298 disp('svals'''), disp(svals'), disp('isv''') 299 disp(isv'), disp('il1'''), disp(il1') 300 disp('svalp'''), disp(svalp'), disp('svalp_err''') 301 disp(svalp_err'), disp('isvalp'), disp(isvalp) 302 disp('nsvalp'), disp(nsvalp), disp('isvals') 303 disp(isvals) 304 end 305 306 %draw plots 307 %warning('off','MATLAB:Axes:NegativeDataInLogAxis') 308 if ( option == 2) 309 % create svd image as thumbnail 310 warning('off','MATLAB:Axes:NegativeDataInLogAxis') 311 312 semilogy(svalp,'--rs','LineWidth',10,'MarkerEdgeColor','g', ... 313 'MarkerFaceColor','g','MarkerSize',20); 314 xlim([0 length(svalp)]); 315 line([0 length(svalp)], [1 1]*tol, 'Color', [0 0 .5]); 316 % marker for tolerance 317 set(gca,'XTick',[],'YTick',[]); 318 319 axisv1 = axis; % in Matlab 7.6 these 2 commands suppress 320 axis(axisv1); % the printing of an (undesired) warning 321 % reason is unclear -- hmmm 322 323 warning('on','MATLAB:Axes:NegativeDataInLogAxis'); 324 325 else 326 % prepare full-res svd image 327 warning('off','MATLAB:Axes:NegativeDataInLogAxis') 328 329 semilogy(svalp,'--rs','LineWidth',5,'MarkerEdgeColor','g', ... 330 'MarkerFaceColor','g','MarkerSize',8); 331 332 333 xlim([0 length(svalp)]); 334 % finish full-res svd image 335 colorm = [0 .5 0]; 336 line(xlim, [1 1]*tol, 'Color', colorm,'Linewidth',1); 337 % marker for rank with MATLAB tolerance 338 339 % fixes missing marker in Okunbor/aft01 plot 340 axisv1 = axis; 341 axisv1(4) = max( max(svalp), axisv1(4) ); 342 axis(axisv1); 343 344 set(gca, 'FontSize', 17); 345 %define locations to place tick marks 346 xtick_v=[]; 347 for ixl = 1: length(isvalp) 348 xtick_v = [ xtick_v, isvalp(ixl): ... 349 (isvalp(ixl)+nsvalp(ixl)-1) ]; 350 end 351 if ( xtick_v(end) ~= length(svalp) ) 352 xtick_v =[ xtick_v length(svalp)]; 353 end 354 set(gca, 'XTick', xtick_v); 355 356 %create labels for x axis 357 xticklabelv = []; 358 ixt = 0; 359 for ixl = 1:length(isvals) 360 for ixl2 = 1:nsvalp(ixl) 361 ixt = ixt+1; 362 if ( (ixl == 1 && nsvalp(ixl) < 10) | ixl2 == 1 | ... 363 ixl2 == nsvalp(ixl) ) 364 xticklabelv{ixt}=int2str(isvals(ixl)+ixl2-1); 365 else 366 xticklabelv{ixt} =''; 367 end 368 end 369 end 370 ixt = ixt+1; 371 xticklabelv{ixt} = int2str(min(m,n)); 372 set(gca, 'XTickLabel', xticklabelv); 373 374 % finish full-res svd image 375 grid on 376 if default_tol 377 text(1, tol, ' max(size(A)) * eps(max(svals))', 'Color',... 378 colorm,'FontSize', 14, 'VerticalAlignment', 'bottom'); 379 else 380 text(1, tol,['tol = ',num2str(tol,'%0.5g')], 'Color',... 381 colorm,'FontSize', 14, 'VerticalAlignment', 'bottom'); 382 end 383 title(['Partial Singular Value Spectrum for ' fullname],... 384 'FontSize',14, 'FontWeight', 'b', 'Interpreter', 'none') ; 385 axisv = axis; 386 387 % include errbounds in plot 388 isvp=find(svalp >= 0 ); 389 isvp = isvp(:)'; 390 bound_lw = 1; % width of line in error bounds 391 for it = isvp 392 if ( svalp_err(it) > 0 ) 393 lowbound = max( axisv(3),svalp(it)-svalp_err(it)); 394 upbound = min(axisv(4),svalp(it)+svalp_err(it)); 395 line([it,it],[lowbound,upbound],'Color','b',... 396 'Linewidth',bound_lw); 397 boundlength = .15; 398 if ( lowbound ~= axisv(3) ) 399 line([it - boundlength,it + boundlength],... 400 [lowbound,lowbound],... 401 'Color','b','Linewidth',bound_lw); 402 end 403 if (upbound ~= axisv(4) ) 404 line([it - boundlength, it+boundlength],... 405 [upbound, upbound],... 406 'Color','b','Linewidth',bound_lw); 407 end 408 end 409 end 410 411 nsvals_zero = sum( svals == 0 ) ; 412 if ( nsvals_zero == 1 ) 413 xlabel([' singular value ',... 414 'number \newline', ... 415 ' Note: one calculated singular value is 0 ',... 416 'and is not shown'],'FontSize', 14); 417 elseif (nsvals_zero > 1 ) 418 xlabel([' singular value ',... 419 'number \newline ',... 420 ' Note: ', int2str(nsvals_zero), ... 421 ' calculated singular values are 0 ', ... 422 'and are not shown'],'FontSize', 14); 423 else 424 xlabel('singular value number', 'FontSize', 14); 425 end 426 ylabel('calculated singular value', 'FontSize', 14); 427 warning('on','MATLAB:Axes:NegativeDataInLogAxis') 428 end 429 end 430end 431