1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbwrap/linsol.cc,v 1.82 2017/01/12 14:44:25 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati	<masarati@aero.polimi.it>
9  * Paolo Mantegazza	<mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
30  */
31 
32 #include "mbconfig.h"           /* This goes first in every *.c,*.cc file */
33 
34 #include "myassert.h"
35 #include "ac/sys_sysinfo.h"
36 
37 #include "parser.h"
38 
39 #include "spmapmh.h"
40 #include "ccmh.h"
41 #include "dirccmh.h"
42 #include "harwrap.h"
43 #include "mschwrap.h"
44 #include "y12wrap.h"
45 #include "umfpackwrap.h"
46 #include "kluwrap.h"
47 #include "parsuperluwrap.h"
48 #include "superluwrap.h"
49 #include "lapackwrap.h"
50 #include "taucswrap.h"
51 #include "naivewrap.h"
52 #include "parnaivewrap.h"
53 
54 
55 #include "linsol.h"
56 
57 /* solver data */
58 const LinSol::solver_t solver[] = {
59 	{ "Empty", NULL,
60 		LinSol::EMPTY_SOLVER,
61 		LinSol::SOLVER_FLAGS_NONE,
62 		LinSol::SOLVER_FLAGS_NONE,
63        		-1., -1. },
64 	{ "Harwell", NULL,
65 		LinSol::HARWELL_SOLVER,
66 		LinSol::SOLVER_FLAGS_NONE,
67 		LinSol::SOLVER_FLAGS_NONE,
68 		-1., -1. },
69 	{ "Lapack", NULL,
70 		LinSol::LAPACK_SOLVER,
71 		LinSol::SOLVER_FLAGS_NONE,
72 		LinSol::SOLVER_FLAGS_NONE,
73 		-1., -1. },
74 	{ "Meschach", NULL,
75 		LinSol::MESCHACH_SOLVER,
76 		LinSol::SOLVER_FLAGS_NONE,
77 		LinSol::SOLVER_FLAGS_NONE,
78 		-1., -1. },
79 	{ "Naive", NULL,
80 		LinSol::NAIVE_SOLVER,
81 		LinSol::SOLVER_FLAGS_ALLOWS_REVERSE_CUTHILL_MC_KEE |
82 			LinSol::SOLVER_FLAGS_ALLOWS_COLAMD |
83 			LinSol::SOLVER_FLAGS_ALLOWS_MMDATA |
84 			LinSol::SOLVER_FLAGS_ALLOWS_MDAPLUSAT |
85 			LinSol::SOLVER_FLAGS_ALLOWS_REVERSE_CUTHILL_MC_KEE |
86 			LinSol::SOLVER_FLAGS_ALLOWS_KING |
87 			LinSol::SOLVER_FLAGS_ALLOWS_SLOAN |
88 			LinSol::SOLVER_FLAGS_ALLOWS_NESTED_DISSECTION |
89 			LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS |
90 #ifdef USE_NAIVE_MULTITHREAD
91 			LinSol::SOLVER_FLAGS_ALLOWS_MT_FCT |
92 #endif /* USE_NAIVE_MULTITHREAD */
93 			0,
94 		LinSol::SOLVER_FLAGS_NONE,
95 		1.e-8, -1. },
96 	{ "SuperLU", NULL,
97 		LinSol::SUPERLU_SOLVER,
98 		LinSol::SOLVER_FLAGS_ALLOWS_COLAMD |
99 			LinSol::SOLVER_FLAGS_ALLOWS_MMDATA |
100 			LinSol::SOLVER_FLAGS_ALLOWS_MAP |
101 			LinSol::SOLVER_FLAGS_ALLOWS_CC |
102 			LinSol::SOLVER_FLAGS_ALLOWS_DIR |
103 #ifdef USE_SUPERLU_MT
104 			LinSol::SOLVER_FLAGS_ALLOWS_MT_FCT |
105 #endif /* USE_SUPERLU_MT */
106 			0,
107 		LinSol::SOLVER_FLAGS_ALLOWS_MAP|LinSol::SOLVER_FLAGS_ALLOWS_COLAMD,
108 		1., -1. },
109 	{ "Taucs", NULL,
110 		LinSol::TAUCS_SOLVER,
111 		LinSol::SOLVER_FLAGS_ALLOWS_MAP|LinSol::SOLVER_FLAGS_ALLOWS_CC|LinSol::SOLVER_FLAGS_ALLOWS_DIR,
112 		LinSol::SOLVER_FLAGS_ALLOWS_MAP,
113        		-1., -1. },
114 	{ "Umfpack", "umfpack3",
115 		LinSol::UMFPACK_SOLVER,
116 		LinSol::SOLVER_FLAGS_ALLOWS_MAP|LinSol::SOLVER_FLAGS_ALLOWS_CC|LinSol::SOLVER_FLAGS_ALLOWS_DIR|LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS,
117 		LinSol::SOLVER_FLAGS_ALLOWS_MAP|LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS,
118 		.1,
119 #ifdef UMFPACK_DROPTOL
120 		0.
121 #else // ! UMFPACK_DROPTOL
122 		-1.
123 #endif // ! UMFPACK_DROPTOL
124 		},
125 	{ "KLU", NULL,
126 		LinSol::KLU_SOLVER,
127 		LinSol::SOLVER_FLAGS_ALLOWS_MAP|LinSol::SOLVER_FLAGS_ALLOWS_CC|LinSol::SOLVER_FLAGS_ALLOWS_DIR|LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS,
128 		LinSol::SOLVER_FLAGS_ALLOWS_MAP|LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS,
129 		.1 },
130 	{ "Y12", NULL,
131 		LinSol::Y12_SOLVER,
132 		LinSol::SOLVER_FLAGS_ALLOWS_MAP|LinSol::SOLVER_FLAGS_ALLOWS_CC|LinSol::SOLVER_FLAGS_ALLOWS_DIR|LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS,
133 		LinSol::SOLVER_FLAGS_ALLOWS_MAP|LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS,
134 		-1., -1. },
135 	{ NULL, NULL,
136 		LinSol::EMPTY_SOLVER,
137 		LinSol::SOLVER_FLAGS_NONE,
138 		LinSol::SOLVER_FLAGS_NONE,
139 		-1., -1. }
140 };
141 
142 /*
143  * Default solver
144  */
145 LinSol::SolverType LinSol::defaultSolver =
146 #if defined(USE_UMFPACK)
147 	LinSol::UMFPACK_SOLVER
148 #elif /* !USE_UMFPACK */ defined(USE_KLU)
149 	LinSol::KLU_SOLVER
150 #else /* !USE_KLU */
151 	// Naive always present
152 	LinSol::NAIVE_SOLVER
153 #endif
154 	;
155 
LinSol(void)156 LinSol::LinSol(void)
157 : currSolver(LinSol::defaultSolver),
158 solverFlags(0),
159 nThreads(1),
160 iWorkSpaceSize(0),
161 blockSize(0),
162 dPivotFactor(-1.),
163 dDropTolerance(0.),
164 iMaxIter(0) // Restore the original behavior by default
165 {
166 	NO_OP;
167 }
168 
~LinSol(void)169 LinSol::~LinSol(void)
170 {
171 	NO_OP;
172 }
173 
174 LinSol::SolverType
GetSolver(void) const175 LinSol::GetSolver(void) const
176 {
177 	return currSolver;
178 }
179 
180 bool
SetSolver(LinSol::SolverType t,unsigned f)181 LinSol::SetSolver(LinSol::SolverType t, unsigned f)
182 {
183 	/* check flags */
184 	if ((::solver[t].s_flags & f) != f) {
185 		return false;
186 	}
187 
188 	solverFlags = f;
189 
190 	switch (t) {
191 	case LinSol::UMFPACK_SOLVER:
192 #ifdef USE_UMFPACK
193 		currSolver = t;
194 		return true;
195 #endif /* USE_UMFPACK */
196 
197 	case LinSol::KLU_SOLVER:
198 #ifdef USE_KLU
199 		currSolver = t;
200 		return true;
201 #endif /* USE_UMFPACK */
202 
203 	case LinSol::SUPERLU_SOLVER:
204 #ifdef USE_SUPERLU
205 		currSolver = t;
206 		return true;
207 #endif /* USE_SUPERLU */
208 
209 	case LinSol::LAPACK_SOLVER:
210 #ifdef USE_LAPACK
211 		currSolver = t;
212 		return true;
213 #endif /* USE_LAPACK */
214 
215 	case LinSol::TAUCS_SOLVER:
216 #ifdef USE_TAUCS
217 		currSolver = t;
218 		return true;
219 #endif /* USE_TAUCS */
220 
221 	case LinSol::Y12_SOLVER:
222 #ifdef USE_Y12
223 		currSolver = t;
224 		return true;
225 #endif /* USE_Y12 */
226 
227 	case LinSol::HARWELL_SOLVER:
228 #ifdef USE_HARWELL
229 		currSolver = t;
230 		return true;
231 #endif /* USE_HARWELL */
232 
233 	case LinSol::MESCHACH_SOLVER:
234 #ifdef USE_MESCHACH
235 		currSolver = t;
236 		return true;
237 #endif /* USE_MESCHACH */
238 
239 		/* else */
240 		silent_cerr(::solver[t].s_name << " unavailable" << std::endl);
241 		return false;
242 
243 	case LinSol::NAIVE_SOLVER:
244 		currSolver = t;
245 		return true;
246 
247 	case LinSol::EMPTY_SOLVER:
248 		currSolver = t;
249 		return true;
250 
251 	default:
252 		return false;
253 	}
254 }
255 
256 unsigned
GetSolverFlags(void) const257 LinSol::GetSolverFlags(void) const
258 {
259 	return solverFlags;
260 }
261 
262 unsigned
GetSolverFlags(SolverType t) const263 LinSol::GetSolverFlags(SolverType t) const
264 {
265 	return ::solver[t].s_flags;
266 }
267 
268 const char *const
GetSolverName(void) const269 LinSol::GetSolverName(void) const
270 {
271 	return ::solver[currSolver].s_name;
272 }
273 
274 const char *const
GetSolverName(SolverType t) const275 LinSol::GetSolverName(SolverType t) const
276 {
277 	return ::solver[t].s_name;
278 }
279 
280 bool
SetSolverFlags(unsigned f)281 LinSol::SetSolverFlags(unsigned f)
282 {
283 	if ((::solver[currSolver].s_flags & f) == f) {
284 		solverFlags = f;
285 		return true;
286 	}
287 
288 	return false;
289 }
290 
291 bool
AddSolverFlags(unsigned f)292 LinSol::AddSolverFlags(unsigned f)
293 {
294 	if ((::solver[currSolver].s_flags & f) == f) {
295 		solverFlags |= f;
296 		return true;
297 	}
298 
299 	return false;
300 }
301 
302 bool
MaskSolverFlags(unsigned f)303 LinSol::MaskSolverFlags(unsigned f)
304 {
305 	if ((::solver[currSolver].s_flags & f) == f) {
306 		solverFlags &= ~f;
307 		return true;
308 	}
309 
310 	return false;
311 }
312 
313 bool
SetNumThreads(unsigned nt)314 LinSol::SetNumThreads(unsigned nt)
315 {
316 	if (GetSolverFlags(currSolver) & LinSol::SOLVER_FLAGS_ALLOWS_MT_FCT) {
317 		if (nt == 0) {
318 			solverFlags &= ~LinSol::SOLVER_FLAGS_ALLOWS_MT_FCT;
319 
320 		} else {
321 			solverFlags |= LinSol::SOLVER_FLAGS_ALLOWS_MT_FCT;
322 		}
323 
324 		nThreads = nt;
325 
326 		return true;
327 	}
328 
329 	return false;
330 }
331 
332 unsigned
GetNumThreads(void) const333 LinSol::GetNumThreads(void) const
334 {
335 	return nThreads;
336 }
337 
338 integer
iGetWorkSpaceSize(void) const339 LinSol::iGetWorkSpaceSize(void) const
340 {
341 	return iWorkSpaceSize;
342 }
343 
344 const doublereal&
dGetPivotFactor(void) const345 LinSol::dGetPivotFactor(void) const
346 {
347 	return dPivotFactor;
348 }
349 
350 const doublereal&
dGetDropTolerance(void) const351 LinSol::dGetDropTolerance(void) const
352 {
353 	return dDropTolerance;
354 }
355 
356 bool
SetWorkSpaceSize(integer i)357 LinSol::SetWorkSpaceSize(integer i)
358 {
359 	switch (currSolver) {
360 	case LinSol::Y12_SOLVER:
361 		iWorkSpaceSize = i;
362 		break;
363 
364 	default:
365 		return false;
366 	}
367 
368 	return true;
369 }
370 
371 bool
SetPivotFactor(const doublereal & d)372 LinSol::SetPivotFactor(const doublereal& d)
373 {
374 	if (::solver[currSolver].s_pivot_factor == -1.) {
375 		return false;
376 	}
377 
378 	dPivotFactor = d;
379 
380 	return true;
381 }
382 
383 bool
SetDropTolerance(const doublereal & d)384 LinSol::SetDropTolerance(const doublereal& d)
385 {
386 	if (::solver[currSolver].s_drop_tolerance == -1.) {
387 		return false;
388 	}
389 
390 	dDropTolerance = d;
391 
392 	return true;
393 }
394 
395 unsigned
GetBlockSize(void) const396 LinSol::GetBlockSize(void) const
397 {
398 	return blockSize;
399 }
400 
401 bool
SetBlockSize(unsigned bs)402 LinSol::SetBlockSize(unsigned bs)
403 {
404 	switch (currSolver) {
405 	case LinSol::UMFPACK_SOLVER:
406 		blockSize = bs;
407 		break;
408 
409 	default:
410 		return false;
411 	}
412 
413 	return true;
414 }
415 
416 bool
SetScale(const SolutionManager::ScaleOpt & s)417 LinSol::SetScale(const SolutionManager::ScaleOpt& s)
418 {
419 	switch (currSolver) {
420 	case LinSol::NAIVE_SOLVER:
421 	case LinSol::UMFPACK_SOLVER:
422 	case LinSol::KLU_SOLVER:
423 		scale = s;
424 		break;
425 
426 	default:
427 		return false;
428 	}
429 
430 	return true;
431 }
432 
433 integer
GetMaxIterations(void) const434 LinSol::GetMaxIterations(void) const
435 {
436 	return iMaxIter;
437 }
438 
439 bool
SetMaxIterations(integer iMaxIterations)440 LinSol::SetMaxIterations(integer iMaxIterations)
441 {
442 	switch (currSolver) {
443 	case LinSol::UMFPACK_SOLVER:
444 		iMaxIter = iMaxIterations;
445 		break;
446 
447 	default:
448 		return false;
449 	}
450 
451 	return true;
452 }
453 
454 SolutionManager *const
GetSolutionManager(integer iNLD,integer iLWS) const455 LinSol::GetSolutionManager(integer iNLD, integer iLWS) const
456 {
457 	SolutionManager *pCurrSM = NULL;
458 	const unsigned type = (solverFlags & LinSol::SOLVER_FLAGS_TYPE_MASK);
459 	const unsigned perm = (solverFlags & LinSol::SOLVER_FLAGS_PERM_MASK);
460 	const bool mt = (solverFlags & LinSol::SOLVER_FLAGS_ALLOWS_MT_FCT);
461 
462 	/* silence warning */
463 	if (mt) {
464 		NO_OP;
465 	}
466 
467 	ASSERT((::solver[currSolver].s_flags & solverFlags) == solverFlags);
468 
469 	if (iLWS == 0) {
470 		iLWS = iWorkSpaceSize;
471 	}
472 
473    	switch (currSolver) {
474      	case LinSol::Y12_SOLVER:
475 #ifdef USE_Y12
476 		switch (type) {
477 		case LinSol::SOLVER_FLAGS_ALLOWS_DIR: {
478 			typedef Y12SparseCCSolutionManager<DirCColMatrixHandler<1> > CCSM;
479 	      		SAFENEWWITHCONSTRUCTOR(pCurrSM, CCSM,
480 					CCSM(iNLD, iLWS, dPivotFactor));
481 			break;
482 		}
483 
484 		case LinSol::SOLVER_FLAGS_ALLOWS_CC: {
485 			typedef Y12SparseCCSolutionManager<CColMatrixHandler<1> > CCSM;
486 	      		SAFENEWWITHCONSTRUCTOR(pCurrSM, CCSM,
487 					CCSM(iNLD, iLWS, dPivotFactor));
488 			break;
489 		}
490 
491 		default:
492       			SAFENEWWITHCONSTRUCTOR(pCurrSM,
493 				Y12SparseSolutionManager,
494 				Y12SparseSolutionManager(iNLD, iLWS,
495 					dPivotFactor));
496 			break;
497 		}
498       		break;
499 #else /* !USE_Y12 */
500       		silent_cerr("Configure with --with-y12 "
501 			"to enable Y12 solver" << std::endl);
502       		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
503 #endif /* !USE_Y12 */
504 
505      	case LinSol::SUPERLU_SOLVER:
506 #ifdef USE_SUPERLU
507 #ifdef USE_SUPERLU_MT
508 		switch (type) {
509 		case LinSol::SOLVER_FLAGS_ALLOWS_DIR: {
510 			typedef ParSuperLUSparseCCSolutionManager<DirCColMatrixHandler<0> > CCSM;
511 	      		SAFENEWWITHCONSTRUCTOR(pCurrSM, CCSM,
512 					CCSM(nThreads, iNLD, dPivotFactor));
513 			break;
514 		}
515 
516 		case LinSol::SOLVER_FLAGS_ALLOWS_CC: {
517 			typedef ParSuperLUSparseCCSolutionManager<CColMatrixHandler<0> > CCSM;
518 	      		SAFENEWWITHCONSTRUCTOR(pCurrSM, CCSM,
519 					CCSM(nThreads, iNLD, dPivotFactor));
520 			break;
521 		}
522 
523 		default:
524 			SAFENEWWITHCONSTRUCTOR(pCurrSM,
525 				ParSuperLUSparseSolutionManager,
526 				ParSuperLUSparseSolutionManager(nThreads, iNLD,
527 					dPivotFactor));
528 			break;
529 		}
530 		break;
531 
532 		if (nThreads == 1) {
533 			silent_cerr("warning, using multithread SuperLU with only one thread; "
534 				<< std::endl);
535 		}
536 #else /* !USE_SUPERLU_MT */
537 		if (nThreads > 1 && mt) {
538 			silent_cerr("multithread SuperLU solver support not compiled; "
539 				<< std::endl);
540 			throw ErrGeneric(MBDYN_EXCEPT_ARGS);
541 		}
542 		switch (type) {
543 		case LinSol::SOLVER_FLAGS_ALLOWS_DIR: {
544 			typedef SuperLUSparseCCSolutionManager<DirCColMatrixHandler<0> > CCSM;
545 	      		SAFENEWWITHCONSTRUCTOR(pCurrSM, CCSM,
546 					CCSM(iNLD, dPivotFactor));
547 			break;
548 		}
549 
550 		case LinSol::SOLVER_FLAGS_ALLOWS_CC: {
551 			typedef SuperLUSparseCCSolutionManager<CColMatrixHandler<0> > CCSM;
552 	      		SAFENEWWITHCONSTRUCTOR(pCurrSM, CCSM,
553 					CCSM(iNLD, dPivotFactor));
554 			break;
555 		}
556 
557 		default:
558 			SAFENEWWITHCONSTRUCTOR(pCurrSM,
559 				SuperLUSparseSolutionManager,
560 				SuperLUSparseSolutionManager(iNLD,
561 					dPivotFactor));
562 			break;
563 		}
564 #endif /* !USE_SUPERLU_MT */
565 		break;
566 #else /* !USE_SUPERLU */
567       		silent_cerr("Configure with --with-superlu "
568 			"to enable superlu solver" << std::endl);
569       		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
570 #endif /* !USE_SUPERLU */
571 
572 	case LinSol::MESCHACH_SOLVER:
573 #ifdef USE_MESCHACH
574 		SAFENEWWITHCONSTRUCTOR(pCurrSM,
575 			MeschachSparseSolutionManager,
576 			MeschachSparseSolutionManager(iNLD, iLWS,
577 				dPivotFactor));
578 		break;
579 #else /* !USE_MESCHACH */
580 		silent_cerr("Configure with --with-meschach "
581 			"to enable Meschach solver" << std::endl);
582       		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
583 #endif /* !USE_MESCHACH */
584 
585 	case LinSol::LAPACK_SOLVER:
586 #ifdef USE_LAPACK
587 		SAFENEWWITHCONSTRUCTOR(pCurrSM,
588 			LapackSolutionManager,
589 			LapackSolutionManager(iNLD, dPivotFactor));
590 		break;
591 #else /* !USE_LAPACK */
592 		silent_cerr("Configure with --with-lapack "
593 			"to enable Lapack dense solver" << std::endl);
594       		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
595 #endif /* !USE_LAPACK */
596 
597      	case LinSol::TAUCS_SOLVER:
598 #ifdef USE_TAUCS
599 		switch (type) {
600 		case LinSol::SOLVER_FLAGS_ALLOWS_DIR: {
601 			typedef TaucsSparseCCSolutionManager<DirCColMatrixHandler<0> > CCSM;
602 	      		SAFENEWWITHCONSTRUCTOR(pCurrSM, CCSM, CCSM(iNLD));
603 			break;
604 		}
605 
606 		case LinSol::SOLVER_FLAGS_ALLOWS_CC: {
607 			typedef TaucsSparseCCSolutionManager<CColMatrixHandler<0> > CCSM;
608 	      		SAFENEWWITHCONSTRUCTOR(pCurrSM, CCSM, CCSM(iNLD));
609 			break;
610 		}
611 
612 		default:
613       			SAFENEWWITHCONSTRUCTOR(pCurrSM,
614 				TaucsSparseSolutionManager,
615 				TaucsSparseSolutionManager(iNLD));
616 			break;
617 		}
618       		break;
619 #else /* !USE_TAUCS */
620 		silent_cerr("Configure with --with-taucs "
621 			"to enable Taucs sparse solver" << std::endl);
622       		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
623 #endif /* !USE_TAUCS */
624 
625  	case LinSol::HARWELL_SOLVER:
626 #ifdef USE_HARWELL
627 		SAFENEWWITHCONSTRUCTOR(pCurrSM,
628 			HarwellSparseSolutionManager,
629 			HarwellSparseSolutionManager(iNLD, iLWS,
630 				dPivotFactor));
631       		break;
632 #else /* !USE_HARWELL */
633       		silent_cerr("Configure with --with-harwell "
634 			"to enable Harwell solver" << std::endl);
635 		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
636 #endif /* !USE_HARWELL */
637 
638 	case LinSol::UMFPACK_SOLVER:
639 #ifdef USE_UMFPACK
640 		{
641 		switch (type) {
642 		case LinSol::SOLVER_FLAGS_ALLOWS_DIR: {
643 			typedef UmfpackSparseCCSolutionManager<DirCColMatrixHandler<0> > CCSM;
644 	      		SAFENEWWITHCONSTRUCTOR(pCurrSM, CCSM,
645 					CCSM(iNLD, dPivotFactor, dDropTolerance, blockSize, scale, iMaxIter));
646 			break;
647 		}
648 
649 		case LinSol::SOLVER_FLAGS_ALLOWS_CC: {
650 			typedef UmfpackSparseCCSolutionManager<CColMatrixHandler<0> > CCSM;
651 	      		SAFENEWWITHCONSTRUCTOR(pCurrSM, CCSM,
652 					CCSM(iNLD, dPivotFactor, dDropTolerance, blockSize, scale, iMaxIter));
653 			break;
654 		}
655 
656 		default:
657 			SAFENEWWITHCONSTRUCTOR(pCurrSM,
658 				UmfpackSparseSolutionManager,
659 				UmfpackSparseSolutionManager(iNLD,
660 					dPivotFactor, dDropTolerance, blockSize, scale, iMaxIter));
661 			break;
662 		}
663 		} break;
664 #else /* !USE_UMFPACK */
665       		silent_cerr("Configure with --with-umfpack "
666 			"to enable Umfpack solver" << std::endl);
667       		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
668 #endif /* !USE_UMFPACK */
669 
670 	case LinSol::KLU_SOLVER:
671 #ifdef USE_KLU
672 		{
673 		switch (type) {
674 		case LinSol::SOLVER_FLAGS_ALLOWS_DIR: {
675 			typedef KLUSparseCCSolutionManager<DirCColMatrixHandler<0> > CCSM;
676 	      		SAFENEWWITHCONSTRUCTOR(pCurrSM, CCSM,
677 					CCSM(iNLD, dPivotFactor, scale));
678 			break;
679 		}
680 
681 		case LinSol::SOLVER_FLAGS_ALLOWS_CC: {
682 			typedef KLUSparseCCSolutionManager<CColMatrixHandler<0> > CCSM;
683 	      		SAFENEWWITHCONSTRUCTOR(pCurrSM, CCSM,
684 					CCSM(iNLD, dPivotFactor, scale));
685 			break;
686 		}
687 
688 		default:
689 			SAFENEWWITHCONSTRUCTOR(pCurrSM,
690 				KLUSparseSolutionManager,
691 				KLUSparseSolutionManager(iNLD,
692 										 dPivotFactor,
693 										 scale));
694 			break;
695 		}
696 
697 		} break;
698 #else /* !USE_KLU */
699       		silent_cerr("Configure with --with-klu "
700 			"to enable KLU solver" << std::endl);
701       		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
702 #endif /* !USE_KLU */
703 
704 	case LinSol::NAIVE_SOLVER:
705 		if (perm == LinSol::SOLVER_FLAGS_ALLOWS_COLAMD) {
706 			if (nThreads == 1) {
707 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
708 					NaiveSparsePermSolutionManager<Colamd_ordering>,
709 					NaiveSparsePermSolutionManager<Colamd_ordering>(iNLD, dPivotFactor, scale));
710 			} else {
711 #ifdef USE_NAIVE_MULTITHREAD
712 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
713 					ParNaiveSparsePermSolutionManager,
714 					ParNaiveSparsePermSolutionManager(nThreads, iNLD, dPivotFactor));
715 #else /* ! USE_NAIVE_MULTITHREAD */
716 				silent_cerr("multithread naive solver support not compiled; "
717 					"you can configure --enable-multithread-naive "
718 					"on a linux ix86 to get it"
719 					<< std::endl);
720 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
721 #endif /* ! USE_NAIVE_MULTITHREAD */
722 			}
723 #ifdef USE_BOOST
724 #ifdef HAVE_BOOST_GRAPH_CUTHILL_MCKEE_ORDERING_HPP
725 		} else if (perm == LinSol::SOLVER_FLAGS_ALLOWS_REVERSE_CUTHILL_MC_KEE) {
726 			if (nThreads == 1) {
727 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
728 					NaiveSparsePermSolutionManager<rcmk_ordering>,
729 					NaiveSparsePermSolutionManager<rcmk_ordering>(iNLD, dPivotFactor, scale));
730 			} else {
731 #ifdef USE_NAIVE_MULTITHREAD
732 #if 0
733 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
734 					ParNaiveSparsePermSolutionManager,
735 					ParNaiveSparsePermSolutionManager(nThreads, iNLD, dPivotFactor));
736 #endif
737 				silent_cerr("multithread naive solver with"
738 					"reverse Cuthill-McKee permutation not"
739 					"available yet. Patches welcome"
740 					<< std::endl);
741 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
742 #else /* ! USE_NAIVE_MULTITHREAD */
743 				silent_cerr("multithread naive solver support not compiled; "
744 					"you can configure --enable-multithread-naive "
745 					"on a linux ix86 to get it"
746 					<< std::endl);
747 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
748 #endif /* ! USE_NAIVE_MULTITHREAD */
749 			}
750 #endif /* HAVE_BOOST_GRAPH_CUTHILL_MCKEE_ORDERING_HPP */
751 
752 #if 0 /* ?!? */
753  		} else if (perm == LinSol::SOLVER_FLAGS_ALLOWS_MMDATA) {
754  			if (nThreads == 1) {
755  				SAFENEWWITHCONSTRUCTOR(pCurrSM,
756  					NaiveSparsePermSolutionManager<amd_ordering>,
757  					NaiveSparsePermSolutionManager<amd_ordering>(iNLD, dPivotFactor, scale));
758  			} else {
759 #ifdef USE_NAIVE_MULTITHREAD
760 #if 0
761 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
762 					ParNaiveSparsePermSolutionManager,
763 					ParNaiveSparsePermSolutionManager(nThreads, iNLD, dPivotFactor));
764 #endif
765  				silent_cerr("multithread naive solver with"
766  					"approximate minimum degree permutation not"
767  					"available yet. Patches welcome"
768  					<< std::endl);
769  				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
770 #else
771  				silent_cerr("multithread naive solver support not compiled; "
772  					"you can configure --enable-multithread-naive "
773  					"on a linux ix86 to get it"
774  					<< std::endl);
775  				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
776 #endif /* USE_NAIVE_MULTITHREAD */
777  			}
778 #endif
779 
780 #ifdef HAVE_BOOST_GRAPH_KING_ORDERING_HPP
781 		} else if (perm == LinSol::SOLVER_FLAGS_ALLOWS_KING) {
782 			if (nThreads == 1) {
783 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
784 					NaiveSparsePermSolutionManager<king_ordering>,
785 					NaiveSparsePermSolutionManager<king_ordering>(iNLD, dPivotFactor, scale));
786 			} else {
787 #ifdef USE_NAIVE_MULTITHREAD
788 #if 0
789 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
790 					ParNaiveSparsePermSolutionManager,
791 					ParNaiveSparsePermSolutionManager(nThreads, iNLD, dPivotFactor));
792 #endif
793 				silent_cerr("multithread naive solver with"
794 					"king permutation not"
795 					"available yet. Patches welcome"
796 					<< std::endl);
797 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
798 #else /* ! USE_NAIVE_MULTITHREAD */
799 				silent_cerr("multithread naive solver support not compiled; "
800 					"you can configure --enable-multithread-naive "
801 					"on a linux ix86 to get it"
802 					<< std::endl);
803 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
804 #endif /* ! USE_NAIVE_MULTITHREAD */
805 			}
806 #endif /* HAVE_BOOST_GRAPH_KING_ORDERING_HPP */
807 
808 #ifdef HAVE_BOOST_GRAPH_SLOAN_ORDERING_HPP
809 		} else if (perm == LinSol::SOLVER_FLAGS_ALLOWS_SLOAN) {
810 			if (nThreads == 1) {
811 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
812 					NaiveSparsePermSolutionManager<sloan_ordering>,
813 					NaiveSparsePermSolutionManager<sloan_ordering>(iNLD, dPivotFactor, scale));
814 			} else {
815 #ifdef USE_NAIVE_MULTITHREAD
816 #if 0
817 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
818 					ParNaiveSparsePermSolutionManager,
819 					ParNaiveSparsePermSolutionManager(nThreads, iNLD, dPivotFactor));
820 #endif
821 				silent_cerr("multithread naive solver with"
822 					"sloan permutation not"
823 					"available yet. Patches welcome"
824 					<< std::endl);
825 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
826 #else /* ! USE_NAIVE_MULTITHREAD */
827 				silent_cerr("multithread naive solver support not compiled; "
828 					"you can configure --enable-multithread-naive "
829 					"on a linux ix86 to get it"
830 					<< std::endl);
831 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
832 #endif /* ! USE_NAIVE_MULTITHREAD */
833 			}
834 #endif /* HAVE_BOOST_GRAPH_SLOAN_ORDERING_HPP */
835 
836 #ifdef HAVE_BOOST_GRAPH_MINIMUM_DEGREE_ORDERING_HPP
837 		} else if (perm == LinSol::SOLVER_FLAGS_ALLOWS_MDAPLUSAT) {
838 			if (nThreads == 1) {
839 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
840 					NaiveSparsePermSolutionManager<md_ordering>,
841 					NaiveSparsePermSolutionManager<md_ordering>(iNLD, dPivotFactor, scale));
842 			} else {
843 #ifdef USE_NAIVE_MULTITHREAD
844 #if 0
845 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
846 					ParNaiveSparsePermSolutionManager,
847 					ParNaiveSparsePermSolutionManager(nThreads, iNLD, dPivotFactor));
848 #endif
849 				silent_cerr("multithread naive solver with"
850 					"md permutation not"
851 					"available yet. Patches welcome"
852 					<< std::endl);
853 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
854 #else /* ! USE_NAIVE_MULTITHREAD */
855 				silent_cerr("multithread naive solver support not compiled; "
856 					"you can configure --enable-multithread-naive "
857 					"on a linux ix86 to get it"
858 					<< std::endl);
859 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
860 #endif /* ! USE_NAIVE_MULTITHREAD */
861 			}
862 #endif /* HAVE_BOOST_GRAPH_MINIMUM_DEGREE_ORDERING_HPP */
863 #endif /* USE_BOOST */
864 
865 		} else if (perm == LinSol::SOLVER_FLAGS_ALLOWS_NESTED_DISSECTION) {
866 #ifdef USE_METIS
867 			if (nThreads == 1) {
868 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
869 					NaiveSparsePermSolutionManager<metis_ordering>,
870 					NaiveSparsePermSolutionManager<metis_ordering>(iNLD, dPivotFactor, scale));
871 			} else {
872 #ifdef USE_NAIVE_MULTITHREAD
873 #if 0
874 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
875 					ParNaiveSparsePermSolutionManager,
876 					ParNaiveSparsePermSolutionManager(nThreads, iNLD, dPivotFactor));
877 #endif
878 				silent_cerr("multithread naive solver with"
879 					"nested dissection permutation not"
880 					"available yet. Patches welcome"
881 					<< std::endl);
882 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
883 #else /* ! USE_NAIVE_MULTITHREAD */
884 				silent_cerr("multithread naive solver support not compiled; "
885 					"you can configure --enable-multithread-naive "
886 					"on a linux ix86 to get it"
887 					<< std::endl);
888 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
889 #endif /* ! USE_NAIVE_MULTITHREAD */
890 			}
891 #else /* ! USE_METIS */
892 			silent_cerr("you should not get here("<< __FILE__ << ":" <<
893 				__LINE__ << ")" << std::endl);
894 			throw ErrGeneric(MBDYN_EXCEPT_ARGS);
895 #endif /* ! USE_METIS */
896 
897 		} else {
898 			if (nThreads == 1) {
899 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
900 					NaiveSparseSolutionManager,
901 					NaiveSparseSolutionManager(iNLD, dPivotFactor, scale));
902 			} else {
903 #ifdef USE_NAIVE_MULTITHREAD
904 				SAFENEWWITHCONSTRUCTOR(pCurrSM,
905 					ParNaiveSparseSolutionManager,
906 					ParNaiveSparseSolutionManager(nThreads, iNLD, dPivotFactor));
907 #else
908 				silent_cerr("multithread naive solver support not compiled; "
909 					"you can configure --enable-multithread-naive "
910 					"on a linux ix86 to get it"
911 					<< std::endl);
912 				throw ErrGeneric(MBDYN_EXCEPT_ARGS);
913 #endif /* USE_NAIVE_MULTITHREAD */
914 			}
915 		}
916 		break;
917 
918 	case LinSol::EMPTY_SOLVER:
919 		break;
920 
921    	default:
922 		ASSERT(0);
923 		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
924 
925 	}
926 
927 	return pCurrSM;
928 }
929 
930