1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbwrap/mschwrap.cc,v 1.34 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 #ifdef USE_MESCHACH
35
36 #include <iostream>
37 #include <iomanip>
38
39 #include "mschwrap.h"
40
41 /* MeschachVectorHandler - begin */
42
MeschachVectorHandler(integer iSize)43 MeschachVectorHandler::MeschachVectorHandler(integer iSize)
44 : pv(VNULL),
45 pdVecm1(NULL)
46 {
47 /* Note: MeschachVectorHandler owns its workspace memory */
48 if (iSize > 0) {
49 pv = v_get(iSize);
50 if (pv == VNULL) {
51 silent_cerr("out of memory?" << std::endl);
52 throw ErrMemory(MBDYN_EXCEPT_ARGS);
53 }
54 pdVecm1 = pv->ve - 1;
55 }
56 }
57
~MeschachVectorHandler(void)58 MeschachVectorHandler::~MeschachVectorHandler(void)
59 {
60 /* Note: MeschachVectorHandler owns its workspace memory */
61 if (pv != VNULL) {
62 if (v_free(pv) != 0) {
63 /* FIXME: hanlde error */
64 }
65 }
66 }
67
68 #ifdef DEBUG
69 /* Usata per il debug */
70 void
IsValid(void) const71 MeschachVectorHandler::IsValid(void) const
72 {
73 ASSERT(pv != VNULL);
74 ASSERT(pv->max_dim > 0);
75 ASSERT(pv->dim > 0);
76 ASSERT(pv->max_dim >= pv->dim);
77 ASSERT(pv->ve != NULL);
78 ASSERT(pdVecm1 != NULL);
79 ASSERT(pdVecm1 == pv->ve - 1);
80 }
81 #endif /* DEBUG */
82
83 void
Resize(integer iNewSize)84 MeschachVectorHandler::Resize(integer iNewSize)
85 {
86 #ifdef DEBUG
87 IsValid();
88 #endif /* DEBUG */
89
90 /* Note: MeschachVectorHandler owns its workspace memory */
91 VEC* p = v_resize(pv, iNewSize);
92 if (p == VNULL) {
93 silent_cerr("out of memory?" << std::endl);
94 throw ErrMemory(MBDYN_EXCEPT_ARGS);
95 }
96 pv = p;
97 pdVecm1 = pv->ve - 1;
98 }
99
100 void
Reset(void)101 MeschachVectorHandler::Reset(void)
102 {
103 #ifdef DEBUG
104 IsValid();
105 #endif /* DEBUG */
106 v_zero(pv);
107 }
108
109 /* MeschachVectorHandler - end */
110
111 /* MeschachSparseMatrixHandler - begin */
112
MeschachSparseMatrixHandler(integer m,integer n,integer maxlen)113 MeschachSparseMatrixHandler::MeschachSparseMatrixHandler(integer m,
114 integer n,
115 integer maxlen)
116 : mat(SMNULL)
117 {
118 if (m > 0 && n > 0) {
119 Create(m, n, maxlen);
120 }
121 }
122
~MeschachSparseMatrixHandler(void)123 MeschachSparseMatrixHandler::~MeschachSparseMatrixHandler(void)
124 {
125 if (mat != SMNULL) {
126
127 /*
128 * Note: MeschachSparseMatrixHandler owns its workspace memory
129 */
130 if (sp_free(mat)) {
131 /* FIXME: handle error */
132 }
133 }
134 }
135
136 /* costruisce la matrice */
137 void
Create(integer m,integer n,integer maxlen)138 MeschachSparseMatrixHandler::Create(integer m,
139 integer n,
140 integer maxlen)
141 {
142 /*
143 * Note: MeschachSparseMatrixHandler owns its workspace memory
144 */
145 if (mat != SMNULL) {
146 mat = sp_resize(mat, m, n);
147 } else {
148 /* FIXME: in case maxlen == 0, use n */
149 mat = sp_get(m, n, maxlen ? maxlen : n);
150 }
151 }
152
153 #ifdef DEBUG
154 void
IsValid(void) const155 MeschachSparseMatrixHandler::IsValid(void) const
156 {
157 ASSERT(mat != SMNULL);
158 }
159 #endif /* DEBUG */
160
161 void
Reset(void)162 MeschachSparseMatrixHandler::Reset(void)
163 {
164 sp_zero(mat);
165 }
166
167 /* MeschachSparseMatrixHandler - end */
168
169 /* MeschachSparseSolutionManager - begin */
170
171 void
Create(unsigned integer iSize,unsigned integer iMaxSize)172 MeschachSparseSolutionManager::Create(unsigned integer iSize,
173 unsigned integer iMaxSize)
174 {
175 if (prhs == NULL) {
176 SAFENEWWITHCONSTRUCTOR(prhs,
177 MeschachVectorHandler,
178 MeschachVectorHandler(iSize));
179 } else {
180 prhs->Resize(iSize);
181 }
182
183 if (pivot == PNULL || pivot->size < iSize) {
184 PERM* p = px_resize(pivot, iSize);
185 if (p == PNULL) {
186 silent_cerr("out of memory?" << std::endl);
187 throw ErrMemory(MBDYN_EXCEPT_ARGS);
188 }
189 pivot = p;
190 }
191
192 if (pmh != NULL
193 && (pmh->iGetNumRows() < (integer)iSize
194 || pmh->iGetNumCols() < (integer)iSize)) {
195 SAFEDELETE(pmh);
196 }
197
198 if (pmh == NULL) {
199 SAFENEWWITHCONSTRUCTOR(pmh,
200 MeschachSparseMatrixHandler,
201 MeschachSparseMatrixHandler(iSize,
202 iSize,
203 iMaxSize));
204 }
205 }
206
207 void
Factor(void)208 MeschachSparseSolutionManager::Factor(void)
209 {
210 #ifdef DEBUG
211 IsValid();
212 #endif /* DEBUG */
213
214 spLUfactor(pmh->pGetMAT(), pivot, alpha);
215 fStatus = FACTORED;
216 }
217
MeschachSparseSolutionManager(integer iSize,integer iMaxSize,const doublereal & a)218 MeschachSparseSolutionManager::MeschachSparseSolutionManager(integer iSize,
219 integer iMaxSize,
220 const doublereal& a)
221 : prhs(NULL), pivot(PNULL), pmh(NULL), fStatus(RESET), alpha (a)
222 {
223 Create(iSize, iMaxSize);
224 MatrReset();
225 }
226
~MeschachSparseSolutionManager(void)227 MeschachSparseSolutionManager::~MeschachSparseSolutionManager(void)
228 {
229 #ifdef DEBUG
230 IsValid();
231 #endif /* DEBUG */
232
233 if (prhs != NULL) {
234 SAFEDELETE(prhs);
235 }
236
237 if (pivot != PNULL) {
238 px_free(pivot);
239 }
240
241 if (pmh != NULL) {
242 SAFEDELETE(pmh);
243 }
244 }
245
246 #ifdef DEBUG
247 void
IsValid(void) const248 MeschachSparseSolutionManager::IsValid(void) const
249 {
250 ASSERT(prhs != NULL);
251 prhs->IsValid();
252 ASSERT(pivot != PNULL);
253 ASSERT(pmh != NULL);
254 pmh->IsValid();
255 ASSERT(prhs->iGetSize() >= pmh->iGetNumCols());
256 ASSERT(pivot->size >= pmh->iGetNumCols());
257 ASSERT(pivot->size >= pmh->iGetNumRows());
258 }
259 #endif /* DEBUG */
260
261 void
MatrReset(void)262 MeschachSparseSolutionManager::MatrReset(void)
263 {
264 #ifdef DEBUG
265 IsValid();
266 #endif /* DEBUG */
267
268 fStatus = RESET;
269 /* FIXME: TOTALLY UNTESTED */
270 pLS->Reset();
271 }
272
273 void
Solve(void)274 MeschachSparseSolutionManager::Solve(void)
275 {
276 #ifdef DEBUG
277 IsValid();
278 #endif /* DEBUG */
279
280 if (fStatus == RESET) {
281 /* Factor() is in charge of switching fStatus to FACTORED */
282 Factor();
283 }
284
285 spLUsolve(pmh->pGetMAT(),
286 pivot,
287 prhs->pGetMeschachVEC(),
288 prhs->pGetMeschachVEC());
289 }
290
291 /* Rende disponibile l'handler per la matrice */
292 MatrixHandler *
pMatHdl(void) const293 MeschachSparseSolutionManager::pMatHdl(void) const
294 {
295 #ifdef DEBUG
296 IsValid();
297 #endif /* DEBUG */
298
299 return pmh;
300 }
301
302 /* Rende disponibile l'handler per il termine noto */
303 VectorHandler *
pResHdl(void) const304 MeschachSparseSolutionManager::pResHdl(void) const
305 {
306 #ifdef DEBUG
307 IsValid();
308 #endif /* DEBUG */
309
310 return prhs;
311 }
312
313 /*
314 * Rende disponibile l'handler per la soluzione (e' lo stesso
315 * del termine noto, ma concettualmente sono separati)
316 */
317 VectorHandler *
pSolHdl(void) const318 MeschachSparseSolutionManager::pSolHdl(void) const
319 {
320 #ifdef DEBUG
321 IsValid();
322 #endif /* DEBUG */
323
324 return prhs;
325 }
326
327 std::ostream&
operator <<(std::ostream & out,const MeschachSparseMatrixHandler & MH)328 operator << (std::ostream& out, const MeschachSparseMatrixHandler& MH)
329 {
330 SPMAT* p = MH.pGetMAT();
331 for (integer i = 0; i < p->m; i++) {
332 for (integer j = 0; j < p->n; j++) {
333 silent_cout(std::setw(16) << sp_get_val(p, i, j));
334 }
335 silent_cout(std::endl);
336 }
337
338 return out;
339 }
340
341 #endif /* USE_MESCHACH */
342
343