MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
wraptest.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbwrap/wraptest.cc,v 1.81 2017/01/12 14:44:25 masarati Exp $ */
2 /*
3  * This library comes with MBDyn (C), a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  *
10  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
11  * via La Masa, 34 - 20156 Milano, Italy
12  * http://www.aero.polimi.it
13  *
14  * Changing this copyright notice is forbidden.
15  *
16  * This program is free software; you can redistribute it and/or modify
17  * it under the terms of the GNU General Public License as published by
18  * the Free Software Foundation (version 2 of the License).
19  *
20  *
21  * This program is distributed in the hope that it will be useful,
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  * GNU General Public License for more details.
25  *
26  * You should have received a copy of the GNU General Public License
27  * along with this program; if not, write to the Free Software
28  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29  */
30 
31 #include "mbconfig.h"
32 
33 #include <cstring>
34 #include <stdlib.h>
35 #include <unistd.h>
36 #include "ac/getopt.h"
37 
38 extern "C" {
39 #include <time.h>
40 #ifdef HAVE_SYS_TIMES_H
41 #include <sys/times.h>
42 #endif /* HAVE_SYS_TIMES_H */
43 }
44 
45 #include <iostream>
46 #include <fstream>
47 
48 #include "solman.h"
49 #include "submat.h"
50 #include "spmapmh.h"
51 #include "ccmh.h"
52 #include "dirccmh.h"
53 #include "y12wrap.h"
54 #include "harwrap.h"
55 #include "mschwrap.h"
56 #include "umfpackwrap.h"
57 #include "parsuperluwrap.h"
58 #include "superluwrap.h"
59 #include "lapackwrap.h"
60 #include "taucswrap.h"
61 #include "naivewrap.h"
62 #include "parnaivewrap.h"
63 #include "wsmpwrap.h"
64 
65 const char *solvers[] = {
66 #if defined(USE_Y12)
67  "y12",
68 #endif
69 #if defined(USE_UMFPACK)
70  "umfpack",
71 #endif
72 #if defined(USE_SUPERLU)
73  "superlu",
74 #endif
75 #if defined(USE_HARWELL)
76  "harwell",
77 #endif
78 #if defined(USE_MESCHACH)
79  "meschach",
80 #endif
81 #if defined(USE_LAPACK)
82  "lapack",
83 #endif
84 #if defined(USE_TAUCS)
85  "taucs",
86 #endif
87 #if defined(USE_WSMP)
88  "wsmp",
89 #endif
90  "naive",
91  NULL
92 };
93 
94 std::vector<doublereal> x_values;
95 std::vector<integer> row_values;
96 std::vector<integer> col_values;
97 std::vector<integer> acol_values;
98 SparseMatrixHandler *spM = 0;
99 
100 void SetupSystem(
101  const bool random,
102  char *random_args,
103  const bool singular,
104  const char *const matrixfilename,
105  const char *const filename,
106  MatrixHandler *const pM,
107  VectorHandler *const pV
108 ) {
109  VariableSubMatrixHandler SBMH(10, 10);
110  FullSubMatrixHandler& WM = SBMH.SetFull();
111 
112  std::ifstream ifile;
113  std::ofstream ofile;
114  int size = 3;
115 
116  if (filename == 0) {
117  if (random) {
118  int size = (*pM).iGetNumRows();
119  if (spM == 0) {
120  spM = new SpMapMatrixHandler(size);
121 
122  int halfband = 0;
123  if (random_args) {
124  char *next;
125 
126  halfband = (int)strtol(random_args,
127  &next, 10);
128  switch (next[0]) {
129  case '\0':
130  random_args = 0;
131  break;
132 
133  case ':':
134  random_args = &next[1];
135  break;
136 
137  default:
138  std::cerr << "unable to parse "
139  "<halfband> "
140  "from -r args"
141  << std::endl;
142  exit(EXIT_FAILURE);
143  }
144 
145  } else {
146  std::cout << "Halfband?" << std::endl;
147  std::cin >> halfband;
148  }
149 
150  int activcol = 0;
151  if (random_args) {
152  char *next;
153 
154  activcol = (int)strtol(random_args,
155  &next, 10);
156  switch (next[0]) {
157  case '\0':
158  random_args = 0;
159  break;
160 
161  case ':':
162  random_args = &next[1];
163  break;
164 
165  default:
166  std::cerr << "unable to parse "
167  "<activcol> "
168  "from -r args"
169  << std::endl;
170  exit(EXIT_FAILURE);
171  }
172 
173  } else {
174  std::cout << "Activcol?" << std::endl;
175  std::cin >> activcol;
176  }
177 
178  double sprfct = 0.9;
179  if (random_args) {
180  char *next;
181 
182  sprfct = (double)strtod(random_args,
183  &next);
184  switch (next[0]) {
185  case '\0':
186  random_args = 0;
187  break;
188 
189  case ':':
190  random_args = &next[1];
191  break;
192 
193  default:
194  std::cerr << "unable to parse "
195  "<sprfct> "
196  "from -r args"
197  << std::endl;
198  exit(EXIT_FAILURE);
199  }
200 
201  } else {
202  std::cout << "Sprfct (hint: 0.9)?"
203  << std::endl;
204  std::cin >> sprfct;
205  }
206 
207  for (int i = 0; i < size; i++) {
208  for (int k = (i - halfband) < 0 ? 0 : i - halfband; k < ((i + halfband) > size ? size : i + halfband); k++) {
209  if (((doublereal)rand())/RAND_MAX > sprfct) {
210  (*spM)(i+1, k+1) = 2.0*(((doublereal)rand())/RAND_MAX - 0.5);
211  }
212  }
213  }
214  for (int i = size - activcol; i < size; i++) {
215  for (int k = 0; k < size; k++) {
216  if (((doublereal)rand())/RAND_MAX > sprfct) {
217  (*spM)(k+1, i+1) = (*spM)(i+1, k+1) = 2.0*(((doublereal)rand())/RAND_MAX - 0.5);
218  }
219  }
220  }
221  for (int i = 0; i < size; i++) {
222  (*spM)(i+1, i+1) = 1;
223  }
224  spM->MakeIndexForm(x_values, row_values, col_values, acol_values, 1);
225  if (matrixfilename != 0) {
226  ofile.open(matrixfilename);
227  ofile << size << std::endl;
228  int n = x_values.size();
229  for (int i=0; i<n; i++) {
230  ofile << row_values[i] << " " <<
231  col_values[i] << " " <<
232  x_values[i] << std::endl;;
233  }
234  ofile.close();
235  }
236  }
237  int n = x_values.size();
238  for (int i=0; i<n; i++) {
239  (*pM)(row_values[i], col_values[i]) = x_values[i];
240  }
241  for (int i = 1; i <= size; i++) {
242  pV->PutCoef(i, pM->operator()(i, size));
243  }
244  } else {
245 
246  WM.ResizeReset(3, 3);
247  WM.PutRowIndex(1,1);
248  WM.PutRowIndex(2,2);
249  WM.PutRowIndex(3,3);
250  WM.PutColIndex(1,1);
251  WM.PutColIndex(2,2);
252  WM.PutColIndex(3,3);
253  WM.PutCoef(1, 1, 1.);
254  WM.PutCoef(2, 2, 2.);
255  WM.PutCoef(2, 3, -1.);
256  WM.PutCoef(3, 2, 11.);
257  WM.PutCoef(3, 1, 10.);
258  if (singular) {
259  WM.PutCoef(3, 3, 0.);
260  } else {
261  WM.PutCoef(3, 3, 3.);
262  }
263 
264  *pM += SBMH;
265 
266  pV->PutCoef(1, 1.);
267  pV->PutCoef(2, 1.);
268  pV->PutCoef(3, 1.);
269 
270  std::cout << *pM << "\n";
271  }
272  } else {
273  ifile.open(filename);
274  ifile >> size;
275 
276  integer count = 0;
277  integer row, col;
278  doublereal x;
279  while (ifile >> row >> col >> x) {
280  if (row > size || col > size) {
281  std::cerr << "Fatal read error of file" << filename << std::endl;
282  std::cerr << "size: " << size << std::endl;
283  std::cerr << "row: " << row << std::endl;
284  std::cerr << "col: " << col << std::endl;
285  std::cerr << "x: " << x << std::endl;
286  }
287  (*pM)(row, col) = x;
288  count++;
289  }
290 
291  ifile.close();
292 
293  for (integer i = 1; i <= size; i++) {
294  pV->PutCoef(i, pM->operator()(i, size));
295  }
296  std::cout << "\nThe matrix has "
297  << pM->iGetNumRows() << " rows, "
298  << pM->iGetNumCols() << " cols "
299  << "and " << count << " nonzeros\n" << std::endl;
300  }
301 }
302 
303 static inline unsigned long long rd_CPU_ts(void)
304 {
305  // See <http://www-unix.mcs.anl.gov/~kazutomo/rdtsc.html>
306  unsigned long long time = 0;
307 #if defined(__i386__)
308  __asm__ __volatile__( "rdtsc" : "=A" (time));
309 #elif defined(__x86_64__)
310  unsigned hi, lo;
311  __asm__ __volatile__ ("rdtsc" : "=a"(lo), "=d"(hi));
312  time = ( (unsigned long long)lo)|( ((unsigned long long)hi)<<32 );
313 #elif defined(__powerpc__)
314  unsigned long int upper, lower, tmp;
315  __asm__ volatile(
316  "0: \n"
317  "\tmftbu %0 \n"
318  "\tmftb %1 \n"
319  "\tmftbu %2 \n"
320  "\tcmpw %2,%0 \n"
321  "\tbne 0b \n"
322  : "=r"(upper),"=r"(lower),"=r"(tmp)
323  );
324  time = upper;
325  time = result << 32;
326  time = result|lower;
327 #endif
328  return time;
329 }
330 
331 static void
332 usage(int err)
333 {
334  std::cerr << "usage: wraptest "
335  "[-c] "
336  "[-d] "
337  "[-f <filename>] "
338  "[-m <solver>] "
339  << std::endl << "\t\t"
340  "[-o] "
341  "[-O <option>[=<value>]] "
342  "[-p <pivot>] "
343  << std::endl << "\t\t"
344  "[-r[<size>[:<halfband>[:<activcol>[:<sprfct>]]]]] "
345  << std::endl << "\t\t"
346  "[-s] "
347  "[-t <nthreads>] "
348  "[-T] "
349  "[-w <filename>] "
350  << std::endl;
351  std::cerr << " -c : if possible, use compressed column matrix format" << std::endl;
352  std::cerr << "\tfor the naive solver: use colamd" << std::endl;
353  std::cerr << " -d : if possible, use dir matrix format" << std::endl;
354  std::cerr << " -f <filename> : load the matrix from <filename>" << std::endl;
355  std::cerr << "\tIf the matrix is loaded from file the solution should be [0 0 .... 1]" << std::endl;
356  std::cerr << "\tThe file format is: size row col x row col x etc..." << std::endl;
357  std::cerr << " -m <solver> : {" << solvers[0];
358  for (int i = 1; solvers[i]; i++) {
359  std::cerr << "|" << solvers[i];
360  }
361  std::cerr << "}" << std::endl;
362  std::cerr << " -o : output of the solution" << std::endl;
363  std::cerr << " -O <option[=<value>]>" << std::endl
364  << "\tblocksize=<blocksize> (umfpack only)" << std::endl;
365  std::cerr << " -p <pivot> : if meaningful, use <pivot> threshold" << std::endl;
366  std::cerr << " -r[<size>[:<halfband>[:<activcol>[:<sprfct>]]]] :" << std::endl;
367  std::cerr << "\tgenerate a random matrix with <size>, <halfband>, <activcol>" << std::endl;
368  std::cerr << "\tand <sprfct> (prompts for values not provided)" << std::endl;
369  std::cerr << " -s : (singular) with the 3x3 matrix, do not set the element (3,3)" << std::endl;
370  std::cerr << " -t : with multithreaded solvers, use <nthreads> threads" << std::endl;
371  std::cerr << " -T : solve A^T x = b" << std::endl;
372  std::cerr << " -w : write the random matrix to <filename>" << std::endl;
373  exit(err);
374 }
375 
376 
377 enum {
378  COLAMD_PREORD,
379  MMDATA_PREORD
380 };
381 
382 int
383 main(int argc, char *argv[])
384 {
385  SolutionManager *pSM = NULL;
386  const char *solver =
387 #if defined(USE_UMFPACK)
388  "umfpack"
389 #elif defined(USE_Y12)
390  "y12"
391 #elif defined(USE_SUPERLU)
392  "superlu"
393 #elif defined(USE_HARWELL)
394  "harwell"
395 #elif defined(USE_MESCHACH)
396  "meschach"
397 #elif defined(USE_LAPACK)
398  "lapack"
399 #elif defined(USE_TAUCS)
400  "taucs"
401 #elif defined(USE_WSMP)
402  "wsmp"
403 #else
404  "naive"
405 #if 0
406  "no solver!!!"
407 #endif
408 #endif /* NO SOLVER !!! */
409  ;
410  clock_t start, end;
411  double cpu_time_used;
412 
413  char *filename = 0;
414  char *matrixfilename = 0;
415  std::ifstream file;
416  bool random(false);
417  char *random_args = 0;
418  bool cc(false);
419  bool dir(false);
420  unsigned nt = 1;
421  unsigned block_size = 0;
422  double dpivot = -1.;
423  double ddroptol = 0.;
424  bool singular(false);
425  bool output_solution(false);
426  bool transpose(false);
427  int size = 3;
428  long long tf;
429  unsigned preord = COLAMD_PREORD;
430  SolutionManager::ScaleWhen ms = SolutionManager::SCALEW_NEVER;
431 
432  while (1) {
433  int opt = getopt(argc, argv, "cdf:m:oO:p:P:r::st:Tw:");
434 
435  if (opt == EOF) {
436  break;
437  }
438 
439  switch (opt) {
440  case 'c':
441  cc = true;
442  break;
443 
444  case 'd':
445  dir = true;
446  break;
447 
448  case 'm':
449  solver = optarg;
450  break;
451 
452  case 's':
453  singular = true;
454  break;
455 
456  case 't':
457  nt = atoi(optarg);
458  if (nt < 1) {
459  nt = 1;
460  }
461  break;
462 
463  case 'T':
464  transpose = true;
465  break;
466 
467  case 'p':
468  dpivot = atof(optarg);
469  break;
470 
471  case 'P':
472  {
473  if (strncasecmp(optarg, "colamd", STRLENOF("colamd")) == 0) {
474  if ((strcasecmp(solver, "umfpack") != 0) &&
475  (strcasecmp(solver, "naive") != 0) &&
476  (strcasecmp(solver, "superlu") != 0)
477  ) {
478  std::cerr << "colamd preordering meaningful only for umfpack"
479  ", naive and superlu solvers" << std::endl;
480  exit(1);
481  }
482  preord = COLAMD_PREORD;
483  } else if (strncasecmp(optarg, "mmdata", STRLENOF("mmdata")) == 0) {
484  if (strcasecmp(solver, "superlu") != 0) {
485  std::cerr << "mmdata preordering meaningful only for superlu solver" << std::endl;
486  exit(1);
487  }
488  preord = MMDATA_PREORD;
489  } else {
490  std::cerr << "unrecognized option \"" << optarg << "\"" << std::endl;
491  exit(EXIT_FAILURE);
492  }
493 
494  break;
495  }
496 
497  case 'O':
498  {
499  if (strncasecmp(optarg, "blocksize=", STRLENOF("blocksize=")) == 0) {
500  char *next;
501 
502  if (strcasecmp(solver, "umfpack") != 0) {
503  std::cerr << "blocksize only meaningful for umfpack solver" << std::endl;
504  }
505 
506  optarg += STRLENOF("blocksize=");
507  block_size = (int)strtol(optarg, &next, 10);
508  if (next[0] != '\0') {
509  std::cerr << "unable to parse blocksize value" << std::endl;
510  exit(EXIT_FAILURE);
511  }
512  if (block_size < 1) {
513  block_size = 0;
514  }
515 
516  } else {
517  std::cerr << "unrecognized option \"" << optarg << "\"" << std::endl;
518  exit(EXIT_FAILURE);
519  }
520 
521  break;
522  }
523 
524  case 'f':
525  filename = optarg;
526  break;
527 
528  case 'o':
529  output_solution = true;
530  break;
531 
532  case 'r':
533  random = true;
534  random_args = optarg;
535  break;
536 
537  case 'w':
538  matrixfilename = optarg;
539  break;
540 
541 
542  default:
543  usage(EXIT_FAILURE);
544  }
545  }
546 
547  if (filename != 0) {
548  file.open(filename);
549  file >> size;
550  file.close();
551 
552  } else if (random) {
553  if (random_args) {
554  char *next;
555 
556  size = (int)strtol(random_args, &next, 10);
557  switch (next[0]) {
558  case '\0':
559  random_args = 0;
560  break;
561 
562  case ':':
563  random_args = &next[1];
564  break;
565 
566  default:
567  std::cerr << "unable to parse <size> "
568  "from -r args" << std::endl;
569  exit(EXIT_FAILURE);
570  }
571 
572  } else {
573  std::cout << "Matrix size?" << std::endl;
574  std::cin >> size;
575  }
576  }
577 
578  std::cerr << std::endl;
579 
580  if (strcasecmp(solver, "taucs") == 0) {
581 #ifdef USE_TAUCS
582  std::cerr << "Taucs solver"
583  if (dir) {
584  std::cerr << " with dir matrix"
585  typedef TaucsSparseCCSolutionManager<DirCColMatrixHandler<0> > CCMH;
586  SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size));
587 
588  } else if (cc) {
589  std::cerr << " with cc matrix"
590  typedef TaucsSparseCCSolutionManager<CColMatrixHandler<0> > CCMH;
591  SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size));
592 
593  } else {
594  SAFENEWWITHCONSTRUCTOR(pSM, TaucsSparseSolutionManager,
595  TaucsSparseSolutionManager(size));
596  }
597  std::cerr << std::endl;
598 #else /* !USE_TAUCS */
599  std::cerr << "need --with-taucs to use Taucs library sparse solver"
600  << std::endl;
601  usage(EXIT_FAILURE);
602 #endif /* !USE_LAPACK */
603 
604  } else if (strcasecmp(solver, "lapack") == 0) {
605 #ifdef USE_LAPACK
606  std::cerr << "Lapack solver" << std::endl;
607  SAFENEWWITHCONSTRUCTOR(pSM, LapackSolutionManager,
608  LapackSolutionManager(size));
609 #else /* !USE_LAPACK */
610  std::cerr << "need --with-lapack to use Lapack library dense solver"
611  << std::endl;
612  usage(EXIT_FAILURE);
613 #endif /* !USE_LAPACK */
614 
615  } else if (strcasecmp(solver, "superlu") == 0) {
616 #ifdef USE_SUPERLU
617  if (dpivot == -1.) {
618  dpivot = 1.;
619  }
620  unsigned pre = SuperLUSolver::SUPERLU_COLAMD;
621  if (preord == MMDATA_PREORD) {
622  pre = SuperLUSolver::SUPERLU_MMDATA;
623  std::cerr << "Using MMDATA preordering" << std::endl;
624  } else {
625  pre = SuperLUSolver::SUPERLU_COLAMD;
626  std::cerr << "Using colamd preordering" << std::endl;
627  }
628  if (nt > 1) {
629 #ifdef USE_SUPERLU_MT
630  std::cerr << "Multithreaded SuperLU solver";
631  if (pre != SuperLUSolver::SUPERLU_COLAMD) {
632  std::cerr << std::end <<
633  "ERROR, mmdata preordering available only for scalar superlu" <<
634  std::endl;
635  exit(1);
636  }
637  if (dir) {
638  std::cerr << " with dir matrix";
639  typedef ParSuperLUSparseCCSolutionManager<DirCColMatrixHandler<0> > CCMH;
640  SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(nt, size, dpivot));
641  } else if (cc) {
642  std::cerr << " with cc matrix";
643  typedef ParSuperLUSparseCCSolutionManager<CColMatrixHandler<0> > CCMH;
644  SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(nt, size, dpivot));
645  } else {
646  SAFENEWWITHCONSTRUCTOR(pSM, ParSuperLUSparseSolutionManager,
647  ParSuperLUSparseSolutionManager(nt, size, dpivot));
648  }
649  std::cerr << " using " << nt << " threads" << std::endl;
650 #else /* !USE_SUPERLU_MT */
651  silent_cerr("multithread SuperLU solver support not compiled; "
652  << std::endl);
653  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
654  } else {
655  std::cerr << "SuperLU solver";
656  if (dir) {
657  std::cerr << " with dir matrix";
658  typedef SuperLUSparseCCSolutionManager<DirCColMatrixHandler<0> > CCMH;
659  SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size, dpivot, pre));
660  } else if (cc) {
661  std::cerr << " with cc matrix";
662  typedef SuperLUSparseCCSolutionManager<CColMatrixHandler<0> > CCMH;
663  SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size, dpivot, pre));
664  } else {
665  SAFENEWWITHCONSTRUCTOR(pSM, SuperLUSparseSolutionManager,
666  SuperLUSparseSolutionManager(size, dpivot, pre));
667  }
668 #endif /* !USE_SUPERLU_MT */
669  }
670  std::cerr << std::endl;
671 #else /* !USE_SUPERLU */
672  std::cerr << "need --with-superlu to use SuperLU library"
673  << std::endl;
674  usage(EXIT_FAILURE);
675 #endif /* !USE_SUPERLU */
676 
677  } else if (strcasecmp(solver, "y12") == 0) {
678 #ifdef USE_Y12
679  std::cerr << "y12 solver";
680  if (dir) {
681  std::cerr << " with dir matrix";
682  typedef Y12SparseCCSolutionManager<DirCColMatrixHandler<1> > CCMH;
683  SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size));
684 
685  } else if (cc) {
686  std::cerr << " with cc matrix";
687  typedef Y12SparseCCSolutionManager<CColMatrixHandler<1> > CCMH;
688  SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size));
689 
690  } else {
691  SAFENEWWITHCONSTRUCTOR(pSM, Y12SparseSolutionManager,
692  Y12SparseSolutionManager(size));
693  }
694  std::cerr << std::endl;
695 #else /* !USE_Y12 */
696  std::cerr << "need --with-y12 to use y12m library"
697  << std::endl;
698  usage(EXIT_FAILURE);
699 #endif /* !USE_Y12 */
700 
701  } else if (strcasecmp(solver, "harwell") == 0) {
702 #ifdef USE_HARWELL
703  std::cerr << "Harwell solver" << std::endl;
704  SAFENEWWITHCONSTRUCTOR(pSM, HarwellSparseSolutionManager,
705  HarwellSparseSolutionManager(size));
706 #else /* !USE_HARWELL */
707  std::cerr << "need --with-harwell to use HSL library"
708  << std::endl;
709  usage(EXIT_FAILURE);
710 #endif /* !USE_HARWELL */
711 
712  } else if (strcasecmp(solver, "meschach") == 0) {
713 #ifdef USE_MESCHACH
714  std::cerr << "Meschach solver" << std::endl;
715  SAFENEWWITHCONSTRUCTOR(pSM, MeschachSparseSolutionManager,
716  MeschachSparseSolutionManager(size));
717 #else /* !USE_MESCHACH */
718  std::cerr << "need --with-meschach to use Meschach library"
719  << std::endl;
720  usage(EXIT_FAILURE);
721 #endif /* !USE_MESCHACH */
722  } else if (strcasecmp(solver, "umfpack") == 0
723  || strcasecmp(solver, "umfpack3") == 0) {
724 #ifdef USE_UMFPACK
725  std::cerr << "Umfpack solver";
726  if (dir) {
727  std::cerr << " with dir matrix";
728  typedef UmfpackSparseCCSolutionManager<DirCColMatrixHandler<0> > CCMH;
729  SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size, dpivot, ddroptol, block_size));
730 
731  } else if (cc) {
732  std::cerr << " with cc matrix";
733  typedef UmfpackSparseCCSolutionManager<CColMatrixHandler<0> > CCMH;
734  SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size, dpivot, ddroptol, block_size));
735 
736  } else {
737  SAFENEWWITHCONSTRUCTOR(pSM,
738  UmfpackSparseSolutionManager,
739  UmfpackSparseSolutionManager(size, dpivot, ddroptol, block_size));
740  }
741  std::cerr << std::endl;
742 #else /* !USE_UMFPACK */
743  std::cerr << "need --with-umfpack to use Umfpack library"
744  << std::endl;
745  usage(EXIT_FAILURE);
746 #endif /* !USE_UMFPACK */
747 
748  } else if (strcasecmp(solver, "wsmp") == 0) {
749 #ifdef USE_WSMP
750  std::cerr << "Wsmp solver";
751  if (dir) {
752  std::cerr << " with dir matrix";
753  typedef WsmpSparseCCSolutionManager<DirCColMatrixHandler<0> > CCMH;
754  SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size, dpivot, block_size, nt));
755 
756  } else if (cc) {
757  std::cerr << " with cc matrix";
758  typedef WsmpSparseCCSolutionManager<CColMatrixHandler<0> > CCMH;
759  SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size, dpivot, block_size, nt));
760 
761  } else {
762  SAFENEWWITHCONSTRUCTOR(pSM,
763  WsmpSparseSolutionManager,
764  WsmpSparseSolutionManager(size, dpivot, block_size, nt));
765  }
766  std::cerr << " using " << nt << " threads " << std::endl;
767  std::cerr << std::endl;
768 #else /* !USE_WSMP */
769  std::cerr << "need --with-wsmp to use Wsmp library"
770  << std::endl;
771  usage(EXIT_FAILURE);
772 #endif /* !USE_WSMP */
773 
774  } else if (strcasecmp(solver, "naive") == 0) {
775  std::cerr << "Naive solver";
776  if (dpivot == -1.) {
777  dpivot = 1.E-8;
778  }
779  if (cc) {
780  std::cerr << " with Colamd ordering";
781  if (nt > 1) {
782 #ifdef USE_NAIVE_MULTITHREAD
783  SAFENEWWITHCONSTRUCTOR(pSM,
784  ParNaiveSparsePermSolutionManager,
785  ParNaiveSparsePermSolutionManager(nt, size, dpivot));
786 #else
787  silent_cerr("multithread naive solver support not compiled; "
788  "you can configure --enable-multithread-naive "
789  "on a linux ix86 to get it"
790  << std::endl);
791  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
792 #endif /* USE_NAIVE_MULTITHREAD */
793  } else {
794  SAFENEWWITHCONSTRUCTOR(pSM,
795  NaiveSparsePermSolutionManager<Colamd_ordering>,
796  NaiveSparsePermSolutionManager<Colamd_ordering>(size, dpivot, ms));
797  }
798  } else {
799  if (nt > 1) {
800 #ifdef USE_NAIVE_MULTITHREAD
801  SAFENEWWITHCONSTRUCTOR(pSM,
802  ParNaiveSparseSolutionManager,
803  ParNaiveSparseSolutionManager(nt, size, dpivot));
804 #else
805  silent_cerr("multithread naive solver support not compiled; "
806  "you can configure --enable-multithread-naive "
807  "on a linux ix86 to get it"
808  << std::endl);
809  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
810 #endif /* USE_NAIVE_MULTITHREAD */
811  } else {
812  SAFENEWWITHCONSTRUCTOR(pSM,
813  NaiveSparseSolutionManager,
814  NaiveSparseSolutionManager(size, dpivot, ms));
815  }
816  }
817  std::cerr << " using " << nt << " threads " << std::endl;
818  } else {
819  std::cerr << "unknown solver '" << solver << "'" << std::endl;
820  usage(EXIT_FAILURE);
821  }
822 
823  pSM->MatrReset();
824 
825  MatrixHandler *pM = pSM->pMatHdl();
826  VectorHandler *pV = pSM->pResHdl();
827  VectorHandler *px = pSM->pSolHdl();
828 
829  pM->Reset();
830 
831  SetupSystem(random, random_args, singular,
832  matrixfilename, filename, pM, pV);
833 
834 #ifdef HAVE_TIMES
835  clock_t ct = 0;
836  struct tms tmsbuf;
837 #endif // HAVE_TIMES
838  try {
839  start = clock();
840  tf = rd_CPU_ts();
841 #ifdef HAVE_TIMES
842  times(&tmsbuf);
843  ct = tmsbuf.tms_utime + tmsbuf.tms_cutime
844  + tmsbuf.tms_stime + tmsbuf.tms_cstime;
845 #endif // HAVE_TIMES
846  if (transpose) {
847  pSM->SolveT();
848  } else {
849  pSM->Solve();
850  }
851  tf = rd_CPU_ts() - tf;
852  end = clock();
853 #ifdef HAVE_TIMES
854  times(&tmsbuf);
855  ct = tmsbuf.tms_utime + tmsbuf.tms_cutime
856  + tmsbuf.tms_stime + tmsbuf.tms_cstime - ct;
857 #endif // HAVE_TIMES
858  } catch (...) {
859  exit(EXIT_FAILURE);
860  }
861  if (output_solution) {
862  for (int i = 1; i <= size; i++) {
863  std::cout << "\tsol[" << i << "] = " << px->operator()(i)
864  << std::endl;
865  }
866  }
867  cpu_time_used = ((double) (end - start));
868  cpu_time_used = cpu_time_used / CLOCKS_PER_SEC;
869  std::cout << "Clock tics to solve: " << tf << std::endl;
870 #ifdef HAVE_TIMES
871  std::cout << "Clock tics to solve: " << ct << std::endl;
872 #endif // HAVE_TIMES
873  std::cout << "Time to solve: " << cpu_time_used << std::endl;
874 
875 
876  std::cout << "\nSecond solve:\n";
877 
878  pSM->MatrReset();
879  pM = pSM->pMatHdl();
880 
881  pM->Reset();
882 
883  SetupSystem(random, random_args, singular, 0, filename, pM, pV);
884 
885  try {
886  start = clock();
887  tf = rd_CPU_ts();
888 #ifdef HAVE_TIMES
889  times(&tmsbuf);
890  ct = tmsbuf.tms_utime + tmsbuf.tms_cutime
891  + tmsbuf.tms_stime + tmsbuf.tms_cstime;
892 #endif // HAVE_TIMES
893  if (transpose) {
894  pSM->SolveT();
895  } else {
896  pSM->Solve();
897  }
898  tf = rd_CPU_ts() - tf;
899  end = clock();
900 #ifdef HAVE_TIMES
901  times(&tmsbuf);
902  ct = tmsbuf.tms_utime + tmsbuf.tms_cutime
903  + tmsbuf.tms_stime + tmsbuf.tms_cstime - ct;
904 #endif // HAVE_TIMES
905  } catch (...) {
906  exit(EXIT_FAILURE);
907  }
908 
909  if (output_solution) {
910  for (int i = 1; i <= size; i++) {
911  std::cout << "\tsol[" << i << "] = " << px->operator()(i)
912  << std::endl;
913  }
914  }
915  cpu_time_used = ((double) (end - start));
916  cpu_time_used = cpu_time_used / CLOCKS_PER_SEC;
917  std::cout << "Clock tics to solve: " << tf << std::endl;
918 #ifdef HAVE_TIMES
919  std::cout << "Clock tics to solve: " << ct << std::endl;
920 #endif // HAVE_TIMES
921  std::cout << "Time to solve: " << cpu_time_used << std::endl;
922 
923  return 0;
924 }
925 
SolutionManager::pResHdl
virtual VectorHandler * pResHdl(void) const =0
FullSubMatrixHandler::PutColIndex
void PutColIndex(integer iSubCol, integer iCol)
Definition: submat.h:325
SparseMatrixHandler::MakeIndexForm
virtual integer MakeIndexForm(doublereal *const Ax, integer *const Arow, integer *const Acol, integer *const AcolSt, int offset=0) const =0
MatrixHandler::iGetNumCols
virtual integer iGetNumCols(void) const =0
filename
static char * filename
Definition: ann_sf.c:71
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
VariableSubMatrixHandler::SetFull
FullSubMatrixHandler & SetFull(void)
Definition: submat.h:1168
solvers
const char * solvers[]
Definition: wraptest.cc:65
FullSubMatrixHandler::PutCoef
void PutCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:672
x_values
std::vector< doublereal > x_values
Definition: wraptest.cc:94
spM
SparseMatrixHandler * spM
Definition: wraptest.cc:98
MatrixHandler::Reset
virtual void Reset(void)=0
rd_CPU_ts
static unsigned long long rd_CPU_ts(void)
Definition: wraptest.cc:303
solver
const LinSol::solver_t solver[]
Definition: linsol.cc:58
row_values
std::vector< integer > row_values
Definition: wraptest.cc:95
SolutionManager::SolveT
virtual void SolveT(void)
Definition: solman.cc:78
SolutionManager::pMatHdl
virtual MatrixHandler * pMatHdl(void) const =0
main
int main(int argc, char *argv[])
Definition: wraptest.cc:383
acol_values
std::vector< integer > acol_values
Definition: wraptest.cc:97
file
static char * file
Definition: ann_in.c:101
count
static int count
Definition: modules.cc:41
SolutionManager::MatrReset
virtual void MatrReset(void)=0
SetupSystem
void SetupSystem(const bool random, char *random_args, const bool singular, const char *const matrixfilename, const char *const filename, MatrixHandler *const pM, VectorHandler *const pV)
Definition: wraptest.cc:100
SolutionManager::Solve
virtual void Solve(void)=0
VectorHandler::PutCoef
virtual void PutCoef(integer iRow, const doublereal &dCoef)=0
SAFENEWWITHCONSTRUCTOR
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
FullSubMatrixHandler::ResizeReset
virtual void ResizeReset(integer, integer)
Definition: submat.cc:182
getopt
int getopt(int argc, char *const argv[], const char *opts)
Definition: getopt.c:93
STRLENOF
#define STRLENOF(s)
Definition: mbdyn.h:166
col_values
std::vector< integer > col_values
Definition: wraptest.cc:96
FullSubMatrixHandler::PutRowIndex
void PutRowIndex(integer iSubRow, integer iRow)
Definition: submat.h:311
optarg
char * optarg
Definition: getopt.c:74
SolutionManager::pSolHdl
virtual VectorHandler * pSolHdl(void) const =0
doublereal
double doublereal
Definition: colamd.c:52
integer
long int integer
Definition: colamd.c:51
usage
static void usage(int err)
Definition: wraptest.cc:332
MatrixHandler::iGetNumRows
virtual integer iGetNumRows(void) const =0