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 }
rheolef::csr< T, sequential >
Definition: csr.h:323
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 format matrix_market
Definition: iorheo-members.h:103
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 format format format format dump
Definition: iorheo-members.h:119
rheolef::value
rheolef::std value
ml
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
Definition: iorheo-members.h:106
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 hb
Definition: iorheo-members.h:100
usage
void usage()
Definition: csr.cc:261
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::vec
see the vec page for the full documentation
Definition: vec.h:79
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 matlab
Definition: iorheo-members.h:109
p
Definition: sphere.icc:25
main
int main(int argc, char **argv)
Definition: csr.cc:270
print_postscript
void print_postscript(std::ostream &s, const csr< T, sequential > &a)
Definition: csr.cc:99
a
Definition: diffusion_isotropic.h:25
rheolef::din
idiststream din
see the diststream page for the full documentation
Definition: diststream.h:427
rheolef::environment
see the environment page for the full documentation
Definition: environment.h:115
rheolef::csr
see the csr page for the full documentation
Definition: asr.h:31
mkgeo_sector.m
m
Definition: mkgeo_sector.sh:118
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::derr
odiststream derr(cerr)
see the diststream page for the full documentation
Definition: diststream.h:436
mkgeo_contraction.filter
filter
Definition: mkgeo_contraction.sh:299
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
mkgeo_ball.a
a
Definition: mkgeo_ball.sh:151
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::dout
odiststream dout(cout)
see the diststream page for the full documentation
Definition: diststream.h:430
rheolef::std
Definition: vec_expr_v2.h:391
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133
T
Expr1::float_type T
Definition: field_expr.h:218