Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
csr.cc
Go to the documentation of this file.
1
21// the csr unix command
22// author: Pierre.Saramito@imag.fr
23// date: 12 may 1997, updated 9 oct 2012
24
25namespace rheolef {
84} // namespace rheolef
85
86#include "rheolef/linalg.h"
87using namespace rheolef;
88using namespace std;
89
90// ------------------------------------------------------------------------
91// print function
92// ------------------------------------------------------------------------
93template <class T>
94struct log10_abs_or_zero {
95 T operator()(T x) { return ((x == T(0.)) ? T(0.) : log10(xabs(x))); }
96};
97template <class T>
98void
100 std::ostream& s,
101 const csr<T,sequential>& a)
102{
104
105 // scales
106 bool logscale = iorheo::getlogscale(s);
107 bool color = iorheo::getcolor(s);
108
109 // title of the drawing
110 // const char title [256] = "Sparse Matrix";
111 const char title [256] = "";
112
113 // ptitle = position of title; 0 under the drawing, else above
114 const int ptitle = 0;
115
116 // size = size of the drawing (unit =cm)
117 const T size = 15;
118
119 // units cm do dot conversion factor and a4 paper size
120 const T paperx = 21.0;
121 const T conv = 2.54; // cm-per-inch
122 T u2dot = 72.0/conv;
123
124 int maxdim = max(a.nrow(), a.ncol());
125 int m = 1 + maxdim;
126
127 // left and right margins (drawing is centered)
128 T lrmrgn = (paperx-size)/2.0;
129
130 // bottom margin : 2 cm
131 T botmrgn = 2;
132
133 // scaling factor
134 T scfct = size*u2dot/m;
135
136 // matrix frame line witdh
137 const T frlw = 0.25;
138
139 // font size for title (cm)
140 const T fnstit = 0.5;
141 int ltit = strlen(title);
142
143 // position of title : centered horizontally
144 // at 1.0 cm vertically over the drawing
145 const T ytitof = 1.0;
146 T xtit = paperx/2.0;
147
148 T ytit = botmrgn + size*int(a.nrow())/T(m) + T(ytitof);
149
150 // almost exact bounding box
151 T xl = lrmrgn*u2dot - scfct*frlw/2;
152 T xr = (lrmrgn+size)*u2dot + scfct*frlw/2;
153 T yb = botmrgn*u2dot - scfct*frlw/2;
154
155 T yt = (botmrgn+size*int(a.nrow())/T(m))*u2dot + scfct*frlw/2.;
156
157 if (ltit == 0)
158 yt = yt + (ytitof+fnstit*0.70)*u2dot;
159
160 // add some room to bounding box
161 const T delt = 10.0;
162 xl -= delt;
163 xr += delt;
164 yb -= delt;
165 yt += delt;
166
167 // correction for title under the drawing
168 if (ptitle == 0 && ltit > 0) {
169 ytit = botmrgn + fnstit*0.3;
170 botmrgn = botmrgn + ytitof + fnstit*0.7;
171 }
172
173 // print header
174
175 s << "%!PS-Adobe-2.0\n";
176 s << "%%Creator: sparskit++ 1.0, 1997, Computer Modern Artists\n";
177 s << "%%Title: sparskit++ CSR matrix\n";
178 s << "%%BoundingBox: " << xl << " " << yb << " " << xr << " " << yt << std::endl;
179 s << "%%EndComments\n";
180 s << "/cm {72 mul 2.54 div} def\n";
181 s << "/mc {72 div 2.54 mul} def\n";
182 s << "/pnum { 72 div 2.54 mul 20 string\n";
183 s << "cvs print ( ) print} def\n";
184 s << "/Cshow {dup stringwidth pop -2 div 0 rmoveto show} def\n";
185
186 // we leave margins etc. in cm so it is easy to modify them if
187 // needed by editing the output file
188 s << "gsave\n";
189 if (ltit > 0)
190 s << "/Helvetica findfont " << fnstit << " cm scalefont setfont \n";
191 s << xtit << " cm " << ytit << " cm moveto\n";
192 s << "(" << title << ") Cshow\n";
193
194 s << lrmrgn << " cm " << botmrgn << " cm translate\n";
195 s << size << " cm " << m << " div dup scale\n";
196
197 // draw a frame around the matrix
198 s << frlw << " setlinewidth\n";
199 s << "newpath\n";
200 s << "0 0 moveto\n";
201 s << a.ncol()+1 << " " << 0 << " lineto\n";
202 s << a.ncol()+1 << " " << a.nrow()+1 << " lineto\n";
203 s << 0 << " " << a.nrow()+1 << " lineto\n";
204 s << "closepath stroke\n";
205
206 s << "1 1 translate\n";
207 s << "0.8 setlinewidth\n";
208 s << "/p {moveto 0 -.40 rmoveto \n";
209 s << " 0 .80 rlineto stroke} def\n";
210
211#ifdef TODO
212 if (color) {
213
214 // the number of used colors
215 const int ncolor = 100;
216
217 T max_;
218 T min_;
219 T inf = std::numeric_limits<T>::max();
220 log10_abs_or_zero<T> filter;
221
222 if (logscale) {
223 max_ = select_op_value (a, a+a.nnz(), -inf, std::greater<T>(), filter);
224 min_ = select_op_value (a, a+a.nnz(), inf, std::less<T>(), filter);
225 } else {
226 max_ = select_max(a, a+a.nnz(), -inf);
227 min_ = select_min(a, a+a.nnz(), inf);
228 }
229 palette tab = palette(min_, max_, ncolor);
230 tab.postscript_print_color_def (s);
231
232 // color: plotting loop
233 for (Index i = 0; i < a.nrow(); i++)
234 for (Index k = ia [i]; k < ia [i+1]; k++) {
235
236 T value = a[k];
237 if (logscale) value = filter(value);
238 int ic = tab.color_index(value);
239 s << "c" << ic << " "
240 << ja [k] << " "
241 << a.nrow()-i-1 << " p\n";
242 }
243 } else {
244#endif // TODO
245 // black & white: plotting loop
246 typename csr<T,sequential>::const_iterator dia_ia = a.begin();
247 for (size_type i = 0, n = a.nrow(); i < n; i++) {
248 for (typename csr<T,sequential>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
249 size_type j = (*p).first;
250 s << j << " " << a.nrow()-i-1 << " p\n";
251 }
252 }
253#ifdef TODO
254 }
255#endif // TODO
256 s << "showpage\n";
257}
258// ------------------------------------------------------------------------
259// unix command
260// ------------------------------------------------------------------------
261void usage()
262{
263 derr << "csr: usage: csr "
264 << "[-ps|-matlab] "
265 << "[-name string] "
266 << "file[.mtx]"
267 << endl;
268 exit (1);
269}
270int main(int argc, char **argv)
271{
272 environment rheolef (argc, argv);
273 check_macro (communicator().size() == 1, "csr: command may be used as mono-process only");
274 if (argc <= 1) usage();
275 bool have_file = false;
276 typedef enum { mtx, ml, ps} out_format_type;
277 out_format_type out_format = ps;
278 string filename;
279 string name = "a";
280 for (int i = 1; i < argc; i++) {
281 if (strcmp (argv[i], "-ps") == 0) out_format = ps;
282 else if (strcmp (argv[i], "-matlab") == 0) out_format = ml;
283 else if (strcmp (argv[i], "-mtx") == 0) out_format = mtx;
284 else if (strcmp (argv[i], "-name") == 0) {
285 if (i == argc-1) { std::cerr << "csr -name: option argument missing" << std::endl; usage(); }
286 name = argv[++i];
287 }
288 else if (strcmp (argv[i], "-") == 0) filename = "-";
289 else if (argv[i][0] == '-') usage();
290 else {
291 filename = argv[i];
292 }
293 }
295 if (filename == "" || filename == "-") {
296 din >> a;
297 } else {
298 idiststream a_in (filename);
299 a_in >> a;
300 }
301 if (out_format == ml) {
302 dout << matlab << setbasename(name) << a;
303 } else if (out_format == ps) {
304 print_postscript (cout, a);
305 } else if (out_format == mtx) {
306 dout << a;
307 }
308#ifdef TODO
309 //
310 // set default options
311 //
312 cout << hb << setbasename("a");
313
314 bool get_rhs = false;
315
316 //
317 // get input options
318 //
319 for (int i = 1; i < argc; i++) {
320 if (strcmp (argv[i], "-transpose-in") == 0) cin >> transpose;
321 else if (strcmp (argv[i], "-notranspose-in") == 0) cin >> notranspose;
322 else if (strcmp (argv[i], "-input-hb") == 0) cin >> hb;
323 else if (strcmp (argv[i], "-input-mm") == 0) cin >> matrix_market;
324 }
325 //
326 // input matrix
327 //
328 csr<Float> a;
329 cin >> a;
330 //
331 // get options
332 //
333 for (int i = 1; i < argc; i++) {
334
335 if (strcmp (argv[i], "-hb") == 0) cout << hb;
336 else if (strcmp (argv[i], "-mm") == 0) cout << matrix_market;
337 else if (strcmp (argv[i], "-ml") == 0) cout << ml;
338 else if (strcmp (argv[i], "-matlab") == 0) cout << matlab;
339 else if (strcmp (argv[i], "-sparse-matlab") == 0) cout << sparse_matlab;
340 else if (strcmp (argv[i], "-ps") == 0) cout << ps;
341 else if (strcmp (argv[i], "-dump") == 0) cout << dump;
342
343 else if (strcmp (argv[i], "-tril") == 0 ||
344 strcmp (argv[i], "-triu") == 0) {
345 int k = 0;
346 int io = i;
347 if (i+1 < argc) {
348 if (isdigit(argv[i+1][0]))
349 k = atoi(argv[++i]);
350 else if (argv[i+1][0] == '-' &&
351 isdigit(argv[i+1][1]))
352 k = - atoi(argv[++i]+1);
353 }
354 if (strcmp (argv[io], "-tril") == 0) a = tril(a,k);
355 else a = triu(a,k);
356 }
357 else if (strcmp (argv[i], "-gibbs") == 0) {
358 permutation p = gibbs(a);
359 a = perm(a, p, p);
360 }
361 else if (strcmp (argv[i], "-rhs") == 0) {
362 cout << hb;
363 get_rhs = true;
364 }
365 else if (strcmp (argv[i], "-logscale") == 0) cout << logscale;
366 else if (strcmp (argv[i], "-nologscale") == 0) cout << nologscale;
367 else if (strcmp (argv[i], "-color") == 0) cout << color;
368 else if (strcmp (argv[i], "-black-and-white") == 0) cout << black_and_white;
369
370 else if (strcmp (argv[i], "-transpose-out") == 0) cout << transpose;
371 else if (strcmp (argv[i], "-notranspose-out") == 0) cout << notranspose;
372 else if (strcmp (argv[i], "-name") == 0) {
373 if (++i >= argc) {
374 cerr << "-name: missing argument" << endl;
375 usage();
376 }
377 cout << setbasename (argv[i]);
378 }
379 // skit input options
380 else if (strcmp (argv[i], "-transpose-in") == 0) ;
381 else if (strcmp (argv[i], "-notranspose-in") == 0) ;
382 else if (strcmp (argv[i], "-input-hb") == 0) ;
383 else if (strcmp (argv[i], "-input-mm") == 0) ;
384
385 else {
386 cerr << "unexpected command-line argument `" << argv[i] << "'" << endl;
387 usage();
388 }
389 }
390 //
391 // output matrix
392 //
393 if (!get_rhs) {
394
395 // simple output
396 cout << a;
397
398 } else {
399
400 unsigned int n = iorheo::getnrhs(cin);
401
402 // Here, we get the number of right-hand-sides.
403 // Next, get the type: has guess and/or exact solution:
404
405 string rhs_type = iorheo::getrhstype(cin);
406
407 // and start to write.
408 // We specify first in header the number of right-hand-side,
409 // and then output the matrix.
410
411 cout << hb << setnrhs(n) << setrhstype(rhs_type) << notranspose << a ;
412
413 // A first loop for optional rhs:
414
415 for (unsigned int i = 0; i < n; i++) {
416 vec<Float> b;
417 cin >> b;
418 cout << b;
419 }
420 //
421 // Test if guess vectors is supplied:
422 //
423 if (rhs_type[1] == 'G') {
424 vec<Float> x0;
425 for (unsigned int i = 0; i < n; i++) {
426 cin >> x0;
427 cout << x0;
428 }
429 }
430 //
431 // and finally, check if an exact solution vector is supplied
432 //
433 if (rhs_type[2] == 'X') {
434 vec<Float> x;
435 for (unsigned int i = 0; i < n; i++) {
436 cin >> x;
437 cout << x;
438 }
439 }
440 }
441#endif // TODO
442}
field::size_type size_type
Definition branch.cc:430
void print_postscript(std::ostream &s, const csr< T, sequential > &a)
Definition csr.cc:99
void usage()
Definition csr.cc:261
see the communicator page for the full documentation
see the csr page for the full documentation
Definition csr.h:317
see the environment page for the full documentation
idiststream: see the diststream page for the full documentation
Definition diststream.h:336
see the vec page for the full documentation
Definition vec.h:79
void mtx(const Expr &a_expr, string a_name)
int main()
Definition field2bb.cc:58
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format hb
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format matlab
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format format format format dump
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format ml
This file is part of Rheolef.
t operator()(const t &a, const t &b)
Definition space.cc:386
STL namespace.
Definition sphere.icc:25