1 /*===========================================================================*\
2 * *
3 * CoMISo *
4 * Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
5 * www.rwth-graphics.de *
6 * *
7 *---------------------------------------------------------------------------*
8 * This file is part of CoMISo. *
9 * *
10 * CoMISo is free software: you can redistribute it and/or modify *
11 * it under the terms of the GNU General Public License as published by *
12 * the Free Software Foundation, either version 3 of the License, or *
13 * (at your option) any later version. *
14 * *
15 * CoMISo is distributed in the hope that it will be useful, *
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18 * GNU General Public License for more details. *
19 * *
20 * You should have received a copy of the GNU General Public License *
21 * along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
22 * *
23 \*===========================================================================*/
24
25 //== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
26 #include <CoMISo/Config/config.hh>
27 #if COMISO_SUITESPARSE_AVAILABLE
28 //=============================================================================
29
30 #include "CholmodSolver.hh"
31
32
33 namespace COMISO {
34
CholmodSolver()35 CholmodSolver::CholmodSolver()
36 {
37 mp_cholmodCommon = new cholmod_common;
38 cholmod_start( mp_cholmodCommon );
39
40 mp_L = 0;
41
42 show_timings_ = false;
43
44 mp_cholmodCommon->nmethods = 1;
45 // use AMD ordering
46 mp_cholmodCommon->method[0].ordering = CHOLMOD_AMD ;
47
48 // use METIS ordering
49 // mp_cholmodCommon->method[0].ordering = CHOLMOD_METIS ;
50
51 // try all methods
52 // mp_cholmodCommon->nmethods = 9;
53 }
54
55
56 //-----------------------------------------------------------------------------
57
58
~CholmodSolver()59 CholmodSolver::~CholmodSolver()
60 {
61 if( mp_L )
62 {
63 cholmod_free_factor( &mp_L, mp_cholmodCommon );
64 }
65
66 cholmod_finish( mp_cholmodCommon );
67 delete mp_cholmodCommon;
68 mp_cholmodCommon = NULL;
69 }
70
71
72 //-----------------------------------------------------------------------------
73
74
calc_system(const std::vector<int> & _colptr,const std::vector<int> & _rowind,const std::vector<double> & _values)75 bool CholmodSolver::calc_system( const std::vector<int>& _colptr,
76 const std::vector<int>& _rowind,
77 const std::vector<double>& _values)
78 {
79 if(show_timings_) sw_.start();
80
81 colptr_ = _colptr;
82 rowind_ = _rowind;
83 values_ = _values;
84
85 int n = colptr_.size()-1;
86
87 cholmod_sparse matA;
88
89 matA.nrow = n;
90 matA.ncol = n;
91 matA.nzmax = _values.size();
92
93 matA.p = &colptr_[0];
94 matA.i = &rowind_[0];
95 matA.x = &values_[0];
96 matA.nz = 0;
97 matA.z = 0;
98
99 // matA.stype = -1;
100 matA.stype = 1;
101 matA.itype = CHOLMOD_INT;
102 matA.xtype = CHOLMOD_REAL;
103 matA.dtype = CHOLMOD_DOUBLE;
104 matA.sorted = 1;
105 matA.packed = 1;
106
107
108 // clean up
109 if( mp_L )
110 {
111 cholmod_free_factor( &mp_L, mp_cholmodCommon );
112 mp_L = 0;
113 }
114
115 if(show_timings_)
116 {
117 std::cerr << " Cholmod Timing cleanup: " << sw_.stop()/1000.0 << "s\n";
118 sw_.start();
119 }
120
121 if( !(mp_L = cholmod_analyze( &matA, mp_cholmodCommon )) )
122 {
123 std::cout << "cholmod_analyze failed" << std::endl;
124 return false;
125 }
126
127 // // show selected ordering method
128 // std::cerr << "best ordering was: " << mp_cholmodCommon->selected << std::endl;
129 // std::cerr << "current ordering was: " << mp_cholmodCommon->current << std::endl;
130
131
132 if(show_timings_)
133 {
134 std::cerr << " Cholmod Timing analyze: " << sw_.stop()/1000.0 << "s\n";
135 sw_.start();
136 }
137
138 if( !cholmod_factorize( &matA, mp_L, mp_cholmodCommon ) )
139 {
140 std::cout << "cholmod_factorize failed" << std::endl;
141 return false;
142 }
143
144 if(show_timings_)
145 {
146 std::cerr << " Cholmod Timing factorize: " << sw_.stop()/1000.0 << "s\n";
147 sw_.start();
148 }
149
150 return true;
151 }
152
153
154 //-----------------------------------------------------------------------------
155
156
update_system(const std::vector<int> & _colptr,const std::vector<int> & _rowind,const std::vector<double> & _values)157 bool CholmodSolver::update_system( const std::vector<int>& _colptr,
158 const std::vector<int>& _rowind,
159 const std::vector<double>& _values )
160 {
161 if( !mp_L )
162 return false;
163
164 colptr_ = _colptr;
165 rowind_ = _rowind;
166 values_ = _values;
167 int n = colptr_.size()-1;
168
169 cholmod_sparse matA;
170
171 matA.nrow = n;
172 matA.ncol = n;
173 matA.nzmax = _values.size();
174
175 matA.p = &colptr_[0];
176 matA.i = &rowind_[0];
177 matA.x = &values_[0];
178 matA.nz = 0;
179 matA.z = 0;
180
181 // matA.stype = -1;
182 matA.stype = 1;
183 matA.itype = CHOLMOD_INT;
184 matA.xtype = CHOLMOD_REAL;
185 matA.dtype = CHOLMOD_DOUBLE;
186 matA.sorted = 1;
187 matA.packed = 1;
188
189
190 if( !cholmod_factorize( &matA, mp_L, mp_cholmodCommon ) )
191 {
192 std::cout << "cholmod_factorize failed" << std::endl;
193 return false;
194 }
195
196 return true;
197 }
198
199
200 //-----------------------------------------------------------------------------
201
202
solve(double * _x,double * _b)203 bool CholmodSolver::solve( double * _x, double * _b)
204 {
205 const unsigned int n = colptr_.size() - 1;
206
207 cholmod_dense *x, b;
208
209
210 b.nrow = n;
211 b.ncol = 1;
212 b.nzmax = n;
213 b.d = b.nrow;
214 b.x = _b;
215 b.z = 0;
216 b.xtype = CHOLMOD_REAL;
217 b.dtype = CHOLMOD_DOUBLE;
218
219 if( !(x = cholmod_solve( CHOLMOD_A, mp_L, &b, mp_cholmodCommon )) )
220 {
221 std::cout << "cholmod_solve failed" << std::endl;
222 return false;
223 }
224
225 for( unsigned int i = 0; i < n; ++i )
226 _x[i] = ((double*)x->x)[i];
227
228 cholmod_free_dense( &x, mp_cholmodCommon );
229
230 return true;
231 }
232
233
234 //-----------------------------------------------------------------------------
235
dimension()236 int CholmodSolver::dimension()
237 {
238 return std::max(int(0), (int)(colptr_.size()-1));
239 }
240
241 //-----------------------------------------------------------------------------
242
243 bool CholmodSolver::
solve(std::vector<double> & _x0,std::vector<double> & _b)244 solve ( std::vector<double>& _x0, std::vector<double>& _b)
245 {
246 return solve( &(_x0[0]), &(_b[0]));
247 }
248
249 //-----------------------------------------------------------------------------
250
251 bool& CholmodSolver::
show_timings()252 show_timings()
253 {
254 return show_timings_;
255 }
256
257
258 }
259
260 //=============================================================================
261 #endif // COMISO_SUITESPARSE_AVAILABLE
262 //=============================================================================
263