MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
spmh.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbmath/spmh.cc,v 1.33 2017/01/12 14:43:54 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 2003-2017
7  *
8  * This code is a partial merge of HmFe and MBDyn.
9  *
10  * Pierangelo Masarati <masarati@aero.polimi.it>
11  * Paolo Mantegazza <mantegazza@aero.polimi.it>
12  *
13  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
14  * via La Masa, 34 - 20156 Milano, Italy
15  * http://www.aero.polimi.it
16  *
17  * Changing this copyright notice is forbidden.
18  *
19  * This program is free software; you can redistribute it and/or modify
20  * it under the terms of the GNU General Public License as published by
21  * the Free Software Foundation (version 2 of the License).
22  *
23  *
24  * This program is distributed in the hope that it will be useful,
25  * but WITHOUT ANY WARRANTY; without even the implied warranty of
26  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27  * GNU General Public License for more details.
28  *
29  * You should have received a copy of the GNU General Public License
30  * along with this program; if not, write to the Free Software
31  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
32  */
33 
34 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
35 
36 #include "spmh.h"
37 
39 : NRows(n), NCols(nn == 0 ? n : nn), NZ(0)
40 {
41  NO_OP;
42 }
43 
45 {
46  NO_OP;
47 }
48 
50  const integer &nn,
51  std::vector<doublereal>&x,
52  const std::vector<integer>& i,
53  const std::vector<integer>& p)
54 : SparseMatrixHandler(n, nn),
55 bMatDuplicate(false),
56 Ax(x),
57 Ai(i),
58 Ap(p)
59 {
60  NZ = Ax.size();
61 }
62 
64 {
65  if (bMatDuplicate) {
66  delete &Ax;
67  }
68 }
69 
70 /* used to sum CC matrices with identical indices */
71 void
73 {
74  /* FIXME: put in stl-ish form;
75  * see if we can use something from optimized blas,
76  * e.g. ATLAS, goto or so... */
77 
78  /* checks - uncomment to enable */
79 #ifdef DEBUG
80  ASSERT(Ax.size() == m.Ax.size());
81  ASSERT(Ai.size() == m.Ai.size());
82  ASSERT(Ap.size() == m.Ap.size());
83  for (std::vector<doublereal>::size_type i = 0; i < Ai.size(); i++) {
84  if (Ai[i] != m.Ai[i]) {
85  silent_cerr("AddUnchecked: Ai[" << i << "] differs" << std::endl);
87  }
88  }
89 
90  for (std::vector<doublereal>::size_type i = 0; i < Ap.size(); i++) {
91  if (Ap[i] != m.Ap[i]) {
92  silent_cerr("AddUnchecked: Ap[" << i << "] differs" << std::endl);
94  }
95  }
96 #endif /* DEBUG */
97 
98  doublereal *d = &Ax[0], *s = &m.Ax[0];
99  std::vector<doublereal>::size_type n = Ax.size();
100  for (std::vector<doublereal>::size_type i = 0; i < n; i++) {
101  d[i] += s[i];
102  }
103 }
104 
105 void
107 {
108  std::fill(Ax.begin(), Ax.end(), 0.);
109 }
110 
111 integer
113  integer *const Ai, integer *const Ap,
114  int offset) const
115 {
117  return Nz();
118 }
119 
120 integer
122  std::vector<integer>& Ai, std::vector<integer>& Ap,
123  int offset) const
124 {
126  return Nz();
127 }
128 
129 integer
131  integer *const Arow, integer *const Acol,
132  integer *const AcolSt, int offset) const
133 {
135  return Nz();
136 }
137 
138 integer
139 CompactSparseMatrixHandler::MakeIndexForm(std::vector<doublereal>& rAx,
140  std::vector<integer>& Arow, std::vector<integer>& Acol,
141  std::vector<integer>& AcolSt, int offset) const
142 {
144  return Nz();
145 }
146 
147 template <int off>
149  const integer &n,
150  const integer &nn,
151  std::vector<doublereal>&x,
152  const std::vector<integer>& i,
153  const std::vector<integer>& p)
154 : CompactSparseMatrixHandler(n, nn, x, i, p),
155 m_end(*this, true)
156 {
157  NO_OP;
158 }
159 
160 template <int off>
162 {
163  NO_OP;
164 }
165 
166 /* Prodotto Matrice per Matrice */
167 template <int off>
170  void (MatrixHandler::*op)(integer iRow, integer iCol, const doublereal& dCoef),
171  MatrixHandler& out, const MatrixHandler& in) const
172 {
173  ASSERT(in.iGetNumRows() == NCols);
174  ASSERT(out.iGetNumRows() == NRows);
175  ASSERT(in.iGetNumCols() == out.iGetNumCols());
176 
177  // NOTE: out must be zeroed by caller
178 
179  integer ncols_in = in.iGetNumCols();
180 
181  for (integer col_idx = 1; col_idx <= NCols; col_idx++) {
182  integer re = Ap[col_idx] - off;
183  integer ri = Ap[col_idx - 1] - off;
184  for ( ; ri < re; ri++) {
185  integer row_idx = Ai[ri] - off + 1;
186  const doublereal& d = Ax[ri];
187  for (integer col_in = 1; col_in <= ncols_in; col_in++) {
188  (out.*op)(row_idx, col_in, d*in(col_idx, col_in));
189  }
190  }
191  }
192 
193  return out;
194 }
195 
196 /* Prodotto Matrice trasposta per Matrice */
197 template <int off>
200  void (MatrixHandler::*op)(integer iRow, integer iCol, const doublereal& dCoef),
201  MatrixHandler& out, const MatrixHandler& in) const
202 {
203  ASSERT(in.iGetNumRows() == NRows);
204  ASSERT(out.iGetNumRows() == NCols);
205  ASSERT(in.iGetNumCols() == out.iGetNumCols());
206  // NOTE: out must be zeroed by caller
207 
208  integer ncols_in = in.iGetNumCols();
209 
210  for (integer col_idx = 1; col_idx <= NCols; col_idx++) {
211  integer re = Ap[col_idx] - off;
212  integer ri = Ap[col_idx - 1] - off;
213  for ( ; ri < re; ri++) {
214  integer row_idx = Ai[ri] - off + 1;
215  const doublereal& d = Ax[ri];
216  for (integer col_in = 1; col_in <= ncols_in; col_in++) {
217  (out.*op)(col_idx, col_in, d*in(row_idx, col_in));
218  }
219  }
220  }
221 
222  return out;
223 }
224 
225 /* Prodotto Matrice per Vettore */
226 template <int off>
229  void (VectorHandler::*op)(integer iRow, const doublereal& dCoef),
230  VectorHandler& out, const VectorHandler& in) const
231 {
232  ASSERT(in.iGetSize() == NCols);
233  ASSERT(out.iGetSize() == NRows);
234 
235  // NOTE: out must be zeroed by caller
236 
237  for (integer col_idx = 1; col_idx <= NCols; col_idx++) {
238  integer re = Ap[col_idx] - off;
239  integer ri = Ap[col_idx - 1] - off;
240  for ( ; ri < re; ri++) {
241  integer row_idx = Ai[ri] - off + 1;
242  const doublereal& d = Ax[ri];
243  (out.*op)(row_idx, d*in(col_idx));
244  }
245  }
246 
247  return out;
248 }
249 
250 /* Prodotto Matrice trasposta per Vettore */
251 template <int off>
254  void (VectorHandler::*op)(integer iRow,
255  const doublereal& dCoef),
256  VectorHandler& out, const VectorHandler& in) const
257 {
258  ASSERT(in.iGetSize() == NRows);
259  ASSERT(out.iGetSize() == NCols);
260 
261  // NOTE: out must be zeroed by caller
262 
263  for (integer col_idx = 1; col_idx <= NCols; col_idx++) {
264  integer re = Ap[col_idx] - off;
265  integer ri = Ap[col_idx - 1] - off;
266  for ( ; ri < re; ri++) {
267  integer row_idx = Ai[ri] - off + 1;
268  const doublereal& d = Ax[ri];
269  (out.*op)(col_idx, d*in(row_idx));
270  }
271  }
272 
273  return out;
274 }
275 
276 template <int off>
277 void
279 {
280  if (is_end) {
281  elem.iRow = m.iGetNumRows();
282  elem.iCol = m.iGetNumCols();
283 
284  } else {
285  elem.iCol = 0;
286  while (m.Ap[elem.iCol + 1] - off == m.Ap[elem.iCol] - off) {
287  if (++elem.iCol == m.iGetNumCols()) {
288  elem.iRow = m.iGetNumRows();
289  return;
290  }
291  }
292  i_idx = m.Ap[elem.iCol] - off;
293  elem.iRow = m.Ai[i_idx] - off;
294  elem.dCoef = m.Ax[i_idx];
295  }
296 }
297 
298 template <int off>
300 : m(m)
301 {
302  reset(is_end);
303 }
304 
305 template <int off>
307 {
308  NO_OP;
309 }
310 
311 template <int off>
314 {
315  i_idx++;
316  while (i_idx == m.Ap[elem.iCol + 1] - off) {
317  if (++elem.iCol == m.iGetNumCols()) {
318  elem.iRow = m.iGetNumRows();
319  return *this;
320  }
321  i_idx = m.Ap[elem.iCol] - off;
322  }
323  elem.iRow = m.Ai[i_idx] - off;
324  elem.dCoef = m.Ax[i_idx];
325 
326  return *this;
327 }
328 
329 template <int off>
332 {
333  return &elem;
334 }
335 
336 template <int off>
339 {
340  return elem;
341 }
342 
343 template <int off>
344 bool
346 {
347  return elem == op.elem;
348 }
349 
350 template <int off>
351 bool
353 {
354  return elem != op.elem;
355 }
356 
357 template class CompactSparseMatrixHandler_tpl<0>;
358 template class CompactSparseMatrixHandler_tpl<1>;
359 
const_iterator(const CompactSparseMatrixHandler_tpl< off > &m, bool is_end=false)
Definition: spmh.cc:299
virtual VectorHandler & MatTVecMul_base(void(VectorHandler::*op)(integer iRow, const doublereal &dCoef), VectorHandler &out, const VectorHandler &in) const
Definition: spmh.cc:253
virtual integer iGetNumCols(void) const =0
virtual VectorHandler & MatVecMul_base(void(VectorHandler::*op)(integer iRow, const doublereal &dCoef), VectorHandler &out, const VectorHandler &in) const
Definition: spmh.cc:228
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
const CompactSparseMatrixHandler_tpl< off > & m
Definition: spmh.h:220
void AddUnchecked(const CompactSparseMatrixHandler &m)
Definition: spmh.cc:72
const CompactSparseMatrixHandler_tpl< off >::const_iterator & operator++(void) const
Definition: spmh.cc:313
virtual integer MakeCompressedColumnForm(doublereal *const Ax, integer *const Ai, integer *const Ap, int offset=0) const
Definition: spmh.cc:112
const std::vector< integer > & Ap
Definition: spmh.h:154
SparseMatrixHandler(const integer &n, const integer &nn=0)
Definition: spmh.cc:38
const integer Nz() const
Definition: spmh.h:104
#define NO_OP
Definition: myassert.h:74
virtual integer iGetSize(void) const =0
const std::vector< integer > & Ai
Definition: spmh.h:153
bool operator!=(const CompactSparseMatrixHandler_tpl< off >::const_iterator &op) const
Definition: spmh.cc:352
Definition: mbdyn.h:76
std::vector< doublereal > & Ax
Definition: spmh.h:152
virtual ~SparseMatrixHandler(void)
Definition: spmh.cc:44
integer NZ
Definition: spmh.h:47
virtual integer MakeIndexForm(doublereal *const Ax, integer *const Arow, integer *const Acol, integer *const AcolSt, int offset=0) const
Definition: spmh.cc:130
Definition: mbdyn.h:77
const SparseMatrixHandler::SparseMatrixElement * operator->(void)
Definition: spmh.cc:331
virtual ~CompactSparseMatrixHandler_tpl(void)
Definition: spmh.cc:161
SparseMatrixHandler::SparseMatrixElement elem
Definition: spmh.h:222
#define ASSERT(expression)
Definition: colamd.c:977
bool operator==(const CompactSparseMatrixHandler_tpl< off >::const_iterator &op) const
Definition: spmh.cc:345
const SparseMatrixHandler::SparseMatrixElement & operator*(void)
Definition: spmh.cc:338
MatrixHandler & MatMatMul_base(void(MatrixHandler::*op)(integer iRow, integer iCol, const doublereal &dCoef), MatrixHandler &out, const MatrixHandler &in) const
Definition: spmh.cc:169
const_iterator m_end
Definition: spmapmh.h:113
CompactSparseMatrixHandler(const integer &n, const integer &nn, std::vector< doublereal > &x, const std::vector< integer > &i, const std::vector< integer > &p)
Definition: spmh.cc:49
integer NRows
Definition: spmh.h:45
CompactSparseMatrixHandler_tpl(const integer &n, const integer &nn, std::vector< doublereal > &x, const std::vector< integer > &i, const std::vector< integer > &p)
Definition: spmh.cc:148
double doublereal
Definition: colamd.c:52
MatrixHandler & MatTMatMul_base(void(MatrixHandler::*op)(integer iRow, integer iCol, const doublereal &dCoef), MatrixHandler &out, const MatrixHandler &in) const
Definition: spmh.cc:199
long int integer
Definition: colamd.c:51
integer NCols
Definition: spmh.h:46
virtual ~CompactSparseMatrixHandler()
Definition: spmh.cc:63
virtual integer iGetNumRows(void) const =0