1 /* -------------------------------------------------------------------------- *
2 * Simbody(tm): SimTKcommon *
3 * -------------------------------------------------------------------------- *
4 * This is part of the SimTK biosimulation toolkit originating from *
5 * Simbios, the NIH National Center for Physics-Based Simulation of *
6 * Biological Structures at Stanford, funded under the NIH Roadmap for *
7 * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. *
8 * *
9 * Portions copyright (c) 2005-12 Stanford University and the Authors. *
10 * Authors: Michael Sherman *
11 * Contributors: *
12 * *
13 * Licensed under the Apache License, Version 2.0 (the "License"); you may *
14 * not use this file except in compliance with the License. You may obtain a *
15 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
16 * *
17 * Unless required by applicable law or agreed to in writing, software *
18 * distributed under the License is distributed on an "AS IS" BASIS, *
19 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
20 * See the License for the specific language governing permissions and *
21 * limitations under the License. *
22 * -------------------------------------------------------------------------- */
23
24 #include "SimTKcommon/SmallMatrix.h"
25 #include "SimTKcommon/Testing.h"
26
27 #include <cstdio>
28 #include <iostream>
29 #include <iomanip>
30 using std::cout;
31 using std::endl;
32 using std::setw;
33
34 #include <complex>
35 using std::complex;
36 using std::sin;
37 using std::cos;
38
39 using namespace SimTK;
40
41 static void dummy();
42
43 //typedef Vec<3,Complex> CVec3;
44 //static CVec3 f(CVec3 v) {
45 // return CVec3();
46 //}
47
f(Complex x)48 static Complex f(Complex x) {
49 return std::sin(x);
50 }
51
52 // TODO: move lots of tests here and check the answers
testNegator()53 void testNegator() {
54 Vec3 vvv(1,2,3); Row3 rrr(1,2,3);
55 SimTK_TEST(vvv-vvv == Vec3(0)); // these are exact results
56 SimTK_TEST(-vvv-vvv == (-2*vvv));
57 SimTK_TEST(rrr-rrr == Row3(0));
58 SimTK_TEST(-rrr-rrr == (-2*rrr));
59 }
60
testElementwiseOps()61 void testElementwiseOps() {
62 Vec3 vvv(1,2,3); Vec3 www(7,9,2);
63 Row3 rrr(1,2,3); Row3 sss(5,4,10);
64
65 Vec3 mv=vvv.elementwiseMultiply(www);
66 Vec3 dv=vvv.elementwiseDivide(www);
67 SimTK_TEST_EQ(mv, Vec3(7,18,6));
68 SimTK_TEST_EQ(dv, Vec3(Real(1)/7, Real(2)/9, Real(3)/2));
69
70 Row3 mr=rrr.elementwiseMultiply(sss);
71 Row3 dr=rrr.elementwiseDivide(sss);
72 SimTK_TEST_EQ(mr, Row3(5,8,30));
73 SimTK_TEST_EQ(dr, Row3(Real(1)/5, Real(2)/4, Real(3)/10));
74
75 Mat22 mmm(1, 2,
76 3, 4);
77 Mat22 nnn(7, 9,
78 2, 3);
79 SymMat22 yyy(1,
80 2, 3);
81 SymMat22 zzz(5,
82 4, 10);
83
84 Mat22 mm=mmm.elementwiseMultiply(nnn);
85 Mat22 dm=mmm.elementwiseDivide(nnn);
86 SymMat22 my=yyy.elementwiseMultiply(zzz);
87 SymMat22 dy=yyy.elementwiseDivide(zzz);
88
89 SimTK_TEST_EQ(mm, Mat22(7,18,6,12));
90 SimTK_TEST_EQ(dm, Mat22(Real(1)/7,Real(2)/9,Real(3)/2,Real(4)/3));
91 SimTK_TEST_EQ(my, SymMat22(5,8,30));
92 SimTK_TEST_EQ(dy, SymMat22(Real(1)/5,Real(2)/4,Real(3)/10));
93 }
94
testSums()95 void testSums() {
96 Mat22 m(1,2,
97 3,4);
98 SimTK_TEST_EQ(m.colSum(), Row2(4,6));
99 SimTK_TEST_EQ(m.rowSum(), Vec2(3,7));
100 SimTK_TEST(m.sum() == m.colSum()); // should be exact
101
102 SymMat22 y(3, /*4*/
103 4, 5);
104 SimTK_TEST_EQ(y.colSum(), Row2(7,9)); // same for real, sym
105 SimTK_TEST_EQ(y.rowSum(), Vec2(7,9));
106 SimTK_TEST(y.sum() == y.colSum()); // should be exact
107
108 Mat22 sm(y); // create fully populated symmetric matrix
109 SimTK_TEST_EQ(sm.rowSum(), y.rowSum());
110 SimTK_TEST_EQ(sm.colSum(), y.colSum());
111
112 Mat<2,2,Complex> mc(1+2*I, 3+4*I,
113 5+6*I, 7+8*I);
114 typedef Row<2,Complex> CRow2;
115 typedef Vec<2,Complex> CVec2;
116 SimTK_TEST_EQ(mc.colSum(), CRow2(6+8*I, 10+12*I));
117 SimTK_TEST_EQ(mc.rowSum(), CVec2(4+6*I, 12+14*I));
118 SimTK_TEST(mc.sum() == mc.colSum()); // should be exact
119
120 // Row sum and col sum for Hermitian are conjugates; not the same.
121 SymMat<2,Complex> yc( 1, /*3+6*I*/
122 3-6*I, 4);
123 SimTK_TEST_EQ(yc.colSum(), CRow2(4-6*I, 7+6*I));
124 SimTK_TEST_EQ(yc.rowSum(), CVec2(4+6*I, 7-6*I));
125 SimTK_TEST(yc.sum() == yc.colSum()); // should be exact
126
127 Mat<2,2,Complex> smc(yc); // create fully populated symmetric matrix
128 SimTK_TEST_EQ(smc.rowSum(), yc.rowSum());
129 SimTK_TEST_EQ(smc.colSum(), yc.colSum());
130 }
131
testMiscellaneous()132 void testMiscellaneous()
133 {
134 cout << std::setprecision(16);
135 cout << "f(.3)=" << f(Complex(0.3)) << endl;
136 Real h = 1e-20;
137 cout << "f(.3 + i*h)/h=" << f(Complex(0.3,h)) / h << endl;
138
139 cout << CNT< Mat<2,3, Vec<2> > >::getNaN() << endl;
140
141 Mat<2,3, Vec<2, Mat<2,2,Complex> > > isThisNaN;
142 cout << "isThisNan? " << isThisNaN << endl;
143
144 const Complex mdc[] = {
145 Complex(1.,2.), Complex(3.,4.), Complex(5.,6.), Complex(7.,8.),
146 Complex(9.,10.), Complex(10.,11.), Complex(.1,.26), Complex(.3,.45),
147 Complex(.5,.64), Complex(.7,.83), Complex(.9,.102), Complex(.10,.111)
148 };
149
150 cout << "*** TEST COMPLEX DOT PRODUCT ***" << endl;
151 Vec<3,Complex> vdot(mdc), wdot(&mdc[3]);
152 Row<3,Complex> rdot(mdc), sdot(&mdc[3]);
153 cout << "v=" << vdot << " w=" << wdot << endl;
154 cout << "r=" << rdot << " s=" << sdot << endl;
155 cout << "v.normalize()=" << vdot.normalize() << endl;
156 cout << "r.normalize()=" << rdot.normalize() << endl;
157
158
159 cout << "--- dot() global function:dot(v,w), rw, vs, rs should be the same" << endl;
160 cout << "vw=" << dot(vdot,wdot) << " rw" << dot(rdot,wdot)
161 << " vs" << dot(vdot,sdot) << " rs" << dot(rdot,sdot) << endl;
162 cout << "--- dot operator* requires row*col meaning Hermitian transpose with sign changes" << endl;
163 cout << "vw=" << ~vdot*wdot << " rw" << rdot*wdot
164 << " vs" << ~vdot*~sdot << " rs" << rdot*~sdot << endl;
165
166 cout << endl << "*** TEST COMPLEX OUTER PRODUCT ***" << endl;
167 cout << "--- outer() global function:dot(v,w), rw, vs, rs should be the same" << endl;
168 cout << "vw=" << outer(vdot,wdot) << " rw" << outer(rdot,wdot)
169 << " vs" << outer(vdot,sdot) << " rs" << outer(rdot,sdot) << endl;
170 cout << "--- outer operator* requires col*row meaning Hermitian transpose with sign changes" << endl;
171 cout << "vw=" << vdot*~wdot << " rw" << ~rdot*~wdot
172 << " vs" << vdot*sdot << " rs" << ~rdot*sdot << endl;
173
174 cout << "*** TEST COMPLEX CROSS PRODUCT ***" << endl;
175 cout << "--- cross() global function:dot(v,w), rw, vs, rs should be the same" << endl;
176 cout << "vw=" << cross(vdot,wdot) << " rw" << cross(rdot,wdot)
177 << " vs" << cross(vdot,sdot) << " rs" << cross(rdot,sdot) << endl;
178 cout << "--- cross operator% involves NO sign changes, but returns row if either arg is a row" << endl;
179 cout << "vw=" << vdot%wdot << " rw" << rdot%wdot
180 << " vs" << vdot%sdot << " rs" << rdot%sdot << endl;
181
182 cout << "*** TEST crossMat() ***" << endl;
183 Mat<3,3,Complex> vcross(crossMat(vdot));
184 Mat<3,3,Complex> rcross(crossMat(rdot));
185 cout << "--- crossMat 3d should be same whether made from row or vec" << endl;
186 cout << "vdot%wdot=" << vdot%wdot << endl;
187 cout << "crossMat(v)=" << vcross << "crossMat(r)=" << rcross;
188 cout << "crossMat(v)*w=" << crossMat(vdot)*wdot << " vcross*w=" << vcross*wdot << endl;
189
190
191 Vec<2,Complex> vdot2 = vdot.getSubVec<2>(0);
192 Vec<2,Complex> wdot2 = wdot.getSubVec<2>(0);
193 Row<2,Complex> rdot2 = rdot.getSubRow<2>(0);
194 Row<2,Complex> vcross2(crossMat(vdot2));
195 Row<2,Complex> rcross2(crossMat(rdot2));
196
197 cout << "--- crossMat 2d should be same whether made from row or vec" << endl;
198 cout << "vdot2, wdot2=" << vdot2 << ", " << wdot2 << " vdot2%wdot2=" << vdot2%wdot2 << endl;
199 cout << "crossMat(v2)=" << vcross2 << "crossMat(r2)=" << rcross2;
200 cout << "crossMat(v2)*w2=" << crossMat(vdot2)*wdot2 << " vcross2*w2=" << vcross2*wdot2 << endl;
201
202 cout << "*********\n";
203
204
205
206 Mat<2,5,float> m25f( 1, 2, 3, 4, 5,
207 6, 7, 8, 9, 10 );
208 cout << "Mat<2,5,float>=" << m25f;
209 cout << "Mat<2,5,float>.normalize()=" << m25f.normalize();
210 cout << "Mat<2,5,float>.sqrt()=" << m25f.sqrt();
211
212 const Mat<1,5,Vec<2,float> >& m15v2f =
213 *reinterpret_cast<const Mat<1,5,Vec<2,float> >*>(&m25f);
214 cout << " m25f@" << &m25f << " m15v2f@" << &m15v2f << endl;
215 cout << "Mat<1,5,Vec<2,float> >=" << m15v2f;;
216 cout << "Mat<1,5,Vec<2,float> >.normalize()=" << m15v2f.normalize();
217
218 const Real twoXthree[] = { 1, 2, 3,
219 4, 5, 6 };
220 const Real threeXone[] = { .1, .001, .00001 };
221 const Real sym33[] = { 1,
222 2, 3,
223 4, -5, 6 };
224 SymMat<3> sm3(sym33);
225 cout << "SymMat<3> sm3=" << sm3;
226 Mat<3,3> m33sm3;
227 for (int i=0; i<3; ++i)
228 for (int j=0; j<=i; ++j)
229 m33sm3(i,j) = m33sm3(j,i) = sm3(i,j);
230 cout << "Mat33(sm3)=" << m33sm3;
231 cout << "sm3*3=" << sm3*3;
232 cout << "sm3+100=" << sm3+100;
233 cout << "m33sm3+100=" << m33sm3+100;
234 cout << "sm3.normalize()=" << sm3.normalize();
235 cout << "m33sm3.normalize()=" << m33sm3.normalize();
236 cout << "sm3+=100:" << (sm3+=100.);
237 cout << "m33sm3+=100:" << (m33sm3+=100.);
238
239 Mat<3,3,Complex> whole(mdc);
240 SymMat<3,Complex,9> sym = SymMat<3,Complex,9>().setFromLower(whole);
241 cout << "whole=" << whole << endl;
242 cout << "sym =" << sym << "(pos~)sym =" << sym.positionalTranspose() << endl;
243
244 cout << "whole.real()=" << whole.real();
245 cout << "whole.imag()=" << whole.imag();
246 cout << "sym.real()=" << sym.real();
247 cout << "sym.imag()=" << sym.imag();
248
249
250
251 Mat<3,4,Complex> mdcp(mdc); cout << "*** Data looks like this: " << mdcp;
252 SymMat<4,negator<Complex> > symp(reinterpret_cast<const negator<conjugate<double> >*>(mdc));
253 cout << " 4x4 Sym<Neg<cmplx>> from (negator<conj>)pointer to data gives this:" << symp;
254 cout << " sym.real()=" << symp.real();
255 cout << " sym.imag()=" << symp.imag();
256 cout << " ~sym.imag()=" << ~symp.imag();
257 cout << "pos~(sym.imag())=" << symp.imag().positionalTranspose();
258 cout << "(pos~sym).imag()=" << symp.positionalTranspose().imag();
259 cout << " -sym.imag()=" << -symp.imag();
260
261 symp(2,1).real() = 99.;
262 cout << "after sym(2,1).real=99, sym=" << symp;
263
264 symp.updPositionalTranspose().imag()(3,1)=123.;
265 cout << "after (pos~sym).imag()(3,1)=123, (pos~sym).imag()=" << symp.positionalTranspose().imag();
266 cout << " ... sym=" << symp;
267
268 Mat<2,3, SymMat<3,Complex> > weird(Row<3,SymMat<3,Complex> >( sym, -sym, sym ),
269 Row<3,SymMat<3,Complex> >( sym, sym, sym ));
270 cout << "weird=" << weird;
271 weird *= 2.;
272 cout << " weird*=2: " << weird;
273 cout << " weird(1)=" << weird(1) << endl;
274 cout << " weird(0,1)=" << weird(0,1) << " [0][1]=" << weird[0][1] << endl;
275
276 cout << " typename(weird)=" << typeid(weird).name() << endl;
277 cout << " typename(weird.real)=" << typeid(weird.real()).name() << endl;
278 cout << " typename(weird.imag)=" << typeid(weird.imag()).name() << endl;
279
280 cout << " weird.real()=" << weird.real();
281 cout << " weird.imag()=" << weird.imag();
282
283 cout << "sizeof(sym<3,cplx>)=" << sizeof(sym) << " sizeof(mat<2,3,sym>=" << sizeof(weird) << endl;
284
285 Mat<2,3> m23(twoXthree);
286 Mat<3,1> m31(threeXone);
287 cout << "m23=" << m23 << endl;
288 cout << "m31=" << m31 << endl;
289 cout << "m23*-m31=" << m23*-m31 << endl;
290 cout << "~ ~m31 * ~-m23=" << ~((~m31)*(~-m23)) << endl;
291
292 Mat<2,3,Complex> c23(m23);
293 Mat<3,1,Complex> c31(m31);
294 cout << "c23=" << c23 << endl;
295 cout << "c31=" << c31 << endl;
296 cout << "c23*c31=" << c23*-c31 << endl;
297 cout << " ~c31 * ~-c23=" << (~c31)*(~-c23) << endl;
298 cout << "~ ~-c31 * ~c23=" << ~((~-c31)*(~c23)) << endl;
299
300
301 Mat<3,4> m34;
302 Mat<3,4,Complex> cm34;
303
304
305 cm34 = mdc;
306 m34 = cm34.real();
307
308 cout << "Mat<3,4,Complex> cm34=" << cm34 << endl;
309 cout << "cm34.diag()=" << cm34.diag() << endl;
310
311 cout << "cm34 + cm34=" << cm34+cm34 << endl; //INTERNAL COMPILER ERROR IN Release MODE
312 cout << "~cm34 * 1000=" << ~cm34 * 1000. << endl;
313
314 cout << "m34=" << m34 << endl;
315 m34 =19.123;
316 cout << "after m34=19.123, m34=" << m34 << endl;
317
318 const double ddd[] = { 11, 12, 13, 14, 15, 16 };
319 const complex<float> ccc[] = { complex<float>(1.,2.),
320 complex<float>(3.,4.),
321 complex<float>(5.,6.),
322 complex<float>(7.,8.) };
323 Vec<2,complex<float>,1> cv2(ccc);
324 cout << "cv2 from array=" << cv2 << endl;
325 cv2 = Vec<2,complex<float> >(complex<float>(1.,2.), complex<float>(3.,4.));
326 cout << "cv2 after assignment=" << cv2 << endl;
327
328 cout << "cv2.real()=" << cv2.real() << " cv2.imag()=" << cv2.imag() << endl;
329
330 Vec<2,negator<complex<float> >,1>& negCv2 = (Vec<2,negator<complex<float> >,1>&)cv2;
331 Vec<2,conjugate<float>,1>& conjCv2 = (Vec<2,conjugate<float>,1>&)cv2;
332 Vec<2,negator<conjugate<float> >,1>& negConjCv2 = (Vec<2,negator<conjugate<float> >,1>&)cv2;
333
334
335
336 Vec<2,complex<float> > testMe = cv2;
337 cout << "testMe=cv2 (init)=" << testMe << endl;
338 testMe = cv2;
339 cout << "testMe=cv2 (assign)=" << testMe << endl;
340
341
342 cout << "(cv2+cv2)/complex<float>(1000,0):" << (cv2 + cv2) / complex<float>(1000,0) << endl;
343 cout << "(cv2+cv2)/1000.f:" << (cv2 + cv2) / 1000.f << endl;
344 cout << "(cv2+cv2)/1000.:" << (cv2 + cv2) / 1000. << endl;
345 cout << "(cv2+cv2)/1000:" << (cv2 + cv2) / 1000 << endl;
346
347 cout << "negCv2=" << negCv2 << endl;
348 cout << "conjCv2=" << conjCv2 << endl;
349 cout << "negConjCv2=" << negConjCv2 << endl;
350 cout << "cv2+negCv2=" << cv2+negCv2 << endl;
351
352 negConjCv2 = complex<float>(8,9);
353 cout << "AFTER negConjCv2 = (8,9):" << endl;
354 cout << " cv2=" << cv2 << endl;
355 cout << " negCv2=" << negCv2 << endl;
356 cout << " conjCv2=" << conjCv2 << endl;
357 cout << " negConjCv2=" << negConjCv2 << endl;
358
359 cout << "cv2: " << cv2 << endl;
360 cout << "cv2T: " << cv2.transpose() << endl;
361 cout << "-cv2: " << -cv2 << endl;
362 cout << "~cv2: " << ~cv2 << endl;
363 cout << "-~cv2: " << -(~cv2) << endl;
364 cout << "~-cv2: " << ~(-cv2) << endl;
365 cout << "~-cv2*10000: " << (~(-cv2))*10000.f << endl;
366
367 (~cv2)[1]=complex<float>(101.1f,202.3f);
368 cout << "after ~cv2[1]=(101.1f,202.3f), cv2= " << cv2 << endl;
369 (-(~cv2))[1]=complex<float>(11.1f,22.3f);
370 cout << "after -~cv2[1]=(11.1f,22.3f), cv2= " << cv2 << endl;
371
372 Vec<3> dv3(ddd), ddv3(ddd+3);
373 dv3[2] = 1000;
374 cout << "dv3=" << dv3 << " ddv3=" << ddv3 << endl;
375 cout << "100(ddv3-dv3)/1000=" << 100.* (ddv3 - dv3) / 1000. << endl;
376
377 Vec<3> xxx(dv3); cout << "copy of dv3 xxx=" << xxx << endl;
378 Vec<3> yyy(*ddd);cout << "copy of *ddd yyy=" << yyy << endl;
379
380 cout << "dv3.norm()=" << dv3.norm() << endl;
381 cout << "cv2=" << cv2 << " cv2.norm()=" << cv2.norm() << endl;
382
383 const Vec<2> v2c[] = {Vec<2>(ddd),Vec<2>(ddd+1)};
384 Vec<2, Vec<2> > vflt(v2c);
385 cout << "vflt 2xvec2=" << vflt << endl;
386 cout << "10.*vflt=" << 10.*vflt << endl;
387 cout << "vflt*10.=" << vflt*10. << endl;
388
389 int ivals[] = {0x10, 0x20, 0x30, 0x40};
390 Vec<4> iv(ivals);
391 cout << "iv=" << iv << endl;
392
393 Vec<2, Vec<2> > v22;
394 v22 = Vec<2>(&ivals[2]);
395 cout << "v22=" << v22 << endl;
396
397
398 // Test dot product
399 {
400 double d[] = {1,2,3,4,5,6,7,8};
401
402
403 Vec<2> v1(&d[0]), v2(&d[2]);
404 Row<2> r1(&d[4]), r2(&d[6]);
405 Vec<2>::TNeg& nv1 = (Vec<2>::TNeg&)v1;
406
407 negator<double> nd(100); cout << endl << "nd=" << nd << endl;
408 cout << "nv1=" << nv1 << endl;
409 cout << "nv1*nd=" << nv1*nd << endl;
410 cout << "nd*nv1=" << nd*nv1 << endl;
411 cout << "nv1/nd=" << nv1/nd << endl << endl;
412
413
414 cout << "v1,v2=" << v1 << v2 << endl;
415 cout << "r1,r2=" << r1 << r2 << endl;
416 cout << "dot r1*v1 =" << dot(r1,v1) << endl;
417 cout << "dot r1*nv1=" << dot(r1,nv1) << endl;
418 cout << "r1*v1 =" << r1*v1 << endl;
419 cout << "r1*nv1=" << r1*nv1 << endl;
420
421 // outer product
422 cout << " outer v1*r1=" << v1*r1 << endl;
423 cout << " outer nv1*r1=" << nv1*r1 << endl;
424
425 // cross product (2d)
426 cout << "cross(v1,v2)=" << cross(v1,v2) << endl;
427 cout << "v1 % v2=" << v1 % v2 << endl;
428 cout << "cross(r1,v2)=" << cross(r1,v2) << endl;
429 cout << "r1 % v2=" << r1 % v2 << endl;
430 cout << "cross(v1,r2)=" << cross(v1,r2) << endl;
431 cout << "v1 % r2=" << v1 % r2 << endl;
432 cout << "cross(r1,r2)=" << cross(r1,r2) << endl;
433 cout << "r1 % r2=" << r1 % r2 << endl;
434
435 // do the cross products with 3d routines
436 Vec3 v13(v1[0],v1[1],0), v23(v2[0],v2[1],0);
437 Row3 r13(r1[0],r1[1],0), r23(r2[0],r2[1],0);
438 cout << "cross(v13,v23)=" << cross(v13,v23) << endl;
439 cout << "v13 % v23=" << v13 % v23 << endl;
440 cout << "cross(r13,v23)=" << cross(r13,v23) << endl;
441 cout << "r13 % v23=" << r13 % v23 << endl;
442 cout << "cross(v13,r23)=" << cross(v13,r23) << endl;
443 cout << "v13 % r23=" << v13 % r23 << endl;
444 cout << "cross(r13,r23)=" << cross(r13,r23) << endl;
445 cout << "r13 % r23=" << r13 % r23 << endl;
446
447 v13[2]=7;
448 cout << "v13=" << v13 << " 2*v13+.001=" << 2*v13+.001 << endl;
449 cout << "v13 % (2*v13)=" << v13 % (2*v13) << endl;
450 cout << "v13 % (2*v13+.001)=" << v13 % (2*v13+.001) << endl;
451
452 cout << endl;
453
454 // test constructors
455 Mat<2,3> mvcols( v1, ~r1, v2 );
456 Mat<3,2> mvrows( ~v1,
457 r1,
458 r2 );
459
460 cout << "mvcols=" << mvcols << endl;
461 cout << "mvrows=" << mvrows << endl;
462
463 Vec<3,float> v2f(39.f, 40.f, 50.L);
464 cout << "v2f=" << v2f << endl;
465
466 cout << "v2f.drop1(0)=" << v2f.drop1(0) << endl;
467 cout << "v2f.drop1(1)=" << v2f.drop1(1) << endl;
468 cout << "v2f.drop1(2)=" << v2f.drop1(2) << endl;
469
470 cout << "v2f.append1(3.3f)=" << v2f.append1(3.3f) << endl;
471 cout << "v2f.append1((short)1)=" << v2f.append1((short)1) << endl;
472 cout << "v2f.drop1(1).append1((unsigned short)0x10)=" << v2f.drop1(1).append1((unsigned short)0x10) << endl;
473 cout << "(~v2f).drop1(1).append1((unsigned short)0x10)=" << (~v2f).drop1(1).append1((unsigned short)0x10) << endl;
474
475 cout << "v2f.insert1(0, 23.f)=" << v2f.insert1(0, 23.f) << endl;
476 cout << "v2f.insert1(1, 23.f)=" << v2f.insert1(1, 23.f) << endl;
477 cout << "v2f.insert1(2, 23.f)=" << v2f.insert1(2, 23.f) << endl;
478 cout << "v2f.insert1(3, 23.f)=" << v2f.insert1(3, 23.f) << endl;
479 cout << "v2f.insert1(2, 23.f).drop1(2)=" << v2f.insert1(2, 23.f).drop1(2) << endl;
480
481 cout << "(~v2f).insert1(2, 23.f).drop1(2)=" << (~v2f).insert1(2, 23.f).drop1(2) << endl;
482
483 Mat33 m33( Row3(1, 2, 3),
484 Row3(4, 5, 6),
485 Row3(.003f, 9.62L, 41.1) );
486 cout << "m33=" << m33 << endl;
487 cout << "v13=" << v13 << endl;
488 cout << "m33*v13=" << m33*v13 << endl;
489 cout << "~m33*v13=" << (~m33)*v13 << endl;
490 cout << "~v13*~m33=" << ~v13*(~m33) << endl;
491 }
492
493 Mat<4,3> ident43;
494 ident43 = 1;
495 cout << "ident43=" << ident43 << endl;
496
497 Mat<4,3> negid43 = -ident43; // requires implicit conversion from Mat<negator<Real>>
498 cout << "negid43=" << negid43 << endl;
499
500 // Absolute value
501 const Vec<3> vorig(1.,-2.,-3.);
502 cout << "vorig=" << vorig << " vorig.abs()=" << vorig.abs() << endl;
503 const Vec<2, Vec<3> > v2orig(vorig,-vorig);
504 cout << "v2orig=" << v2orig << " v2orig.abs()=" << v2orig.abs() << endl;
505
506 Vec<3> nvorig = -vorig;
507 Vec<2, Vec<3> > nv2orig = -v2orig;
508
509 Mat<3,2> morig(vorig,-vorig);
510 cout << "morig=" << morig << " morig.abs()=" << morig.abs() << endl;
511
512 //Mat<2,2, Vec<3> > m2orig(v2orig, (Vec<2, Vec<3> >)-v2orig);
513 //cout << "m2orig=" << m2orig << " m2orig.abs()=" << m2orig.abs() << endl;
514
515 cout << "vorig.getSubVec<2>(1)=" << vorig.getSubVec<2>(1) << endl;
516
517 negid43.updSubMat<2,2>(2,1) = -27.;
518 cout << "after negid43.updSubMat<2,2>(2,1) = -27., negid43=" << negid43;
519
520 cout << "negid43[2].getSubRow<2>(1)=" << negid43[2].getSubRow<2>(1) << endl;
521
522 cout << "CHECK DIAGONAL LENGTH FOR RECTANGULAR MATRICES" << endl;
523 Mat<3,2, Row3, 1, 2> H;
524 cout << "H[" << H.nrow() << "," << H.ncol() << "]" << endl;
525 cout << "H.diag()[" << H.diag().nrow() << "," << H.diag().ncol() << "]" << endl;
526 Mat<3,2, Row3, 1, 2>::TransposeType Ht;
527 cout << "Ht[" << Ht.nrow() << "," << Ht.ncol() << "]" << endl;
528 cout << "Ht.diag()[" << Ht.diag().nrow() << "," << Ht.diag().ncol() << "]" << endl;
529 }
530
testMatInverse()531 void testMatInverse() {
532 Matrix m = Test::randMatrix(20,20);
533 Matrix mi = m.invert();
534 Matrix id(20,20); id=1; // identity
535 SimTK_TEST_EQ_SIZE(m*mi, id, 20);
536 }
537
538
main()539 int main() {
540 SimTK_START_TEST("MatVecTest");
541
542 SimTK_SUBTEST(testSums);
543 SimTK_SUBTEST(testNegator);
544 SimTK_SUBTEST(testElementwiseOps);
545 SimTK_SUBTEST(testMiscellaneous);
546 SimTK_SUBTEST(testMatInverse);
547
548 SimTK_END_TEST();
549 }
550
551
552
553