1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbwrap/naivewrap.cc,v 1.65 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 /*
33 * The Naive Solver is copyright (C) 2004 by
34 * Paolo Mantegazza <mantegazza@aero.polimi.it>
35 */
36
37 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
38
39 #include <sys/types.h>
40 #include <unistd.h>
41 #include <signal.h>
42
43 #include <set>
44
45 #include "spmh.h"
46 #include "spmapmh.h"
47 #include "naivewrap.h"
48 #include "mthrdslv.h"
49 #include "dgeequ.h"
50
51 /* NaiveSolver - begin */
NaiveSolver(const integer & size,const doublereal & dMP,NaiveMatrixHandler * const a)52 NaiveSolver::NaiveSolver(const integer &size, const doublereal& dMP,
53 NaiveMatrixHandler *const a)
54 : LinearSolver(0),
55 iSize(size),
56 dMinPiv(dMP < 0 ? 0 : dMP),
57 piv(size),
58 A(a)
59 {
60 NO_OP;
61 }
62
~NaiveSolver(void)63 NaiveSolver::~NaiveSolver(void)
64 {
65 NO_OP;
66 }
67
68 void
SetMat(NaiveMatrixHandler * const a)69 NaiveSolver::SetMat(NaiveMatrixHandler *const a)
70 {
71 A = a;
72 }
73
74 void
Reset(void)75 NaiveSolver::Reset(void)
76 {
77 bHasBeenReset = true;
78 }
79
80 void
Solve(void) const81 NaiveSolver::Solve(void) const
82 {
83 if (bHasBeenReset) {
84 const_cast<NaiveSolver *>(this)->Factor();
85 bHasBeenReset = false;
86 }
87
88 integer rc = naivslv(A->ppdRows, iSize, A->piNzc, A->ppiCols,
89 LinearSolver::pdRhs, LinearSolver::pdSol, &piv[0]);
90 integer err = (rc & NAIVE_MASK);
91 if (err) {
92 switch (err) {
93 case NAIVE_ERANGE:
94 silent_cerr("NaiveSolver: ERANGE"
95 << std::endl);
96 break;
97
98 default:
99 silent_cerr("NaiveSolver: (" << rc << ")"
100 << std::endl);
101 break;
102 }
103
104 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
105 }
106 }
107
108 void
Factor(void)109 NaiveSolver::Factor(void)
110 throw(LinearSolver::ErrFactor)
111 {
112 integer rc = naivfct(A->ppdRows, iSize,
113 A->piNzr, A->ppiRows,
114 A->piNzc, A->ppiCols,
115 A->ppnonzero,
116 &piv[0], dMinPiv);
117
118 integer err = (rc & NAIVE_MASK);
119 if (err) {
120 integer idx = (rc & NAIVE_MAX);
121 switch (err) {
122 case NAIVE_ENULCOL:
123 silent_cerr("NaiveSolver: ENULCOL(" << idx << ")"
124 << std::endl);
125 throw LinearSolver::ErrNullColumn(idx, MBDYN_EXCEPT_ARGS);
126
127 case NAIVE_ENOPIV:
128 silent_cerr("NaiveSolver: ENOPIV(" << idx << ")"
129 << std::endl);
130 throw LinearSolver::ErrNoPivot(idx + 1, MBDYN_EXCEPT_ARGS);
131
132 case NAIVE_ERANGE:
133 silent_cerr("NaiveSolver: ERANGE"
134 << std::endl);
135 break;
136
137 default:
138 silent_cerr("NaiveSolver: (" << rc << ")"
139 << std::endl);
140 break;
141 }
142
143 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
144 }
145 }
146
147 /* NaiveSolver - end */
148
149 /* NaiveSparseSolutionManager - begin */
150
NaiveSparseSolutionManager(const integer Dim,const doublereal dMP,const ScaleOpt & s)151 NaiveSparseSolutionManager::NaiveSparseSolutionManager(const integer Dim,
152 const doublereal dMP, const ScaleOpt& s)
153 : A(0),
154 VH(Dim),
155 scale(s),
156 pMatScale(0)
157 {
158 switch (scale.algorithm) {
159 case SCALEA_NONE:
160 case SCALEA_UNDEF:
161 scale.when = SCALEW_NEVER; // Default scaling is not supported for this solver
162 break;
163
164 default:
165 // Allocate MatrixScale<T> on demand
166 ;
167 }
168
169 SAFENEWWITHCONSTRUCTOR(A, NaiveMatrixHandler, NaiveMatrixHandler(Dim));
170 SAFENEWWITHCONSTRUCTOR(pLS, NaiveSolver, NaiveSolver(Dim, dMP, A));
171
172 pLS->pdSetResVec(VH.pdGetVec());
173 pLS->pdSetSolVec(VH.pdGetVec());
174
175 pLS->SetSolutionManager(this);
176 }
177
~NaiveSparseSolutionManager(void)178 NaiveSparseSolutionManager::~NaiveSparseSolutionManager(void)
179 {
180 if (A != 0) {
181 SAFEDELETE(A);
182 A = 0;
183 }
184
185 if (pMatScale != 0) {
186 SAFEDELETE(pMatScale);
187 }
188 }
189
190 void
MatrReset()191 NaiveSparseSolutionManager::MatrReset()
192 {
193 pLS->Reset();
194 }
195
196
197 template <class MH>
198 void
ScaleMatrixAndRightHandSide(MH & mh)199 NaiveSparseSolutionManager::ScaleMatrixAndRightHandSide(MH& mh)
200 {
201 if (scale.when != SCALEW_NEVER) {
202 MatrixScale<MH>& rMatScale = GetMatrixScale<MH>();
203
204 if (pLS->bReset()) {
205 if (!rMatScale.bGetInitialized()
206 || scale.when == SolutionManager::SCALEW_ALWAYS) {
207 // (re)compute
208 rMatScale.ComputeScaleFactors(mh);
209 }
210 // in any case scale matrix and right-hand-side
211 rMatScale.ScaleMatrix(mh);
212
213 if (silent_err) {
214 rMatScale.Report(std::cerr);
215 }
216 }
217
218 rMatScale.ScaleRightHandSide(VH);
219 }
220 }
221
222 template <typename MH>
GetMatrixScale()223 MatrixScale<MH>& NaiveSparseSolutionManager::GetMatrixScale()
224 {
225 if (pMatScale == 0) {
226 pMatScale = MatrixScale<MH>::Allocate(scale);
227 }
228
229 // Will throw std::bad_cast if the type does not match
230 return dynamic_cast<MatrixScale<MH>&>(*pMatScale);
231 }
232
233 void
ScaleSolution(void)234 NaiveSparseSolutionManager::ScaleSolution(void)
235 {
236 if (scale.when != SCALEW_NEVER) {
237 ASSERT(pMatScale != 0);
238 // scale solution
239 pMatScale->ScaleSolution(VH);
240 }
241 }
242
243 /* Risolve il sistema Fattorizzazione + Backward Substitution */
244 void
Solve(void)245 NaiveSparseSolutionManager::Solve(void)
246 {
247 ScaleMatrixAndRightHandSide(*A);
248
249 pLS->Solve();
250
251 ScaleSolution();
252 }
253
254 /* Rende disponibile l'handler per la matrice */
255 MatrixHandler*
pMatHdl(void) const256 NaiveSparseSolutionManager::pMatHdl(void) const
257 {
258 return A;
259 }
260
261 /* Rende disponibile l'handler per il termine noto */
262 MyVectorHandler*
pResHdl(void) const263 NaiveSparseSolutionManager::pResHdl(void) const
264 {
265 return &VH;
266 }
267
268 /* Rende disponibile l'handler per la soluzione */
269 MyVectorHandler*
pSolHdl(void) const270 NaiveSparseSolutionManager::pSolHdl(void) const
271 {
272 return &VH;
273 }
274
275 /* NaiveSparseSolutionManager - end */
276
277 /* NaivePermSparseSolutionManager - begin */
278
279
280 template<class T>
NaiveSparsePermSolutionManager(const integer Dim,const doublereal dMP,const ScaleOpt & scale)281 NaiveSparsePermSolutionManager<T>::NaiveSparsePermSolutionManager(
282 const integer Dim,
283 const doublereal dMP,
284 const ScaleOpt& scale)
285 : NaiveSparseSolutionManager(Dim, dMP, scale),
286 dMinPiv(dMP < 0 ? 0 : dMP),
287 TmpH(Dim),
288 ePermState(PERM_NO)
289 {
290 perm.resize(Dim, 0);
291 invperm.resize(Dim, 0);
292
293 // replace matrix handler
294 SAFEDELETE(A);
295 A = 0;
296 SAFENEWWITHCONSTRUCTOR(A, NaivePermMatrixHandler,
297 NaivePermMatrixHandler(Dim, perm, invperm));
298
299 dynamic_cast<NaiveSolver *>(pLS)->SetMat(A);
300
301 MatrInitialize();
302 }
303
304 template<class T>
~NaiveSparsePermSolutionManager(void)305 NaiveSparsePermSolutionManager<T>::~NaiveSparsePermSolutionManager(void)
306 {
307 NO_OP;
308 }
309
310 template<class T>
311 void
MatrReset(void)312 NaiveSparsePermSolutionManager<T>::MatrReset(void)
313 {
314 if (ePermState == PERM_INTERMEDIATE) {
315 ePermState = PERM_READY;
316 }
317
318 NaiveSparseSolutionManager::MatrReset();
319 }
320
321 template<class T>
322 void
BackPerm(void)323 NaiveSparsePermSolutionManager<T>::BackPerm(void)
324 {
325 /* NOTE: use whatever is stored in pLS - someone could
326 * trick us into using its memory */
327 doublereal *pd = pLS->pdGetResVec();
328
329 ASSERT(pd != TmpH.pdGetVec());
330
331 for (integer i = 0; i < A->iGetNumCols(); i++) {
332 pd[invperm[i]] = TmpH(i + 1);
333 }
334 }
335
336
337 /* Risolve il sistema: Fattorizzazione + Backward Substitution */
338 template<class T>
339 void
Solve(void)340 NaiveSparsePermSolutionManager<T>::Solve(void)
341 {
342 doublereal *pd = 0;
343
344 // Will throw std::bad_cast if the data type does not match
345 ScaleMatrixAndRightHandSide(dynamic_cast<NaivePermMatrixHandler&>(*A));
346
347 if (ePermState == PERM_NO) {
348 ComputePermutation();
349
350 } else if (ePermState == PERM_READY) {
351 /* We need to use local storage to allow BackPerm();
352 * save and restore original pointer */
353 pd = pLS->pdSetSolVec(TmpH.pdGetVec());
354 }
355
356 pLS->Solve();
357
358 if (ePermState == PERM_READY) {
359 BackPerm();
360
361 ASSERT(pd != 0);
362 pLS->pdSetSolVec(pd);
363 }
364
365 ScaleSolution();
366 }
367
368 /* Inizializzatore "speciale" */
369 template<class T>
370 void
MatrInitialize()371 NaiveSparsePermSolutionManager<T>::MatrInitialize()
372 {
373 ePermState = PERM_NO;
374 for (integer i = 0; i < A->iGetNumRows(); i++) {
375 perm[i] = i;
376 invperm[i] = i;
377 }
378
379 MatrReset();
380 }
381
382 // explicit specializations
383
384 extern "C" {
385 #include "colamd.h"
386 }
387
388 template<>
389 void
ComputePermutation(void)390 NaiveSparsePermSolutionManager<Colamd_ordering>::ComputePermutation(void)
391 {
392 std::vector<integer> Ai;
393 A->MakeCCStructure(Ai, invperm);
394 doublereal knobs[COLAMD_KNOBS];
395 integer stats[COLAMD_STATS];
396 integer Alen = mbdyn_colamd_recommended(Ai.size(), A->iGetNumRows(),
397 A->iGetNumCols());
398 Ai.resize(Alen);
399 mbdyn_colamd_set_defaults(knobs);
400 if (!mbdyn_colamd(A->iGetNumRows(), A->iGetNumCols(), Alen,
401 &Ai[0], &invperm[0], knobs, stats))
402 {
403 silent_cerr("colamd permutation failed" << std::endl);
404 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
405 }
406 for (integer i = 0; i < A->iGetNumRows(); i++) {
407 perm[invperm[i]] = i;
408 }
409 ePermState = PERM_INTERMEDIATE;
410 }
411
412
413 #ifdef USE_BOOST
414 /* NaivePermSparseSolutionManager - end */
415 #include "boost/config.hpp"
416 #include "boost/graph/adjacency_list.hpp"
417 #include "boost/graph/properties.hpp"
418 #include "boost/graph/bandwidth.hpp"
419 #include <boost/graph/wavefront.hpp>
420
421 #ifdef HAVE_BOOST_GRAPH_CUTHILL_MCKEE_ORDERING_HPP
422 #include "boost/graph/cuthill_mckee_ordering.hpp"
423 #endif /* HAVE_BOOST_GRAPH_CUTHILL_MCKEE_ORDERING_HPP */
424 #ifdef HAVE_BOOST_GRAPH_KING_ORDERING_HPP
425 #include "boost/graph/king_ordering.hpp"
426 #endif /* HAVE_BOOST_GRAPH_KING_ORDERING_HPP */
427 #ifdef HAVE_BOOST_GRAPH_MINIMUM_DEGREE_ORDERING_HPP
428 #include <boost/graph/sloan_ordering.hpp>
429 #endif /* HAVE_BOOST_GRAPH_MINIMUM_DEGREE_ORDERING_HPP */
430 #ifdef HAVE_BOOST_GRAPH_SLOAN_ORDERING_HPP
431 #include "boost/graph/minimum_degree_ordering.hpp"
432 #endif /* HAVE_BOOST_GRAPH_SLOAN_ORDERING_HPP */
433
434 #ifdef HAVE_BOOST_GRAPH_CUTHILL_MCKEE_ORDERING_HPP
435 template<>
436 void
ComputePermutation(void)437 NaiveSparsePermSolutionManager<rcmk_ordering>::ComputePermutation(void)
438 {
439 std::vector<integer> Ai;
440 std::vector<integer> Ac;
441
442
443 invperm.resize(A->iGetNumCols());
444 A->MakeCCStructure(Ai, Ac);
445
446
447 /* boost */
448
449 typedef boost::adjacency_list<
450 boost::setS,
451 boost::vecS,
452 boost::undirectedS,
453 boost::property<
454 boost::vertex_color_t,
455 boost::default_color_type,
456 boost::property<
457 boost::vertex_degree_t,
458 integer
459 // ,
460 // boost::property<
461 // boost::vertex_priority_t,
462 // double
463 // >
464 >
465 >
466 > Graph;
467 typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;
468 typedef boost::graph_traits<Graph>::vertices_size_type size_type;
469
470
471 Graph G(A->iGetNumRows());
472 for (int col=0; col<A->iGetNumCols(); col++) {
473 for (int i = Ac[col]; i < Ac[col + 1]; i++) {
474 int row = Ai[i];
475 if (row != col) {
476 boost::add_edge(row, col, G);
477 boost::add_edge(col, row, G);
478 }
479 }
480 }
481
482
483 std::vector<Vertex> inv_perm(num_vertices(G));
484 boost::cuthill_mckee_ordering(G, inv_perm.rbegin());
485 #if 0
486 boost::sloan_ordering(G, inv_perm.rbegin(),
487 boost::get(boost::vertex_color, G),
488 boost::make_degree_map(G),
489 boost::get(boost::vertex_priority, G));
490 #endif
491
492 for (integer i = 0; i < A->iGetNumRows(); i++) {
493 invperm[i] = inv_perm[i];
494 perm[invperm[i]] = i;
495 }
496 ePermState = PERM_INTERMEDIATE;
497 }
498 #endif /* HAVE_BOOST_GRAPH_CUTHILL_MCKEE_ORDERING_HPP */
499
500 #ifdef HAVE_BOOST_GRAPH_SLOAN_ORDERING_HPP
501 template<>
502 void
ComputePermutation(void)503 NaiveSparsePermSolutionManager<sloan_ordering>::ComputePermutation(void)
504 {
505 std::vector<integer> Ai;
506 std::vector<integer> Ac;
507
508
509 invperm.resize(A->iGetNumCols());
510 A->MakeCCStructure(Ai, Ac);
511
512
513 /* boost */
514
515 typedef boost::adjacency_list<
516 boost::setS,
517 boost::vecS,
518 boost::undirectedS,
519 boost::property<
520 boost::vertex_color_t,
521 boost::default_color_type,
522 boost::property<
523 boost::vertex_degree_t,
524 integer,
525 boost::property<
526 boost::vertex_priority_t,
527 double
528 >
529 >
530 >
531 > Graph;
532 typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;
533 typedef boost::graph_traits<Graph>::vertices_size_type size_type;
534
535
536 Graph G(A->iGetNumRows());
537 for (int col=0; col<A->iGetNumCols(); col++) {
538 for (int i = Ac[col]; i < Ac[col + 1]; i++) {
539 int row = Ai[i];
540 if (row != col) {
541 boost::add_edge(row, col, G);
542 boost::add_edge(col, row, G);
543 }
544 }
545 }
546
547
548 std::vector<Vertex> inv_perm(num_vertices(G));
549 std::vector<Vertex> _perm(num_vertices(G));
550 boost::sloan_ordering(G, _perm.begin(),
551 boost::get(boost::vertex_color, G),
552 boost::make_degree_map(G),
553 boost::get(boost::vertex_priority, G));
554
555 for (integer i = 0; i < A->iGetNumRows(); i++) {
556 // invperm[i] = inv_perm[i];
557 // perm[invperm[i]] = i;
558 perm[i] = _perm[i];
559 invperm[perm[i]] = i;
560 }
561 ePermState = PERM_INTERMEDIATE;
562 }
563 #endif /* HAVE_BOOST_GRAPH_SLOAN_ORDERING_HPP */
564
565 #ifdef HAVE_BOOST_GRAPH_KING_ORDERING_HPP
566 template<>
567 void
ComputePermutation(void)568 NaiveSparsePermSolutionManager<king_ordering>::ComputePermutation(void)
569 {
570 std::vector<integer> Ai;
571 std::vector<integer> Ac;
572
573
574 invperm.resize(A->iGetNumCols());
575 A->MakeCCStructure(Ai, Ac);
576
577
578 /* boost */
579
580 typedef boost::adjacency_list<
581 boost::setS,
582 boost::vecS,
583 boost::undirectedS,
584 boost::property<
585 boost::vertex_color_t,
586 boost::default_color_type,
587 boost::property<
588 boost::vertex_degree_t,
589 integer
590 >
591 >
592 > Graph;
593 typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;
594 typedef boost::graph_traits<Graph>::vertices_size_type size_type;
595
596
597 Graph G(A->iGetNumRows());
598 for (int col=0; col<A->iGetNumCols(); col++) {
599 for (int i = Ac[col]; i < Ac[col + 1]; i++) {
600 int row = Ai[i];
601 if (row != col) {
602 boost::add_edge(row, col, G);
603 boost::add_edge(col, row, G);
604 }
605 }
606 }
607
608 std::vector<Vertex> inv_perm(num_vertices(G));
609 boost::king_ordering(G, inv_perm.rbegin());
610
611 for (integer i = 0; i < A->iGetNumRows(); i++) {
612 invperm[i] = inv_perm[i];
613 perm[invperm[i]] = i;
614 }
615 ePermState = PERM_INTERMEDIATE;
616 }
617 #endif /* HAVE_BOOST_GRAPH_KING_ORDERING_HPP */
618
619 #ifdef HAVE_BOOST_GRAPH_MINIMUM_DEGREE_ORDERING_HPP
620 template<>
621 void
ComputePermutation(void)622 NaiveSparsePermSolutionManager<md_ordering>::ComputePermutation(void)
623 {
624 std::vector<integer> Ai;
625 std::vector<integer> Ac;
626
627
628 invperm.resize(A->iGetNumCols());
629 A->MakeCCStructure(Ai, Ac);
630
631 typedef boost::adjacency_list<
632 boost::setS,
633 boost::vecS,
634 boost::directedS
635 > Graph;
636 typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;
637 typedef boost::graph_traits<Graph>::vertices_size_type size_type;
638
639
640 Graph G(A->iGetNumRows());
641 for (int col=0; col<A->iGetNumCols(); col++) {
642 for (int i = Ac[col]; i < Ac[col + 1]; i++) {
643 int row = Ai[i];
644 if (row != col) {
645 boost::add_edge(row, col, G);
646 boost::add_edge(col, row, G);
647 }
648 }
649 }
650
651 boost::property_map<Graph, boost::vertex_index_t>::type
652 id = boost::get(boost::vertex_index, G);
653 std::vector<Vertex> inv_perm(num_vertices(G));
654 std::vector<Vertex> _perm(num_vertices(G));
655 std::vector<integer> degree(A->iGetNumRows(), 0);
656 std::vector<integer> supernode_sizes(A->iGetNumRows(), 1); // init has to be 1
657 int delta = 0;
658
659 boost::minimum_degree_ordering(
660 G,
661 make_iterator_property_map(°ree[0], id, degree[0]),
662 &invperm[0],
663 &perm[0],
664 boost::make_iterator_property_map(
665 &supernode_sizes[0],
666 id,
667 supernode_sizes[0]
668 ),
669 delta,
670 id
671 );
672
673 // for (integer i = 0; i < A->iGetNumRows(); i++) {
674 // invperm[i] = inv_perm[i];
675 // perm[invperm[i]] = i;
676 // }
677 ePermState = PERM_INTERMEDIATE;
678
679
680
681 }
682 #endif /* HAVE_BOOST_GRAPH_MINIMUM_DEGREE_ORDERING_HPP */
683 #endif /* USE_BOOST */
684
685 #ifdef USE_METIS
686 extern "C" {
687 #include "metiswrap.h"
688 }
689
690 template<>
691 void
ComputePermutation(void)692 NaiveSparsePermSolutionManager<metis_ordering>::ComputePermutation(void)
693 {
694 std::vector<integer> Ai;
695 std::vector<integer> Ac;
696
697
698 invperm.resize(A->iGetNumCols());
699 A->MakeCCStructure(Ai, Ac);
700
701
702 /* metis */
703
704 typedef std::set<integer> row_cont_type;
705 std::vector<row_cont_type> col_indices(A->iGetNumCols());
706
707 integer NZ = 0;
708 {
709 integer ncols = A->iGetNumCols();
710 for (integer col_id=0; col_id<ncols; col_id++) {
711 for (integer i=Ac[col_id]; i<Ac[col_id+1]; i++) {
712 int row_id = Ai[i];
713 if (row_id != col_id) {
714 //A
715 row_cont_type& row1 = col_indices[col_id];
716 if (row1.find(row_id) == row1.end()) {
717 NZ++;
718 row1.insert(row_id);
719 }
720 // + A^T
721 row_cont_type& row2 = col_indices[row_id];
722 if (row2.find(col_id) == row2.end()) {
723 NZ++;
724 row2.insert(col_id);
725 }
726 }
727 }
728 }
729 }
730 Ai.resize(NZ);
731 {
732 integer x_ptr = 0;
733 row_cont_type::iterator ri;
734 row_cont_type::const_iterator re;
735 integer NCols = A->iGetNumCols();
736 for (integer col = 0; col < NCols; col++) {
737 Ac[col] = x_ptr;
738 re = col_indices[col].end();
739 for (ri = col_indices[col].begin(); ri != re; ri++) {
740 Ai[x_ptr] = *ri;
741 x_ptr++;
742 }
743 }
744 Ac[NCols] = x_ptr;
745 }
746 int numflag = 0;
747 int n = A->iGetNumCols();
748 int options[8] = {1, 3, 1, 3, 0, 1, 200, 3};
749 // METIS_EdgeND(&n, &(Ac[0]), &Ai[0], &numflag, options, &(perm[0]), &(invperm[0]));
750 METIS_NodeND(&n, &(Ac[0]), &Ai[0], &numflag, options, &perm[0], &invperm[0]);
751
752 ePermState = PERM_INTERMEDIATE;
753
754 }
755
756 #endif //USE_METIS
757
758 // #ifdef HAVE_UMFPACK
759 // extern "C" {
760 // #include "amd.h"
761 // }
762 // template<>
763 // void
764 // NaiveSparsePermSolutionManager<amd_ordering>::ComputePermutation(void)
765 // {
766 // std::vector<integer> Ai;
767 // std::vector<integer> Ac;
768 //
769 //
770 // invperm.resize(A->iGetNumCols());
771 // A->MakeCCStructure(Ai, Ac);
772 //
773 // /* amd */
774 // double Control[AMD_CONTROL], Info[AMD_INFO];
775 // amd_defaults(Control);
776 // amd_order(A->iGetNumCols(), &(Ac[0]), &(Ai[0]), &(invperm[0]), Control, Info);
777 // for (integer i = 0; i < A->iGetNumRows(); i++) {
778 // perm[invperm[i]] = i;
779 // }
780 //
781 //
782 //
783 // }
784 // #endif //HAVE_UMFPACK
785
786
787 //explicit instantiations:
788 template class NaiveSparsePermSolutionManager<Colamd_ordering>;
789 #ifdef USE_BOOST
790 #ifdef HAVE_BOOST_GRAPH_CUTHILL_MCKEE_ORDERING_HPP
791 template class NaiveSparsePermSolutionManager<rcmk_ordering>;
792 #endif /* HAVE_BOOST_GRAPH_CUTHILL_MCKEE_ORDERING_HPP */
793 #ifdef HAVE_BOOST_GRAPH_KING_ORDERING_HPP
794 template class NaiveSparsePermSolutionManager<king_ordering>;
795 #endif /* HAVE_BOOST_GRAPH_KING_ORDERING_HPP */
796 #ifdef HAVE_BOOST_GRAPH_MINIMUM_DEGREE_ORDERING_HPP
797 template class NaiveSparsePermSolutionManager<md_ordering>;
798 #endif /* HAVE_BOOST_GRAPH_MINIMUM_DEGREE_ORDERING_HPP */
799 #ifdef HAVE_BOOST_GRAPH_SLOAN_ORDERING_HPP
800 template class NaiveSparsePermSolutionManager<sloan_ordering>;
801 #endif /* HAVE_BOOST_GRAPH_SLOAN_ORDERING_HPP */
802 #endif /* USE_BOOST */
803 #ifdef USE_METIS
804 template class NaiveSparsePermSolutionManager<metis_ordering>;
805 #endif //USE_METIS
806 // #ifdef HAVE_UMFPACK
807 // template class NaiveSparsePermSolutionManager<amd_ordering>;
808 // #endif //HAVE_UMFPACK
809
810