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
SetupSystem(const bool random,char * random_args,const bool singular,const char * const matrixfilename,const char * const filename,MatrixHandler * const pM,VectorHandler * const pV)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
rd_CPU_ts(void)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 << 32) | lower;
325 #endif
326 return time;
327 }
328
329 static void
usage(int err)330 usage(int err)
331 {
332 std::cerr << "usage: wraptest "
333 "[-c] "
334 "[-d] "
335 "[-f <filename>] "
336 "[-m <solver>] "
337 << std::endl << "\t\t"
338 "[-o] "
339 "[-O <option>[=<value>]] "
340 "[-p <pivot>] "
341 << std::endl << "\t\t"
342 "[-r[<size>[:<halfband>[:<activcol>[:<sprfct>]]]]] "
343 << std::endl << "\t\t"
344 "[-s] "
345 "[-t <nthreads>] "
346 "[-T] "
347 "[-w <filename>] "
348 << std::endl;
349 std::cerr << " -c : if possible, use compressed column matrix format" << std::endl;
350 std::cerr << "\tfor the naive solver: use colamd" << std::endl;
351 std::cerr << " -d : if possible, use dir matrix format" << std::endl;
352 std::cerr << " -f <filename> : load the matrix from <filename>" << std::endl;
353 std::cerr << "\tIf the matrix is loaded from file the solution should be [0 0 .... 1]" << std::endl;
354 std::cerr << "\tThe file format is: size row col x row col x etc..." << std::endl;
355 std::cerr << " -m <solver> : {" << solvers[0];
356 for (int i = 1; solvers[i]; i++) {
357 std::cerr << "|" << solvers[i];
358 }
359 std::cerr << "}" << std::endl;
360 std::cerr << " -o : output of the solution" << std::endl;
361 std::cerr << " -O <option[=<value>]>" << std::endl
362 << "\tblocksize=<blocksize> (umfpack only)" << std::endl;
363 std::cerr << " -p <pivot> : if meaningful, use <pivot> threshold" << std::endl;
364 std::cerr << " -r[<size>[:<halfband>[:<activcol>[:<sprfct>]]]] :" << std::endl;
365 std::cerr << "\tgenerate a random matrix with <size>, <halfband>, <activcol>" << std::endl;
366 std::cerr << "\tand <sprfct> (prompts for values not provided)" << std::endl;
367 std::cerr << " -s : (singular) with the 3x3 matrix, do not set the element (3,3)" << std::endl;
368 std::cerr << " -t : with multithreaded solvers, use <nthreads> threads" << std::endl;
369 std::cerr << " -T : solve A^T x = b" << std::endl;
370 std::cerr << " -w : write the random matrix to <filename>" << std::endl;
371 exit(err);
372 }
373
374
375 enum {
376 COLAMD_PREORD,
377 MMDATA_PREORD
378 };
379
380 int
main(int argc,char * argv[])381 main(int argc, char *argv[])
382 {
383 SolutionManager *pSM = NULL;
384 const char *solver =
385 #if defined(USE_UMFPACK)
386 "umfpack"
387 #elif defined(USE_Y12)
388 "y12"
389 #elif defined(USE_SUPERLU)
390 "superlu"
391 #elif defined(USE_HARWELL)
392 "harwell"
393 #elif defined(USE_MESCHACH)
394 "meschach"
395 #elif defined(USE_LAPACK)
396 "lapack"
397 #elif defined(USE_TAUCS)
398 "taucs"
399 #elif defined(USE_WSMP)
400 "wsmp"
401 #else
402 "naive"
403 #if 0
404 "no solver!!!"
405 #endif
406 #endif /* NO SOLVER !!! */
407 ;
408 clock_t start, end;
409 double cpu_time_used;
410
411 char *filename = 0;
412 char *matrixfilename = 0;
413 std::ifstream file;
414 bool random(false);
415 char *random_args = 0;
416 bool cc(false);
417 bool dir(false);
418 unsigned nt = 1;
419 unsigned block_size = 0;
420 double dpivot = -1.;
421 double ddroptol = 0.;
422 bool singular(false);
423 bool output_solution(false);
424 bool transpose(false);
425 int size = 3;
426 long long tf;
427 unsigned preord = COLAMD_PREORD;
428 SolutionManager::ScaleWhen ms = SolutionManager::SCALEW_NEVER;
429
430 while (1) {
431 int opt = getopt(argc, argv, "cdf:m:oO:p:P:r::st:Tw:");
432
433 if (opt == EOF) {
434 break;
435 }
436
437 switch (opt) {
438 case 'c':
439 cc = true;
440 break;
441
442 case 'd':
443 dir = true;
444 break;
445
446 case 'm':
447 solver = optarg;
448 break;
449
450 case 's':
451 singular = true;
452 break;
453
454 case 't':
455 nt = atoi(optarg);
456 if (nt < 1) {
457 nt = 1;
458 }
459 break;
460
461 case 'T':
462 transpose = true;
463 break;
464
465 case 'p':
466 dpivot = atof(optarg);
467 break;
468
469 case 'P':
470 {
471 if (strncasecmp(optarg, "colamd", STRLENOF("colamd")) == 0) {
472 if ((strcasecmp(solver, "umfpack") != 0) &&
473 (strcasecmp(solver, "naive") != 0) &&
474 (strcasecmp(solver, "superlu") != 0)
475 ) {
476 std::cerr << "colamd preordering meaningful only for umfpack"
477 ", naive and superlu solvers" << std::endl;
478 exit(1);
479 }
480 preord = COLAMD_PREORD;
481 } else if (strncasecmp(optarg, "mmdata", STRLENOF("mmdata")) == 0) {
482 if (strcasecmp(solver, "superlu") != 0) {
483 std::cerr << "mmdata preordering meaningful only for superlu solver" << std::endl;
484 exit(1);
485 }
486 preord = MMDATA_PREORD;
487 } else {
488 std::cerr << "unrecognized option \"" << optarg << "\"" << std::endl;
489 exit(EXIT_FAILURE);
490 }
491
492 break;
493 }
494
495 case 'O':
496 {
497 if (strncasecmp(optarg, "blocksize=", STRLENOF("blocksize=")) == 0) {
498 char *next;
499
500 if (strcasecmp(solver, "umfpack") != 0) {
501 std::cerr << "blocksize only meaningful for umfpack solver" << std::endl;
502 }
503
504 optarg += STRLENOF("blocksize=");
505 block_size = (int)strtol(optarg, &next, 10);
506 if (next[0] != '\0') {
507 std::cerr << "unable to parse blocksize value" << std::endl;
508 exit(EXIT_FAILURE);
509 }
510 if (block_size < 1) {
511 block_size = 0;
512 }
513
514 } else {
515 std::cerr << "unrecognized option \"" << optarg << "\"" << std::endl;
516 exit(EXIT_FAILURE);
517 }
518
519 break;
520 }
521
522 case 'f':
523 filename = optarg;
524 break;
525
526 case 'o':
527 output_solution = true;
528 break;
529
530 case 'r':
531 random = true;
532 random_args = optarg;
533 break;
534
535 case 'w':
536 matrixfilename = optarg;
537 break;
538
539
540 default:
541 usage(EXIT_FAILURE);
542 }
543 }
544
545 if (filename != 0) {
546 file.open(filename);
547 file >> size;
548 file.close();
549
550 } else if (random) {
551 if (random_args) {
552 char *next;
553
554 size = (int)strtol(random_args, &next, 10);
555 switch (next[0]) {
556 case '\0':
557 random_args = 0;
558 break;
559
560 case ':':
561 random_args = &next[1];
562 break;
563
564 default:
565 std::cerr << "unable to parse <size> "
566 "from -r args" << std::endl;
567 exit(EXIT_FAILURE);
568 }
569
570 } else {
571 std::cout << "Matrix size?" << std::endl;
572 std::cin >> size;
573 }
574 }
575
576 std::cerr << std::endl;
577
578 if (strcasecmp(solver, "taucs") == 0) {
579 #ifdef USE_TAUCS
580 std::cerr << "Taucs solver"
581 if (dir) {
582 std::cerr << " with dir matrix"
583 typedef TaucsSparseCCSolutionManager<DirCColMatrixHandler<0> > CCMH;
584 SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size));
585
586 } else if (cc) {
587 std::cerr << " with cc matrix"
588 typedef TaucsSparseCCSolutionManager<CColMatrixHandler<0> > CCMH;
589 SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size));
590
591 } else {
592 SAFENEWWITHCONSTRUCTOR(pSM, TaucsSparseSolutionManager,
593 TaucsSparseSolutionManager(size));
594 }
595 std::cerr << std::endl;
596 #else /* !USE_TAUCS */
597 std::cerr << "need --with-taucs to use Taucs library sparse solver"
598 << std::endl;
599 usage(EXIT_FAILURE);
600 #endif /* !USE_LAPACK */
601
602 } else if (strcasecmp(solver, "lapack") == 0) {
603 #ifdef USE_LAPACK
604 std::cerr << "Lapack solver" << std::endl;
605 SAFENEWWITHCONSTRUCTOR(pSM, LapackSolutionManager,
606 LapackSolutionManager(size));
607 #else /* !USE_LAPACK */
608 std::cerr << "need --with-lapack to use Lapack library dense solver"
609 << std::endl;
610 usage(EXIT_FAILURE);
611 #endif /* !USE_LAPACK */
612
613 } else if (strcasecmp(solver, "superlu") == 0) {
614 #ifdef USE_SUPERLU
615 if (dpivot == -1.) {
616 dpivot = 1.;
617 }
618 unsigned pre = SuperLUSolver::SUPERLU_COLAMD;
619 if (preord == MMDATA_PREORD) {
620 pre = SuperLUSolver::SUPERLU_MMDATA;
621 std::cerr << "Using MMDATA preordering" << std::endl;
622 } else {
623 pre = SuperLUSolver::SUPERLU_COLAMD;
624 std::cerr << "Using colamd preordering" << std::endl;
625 }
626 if (nt > 1) {
627 #ifdef USE_SUPERLU_MT
628 std::cerr << "Multithreaded SuperLU solver";
629 if (pre != SuperLUSolver::SUPERLU_COLAMD) {
630 std::cerr << std::end <<
631 "ERROR, mmdata preordering available only for scalar superlu" <<
632 std::endl;
633 exit(1);
634 }
635 if (dir) {
636 std::cerr << " with dir matrix";
637 typedef ParSuperLUSparseCCSolutionManager<DirCColMatrixHandler<0> > CCMH;
638 SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(nt, size, dpivot));
639 } else if (cc) {
640 std::cerr << " with cc matrix";
641 typedef ParSuperLUSparseCCSolutionManager<CColMatrixHandler<0> > CCMH;
642 SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(nt, size, dpivot));
643 } else {
644 SAFENEWWITHCONSTRUCTOR(pSM, ParSuperLUSparseSolutionManager,
645 ParSuperLUSparseSolutionManager(nt, size, dpivot));
646 }
647 std::cerr << " using " << nt << " threads" << std::endl;
648 #else /* !USE_SUPERLU_MT */
649 silent_cerr("multithread SuperLU solver support not compiled; "
650 << std::endl);
651 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
652 } else {
653 std::cerr << "SuperLU solver";
654 if (dir) {
655 std::cerr << " with dir matrix";
656 typedef SuperLUSparseCCSolutionManager<DirCColMatrixHandler<0> > CCMH;
657 SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size, dpivot, pre));
658 } else if (cc) {
659 std::cerr << " with cc matrix";
660 typedef SuperLUSparseCCSolutionManager<CColMatrixHandler<0> > CCMH;
661 SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size, dpivot, pre));
662 } else {
663 SAFENEWWITHCONSTRUCTOR(pSM, SuperLUSparseSolutionManager,
664 SuperLUSparseSolutionManager(size, dpivot, pre));
665 }
666 #endif /* !USE_SUPERLU_MT */
667 }
668 std::cerr << std::endl;
669 #else /* !USE_SUPERLU */
670 std::cerr << "need --with-superlu to use SuperLU library"
671 << std::endl;
672 usage(EXIT_FAILURE);
673 #endif /* !USE_SUPERLU */
674
675 } else if (strcasecmp(solver, "y12") == 0) {
676 #ifdef USE_Y12
677 std::cerr << "y12 solver";
678 if (dir) {
679 std::cerr << " with dir matrix";
680 typedef Y12SparseCCSolutionManager<DirCColMatrixHandler<1> > CCMH;
681 SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size));
682
683 } else if (cc) {
684 std::cerr << " with cc matrix";
685 typedef Y12SparseCCSolutionManager<CColMatrixHandler<1> > CCMH;
686 SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size));
687
688 } else {
689 SAFENEWWITHCONSTRUCTOR(pSM, Y12SparseSolutionManager,
690 Y12SparseSolutionManager(size));
691 }
692 std::cerr << std::endl;
693 #else /* !USE_Y12 */
694 std::cerr << "need --with-y12 to use y12m library"
695 << std::endl;
696 usage(EXIT_FAILURE);
697 #endif /* !USE_Y12 */
698
699 } else if (strcasecmp(solver, "harwell") == 0) {
700 #ifdef USE_HARWELL
701 std::cerr << "Harwell solver" << std::endl;
702 SAFENEWWITHCONSTRUCTOR(pSM, HarwellSparseSolutionManager,
703 HarwellSparseSolutionManager(size));
704 #else /* !USE_HARWELL */
705 std::cerr << "need --with-harwell to use HSL library"
706 << std::endl;
707 usage(EXIT_FAILURE);
708 #endif /* !USE_HARWELL */
709
710 } else if (strcasecmp(solver, "meschach") == 0) {
711 #ifdef USE_MESCHACH
712 std::cerr << "Meschach solver" << std::endl;
713 SAFENEWWITHCONSTRUCTOR(pSM, MeschachSparseSolutionManager,
714 MeschachSparseSolutionManager(size));
715 #else /* !USE_MESCHACH */
716 std::cerr << "need --with-meschach to use Meschach library"
717 << std::endl;
718 usage(EXIT_FAILURE);
719 #endif /* !USE_MESCHACH */
720 } else if (strcasecmp(solver, "umfpack") == 0
721 || strcasecmp(solver, "umfpack3") == 0) {
722 #ifdef USE_UMFPACK
723 std::cerr << "Umfpack solver";
724 if (dir) {
725 std::cerr << " with dir matrix";
726 typedef UmfpackSparseCCSolutionManager<DirCColMatrixHandler<0> > CCMH;
727 SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size, dpivot, ddroptol, block_size));
728
729 } else if (cc) {
730 std::cerr << " with cc matrix";
731 typedef UmfpackSparseCCSolutionManager<CColMatrixHandler<0> > CCMH;
732 SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size, dpivot, ddroptol, block_size));
733
734 } else {
735 SAFENEWWITHCONSTRUCTOR(pSM,
736 UmfpackSparseSolutionManager,
737 UmfpackSparseSolutionManager(size, dpivot, ddroptol, block_size));
738 }
739 std::cerr << std::endl;
740 #else /* !USE_UMFPACK */
741 std::cerr << "need --with-umfpack to use Umfpack library"
742 << std::endl;
743 usage(EXIT_FAILURE);
744 #endif /* !USE_UMFPACK */
745
746 } else if (strcasecmp(solver, "wsmp") == 0) {
747 #ifdef USE_WSMP
748 std::cerr << "Wsmp solver";
749 if (dir) {
750 std::cerr << " with dir matrix";
751 typedef WsmpSparseCCSolutionManager<DirCColMatrixHandler<0> > CCMH;
752 SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size, dpivot, block_size, nt));
753
754 } else if (cc) {
755 std::cerr << " with cc matrix";
756 typedef WsmpSparseCCSolutionManager<CColMatrixHandler<0> > CCMH;
757 SAFENEWWITHCONSTRUCTOR(pSM, CCMH, CCMH(size, dpivot, block_size, nt));
758
759 } else {
760 SAFENEWWITHCONSTRUCTOR(pSM,
761 WsmpSparseSolutionManager,
762 WsmpSparseSolutionManager(size, dpivot, block_size, nt));
763 }
764 std::cerr << " using " << nt << " threads " << std::endl;
765 std::cerr << std::endl;
766 #else /* !USE_WSMP */
767 std::cerr << "need --with-wsmp to use Wsmp library"
768 << std::endl;
769 usage(EXIT_FAILURE);
770 #endif /* !USE_WSMP */
771
772 } else if (strcasecmp(solver, "naive") == 0) {
773 std::cerr << "Naive solver";
774 if (dpivot == -1.) {
775 dpivot = 1.E-8;
776 }
777 if (cc) {
778 std::cerr << " with Colamd ordering";
779 if (nt > 1) {
780 #ifdef USE_NAIVE_MULTITHREAD
781 SAFENEWWITHCONSTRUCTOR(pSM,
782 ParNaiveSparsePermSolutionManager,
783 ParNaiveSparsePermSolutionManager(nt, size, dpivot));
784 #else
785 silent_cerr("multithread naive solver support not compiled; "
786 "you can configure --enable-multithread-naive "
787 "on a linux ix86 to get it"
788 << std::endl);
789 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
790 #endif /* USE_NAIVE_MULTITHREAD */
791 } else {
792 SAFENEWWITHCONSTRUCTOR(pSM,
793 NaiveSparsePermSolutionManager<Colamd_ordering>,
794 NaiveSparsePermSolutionManager<Colamd_ordering>(size, dpivot, ms));
795 }
796 } else {
797 if (nt > 1) {
798 #ifdef USE_NAIVE_MULTITHREAD
799 SAFENEWWITHCONSTRUCTOR(pSM,
800 ParNaiveSparseSolutionManager,
801 ParNaiveSparseSolutionManager(nt, size, dpivot));
802 #else
803 silent_cerr("multithread naive solver support not compiled; "
804 "you can configure --enable-multithread-naive "
805 "on a linux ix86 to get it"
806 << std::endl);
807 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
808 #endif /* USE_NAIVE_MULTITHREAD */
809 } else {
810 SAFENEWWITHCONSTRUCTOR(pSM,
811 NaiveSparseSolutionManager,
812 NaiveSparseSolutionManager(size, dpivot, ms));
813 }
814 }
815 std::cerr << " using " << nt << " threads " << std::endl;
816 } else {
817 std::cerr << "unknown solver '" << solver << "'" << std::endl;
818 usage(EXIT_FAILURE);
819 }
820
821 pSM->MatrReset();
822
823 MatrixHandler *pM = pSM->pMatHdl();
824 VectorHandler *pV = pSM->pResHdl();
825 VectorHandler *px = pSM->pSolHdl();
826
827 pM->Reset();
828
829 SetupSystem(random, random_args, singular,
830 matrixfilename, filename, pM, pV);
831
832 #ifdef HAVE_TIMES
833 clock_t ct = 0;
834 struct tms tmsbuf;
835 #endif // HAVE_TIMES
836 try {
837 start = clock();
838 tf = rd_CPU_ts();
839 #ifdef HAVE_TIMES
840 times(&tmsbuf);
841 ct = tmsbuf.tms_utime + tmsbuf.tms_cutime
842 + tmsbuf.tms_stime + tmsbuf.tms_cstime;
843 #endif // HAVE_TIMES
844 if (transpose) {
845 pSM->SolveT();
846 } else {
847 pSM->Solve();
848 }
849 tf = rd_CPU_ts() - tf;
850 end = clock();
851 #ifdef HAVE_TIMES
852 times(&tmsbuf);
853 ct = tmsbuf.tms_utime + tmsbuf.tms_cutime
854 + tmsbuf.tms_stime + tmsbuf.tms_cstime - ct;
855 #endif // HAVE_TIMES
856 } catch (...) {
857 exit(EXIT_FAILURE);
858 }
859 if (output_solution) {
860 for (int i = 1; i <= size; i++) {
861 std::cout << "\tsol[" << i << "] = " << px->operator()(i)
862 << std::endl;
863 }
864 }
865 cpu_time_used = ((double) (end - start));
866 cpu_time_used = cpu_time_used / CLOCKS_PER_SEC;
867 std::cout << "Clock tics to solve: " << tf << std::endl;
868 #ifdef HAVE_TIMES
869 std::cout << "Clock tics to solve: " << ct << std::endl;
870 #endif // HAVE_TIMES
871 std::cout << "Time to solve: " << cpu_time_used << std::endl;
872
873
874 std::cout << "\nSecond solve:\n";
875
876 pSM->MatrReset();
877 pM = pSM->pMatHdl();
878
879 pM->Reset();
880
881 SetupSystem(random, random_args, singular, 0, filename, pM, pV);
882
883 try {
884 start = clock();
885 tf = rd_CPU_ts();
886 #ifdef HAVE_TIMES
887 times(&tmsbuf);
888 ct = tmsbuf.tms_utime + tmsbuf.tms_cutime
889 + tmsbuf.tms_stime + tmsbuf.tms_cstime;
890 #endif // HAVE_TIMES
891 if (transpose) {
892 pSM->SolveT();
893 } else {
894 pSM->Solve();
895 }
896 tf = rd_CPU_ts() - tf;
897 end = clock();
898 #ifdef HAVE_TIMES
899 times(&tmsbuf);
900 ct = tmsbuf.tms_utime + tmsbuf.tms_cutime
901 + tmsbuf.tms_stime + tmsbuf.tms_cstime - ct;
902 #endif // HAVE_TIMES
903 } catch (...) {
904 exit(EXIT_FAILURE);
905 }
906
907 if (output_solution) {
908 for (int i = 1; i <= size; i++) {
909 std::cout << "\tsol[" << i << "] = " << px->operator()(i)
910 << std::endl;
911 }
912 }
913 cpu_time_used = ((double) (end - start));
914 cpu_time_used = cpu_time_used / CLOCKS_PER_SEC;
915 std::cout << "Clock tics to solve: " << tf << std::endl;
916 #ifdef HAVE_TIMES
917 std::cout << "Clock tics to solve: " << ct << std::endl;
918 #endif // HAVE_TIMES
919 std::cout << "Time to solve: " << cpu_time_used << std::endl;
920
921 return 0;
922 }
923
924