1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/readlinsol.cc,v 1.36 2017/01/12 14:46:10 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 extern "C" { int get_nprocs_conf(void); int get_nprocs(void); }
35 #include "dataman.h"
36 #include "readlinsol.h"
37 
38 extern "C" {int get_nprocs(void);}
39 
40 void
ReadLinSol(LinSol & cs,HighParser & HP,bool bAllowEmpty)41 ReadLinSol(LinSol& cs, HighParser &HP, bool bAllowEmpty)
42 {
43    	/* parole chiave */
44    	const char* sKeyWords[] = {
45 		::solver[LinSol::EMPTY_SOLVER].s_name,
46 		::solver[LinSol::HARWELL_SOLVER].s_name,
47 		::solver[LinSol::LAPACK_SOLVER].s_name,
48 		::solver[LinSol::MESCHACH_SOLVER].s_name,
49 		::solver[LinSol::NAIVE_SOLVER].s_name,
50 		::solver[LinSol::SUPERLU_SOLVER].s_name,
51 		::solver[LinSol::TAUCS_SOLVER].s_name,
52 		::solver[LinSol::UMFPACK_SOLVER].s_name,
53 		::solver[LinSol::UMFPACK_SOLVER].s_alias,
54 		::solver[LinSol::KLU_SOLVER].s_name,
55 		::solver[LinSol::Y12_SOLVER].s_name,
56 		NULL
57 	};
58 
59 	enum KeyWords {
60 		EMPTY,
61 		HARWELL,
62 		LAPACK,
63 		MESCHACH,
64 		NAIVE,
65 		SUPERLU,
66 		TAUCS,
67 		UMFPACK,
68 		UMFPACK3,
69 		KLU,
70 		Y12,
71 
72 		LASTKEYWORD
73 	};
74 
75    	/* tabella delle parole chiave */
76    	KeyTable K(HP, sKeyWords);
77 
78 	bool bGotIt = false;
79 	switch (KeyWords(HP.GetWord())) {
80 	case EMPTY:
81 		if (!bAllowEmpty) {
82 			silent_cerr("empty solver not allowed at line "
83 				<< HP.GetLineData() << std::endl);
84 			throw ErrGeneric(MBDYN_EXCEPT_ARGS);
85 		}
86 
87 		cs.SetSolver(LinSol::EMPTY_SOLVER);
88 		DEBUGLCOUT(MYDEBUG_INPUT,
89 				"No LU solver" << std::endl);
90 		bGotIt = true;
91 		break;
92 
93 	case HARWELL:
94 		cs.SetSolver(LinSol::HARWELL_SOLVER);
95 		DEBUGLCOUT(MYDEBUG_INPUT,
96 				"Using harwell sparse LU solver" << std::endl);
97 #ifdef USE_HARWELL
98 		bGotIt = true;
99 #endif /* USE_HARWELL */
100 		break;
101 
102 	case LAPACK:
103 		cs.SetSolver(LinSol::LAPACK_SOLVER);
104 		DEBUGLCOUT(MYDEBUG_INPUT,
105 				"Using Lapack dense LU solver" << std::endl);
106 #ifdef USE_LAPACK
107 		bGotIt = true;
108 #endif /* USE_LAPACK */
109 		break;
110 
111 	case MESCHACH:
112 		cs.SetSolver(LinSol::MESCHACH_SOLVER);
113 		DEBUGLCOUT(MYDEBUG_INPUT,
114 				"Using meschach sparse LU solver"
115 				<< std::endl);
116 #ifdef USE_MESCHACH
117 		bGotIt = true;
118 #endif /* USE_MESCHACH */
119 		break;
120 
121 	case NAIVE:
122 		cs.SetSolver(LinSol::NAIVE_SOLVER);
123 		bGotIt = true;
124 		break;
125 
126 	case SUPERLU:
127 		/*
128 		 * FIXME: use CC as default???
129 		 */
130 		cs.SetSolver(LinSol::SUPERLU_SOLVER);
131 		DEBUGLCOUT(MYDEBUG_INPUT,
132 				"Using SuperLU sparse LU solver" << std::endl);
133 #ifdef USE_SUPERLU
134 		bGotIt = true;
135 #endif /* USE_SUPERLU */
136 		break;
137 
138 	case TAUCS:
139 		/*
140 		 * FIXME: use CC as default???
141 		 */
142 		cs.SetSolver(LinSol::TAUCS_SOLVER);
143 		DEBUGLCOUT(MYDEBUG_INPUT,
144 				"Using Taucs sparse solver" << std::endl);
145 #ifdef USE_TAUCS
146 		bGotIt = true;
147 #endif /* USE_TAUCS */
148 		break;
149 
150 	case UMFPACK3:
151 		pedantic_cerr("\"umfpack3\" is deprecated; "
152 				"use \"umfpack\" instead" << std::endl);
153 	case UMFPACK:
154 		/*
155 		 * FIXME: use CC as default???
156 		 */
157 		cs.SetSolver(LinSol::UMFPACK_SOLVER);
158 		DEBUGLCOUT(MYDEBUG_INPUT,
159 				"Using umfpack sparse LU solver" << std::endl);
160 #ifdef USE_UMFPACK
161 		bGotIt = true;
162 #endif /* USE_UMFPACK */
163 		break;
164 
165 	case KLU:
166 		/*
167 		 * FIXME: use CC as default???
168 		 */
169 		cs.SetSolver(LinSol::KLU_SOLVER);
170 		DEBUGLCOUT(MYDEBUG_INPUT,
171 				"Using KLU sparse LU solver" << std::endl);
172 #ifdef USE_KLU
173 		bGotIt = true;
174 #endif /* USE_KLU */
175 		break;
176 
177 	case Y12:
178 		/*
179 		 * FIXME: use CC as default???
180 		 */
181 		cs.SetSolver(LinSol::Y12_SOLVER);
182 		DEBUGLCOUT(MYDEBUG_INPUT,
183 				"Using y12 sparse LU solver" << std::endl);
184 #ifdef USE_Y12
185 		bGotIt = true;
186 #endif /* USE_Y12 */
187 		break;
188 
189 	default:
190 		silent_cerr("unknown solver" << std::endl);
191 		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
192 	}
193 
194 	const LinSol::solver_t	currSolver = ::solver[cs.GetSolver()];
195 
196 	if (!bGotIt) {
197 		silent_cerr(currSolver.s_name << " solver "
198 			"not available; requested at line " << HP.GetLineData()
199 			<< std::endl);
200 		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
201 	}
202 
203 	cs.SetSolverFlags(currSolver.s_default_flags);
204 
205 	/* map? */
206 	if (HP.IsKeyWord("map")) {
207 		if (currSolver.s_flags & LinSol::SOLVER_FLAGS_ALLOWS_MAP) {
208 			cs.MaskSolverFlags(LinSol::SOLVER_FLAGS_TYPE_MASK);
209 			cs.AddSolverFlags(LinSol::SOLVER_FLAGS_ALLOWS_MAP);
210 			pedantic_cout("using map matrix handling for "
211 					<< currSolver.s_name
212 					<< " solver" << std::endl);
213 		} else {
214 			pedantic_cerr("map is meaningless for "
215 					<< currSolver.s_name
216 					<< " solver" << std::endl);
217 		}
218 
219 	/* CC? */
220 	} else if (HP.IsKeyWord("column" "compressed") || HP.IsKeyWord("cc")) {
221 		if (currSolver.s_flags & LinSol::SOLVER_FLAGS_ALLOWS_CC) {
222 			cs.MaskSolverFlags(LinSol::SOLVER_FLAGS_TYPE_MASK);
223 			cs.AddSolverFlags(LinSol::SOLVER_FLAGS_ALLOWS_CC);
224 			pedantic_cout("using column compressed matrix handling for "
225 					<< currSolver.s_name
226 					<< " solver" << std::endl);
227 
228 		} else {
229 			pedantic_cerr("column compressed is meaningless for "
230 					<< currSolver.s_name
231 					<< " solver" << std::endl);
232 		}
233 
234 	/* direct? */
235 	} else if (HP.IsKeyWord("direct" "access") || HP.IsKeyWord("dir")) {
236 		if (currSolver.s_flags & LinSol::SOLVER_FLAGS_ALLOWS_DIR) {
237 			cs.MaskSolverFlags(LinSol::SOLVER_FLAGS_TYPE_MASK);
238 			cs.AddSolverFlags(LinSol::SOLVER_FLAGS_ALLOWS_DIR);
239 			pedantic_cout("using direct access matrix handling for "
240 					<< currSolver.s_name
241 					<< " solver" << std::endl);
242 		} else {
243 			pedantic_cerr("direct is meaningless for "
244 					<< currSolver.s_name
245 					<< " solver" << std::endl);
246 		}
247 	}
248 
249 	/* colamd? */
250 	if (HP.IsKeyWord("colamd")) {
251 		if (currSolver.s_flags & LinSol::SOLVER_FLAGS_ALLOWS_COLAMD) {
252 			cs.AddSolverFlags(LinSol::SOLVER_FLAGS_ALLOWS_COLAMD);
253 			pedantic_cout("using colamd symmetric preordering for "
254 					<< currSolver.s_name
255 					<< " solver" << std::endl);
256 		} else {
257 			pedantic_cerr("colamd preordering is meaningless for "
258 					<< currSolver.s_name
259 					<< " solver" << std::endl);
260 		}
261 	/* amd ata? */
262 	} else if (HP.IsKeyWord("mmdata")) {
263 		silent_cerr("approximate minimum degree solver support is still TODO"
264 			"task: detect (or import) the MD library;"
265 			"uncomment the relevant bits in naivewrap;"
266 			"remove this check (readlinsol.cc)."
267 			"Patches welcome"
268 			<< std::endl);
269 		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
270 		if (currSolver.s_flags & LinSol::SOLVER_FLAGS_ALLOWS_MMDATA) {
271 			cs.AddSolverFlags(LinSol::SOLVER_FLAGS_ALLOWS_MMDATA);
272 			pedantic_cout("using mmd symmetric preordering for "
273 					<< currSolver.s_name
274 					<< " solver" << std::endl);
275 		} else {
276 			pedantic_cerr("mmdata preordering is meaningless for "
277 					<< currSolver.s_name
278 					<< " solver" << std::endl);
279 		}
280 	/* minimum degree ?*/
281 	} else if (HP.IsKeyWord("minimum" "degree")) {
282 		if (currSolver.s_flags & LinSol::SOLVER_FLAGS_ALLOWS_MDAPLUSAT) {
283 			cs.AddSolverFlags(LinSol::SOLVER_FLAGS_ALLOWS_MDAPLUSAT);
284 			pedantic_cout("using minimum degree symmetric preordering of A+A^T for "
285 					<< currSolver.s_name
286 					<< " solver" << std::endl);
287 		} else {
288 			pedantic_cerr("md preordering is meaningless for "
289 					<< currSolver.s_name
290 					<< " solver" << std::endl);
291 		}
292 	/* Reverse Kuthill McKee? */
293 	} else if (HP.IsKeyWord("rcmk")) {
294 		if (currSolver.s_flags & LinSol::SOLVER_FLAGS_ALLOWS_REVERSE_CUTHILL_MC_KEE) {
295 			cs.AddSolverFlags(LinSol::SOLVER_FLAGS_ALLOWS_REVERSE_CUTHILL_MC_KEE);
296 			pedantic_cout("using rcmk symmetric preordering for "
297 					<< currSolver.s_name
298 					<< " solver" << std::endl);
299 
300 		} else {
301 			pedantic_cerr("rcmk preordering is meaningless for "
302 					<< currSolver.s_name
303 					<< " solver" << std::endl);
304 		}
305 	/* king ?*/
306 	} else if (HP.IsKeyWord("king")) {
307 		if (currSolver.s_flags & LinSol::SOLVER_FLAGS_ALLOWS_KING) {
308 			cs.AddSolverFlags(LinSol::SOLVER_FLAGS_ALLOWS_KING);
309 			pedantic_cout("using king symmetric preordering for "
310 					<< currSolver.s_name
311 					<< " solver" << std::endl);
312 
313 		} else {
314 			pedantic_cerr("king preordering is meaningless for "
315 					<< currSolver.s_name
316 					<< " solver" << std::endl);
317 		}
318 	/* sloan ? */
319 	} else if (HP.IsKeyWord("sloan")) {
320 		if (currSolver.s_flags & LinSol::SOLVER_FLAGS_ALLOWS_KING) {
321 			cs.AddSolverFlags(LinSol::SOLVER_FLAGS_ALLOWS_KING);
322 			pedantic_cout("using sloan symmetric preordering for "
323 					<< currSolver.s_name
324 					<< " solver" << std::endl);
325 
326 		} else {
327 			pedantic_cerr("sloan preordering is meaningless for "
328 					<< currSolver.s_name
329 					<< " solver" << std::endl);
330 		}
331 	/* nested dissection ? */
332 	} else if (HP.IsKeyWord("nested" "dissection")) {
333 #ifdef USE_METIS
334 		if (currSolver.s_flags & LinSol::SOLVER_FLAGS_ALLOWS_NESTED_DISSECTION) {
335 			cs.AddSolverFlags(LinSol::SOLVER_FLAGS_ALLOWS_NESTED_DISSECTION);
336 			pedantic_cout("using nested dissection symmetric preordering for "
337 					<< currSolver.s_name
338 					<< " solver" << std::endl);
339 
340 		} else {
341 			pedantic_cerr("nested dissection preordering is meaningless for "
342 					<< currSolver.s_name
343 					<< " solver" << std::endl);
344 		}
345 #else //!USE_METIS
346 		silent_cerr("nested dissection permutation not built in;"
347 			"please configure --with-metis to get it"
348 			<< std::endl);
349 		throw ErrGeneric(MBDYN_EXCEPT_ARGS);
350 #endif //USE_METIS
351 	}
352 
353 	/* multithread? */
354 	if (HP.IsKeyWord("multi" "thread") || HP.IsKeyWord("mt")) {
355 		int nThreads = HP.GetInt();
356 
357 		if (currSolver.s_flags & LinSol::SOLVER_FLAGS_ALLOWS_MT_FCT) {
358 			cs.AddSolverFlags(LinSol::SOLVER_FLAGS_ALLOWS_MT_FCT);
359 			if (nThreads < 1) {
360 				silent_cerr("illegal thread number, using 1" << std::endl);
361 				nThreads = 1;
362 			}
363 			cs.SetNumThreads(nThreads);
364 
365 		} else if (nThreads != 1) {
366 			pedantic_cerr("multithread is meaningless for "
367 					<< currSolver.s_name
368 					<< " solver" << std::endl);
369 		}
370 
371 
372 	} else {
373 		if (currSolver.s_flags & LinSol::SOLVER_FLAGS_ALLOWS_MT_FCT) {
374 			int nThreads = get_nprocs();
375 
376 			if (nThreads > 1) {
377 				silent_cout("no multithread requested "
378 						"with a potential "
379 						"of " << nThreads << " CPUs"
380 						<< std::endl);
381 			}
382 
383 			cs.SetNumThreads(nThreads);
384 		}
385 	}
386 
387 	if (HP.IsKeyWord("workspace" "size")) {
388 		integer iWorkSpaceSize = HP.GetInt();
389 		if (iWorkSpaceSize < 0) {
390 			iWorkSpaceSize = 0;
391 		}
392 
393 		switch (cs.GetSolver()) {
394 		case LinSol::Y12_SOLVER:
395 			cs.SetWorkSpaceSize(iWorkSpaceSize);
396 			break;
397 
398 		default:
399 			pedantic_cerr("workspace size is meaningless for "
400 					<< currSolver.s_name
401 					<< " solver" << std::endl);
402 			break;
403 		}
404 	}
405 
406 	if (HP.IsKeyWord("pivot" "factor")) {
407 		doublereal dPivotFactor = HP.GetReal();
408 
409 		if (currSolver.s_pivot_factor == -1.) {
410 			pedantic_cerr("pivot factor is meaningless for "
411 					<< currSolver.s_name
412 					<< " solver" << std::endl);
413 
414 		} else {
415 			if (dPivotFactor <= 0. || dPivotFactor > 1.) {
416 				silent_cerr("pivot factor " << dPivotFactor
417 						<< " is out of bounds; "
418 						"using default "
419 						"(" << currSolver.s_pivot_factor << ")"
420 						<< std::endl);
421 				dPivotFactor = currSolver.s_pivot_factor;
422 			}
423 			cs.SetPivotFactor(dPivotFactor);
424 		}
425 
426 	} else {
427 		if (currSolver.s_pivot_factor != -1.) {
428 			cs.SetPivotFactor(currSolver.s_pivot_factor);
429 		}
430 	}
431 
432 	if (HP.IsKeyWord("drop" "tolerance")) {
433 		doublereal dDropTolerance = HP.GetReal();
434 
435 		if (currSolver.s_drop_tolerance == -1.) {
436 			pedantic_cerr("\"drop tolerance\" is meaningless for "
437 					<< currSolver.s_name
438 					<< " solver" << std::endl);
439 
440 		} else {
441 			if (dDropTolerance < 0.) {
442 				silent_cerr("drop tolerance " << dDropTolerance
443 						<< " is out of bounds; "
444 						"using default "
445 						"(" << currSolver.s_drop_tolerance << ")"
446 						<< std::endl);
447 				dDropTolerance = currSolver.s_drop_tolerance;
448 			}
449 			cs.SetDropTolerance(dDropTolerance);
450 		}
451 
452 	} else {
453 		if (currSolver.s_drop_tolerance != -1.) {
454 			cs.SetDropTolerance(currSolver.s_drop_tolerance);
455 		}
456 	}
457 
458 	if (HP.IsKeyWord("block" "size")) {
459 		integer blockSize = HP.GetInt();
460 		if (blockSize < 1) {
461 			silent_cerr("illegal negative block size; "
462 					"using default" << std::endl);
463 			blockSize = 0;
464 		}
465 
466 		switch (cs.GetSolver()) {
467 		case LinSol::UMFPACK_SOLVER:
468 			cs.SetBlockSize(blockSize);
469 			break;
470 
471 		default:
472 			pedantic_cerr("block size is meaningless for "
473 					<< currSolver.s_name
474 					<< " solver" << std::endl);
475 			break;
476 		}
477 	}
478 
479 	if (HP.IsKeyWord("scale")) {
480 		switch (cs.GetSolver()) {
481 		case LinSol::NAIVE_SOLVER:
482 		case LinSol::KLU_SOLVER:
483 		case LinSol::UMFPACK_SOLVER: {
484 			SolutionManager::ScaleOpt scale;
485 
486 			if (HP.IsKeyWord("no")) {
487 				scale.algorithm = SolutionManager::SCALEA_NONE;
488 
489 			} else if (HP.IsKeyWord("max") || HP.IsKeyWord("row" "max")) {
490 				scale.algorithm = SolutionManager::SCALEA_ROW_MAX;
491 
492 			} else if (HP.IsKeyWord("sum") || HP.IsKeyWord("row" "sum")) {
493 				scale.algorithm = SolutionManager::SCALEA_ROW_SUM;
494 			} else if (HP.IsKeyWord("column" "max")) {
495 				scale.algorithm = SolutionManager::SCALEA_COL_MAX;
496 			} else if (HP.IsKeyWord("column" "sum")) {
497 				scale.algorithm = SolutionManager::SCALEA_COL_SUM;
498 			} else if (HP.IsKeyWord("lapack")) {
499 				scale.algorithm = SolutionManager::SCALEA_LAPACK;
500 			} else if (HP.IsKeyWord("iterative")) {
501 				scale.algorithm = SolutionManager::SCALEA_ITERATIVE;
502 			} else if (HP.IsKeyWord("row" "max" "column" "max") || HP.IsKeyWord("sinkhorn" "knopp")) {
503 				scale.algorithm = SolutionManager::SCALEA_ROW_MAX_COL_MAX;
504 			}
505 
506 			if (HP.IsKeyWord("scale" "tolerance")) {
507 				scale.dTol = HP.GetReal();
508 			}
509 
510 			if (HP.IsKeyWord("scale" "iterations")) {
511 				scale.iMaxIter = HP.GetInt();
512 			}
513 
514 			if (HP.IsKeyWord("once")) {
515 				scale.when = SolutionManager::SCALEW_ONCE;
516 
517 			} else if (HP.IsKeyWord("always")) {
518 				scale.when = SolutionManager::SCALEW_ALWAYS;
519 
520 			} else if (HP.IsKeyWord("never")) {
521 				scale.when = SolutionManager::SCALEW_NEVER;
522 			}
523 
524 			switch (scale.when) {
525 			case SolutionManager::SCALEW_ONCE:
526 			case SolutionManager::SCALEW_ALWAYS:
527 				switch (scale.algorithm) {
528 				case SolutionManager::SCALEA_UNDEF:
529 					scale.algorithm = SolutionManager::SCALEA_LAPACK; // Restore the original behavior for Naive
530 					break;
531 
532 				default:
533 					// Use the value provided by the input file
534 					;
535 				}
536 				break;
537 
538 			case SolutionManager::SCALEW_NEVER:
539 				scale.algorithm = SolutionManager::SCALEA_NONE;
540 				break;
541 			}
542 
543 			if (HP.IsKeyWord("verbose") && HP.GetYesNoOrBool()) {
544 				scale.uFlags |= SolutionManager::SCALEF_VERBOSE;
545 			}
546 
547 			if (HP.IsKeyWord("warnings") && HP.GetYesNoOrBool()) {
548 				scale.uFlags |= SolutionManager::SCALEF_WARN;
549 			}
550 
551 			unsigned uCondFlag = 0;
552 
553 			if (HP.IsKeyWord("print" "condition" "number")) {
554 				if (HP.IsKeyWord("norm")) {
555 					if (HP.IsKeyWord("inf")) {
556 						uCondFlag = SolutionManager::SCALEF_COND_NUM_INF;
557 					} else {
558 						const doublereal dNorm = HP.GetReal();
559 
560 						if (dNorm != 1.) {
561 							silent_cerr("Only one norm or infinity norm are supported for condition numbers at line " << HP.GetLineData() << std::endl);
562 							throw ErrGeneric(MBDYN_EXCEPT_ARGS);
563 						}
564 
565 						uCondFlag = SolutionManager::SCALEF_COND_NUM_1;
566 					}
567 				} else {
568 					uCondFlag = SolutionManager::SCALEF_COND_NUM_1;
569 				}
570 			}
571 
572 			if (uCondFlag != 0 && HP.GetYesNoOrBool()) {
573 				scale.uFlags |= uCondFlag;
574 			}
575 
576 			if (!cs.SetScale(scale)) {
577 				silent_cerr("Warning: Scale options are not available for "
578 						    << cs.GetSolverName()
579 						    << " at line "
580 						    << HP.GetLineData() << std::endl);
581 			}
582 
583 			} break;
584 
585 		default:
586 			pedantic_cerr("scale is meaningless for "
587 					<< currSolver.s_name
588 					<< " solver" << std::endl);
589 			break;
590 		}
591 	}
592 
593 	if (HP.IsKeyWord("max" "iterations")) {
594 		if (!cs.SetMaxIterations(HP.GetInt())) {
595 			silent_cerr("Warning: iterative refinement is not supported by " << cs.GetSolverName() << " at line " << HP.GetLineData() << std::endl);
596 		}
597 	}
598 
599 	switch (cs.GetSolver()) {
600 	case LinSol::NAIVE_SOLVER:
601 		if (!(cs.GetSolverFlags() & LinSol::SOLVER_FLAGS_ALLOWS_COLAMD)) {
602 			silent_cout("warning: \"naive\" solver should be used with \"colamd\"" << std::endl);
603 		}
604 		break;
605 
606 	// add more warnings...
607 
608 	default:
609 		break;
610 	}
611 }
612 
RestartLinSol(std::ostream & out,const LinSol & cs)613 std::ostream & RestartLinSol(std::ostream& out, const LinSol& cs)
614 {
615 	out << cs.GetSolverName();
616 
617 	const LinSol::solver_t	currSolver = ::solver[cs.GetSolver()];
618 
619 	if (cs.GetSolverFlags() != currSolver.s_default_flags) {
620 		unsigned f = cs.GetSolverFlags();
621 		if((f & LinSol::SOLVER_FLAGS_ALLOWS_MAP) ==
622 			LinSol::SOLVER_FLAGS_ALLOWS_MAP) {
623 			/*Map*/
624 			out << ", map ";
625 		}
626 		if((f & LinSol::SOLVER_FLAGS_ALLOWS_CC) ==
627 			LinSol::SOLVER_FLAGS_ALLOWS_CC) {
628 			/*column compressed*/
629 			out << ", cc ";
630 		}
631 		if((f & LinSol::SOLVER_FLAGS_ALLOWS_DIR) ==
632 			LinSol::SOLVER_FLAGS_ALLOWS_DIR) {
633 			/*direct access*/
634 			out << ", dir ";
635 		}
636 		if((f & LinSol::SOLVER_FLAGS_ALLOWS_COLAMD) ==
637 			LinSol::SOLVER_FLAGS_ALLOWS_COLAMD) {
638 			/*colamd*/
639 			out << ", colamd ";
640 		}
641 		unsigned nt = cs.GetNumThreads();
642 		if(((f & LinSol::SOLVER_FLAGS_ALLOWS_MT_FCT) ==
643 			 LinSol::SOLVER_FLAGS_ALLOWS_MT_FCT) &&
644 			(nt != 1)) {
645 			/*multi thread*/
646 			out << ", mt , " << nt;
647 		}
648 	}
649 	integer ws = cs.iGetWorkSpaceSize();
650 	if(ws > 0) {
651 		/*workspace size*/
652 		out << ", workspace size, " << ws;
653 	}
654 	doublereal pf = cs.dGetPivotFactor();
655 	if(pf != -1.) {
656 		/*pivot factor*/
657 		out << ", pivot factor, " << pf;
658 	}
659 	unsigned bs = cs.GetBlockSize();
660 	if(bs != 0) {
661 		/*block size*/
662 		out << ", block size, " << bs ;
663 	}
664 	out << ";" << std::endl;
665 	return out;
666 }
667