MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
matmultest.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbmath/matmultest.cc,v 1.19 2017/01/12 14:43:53 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 #include <iostream>
35 #include <iomanip>
36 
37 #include "fullmh.h"
38 #include "spmapmh.h"
39 #include "ccmh.h"
40 #include "dirccmh.h"
41 #include "naivemh.h"
42 
43 static doublereal mat[5][5] = {
44  { 11., 0., 13., 0., 15. },
45  { 0., 22., 0., 24., 0. },
46  { 31., 0., 33., 0., 35. },
47  { 0., 42., 0., 44., 0. },
48  { 51., 0., 53., 0., 55. }
49 };
50 
51 static int
53 {
54  for (int r = 0; r < 5; r++) {
55  for (int c = 0; c < 5; c++) {
56  if (mat[r][c] != mh(r + 1, c + 1)) {
57  return -1;
58  }
59  }
60  }
61 
62  return 0;
63 }
64 
65 static int
67 {
68  for (int r = 0; r < 5; r++) {
69  for (int c = 0; c < 5; c++) {
70  if (mat[c][r] != mh(r + 1, c + 1)) {
71  return -1;
72  }
73  }
74  }
75 
76  return 0;
77 }
78 
79 static int
80 check_vec(VectorHandler& vh, unsigned c)
81 {
82  c--;
83  for (int r = 0; r < 5; r++) {
84  if (mat[r][c] != vh(r + 1)) {
85  return -1;
86  }
87  }
88 
89  return 0;
90 }
91 
92 static int
94 {
95  r--;
96  for (int c = 0; c < 5; c++) {
97  if (mat[r][c] != vh(c + 1)) {
98  return -1;
99  }
100  }
101 
102  return 0;
103 }
104 
105 int
106 main(void)
107 {
108  std::vector<integer> perm(5), invperm(5);
109  perm[0] = 4;
110  perm[1] = 3;
111  perm[2] = 2;
112  perm[3] = 1;
113  perm[4] = 0;
114  for (int i = 0; i < 5; i++) {
115  invperm[perm[i]] = i;
116  }
117 
118  SpMapMatrixHandler spm(5, 5);
119  NaiveMatrixHandler nm(5);
120  NaivePermMatrixHandler npm(5, perm, invperm);
121 
122  for (int r = 0; r < 5; r++) {
123  for (int c = 0; c < 5; c++) {
124  if (mat[r][c] != 0.) {
125  spm(r + 1, c + 1) = mat[r][c];
126  nm(r + 1, c + 1) = mat[r][c];
127  npm(r + 1, c + 1) = mat[r][c];
128  }
129  }
130  }
131 
132  std::cout << "matrix in sparse form: " << std::endl
133  << spm << std::endl;
134 
135  std::cout << "matrix in naive form: " << std::endl
136  << nm << std::endl;
137 
138  std::cout << "matrix in naive permuted form: " << std::endl
139  << npm << std::endl;
140 
141  std::vector<doublereal> Ax0;
142  std::vector<integer> Ai0, Ap0;
143  spm.MakeCompressedColumnForm(Ax0, Ai0, Ap0, 0);
144 
145  CColMatrixHandler<0> ccm0(Ax0, Ai0, Ap0);
146 
147  std::cout << "matrix in cc<0> form: " << std::endl
148  << ccm0 << std::endl;
149 
150  std::cout << "matrix in cc<0> form again: " << std::endl;
151  const CColMatrixHandler<0>& const_ccm0 = ccm0;
152  for (int ir = 1; ir <= 5; ir++) {
153  for (int ic = 1; ic <= 5; ic++) {
154  std::cout << std::setw(16) << const_ccm0(ir, ic);
155  }
156  std::cout << std::endl;
157  }
158 
159  std::vector<doublereal> Ax1;
160  std::vector<integer> Ai1, Ap1;
161  spm.MakeCompressedColumnForm(Ax1, Ai1, Ap1, 1);
162 
163  CColMatrixHandler<1> ccm1(Ax1, Ai1, Ap1);
164 
165  std::cout << "matrix in cc<1> form: " << std::endl
166  << ccm1 << std::endl;
167 
168  std::cout << "matrix in cc<1> form again: " << std::endl;
169  const CColMatrixHandler<1>& const_ccm1 = ccm1;
170  for (int ir = 1; ir <= 5; ir++) {
171  for (int ic = 1; ic <= 5; ic++) {
172  std::cout << std::setw(16) << const_ccm1(ir, ic);
173  }
174  std::cout << std::endl;
175  }
176 
177  DirCColMatrixHandler<0> dirm0(Ax0, Ai0, Ap0);
178  std::cout << "matrix in dir<0> form: " << std::endl
179  << dirm0 << std::endl;
180 
181  std::cout << "matrix in dir<0> form again: " << std::endl;
182  const DirCColMatrixHandler<0>& const_dirm0 = dirm0;
183  for (int ir = 1; ir <= 5; ir++) {
184  for (int ic = 1; ic <= 5; ic++) {
185  std::cout << std::setw(16) << const_dirm0(ir, ic);
186  }
187  std::cout << std::endl;
188  }
189 
190  DirCColMatrixHandler<1> dirm1(Ax1, Ai1, Ap1);
191  std::cout << "matrix in dir<1> form: " << std::endl
192  << dirm1 << std::endl;
193 
194  std::cout << "matrix in dir<1> form again: " << std::endl;
195  const DirCColMatrixHandler<1>& const_dirm1 = dirm1;
196  for (int ir = 1; ir <= 5; ir++) {
197  for (int ic = 1; ic <= 5; ic++) {
198  std::cout << std::setw(16) << const_dirm1(ir, ic);
199  }
200  std::cout << std::endl;
201  }
202 
203  MyVectorHandler v(5), out(5);
204  v.Reset();
205  for (int i = 1; i <= 5; i++) {
206  v(i) = 1.;
207 
208  spm.MatVecMul(out, v);
209  std::cout << "sp*v(" << i << ")=" << std::endl
210  << out << std::endl;
211  if (check_vec(out, i)) {
212  std::cerr << "*** failed!" << std::endl;
213  }
214 
215  spm.MatTVecMul(out, v);
216  std::cout << "sp^T*v(" << i << ")=" << std::endl
217  << out << std::endl;
218  if (check_vec_transpose(out, i)) {
219  std::cerr << "*** failed!" << std::endl;
220  }
221 
222  nm.MatVecMul(out, v);
223  std::cout << "naive*v(" << i << ")=" << std::endl
224  << out << std::endl;
225  if (check_vec(out, i)) {
226  std::cerr << "*** failed!" << std::endl;
227  }
228 
229  nm.MatTVecMul(out, v);
230  std::cout << "naive^T*v(" << i << ")=" << std::endl
231  << out << std::endl;
232  if (check_vec_transpose(out, i)) {
233  std::cerr << "*** failed!" << std::endl;
234  }
235 
236  npm.MatVecMul(out, v);
237  std::cout << "naiveperm*v(" << i << ")=" << std::endl
238  << out << std::endl;
239  if (check_vec(out, i)) {
240  std::cerr << "*** failed!" << std::endl;
241  }
242 
243  npm.MatTVecMul(out, v);
244  std::cout << "naiveperm^T*v(" << i << ")=" << std::endl
245  << out << std::endl;
246  if (check_vec_transpose(out, i)) {
247  std::cerr << "*** failed!" << std::endl;
248  }
249 
250  ccm0.MatVecMul(out, v);
251  std::cout << "cc<0>*v(" << i << ")=" << std::endl
252  << out << std::endl;
253  if (check_vec(out, i)) {
254  std::cerr << "*** failed!" << std::endl;
255  }
256 
257  ccm0.MatTVecMul(out, v);
258  std::cout << "cc<0>^T*v(" << i << ")=" << std::endl
259  << out << std::endl;
260  if (check_vec_transpose(out, i)) {
261  std::cerr << "*** failed!" << std::endl;
262  }
263 
264  ccm1.MatVecMul(out, v);
265  std::cout << "cc<1>*v(" << i << ")=" << std::endl
266  << out << std::endl;
267  if (check_vec(out, i)) {
268  std::cerr << "*** failed!" << std::endl;
269  }
270 
271  ccm1.MatTVecMul(out, v);
272  std::cout << "cc<1>^T*v(" << i << ")=" << std::endl
273  << out << std::endl;
274  if (check_vec_transpose(out, i)) {
275  std::cerr << "*** failed!" << std::endl;
276  }
277 
278  dirm0.MatVecMul(out, v);
279  std::cout << "dir<0>*v(" << i << ")=" << std::endl
280  << out << std::endl;
281  if (check_vec(out, i)) {
282  std::cerr << "*** failed!" << std::endl;
283  }
284 
285  dirm0.MatTVecMul(out, v);
286  std::cout << "dir<0>^T*v(" << i << ")=" << std::endl
287  << out << std::endl;
288  if (check_vec_transpose(out, i)) {
289  std::cerr << "*** failed!" << std::endl;
290  }
291 
292  dirm1.MatVecMul(out, v);
293  std::cout << "dir<1>*v(" << i << ")=" << std::endl
294  << out << std::endl;
295  if (check_vec(out, i)) {
296  std::cerr << "*** failed!" << std::endl;
297  }
298 
299  dirm1.MatTVecMul(out, v);
300  std::cout << "dir<1>^T*v(" << i << ")=" << std::endl
301  << out << std::endl;
302  if (check_vec_transpose(out, i)) {
303  std::cerr << "*** failed!" << std::endl;
304  }
305 
306  v(i) = 0.;
307  }
308 
309  FullMatrixHandler fmin(5, 5), fmout(5, 5);
310  fmin.Reset();
311 
312  fmin(1, 1) = 1.;
313  fmin(2, 2) = 1.;
314  fmin(3, 3) = 1.;
315  fmin(4, 4) = 1.;
316  fmin(5, 5) = 1.;
317 
318  spm.MatMatMul(fmout, fmin);
319  std::cout << "sp*eye=" << std::endl
320  << fmout << std::endl;
321  if (check_mat(fmout)) {
322  std::cerr << "*** failed!" << std::endl;
323  }
324 
325  spm.MatTMatMul(fmout, fmin);
326  std::cout << "sp^T*eye=" << std::endl
327  << fmout << std::endl;
328  if (check_mat_transpose(fmout)) {
329  std::cerr << "*** failed!" << std::endl;
330  }
331 
332  nm.MatMatMul(fmout, fmin);
333  std::cout << "naive*eye=" << std::endl
334  << fmout << std::endl;
335  if (check_mat(fmout)) {
336  std::cerr << "*** failed!" << std::endl;
337  }
338 
339  nm.MatTMatMul(fmout, fmin);
340  std::cout << "naive^T*eye=" << std::endl
341  << fmout << std::endl;
342  if (check_mat_transpose(fmout)) {
343  std::cerr << "*** failed!" << std::endl;
344  }
345 
346  npm.MatMatMul(fmout, fmin);
347  std::cout << "naiveperm*eye=" << std::endl
348  << fmout << std::endl;
349  if (check_mat(fmout)) {
350  std::cerr << "*** failed!" << std::endl;
351  }
352 
353  npm.MatTMatMul(fmout, fmin);
354  std::cout << "naiveperm^T*eye=" << std::endl
355  << fmout << std::endl;
356  if (check_mat_transpose(fmout)) {
357  std::cerr << "*** failed!" << std::endl;
358  }
359 
360  ccm0.MatMatMul(fmout, fmin);
361  std::cout << "cc<0>*eye=" << std::endl
362  << fmout << std::endl;
363  if (check_mat(fmout)) {
364  std::cerr << "*** failed!" << std::endl;
365  }
366 
367  ccm0.MatTMatMul(fmout, fmin);
368  std::cout << "cc<0>^T*eye=" << std::endl
369  << fmout << std::endl;
370  if (check_mat_transpose(fmout)) {
371  std::cerr << "*** failed!" << std::endl;
372  }
373 
374  ccm1.MatMatMul(fmout, fmin);
375  std::cout << "cc<1>*eye=" << std::endl
376  << fmout << std::endl;
377  if (check_mat(fmout)) {
378  std::cerr << "*** failed!" << std::endl;
379  }
380 
381  ccm1.MatTMatMul(fmout, fmin);
382  std::cout << "cc<1>^T*eye=" << std::endl
383  << fmout << std::endl;
384  if (check_mat_transpose(fmout)) {
385  std::cerr << "*** failed!" << std::endl;
386  }
387 
388  dirm0.MatMatMul(fmout, fmin);
389  std::cout << "dir<0>*eye=" << std::endl
390  << fmout << std::endl;
391  if (check_mat(fmout)) {
392  std::cerr << "*** failed!" << std::endl;
393  }
394 
395  dirm0.MatTMatMul(fmout, fmin);
396  std::cout << "dir<0>^T*eye=" << std::endl
397  << fmout << std::endl;
398  if (check_mat_transpose(fmout)) {
399  std::cerr << "*** failed!" << std::endl;
400  }
401 
402  dirm1.MatMatMul(fmout, fmin);
403  std::cout << "dir<1>*eye=" << std::endl
404  << fmout << std::endl;
405  if (check_mat(fmout)) {
406  std::cerr << "*** failed!" << std::endl;
407  }
408 
409  dirm1.MatTMatMul(fmout, fmin);
410  std::cout << "dir<1>^T*eye=" << std::endl
411  << fmout << std::endl;
412  if (check_mat_transpose(fmout)) {
413  std::cerr << "*** failed!" << std::endl;
414  }
415 
416  return 0;
417 }
418 
static int check_vec_transpose(VectorHandler &vh, unsigned r)
Definition: matmultest.cc:93
integer MakeCompressedColumnForm(doublereal *const Ax, integer *const Ai, integer *const Ap, int offset=0) const
Definition: spmapmh.cc:80
static int check_mat(MatrixHandler &mh)
Definition: matmultest.cc:52
static doublereal mat[5][5]
Definition: matmultest.cc:43
int main(void)
Definition: matmultest.cc:106
static std::stack< cleanup * > c
Definition: cleanup.cc:59
static int check_mat_transpose(MatrixHandler &mh)
Definition: matmultest.cc:66
virtual VectorHandler & MatTVecMul(VectorHandler &out, const VectorHandler &in) const
Definition: mh.cc:341
double doublereal
Definition: colamd.c:52
virtual MatrixHandler & MatMatMul(MatrixHandler &out, const MatrixHandler &in) const
Definition: mh.cc:232
virtual MatrixHandler & MatTMatMul(MatrixHandler &out, const MatrixHandler &in) const
Definition: mh.cc:241
virtual VectorHandler & MatVecMul(VectorHandler &out, const VectorHandler &in) const
Definition: mh.cc:332
static int check_vec(VectorHandler &vh, unsigned c)
Definition: matmultest.cc:80