Mercurial > hg > octave-lyh

comparison doc/interpreter/sparseimages.m @ 16994:333243133364

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use matrix concatenation for strings, rather than cstrcat(), for clarity and performance. matrix concatenation is ~80% faster than cstrcat(). * doc/interpreter/strings.txi: Document preference for matrix concatenation as opposed to alternatives. * doc/interpreter/geometryimages.m, doc/interpreter/interpimages.m, doc/interpreter/plotimages.m, doc/interpreter/sparseimages.m, doc/interpreter/splineimages.m, scripts/general/genvarname.m, scripts/general/int2str.m, scripts/general/num2str.m, scripts/help/__makeinfo__.m, scripts/help/help.m, scripts/miscellaneous/copyfile.m, scripts/miscellaneous/dir.m, scripts/miscellaneous/edit.m, scripts/miscellaneous/fact.m, scripts/miscellaneous/fullfile.m, scripts/miscellaneous/mkoctfile.m, scripts/miscellaneous/movefile.m, scripts/miscellaneous/perl.m, scripts/miscellaneous/python.m, scripts/miscellaneous/run.m, scripts/miscellaneous/tempdir.m, scripts/miscellaneous/unpack.m, scripts/pkg/private/configure_make.m, scripts/pkg/private/create_pkgadddel.m, scripts/pkg/private/extract_pkg.m, scripts/pkg/private/get_description.m, scripts/pkg/private/get_forge_pkg.m, scripts/pkg/private/getarch.m, scripts/pkg/private/getarchprefix.m, scripts/pkg/private/install.m, scripts/pkg/private/installed_packages.m, scripts/pkg/private/load_packages_and_dependencies.m, scripts/pkg/private/rebuild.m, scripts/pkg/private/repackage.m, scripts/pkg/private/shell.m, scripts/pkg/private/uninstall.m, scripts/plot/private/__go_draw_axes__.m, scripts/signal/spectral_adf.m, scripts/signal/spectral_xdf.m, scripts/statistics/tests/z_test.m, scripts/statistics/tests/z_test_2.m, scripts/strings/mat2str.m, scripts/strings/strtok.m, scripts/testfun/__run_test_suite__.m, scripts/testfun/assert.m, scripts/testfun/demo.m, scripts/testfun/speed.m, scripts/testfun/test.m, test/eval-catch.tst, test/io.tst, test/try.tst: Replace cstrcat() with matrix concatenation where possible.
author Rik <rik@octave.org>
date Wed, 17 Jul 2013 14:02:32 -0700
parents 2a4f83826024
children
comparison
equal deleted inserted replaced
16993:78f57b14535c 16994:333243133364
24 endif 24 endif
25 25
26 if (__have_feature__ ("COLAMD") 26 if (__have_feature__ ("COLAMD")
27 && __have_feature__ ("CHOLMOD") 27 && __have_feature__ ("CHOLMOD")
28 && __have_feature__ ("UMFPACK")) 28 && __have_feature__ ("UMFPACK"))
29 if (strcmp(typ,"txt")) 29 if (strcmp (typ,"txt"))
30 txtimages (nm, 15, typ); 30 txtimages (nm, 15, typ);
31 else 31 else
32 if (strcmp (nm, "gplot")) 32 if (strcmp (nm, "gplot"))
33 gplotimages ("gplot", typ); 33 gplotimages ("gplot", typ);
34 elseif (strcmp (nm, "grid")) 34 elseif (strcmp (nm, "grid"))
60 function gplotimages (nm, typ) 60 function gplotimages (nm, typ)
61 hide_output (); 61 hide_output ();
62 if (strcmp (typ, "eps")) 62 if (strcmp (typ, "eps"))
63 d_typ = "-depsc2"; 63 d_typ = "-depsc2";
64 else 64 else
65 d_typ = cstrcat ("-d", typ); 65 d_typ = ["-d" typ];
66 endif 66 endif
67 67
68 A = sparse ([2,6,1,3,2,4,3,5,4,6,1,5], 68 A = sparse ([2,6,1,3,2,4,3,5,4,6,1,5],
69 [1,1,2,2,3,3,4,4,5,5,6,6], 1, 6, 6); 69 [1,1,2,2,3,3,4,4,5,5,6,6], 1, 6, 6);
70 xy = [0,4,8,6,4,2;5,0,5,7,5,7]'; 70 xy = [0,4,8,6,4,2;5,0,5,7,5,7]';
71 gplot (A, xy) 71 gplot (A, xy);
72 print (cstrcat (nm, ".", typ), d_typ) 72 print ([nm "." typ], d_typ);
73 hide_output (); 73 hide_output ();
74 endfunction 74 endfunction
75 75
76 function txtimages(nm, n, typ) 76 function txtimages(nm, n, typ)
77 a = 10*speye(n) + sparse(1:n,ceil([1:n]/2),1,n,n) + ... 77 a = 10*speye(n) + sparse(1:n,ceil([1:n]/2),1,n,n) + ...
82 fputs (fid, "+---------------------------------+\n"); 82 fputs (fid, "+---------------------------------+\n");
83 fputs (fid, "| Image unavailable in text mode. |\n"); 83 fputs (fid, "| Image unavailable in text mode. |\n");
84 fputs (fid, "+---------------------------------+\n"); 84 fputs (fid, "+---------------------------------+\n");
85 fclose (fid); 85 fclose (fid);
86 elseif (strcmp (nm, "spmatrix")) 86 elseif (strcmp (nm, "spmatrix"))
87 printsparse(a,cstrcat("spmatrix.",typ)); 87 printsparse (a, ["spmatrix." typ]);
88 else 88 else
89 if (__have_feature__ ("COLAMD") 89 if (__have_feature__ ("COLAMD")
90 && __have_feature__ ("CHOLMOD")) 90 && __have_feature__ ("CHOLMOD"))
91 if (strcmp (nm, "spchol")) 91 if (strcmp (nm, "spchol"))
92 r1 = chol(a); 92 r1 = chol (a);
93 printsparse(r1,cstrcat("spchol.",typ)); 93 printsparse (r1, ["spchol." typ]);
94 elseif (strcmp (nm, "spcholperm")) 94 elseif (strcmp (nm, "spcholperm"))
95 [r2,p2,q2]=chol(a); 95 [r2,p2,q2] = chol (a);
96 printsparse(r2,cstrcat("spcholperm.",typ)); 96 printsparse(r2, ["spcholperm." typ]);
97 endif 97 endif
98 ## printf("Text NNZ: Matrix %d, Chol %d, PermChol %d\n",nnz(a),nnz(r1),nnz(r2)); 98 ## printf("Text NNZ: Matrix %d, Chol %d, PermChol %d\n",nnz(a),nnz(r1),nnz(r2));
99 endif 99 endif
100 endif 100 endif
101 endfunction 101 endfunction
103 function otherimages(nm, n, typ) 103 function otherimages(nm, n, typ)
104 hide_output (); 104 hide_output ();
105 if (strcmp (typ, "eps")) 105 if (strcmp (typ, "eps"))
106 d_typ = "-depsc2"; 106 d_typ = "-depsc2";
107 else 107 else
108 d_typ = cstrcat ("-d", typ); 108 d_typ = ["-d" typ];
109 endif 109 endif
110 110
111 a = 10*speye(n) + sparse(1:n,ceil([1:n]/2),1,n,n) + ... 111 a = 10*speye(n) + sparse(1:n,ceil([1:n]/2),1,n,n) + ...
112 sparse(ceil([1:n]/2),1:n,1,n,n); 112 sparse(ceil([1:n]/2),1:n,1,n,n);
113 if (strcmp (nm, "spmatrix")) 113 if (strcmp (nm, "spmatrix"))
114 spy(a); 114 spy (a);
115 axis("ij") 115 axis ("ij");
116 print(cstrcat("spmatrix.",typ), d_typ) 116 print (["spmatrix." typ], d_typ);
117 hide_output (); 117 hide_output ();
118 else 118 else
119 if (__have_feature__ ("COLAMD") 119 if (__have_feature__ ("COLAMD")
120 && __have_feature__ ("CHOLMOD")) 120 && __have_feature__ ("CHOLMOD"))
121 if (strcmp (nm, "spchol")) 121 if (strcmp (nm, "spchol"))
122 r1 = chol(a); 122 r1 = chol (a);
123 spy(r1); 123 spy (r1);
124 axis("ij") 124 axis ("ij");
125 print(cstrcat("spchol.",typ), d_typ) 125 print (["spchol." typ], d_typ);
126 hide_output (); 126 hide_output ();
127 elseif (strcmp (nm, "spcholperm")) 127 elseif (strcmp (nm, "spcholperm"))
128 [r2,p2,q2]=chol(a); 128 [r2,p2,q2] = chol (a);
129 spy(r2); 129 spy (r2);
130 axis("ij") 130 axis ("ij");
131 print(cstrcat("spcholperm.",typ), d_typ) 131 print (["spcholperm." typ], d_typ);
132 hide_output (); 132 hide_output ();
133 endif 133 endif
134 ## printf("Image NNZ: Matrix %d, Chol %d, PermChol %d\n",nnz(a),nnz(r1),nnz(r2)); 134 ## printf("Image NNZ: Matrix %d, Chol %d, PermChol %d\n",nnz(a),nnz(r1),nnz(r2));
135 endif 135 endif
136 endif 136 endif
137 endfunction 137 endfunction
138 138
139 function printsparse(a, nm) 139 function printsparse (a, nm)
140 fid = fopen (nm,"wt"); 140 fid = fopen (nm,"wt");
141 fputs (fid, "\n"); 141 fputs (fid, "\n");
142 for i = 1:size(a,1) 142 for i = 1:rows (a)
143 if (rem(i,5) == 0) 143 if (rem (i,5) == 0)
144 fprintf (fid," %2d - ", i); 144 fprintf (fid," %2d - ", i);
145 else 145 else
146 fprintf (fid," | "); 146 fprintf (fid," | ");
147 endif 147 endif
148 for j = 1:size(a,2) 148 for j = 1:columns (a)
149 if (a(i,j) == 0) 149 if (a(i,j) == 0)
150 fprintf(fid," ") 150 fprintf (fid," ");
151 else 151 else
152 fprintf(fid," *") 152 fprintf (fid," *");
153 endif 153 endif
154 endfor 154 endfor
155 fprintf(fid,"\n") 155 fprintf (fid,"\n");
156 endfor 156 endfor
157 fprintf(fid," |-"); 157 fprintf (fid," |-");
158 for j=1:size(a,2) 158 for j=1:columns (a)
159 if (rem(j,5)==0) 159 if (rem (j,5) == 0)
160 fprintf(fid,"-|"); 160 fprintf (fid,"-|");
161 else 161 else
162 fprintf(fid,"--"); 162 fprintf (fid,"--");
163 endif 163 endif
164 endfor 164 endfor
165 fprintf(fid,"\n") 165 fprintf (fid,"\n");
166 fprintf(fid," "); 166 fprintf (fid," ");
167 for j=1:size(a,2) 167 for j=1:columns (a)
168 if (rem(j,5)==0) 168 if (rem (j,5) == 0)
169 fprintf(fid,"%2d",j); 169 fprintf (fid,"%2d",j);
170 else 170 else
171 fprintf(fid," "); 171 fprintf (fid," ");
172 endif 172 endif
173 endfor 173 endfor
174 fclose(fid); 174 fclose (fid);
175 endfunction 175 endfunction
176 176
177 function femimages (nm, typ) 177 function femimages (nm, typ)
178 hide_output (); 178 hide_output ();
179 if (strcmp (typ, "eps")) 179 if (strcmp (typ, "eps"))
180 d_typ = "-depsc2"; 180 d_typ = "-depsc2";
181 else 181 else
182 d_typ = cstrcat ("-d", typ); 182 d_typ = ["-d" typ];
183 endif 183 endif
184 184
185 if (__have_feature__ ("COLAMD") 185 if (__have_feature__ ("COLAMD")
186 && __have_feature__ ("CHOLMOD") 186 && __have_feature__ ("CHOLMOD")
187 && __have_feature__ ("UMFPACK")) 187 && __have_feature__ ("UMFPACK"))
188 ## build a rectangle 188 ## build a rectangle
189 node_y = [1;1.2;1.5;1.8;2]*ones(1,11); 189 node_y = [1;1.2;1.5;1.8;2]*ones(1,11);
190 node_x = ones(5,1)*[1,1.05,1.1,1.2,1.3,1.5,1.7,1.8,1.9,1.95,2]; 190 node_x = ones(5,1)*[1,1.05,1.1,1.2,1.3,1.5,1.7,1.8,1.9,1.95,2];
191 nodes = [node_x(:), node_y(:)]; 191 nodes = [node_x(:), node_y(:)];
192 192
193 [h,w] = size(node_x); 193 [h,w] = size (node_x);
194 elems = []; 194 elems = [];
195 for idx = 1:w-1 195 for idx = 1:w-1
196 widx = (idx-1)*h; 196 widx = (idx-1)*h;
197 elems = [elems; widx+[(1:h-1);(2:h);h+(1:h-1)]']; 197 elems = [elems; widx+[(1:h-1);(2:h);h+(1:h-1)]'];
198 elems = [elems; widx+[(2:h);h+(2:h);h+(1:h-1)]']; 198 elems = [elems; widx+[(2:h);h+(2:h);h+(1:h-1)]'];
199 endfor 199 endfor
200 200
201 E = size(elems,1); #No. of elements 201 E = size (elems,1); #No. of elements
202 N = size(nodes,1); #No. of elements 202 N = size (nodes,1); #No. of elements
203 D = size(elems,2); #dimensions+1 203 D = size (elems,2); #dimensions+1
204 204
205 ## Plot FEM Geometry 205 ## Plot FEM Geometry
206 elemx = elems(:,[1,2,3,1])'; 206 elemx = elems(:,[1,2,3,1])';
207 xelems = reshape( nodes(elemx, 1), 4, E); 207 xelems = reshape (nodes(elemx, 1), 4, E);
208 yelems = reshape( nodes(elemx, 2), 4, E); 208 yelems = reshape (nodes(elemx, 2), 4, E);
209 209
210 ## Set element conductivity 210 ## Set element conductivity
211 conductivity = [1*ones(1,16),2*ones(1,48),1*ones(1,16)]; 211 conductivity = [1*ones(1,16),2*ones(1,48),1*ones(1,16)];
212 212
213 ## Dirichlet boundary conditions 213 ## Dirichlet boundary conditions
219 ## length and element conductivity 219 ## length and element conductivity
220 N_nodes = []; 220 N_nodes = [];
221 N_value = []; 221 N_value = [];
222 222
223 ## Calculate connectivity matrix 223 ## Calculate connectivity matrix
224 C = sparse((1:D*E), reshape(elems',D*E,1),1, D*E, N); 224 C = sparse ((1:D*E), reshape (elems',D*E,1),1, D*E, N);
225 225
226 ## Calculate stiffness matrix 226 ## Calculate stiffness matrix
227 Siidx = floor([0:D*E-1]'/D)*D*ones(1,D) + ones(D*E,1)*(1:D) ; 227 Siidx = floor ([0:D*E-1]'/D)*D*ones(1,D) + ones(D*E,1)*(1:D) ;
228 Sjidx = [1:D*E]'*ones(1,D); 228 Sjidx = [1:D*E]'*ones (1,D);
229 Sdata = zeros(D*E,D); 229 Sdata = zeros (D*E,D);
230 dfact = prod(2:(D-1)); 230 dfact = prod (2:(D-1));
231 for j = 1:E 231 for j = 1:E
232 a = inv([ ones(D,1), nodes( elems(j,:), : ) ]); 232 a = inv ([ ones(D,1), nodes( elems(j,:), : ) ]);
233 const = conductivity(j)*2/dfact/abs(det(a)); 233 const = conductivity(j)*2/dfact/abs(det(a));
234 Sdata(D*(j-1)+(1:D),:)= const * a(2:D,:)'*a(2:D,:); 234 Sdata(D*(j-1)+(1:D),:)= const * a(2:D,:)'*a(2:D,:);
235 endfor 235 endfor
236 236
237 ## Element-wise system matrix 237 ## Element-wise system matrix
238 SE = sparse(Siidx,Sjidx,Sdata); 238 SE = sparse (Siidx,Sjidx,Sdata);
239 ## Global system matrix 239 ## Global system matrix
240 S = C'* SE *C; 240 S = C'* SE *C;
241 241
242 ## Set Dirichlet boundary 242 ## Set Dirichlet boundary
243 V = zeros(N,1); 243 V = zeros (N,1);
244 V(D_nodes) = D_value; 244 V(D_nodes) = D_value;
245 idx = 1:N; 245 idx = 1:N;
246 idx(D_nodes) = []; 246 idx(D_nodes) = [];
247 247
248 ## Set Neumann boundary 248 ## Set Neumann boundary
249 Q = zeros(N,1); 249 Q = zeros (N,1);
250 Q(N_nodes) = N_value; # FIXME 250 Q(N_nodes) = N_value; # FIXME
251 251
252 V(idx) = S(idx,idx)\( Q(idx) - S(idx,D_nodes)*V(D_nodes) ); 252 V(idx) = S(idx,idx)\( Q(idx) - S(idx,D_nodes)*V(D_nodes) );
253 253
254 velems = reshape( V(elemx), 4, E); 254 velems = reshape ( V(elemx), 4, E);
255 255
256 sz = size(xelems,2); 256 sz = size (xelems,2);
257 257
258 plot3 (xelems, yelems, velems); 258 plot3 (xelems, yelems, velems);
259 view (10, 10); 259 view (10, 10);
260 print(cstrcat(nm,".",typ), d_typ) 260 print ([nm "." typ], d_typ);
261 hide_output (); 261 hide_output ();
262 endif 262 endif
263 endfunction 263 endfunction
264 264
265 ## There is no sparse matrix implementation available because of missing 265 ## There is no sparse matrix implementation available because of missing
279 279
280 hide_output (); 280 hide_output ();
281 if (strcmp (typ, "eps")) 281 if (strcmp (typ, "eps"))
282 d_typ = "-depsc2"; 282 d_typ = "-depsc2";
283 else 283 else
284 d_typ = cstrcat ("-d", typ); 284 d_typ = ["-d" typ];
285 endif 285 endif
286 286
287 x = y = linspace (-8, 8, 41)'; 287 x = y = linspace (-8, 8, 41)';
288 [xx, yy] = meshgrid (x, y); 288 [xx, yy] = meshgrid (x, y);
289 r = sqrt (xx .^ 2 + yy .^ 2) + eps; 289 r = sqrt (xx .^ 2 + yy .^ 2) + eps;
290 z = sin (r) ./ r; 290 z = sin (r) ./ r;
291 unwind_protect 291 unwind_protect
292 mesh (x, y, z); 292 mesh (x, y, z);
293 title ("Sorry, graphics are unavailable because Octave was\ncompiled without a sparse matrix implementation."); 293 title ("Sorry, graphics are unavailable because Octave was\ncompiled without a sparse matrix implementation.");
294 unwind_protect_cleanup 294 unwind_protect_cleanup
295 print (cstrcat (nm, ".", typ), d_typ); 295 print ([nm "." typ], d_typ);
296 hide_output (); 296 hide_output ();
297 end_unwind_protect 297 end_unwind_protect
298 endif 298 endif
299 endfunction 299 endfunction
300 300