Rheolef  7.1
an efficient C++ finite element environment
csr.cc
Go to the documentation of this file.
1 // the csr unix command
22 // author: Pierre.Saramito@imag.fr
23 // date: 12 may 1997, updated 9 oct 2012
24 
25 namespace rheolef {
84 } // namespace rheolef
85 
86 #include "rheolef/linalg.h"
87 using namespace rheolef;
88 using namespace std;
89 
90 // ------------------------------------------------------------------------
91 // print function
92 // ------------------------------------------------------------------------
93 template <class T>
94 struct log10_abs_or_zero {
95  T operator()(T x) { return ((x == T(0.)) ? T(0.) : log10(xabs(x))); }
96 };
97 template <class T>
98 void
100  std::ostream& s,
101  const csr<T,sequential>& a)
102 {
103  typedef typename csr<T,sequential>::size_type size_type;
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 // ------------------------------------------------------------------------
261 void usage()
262 {
263  derr << "csr: usage: csr "
264  << "[-ps|-matlab] "
265  << "[-name string] "
266  << "file[.mtx]"
267  << endl;
268  exit (1);
269 }
270 int 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:425
void print_postscript(std::ostream &s, const csr< T, sequential > &a)
Definition: csr.cc:99
void usage()
Definition: csr.cc:261
int main(int argc, char **argv)
Definition: csr.cc:270
see the csr page for the full documentation
Definition: csr.h:317
see the environment page for the full documentation
Definition: environment.h:104
see the vec page for the full documentation
Definition: vec.h:79
size_t size_type
Definition: basis_get.cc:76
idiststream din
see the diststream page for the full documentation
Definition: diststream.h:427
rheolef::std value
odiststream dout(cout)
see the diststream page for the full documentation
Definition: diststream.h:430
odiststream derr(cerr)
see the diststream page for the full documentation
Definition: diststream.h:436
Expr1::float_type T
Definition: field_expr.h:261
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 matrix_market
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.
Definition: sphere.icc:25