1 /* -------------------------------------------------------------------------- *
2 * SimTK Simbody: 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.h"
25
26 #include <cstdlib> // for rand()
27 #include <cstdio>
28 #include <ctime>
29 #include <string>
30 #include <complex>
31 #include <iostream>
32 using std::cout;
33 using std::endl;
34 using std::complex;
35
36 using namespace SimTK;
37
38 // Numerical Recipes declarations.
39 namespace NR {
40 template <class DP>
41 void lubksb(const int N, const DP* a/*N,N*/, const int* indx/*N*/,
42 DP* b/*N*/);
43
44 template <class DP>
45 void ludcmp(const int N, DP* a/*N,N*/, int* indx/*N*/, DP &d);
46
47 template <class DP>
48 void luinvert(int N, DP* a/*N,N*/, DP* y/*N,N*/);
49 }
50
51 // Some explicit instantiations just to make sure everything's there.
52 namespace SimTK {
53 template class Matrix_<Real>;
54 template class Vector_<Complex>;
55 template class RowVector_< conjugate<float> >;
56 //template class MatrixBase< Mat<3,4,Vec2> >;
57
58 template class MatrixView_< complex<double> >;
59 template class VectorView_< negator<float> >;
60 template class RowVectorView_< negator< conjugate<float> > >;
61 }
62
63 template <class NT>
dump(const String & s,const Matrix_<NT> & mm)64 void dump(const String& s, const Matrix_<NT>& mm) {
65 cout << s << ":" << endl;
66 for (int i=0; i<mm.nrow(); ++i) {
67 for (int j=0; j<mm.ncol(); ++j)
68 cout << mm(i,j) << " ";
69 cout << endl;
70 }
71 }
72
73 template <class NT>
dump(const String & s,const MatrixView_<NT> & mv)74 void dump(const String& s, const MatrixView_<NT>& mv)
75 { dump(s,(const Matrix_<NT>&)mv); }
76
77 extern "C" {
78 void SimTK_version_SimTKlapack(int* major, int* minor, int* build);
79 void SimTK_about_SimTKlapack(const char* key, int maxlen, char* value);
80 }
81
testSums()82 void testSums() {
83 Matrix m(Mat22(1,2,
84 3,4));
85 SimTK_TEST_EQ(m.colSum(), RowVector(Row2(4,6)));
86 SimTK_TEST_EQ(m.rowSum(), Vector(Vec2(3,7)));
87 SimTK_TEST_EQ(m.sum(), m.colSum()); // should be exact
88
89 Matrix_<Complex> mc(Mat<2,2,Complex>(1+2*I, 3+4*I,
90 5+6*I, 7+8*I));
91 typedef Row<2,Complex> CRow2;
92 typedef Vec<2,Complex> CVec2;
93 typedef RowVector_<Complex> CRowVector;
94 typedef Vector_<Complex> CVector;
95 SimTK_TEST_EQ(mc.colSum(), CRowVector(CRow2(6+8*I, 10+12*I)));
96 SimTK_TEST_EQ(mc.rowSum(), CVector(CVec2(4+6*I, 12+14*I)));
97 SimTK_TEST_EQ(mc.sum(), mc.colSum()); // should be exact
98 }
99
testCharacter()100 void testCharacter() {
101 MatrixCommitment mc;
102 // mc = MatrixCommitment::Symmetric();
103 //mc.commitOutline(MatrixOutline::Square);
104 cout << mc;
105 cout << mc.calcDefaultCharacter(0,0);
106 Matrix m(mc);
107 m.resize(2,3);
108 m = 5;
109 m.dump("m");
110 cout << "NOW: " << m.getMatrixCharacter();
111 cout << "m.diag()=" << m.diag() << endl;
112 cout << "m.diag() " << m.diag().getMatrixCharacter();
113 cout << "m[1]=" << m[1] << endl;
114 Matrix mm = m*3;
115 cout << "mm=" << mm;
116
117 MatrixView mmv = mm(0,1,2,2);
118 cout << "mmv(mm(0,1,2,2):" << mmv;
119 cout << "mm(0,0,2,2):" << mm(0,0,2,2);
120 mmv.clear();
121 mmv = mm(0,0,1,2);
122 cout << "mmv.clear(), then mmv=mm(0,0,1,2): " << mmv;
123
124
125 mc = MatrixCommitment( MatrixStructure(MatrixStructure::Triangular, MatrixStructure::Upper) );
126 cout << "mc=" << mc << " actual=" << mc.calcDefaultCharacter(0,0);
127 Matrix t(mc);
128 t.resize(5,3);
129 t.dump("t");
130
131 t.clear();
132 t.commitTo(MatrixStructure(MatrixStructure::Hermitian, MatrixStructure::Upper));
133 t.dump("t now hermitian, prior to resize");
134 cout << "t commit=" << t.getCharacterCommitment() << "t actual=" << t.getMatrixCharacter();
135 t.resize(3,5);
136 t.dump("t after resize");
137 cout << "t commit=" << t.getCharacterCommitment() << "t actual=" << t.getMatrixCharacter();
138 t = 123;
139 t(0,2)=-2;
140 t.dump("t now hermitian");
141 cout << "t commit=" << t.getCharacterCommitment() << "t actual=" << t.getMatrixCharacter();
142
143 Vector v(10);
144 for (int i=0; i<10; ++i) v[i] = i*.1;
145 v[4] = -17.3;
146 cout << "v commitment: " << v.getCharacterCommitment();
147 cout << "v character: " << v.getMatrixCharacter();
148 cout << "v=" << v << endl;
149 Vector w1(10, Real(1));
150 cout << "weights w1=" << w1 << endl;
151 int worst;
152 cout << "|w1*v|_rms=" << v.weightedNormRMS(w1,&worst) << "\n";
153 cout << "(worst=" << worst << ")\n";
154 cout << "|w1*v|_inf=" << v.weightedNormInf(w1,&worst) << "\n";
155 cout << "(worst=" << worst << ")\n";
156 cout << "|v|_rms=" << v.normRMS(&worst) << "\n";
157 cout << "(worst=" << worst << ")\n";
158 cout << "|v|_inf=" << v.normInf(&worst) << "\n";
159 cout << "(worst=" << worst << ")\n";
160 cout << "2nd sig: |v|_rms=" << v.normRMS() << "\n";
161 cout << "2nd sig: |v|_inf=" << v.normInf() << "\n";
162
163 w1(9) = 100;
164 cout << "weights w1=" << w1 << endl;
165 cout << "|w1*v|_rms=" << v.weightedNormRMS(w1,&worst) << "\n";
166 cout << "(worst=" << worst << ")\n";
167 cout << "|w1*v|_inf=" << v.weightedNormInf(w1,&worst) << "\n";
168 cout << "(worst=" << worst << ")\n";
169 cout << "2nd sig: |w1*v|_rms=" << v.weightedNormRMS(w1) << "\n";
170 cout << "2nd sig: |w1*v|_inf=" << v.weightedNormInf(w1) << "\n";
171
172 cout << "rms(zero length)=" << Vector().normRMS(&worst) << "\n";
173 cout << " worst=" << worst << "\n";
174 cout << "inf(zero length)=" << Vector().normInf(&worst) << "\n";
175 cout << " worst=" << worst << "\n";
176
177
178 Array_<int> vx;
179 vx.push_back(2); vx.push_back(5); vx.push_back(7); vx.push_back(8);
180
181 VectorView vxx = v(vx);
182 cout << "vxx character: " << vxx.getMatrixCharacter();
183 cout << "vxx=" << vxx << endl;
184 cout << "vxx(1,2)=" << vxx(1,2) << endl;
185
186 Complex cmplx[] = {Complex(1,2), Complex(3,4), Complex(-.2,.3),
187 Complex(-100,200), Complex(20,40), Complex(-.001,.002)};
188
189 ComplexMatrix cm(2,3,2,cmplx); // shared data is in column order
190 cout << "cm(2,3,lda=2,cmplx)=" << cm;
191
192 ComplexMatrix cm2(2,3,cmplx); // init constructor takes data in row order
193 cout << "cm(2,3,cmplx)=" << cm2;
194
195 // Share data but note that it is row order. Using MatrixBase to get to
196 // a general constructor.
197 MatrixCharacter::LapackFull mchar(2,3);
198 mchar.setStorage(MatrixStorage(MatrixStorage::NoPacking,MatrixStorage::RowOrder));
199 MatrixBase<Complex> cm3(MatrixCommitment(), mchar, 3, cmplx);
200 cout << "cm3(RowOrder,2,3,lda=3,cmplx)=" << cm3;
201 cm3(0,1)=99;
202 cout << "cm3(0,1)=99 ->" << cm3;
203 cout << "cm (should have changed too) ->" << cm;
204 cout << "cm2 (should not have changed) ->" << cm2;
205 }
206
207 // Test "scalar" multiply for Vectors and RowVectors that have CNT types
208 // as elements.
209 // NOTE: the Vector elements and the "scalar" CNT must conform for this
210 // to work; nonconforming will cause obscure compile errors.
testScalarMultiply()211 void testScalarMultiply() {
212 cout << "\n------ TEST SCALAR MULTIPLY ------\n";
213 Mat33 m33( .03, .04, .05,
214 .06, .08, .09,
215 .07, .10, .11 );
216 Vector_< SpatialVec > vs(3, SpatialVec(Vec3(1,2,3), Vec3(4,5,6)));
217 cout << "vs=" << vs << endl;
218 cout << "vs*SpatialRow=" << vs*SpatialRow(Row3(.1)) << endl;
219 cout << "SpatialRow*vs=" << SpatialRow(Row3(.1))*vs << endl;
220 cout << "m33 * SpatialVec=" << m33 * SpatialVec(Vec3(1,2,3), Vec3(4,5,6)) << endl;
221 cout << "m33 * vs=" << SpatialMat(m33) * vs << endl; // note cast to conforming diagonal matrix
222 cout << "SpatialRow * m33=" << ~SpatialVec(Vec3(1,2,3), Vec3(4,5,6)) * m33 << endl;
223 cout << "~vs*m33=" << ~vs * SpatialMat(m33) << endl; // note cast to conforming diagonal matrix
224 RowVector_< SymMat22 > rv(3, SymMat22(1,2,3));
225 cout << "rv=" << rv << endl;
226 cout << "rv*Vec2=" << rv*Vec2(.1,.2) << endl;
227 cout << "Row2*rv=" << Row2(.1,.2)*rv << endl;
228 cout << "------ END TEST SCALAR MULTIPLY ------\n\n";
229 }
230
231
testAjaysBlock()232 void testAjaysBlock() {
233 cout << "\n------ TEST AJAY'S BLOCK ------\n";
234 const int nu =7, nm=4;
235 Matrix J(6,nu);
236 for (int i=0; i<6; ++i)
237 for (int j=0; j<nu; ++j)
238 J(i,j) = 1000*i+j;
239 Matrix t = ~J(0,3,3,nm);
240 cout << "J=" << J << endl;
241 cout << "t=" << t;
242 cout << "\n------ END TEST AJAY'S BLOCK ------\n";
243 }
244
main()245 int main()
246 {
247 try {
248 SimTK_DEBUG("Running BigMatrixTest ...\n");
249
250 testSums();
251
252 Matrix assignToMe(5,4);
253 assignToMe.elementwiseAssign(1.);
254 std::cout << "assignToMe=" << assignToMe;
255 assignToMe.elementwiseAssign(14);
256 std::cout << "assignToMe=" << assignToMe;
257
258 testAjaysBlock();
259 testScalarMultiply();
260 testCharacter();
261
262 int major,minor,build;
263 char out[100];
264 const char* keylist[] = { "version", "library", "type", "debug", "authors", "copyright", "svn_revision", 0 };
265
266 //SimTK_version_SimTKlapack(&major,&minor,&build);
267 //std::printf("SimTKlapack library version: %d.%d.%d\n", major, minor, build);
268
269
270 SimTK_version_SimTKcommon(&major,&minor,&build);
271 std::printf("==> SimTKcommon library version: %d.%d.%d\n", major, minor, build);
272 std::printf(" SimTK_about_SimTKcommon():\n");
273 for (const char** p = keylist; *p; ++p) {
274 SimTK_about_SimTKcommon(*p, 100, out);
275 std::printf(" about(%s)='%s'\n", *p, out);
276 }
277
278 const Real vvvdata[] = {1,2,.1,.2,9,10,-22,-23,-24,25};
279 Vector vvv(10,vvvdata);
280 cout << "vvv=" << vvv << endl;
281 cout << "vvv(2,5)=" << vvv(2,5) << endl;
282 Vector vvv25;
283 vvv25.viewAssign(vvv(2,5));
284 cout << "vvv25=" << vvv25 << endl;
285 cout << "vvv(2,0)=" << vvv(2,0) << endl;
286 Vector vvv20;
287 vvv20.viewAssign(vvv(2,0));
288 cout << "vvv20=" << vvv20 << endl;
289 cout << "vvv(0,0)=" << vvv(0,0) << endl;
290 Vector vvv00;
291 vvv00.viewAssign(vvv(0,0));
292 cout << "vvv00=" << vvv00 << endl;
293 Vector vb;
294 vvv00 = vb;
295 cout << "after vvv00=vb [null], vvv00(" << vvv00.nrow() << "," << vvv00.ncol() << ")=" << vvv00 << endl;
296
297 const Complex mdc[] = { Complex(1.,2.),
298 Complex(3.,4.),
299 Complex(5.,6.),
300 Complex(7.,8.),
301 Complex(9.,10.),
302 Complex(10.,11.),
303 Complex(.1,.26),
304 Complex(.3,.45),
305 Complex(.5,.64),
306 Complex(.7,.83),
307 Complex(.9,.102),
308 Complex(.10,.111)
309 };
310 Matrix_<Complex> md(2,2,mdc);
311 cout << "2x2 complex Matrix md=" << md;
312 Mat<2,2,Complex> md_mat(mdc);
313 cout << "2x2 complex Mat md_mat=" << md_mat;
314 cout << " sizeof(Complex)=" << sizeof(Complex) << " sizeof(md_mat)=" << sizeof(md_mat) << endl;
315 Matrix_<Mat<2,2,Complex> > mm22c;
316 Mat<2,2, Mat<2,2,Complex> > mm22c_mat;
317 cout << " sizeof(mm22c_mat)=" << sizeof(mm22c_mat) << " should be " << 16*sizeof(Complex) << endl;
318 mm22c.resize(2,2); cout << " sizeof(mm22c(2,2))=" << ((char*)(&mm22c(1,1)(1,1)+1)) - ((char*)&mm22c(0,0)(0,0)) << endl;
319
320 cout << "scalar md.normRMS: " << md.normRMS() << endl;
321 try {
322 cout << "nonscalar mm22c.normRMS: " << mm22c.normRMS() << endl;
323 } catch(const std::exception& e) {
324 cout << "SHOULD FAIL FOR NONSCALAR ELEMENTS: " << e.what() << endl;
325 }
326
327
328 mm22c.setToZero(); cout << "mm22c after setToZero():" << mm22c;
329 mm22c.setToNaN(); cout << "mm22c after setToNaN():" << mm22c;
330
331 mm22c.dump("**** mm22c ****");
332
333 cout << "~md=" << ~md;
334 cout << "~md(1)=" << ~md(1) << endl;
335 cout << "(~md)(1)=" << (~md)(1) << endl;
336 cout << "~md[1]=" << ~md[1] << endl;
337 cout << "(~md)[1]=" << (~md)[1] << endl;
338
339 dump("2x2 complex md",md);
340 dump("md(0,1,2,1)",md(0,1,2,1));
341 const ComplexMatrixView& mvc = md(0,1,2,1);
342 dump("mvc=md(0,1,2,1)", mvc);
343 md(1,0) *= 10.;
344 dump("md after md(1,0) *= 10.",md);
345 md(0,1,2,1) *= Complex(10.,100.);
346 dump("md after md(0,1,2,1)*=(10+100i)",md);
347
348 Matrix_<Complex> mm(3,4);
349 mm = 2390.;
350 cout << "after [3x4 complex] mm = 2390, mm: " << mm << endl;
351 for (int i=0; i<mm.nrow(); ++i) {
352 for (int j=0; j<mm.ncol(); ++j)
353 mm(i,j)=(i+1)*(j+1);
354 }
355 cout << "after mm(i,j)=(i+1)*(j+1), mm: " << mm << endl;
356
357 Vector mmColScale(4), mmRowScale(3);
358 mmColScale[0]=1; mmColScale[1]=10; mmColScale[2]=100; mmColScale[3]=1000;
359 mmRowScale[0]=-1000; mmRowScale[1]=-100; mmRowScale[2]=-10;
360 Vector_<double> mmRowScaleR(3); for(int i=0;i<3;++i) mmRowScaleR[i]=(float)mmRowScale[i];
361
362 cout << "mm=" << mm << " mmColScale=" << mmColScale << endl;
363 mm.colScaleInPlace(mmColScale);
364 cout << "after col scale, mm=" << mm;
365
366 mm.colScaleInPlace(mmColScale.elementwiseInvert());
367 cout << "after col UNscale mm=" << mm;
368
369 cout << " mmRowScale=" << mmRowScale << endl;
370 mm.rowScaleInPlace(mmRowScale);
371 cout << "after row scale, mm=" << mm;
372
373 mm.rowScaleInPlace(mmRowScale.elementwiseInvert());
374 cout << "after row UNscale mm=" << mm;
375
376 mm.rowScaleInPlace(mmRowScaleR);
377 cout << "after LONG DOUBLE row scale mm=" << mm;
378
379 cout << "mm.rowScale(double)=" << mm.rowScale(mmRowScaleR);
380 cout << "type(mm.rowScale(double))=" << typeid(mm.rowScale(mmRowScaleR)).name() << endl;
381
382 Vector_<Vec2> mmVCol(4); for (int i=0; i<4; ++i) mmVCol[i] = Vec2(4*i, 4*i+1);
383 cout << "mm Vec2 colScale=" << mmVCol << endl;
384 cout << "mm.colScale(Vec2)=" << mm.colScale(mmVCol);
385 cout << "type(mm.colScale(Vec2))=" << typeid(mm.colScale(mmVCol)).name() << endl;
386
387 mm *= 1000.; dump("mm(3,4) after *=1000", mm);
388 mm.dump("*** mm ***");
389
390 cout << "mm using << op: " << mm;
391 cout << "~mm: " << ~mm;
392 cout << "mm.norm()=" << mm.norm() << endl;
393 cout << "mm(3)=" << mm(3) << endl;
394 cout << "mm[2]=" << mm[2] << endl;
395 cout << "mm + 2*(-mm)=" << mm + 2*(-mm) << endl;
396
397 Matrix_<Complex> mm2 = mm;
398 cout << "mm2=" << mm2;
399 cout << "mm2(1,1,1,2)=" << mm2(1,1,1,2);
400 mm2(1,1,1,2)=7;
401 cout << "after mm2(1112)=7, mm2=" << mm2;
402
403 cout << "mm=" << mm;
404 (~mm)(1,1,1,2)=99;
405 cout << "after (~mm)(1112)=99, mm=" << mm;
406
407 cout << "\n--- DIAGONAL TEST ---\n";
408 cout << "mm2.diag()=" << mm2.diag() << endl;
409
410 mm2 = 99;
411 cout << "after mm2=99: mm2=" << mm2;
412
413 mm2(0,2,2,2) = std::complex<double>(-99,-99);
414 cout << "after mm2(0,2,2,2) = (-99,-99): mm2=" << mm2;
415
416 mm2.updDiag() *= .001;
417 cout << "after mm2.updDiag()*=.001: mm2=" << mm2;
418
419 // scalar assign to a row-shaped matrix should set only
420 // the first element to the scalar and all else zero
421 // (since that's the diagonal), while
422 // scalar assign to a row should set all the elements.
423 mm2(0,0,1,mm2.ncol()) = 1; // row shaped block of first row
424 mm2[mm2.nrow()-1] = 1; // last row
425 cout << "after mm2(0,0,1,n)=1 and mm2[m-1]=1: mm2=" << mm2;
426
427 // ditto for columns
428 mm2(0,1,mm2.nrow(),1) = 2;
429 mm2(2) = 2;
430
431 cout << "after mm2(0,1,m,1)=2 and mm2(2)=2: mm2=" << mm2;
432
433
434 cout << "\n-------- ASSIGN TEST --------\n";
435 Matrix_< Vec<2,Complex> > rr(2,3);
436 for (int i=0;i<2;++i) for (int j=0;j<3;++j)
437 rr(i,j) = Vec<2,Complex>(mdc[i]*(j+1), mdc[i]*(j-1));
438 Matrix_< Vec<2,Complex> > rrAssign;
439 (rrAssign = rr(0,1,2,2)) *= 1000.;
440 cout << "rr=" << rr;
441 cout << "(rrAssign=rr(0,1,2,2)) *= 1000.; rrAssign=" << rrAssign;
442 cout << "rr=" << rr;
443 (rrAssign.viewAssign(rr(0,1,2,2))) *= 100.;
444 cout << "(rrAssign.viewAssign(rr(0,1,2,2))) *= 100; rrAssign=" << rrAssign;
445 cout << "rr=" << rr;
446 cout << "-------- END ASSIGN TEST --------\n\n";
447
448 cout << "\n-------- RESIZE KEEP TEST --------\n";
449 Vector resizeMe(5); for (int i=0; i<5; ++i) resizeMe[i]=i;
450 cout << "resizeMe=" << resizeMe << endl;
451 resizeMe.resize(10);
452 cout << "after resize(10), resizeMe=" << resizeMe << endl;
453 Vector resizeMe2(5); for (int i=0; i<5; ++i) resizeMe2[i]=i;
454 cout << "resizeMe2=" << resizeMe2 << endl;
455 resizeMe2.resizeKeep(10);
456 cout << "after resizeKeep(10), resizeMe2=" << resizeMe2 << endl;
457 Matrix resizeMem(2,3); for (int i=0; i<2; ++i) for (int j=0; j<3; ++j) resizeMem(i,j)=(i+1)*(j+1);
458 cout << "resizeMem(2,3)=" << resizeMem << endl;
459 resizeMem.resizeKeep(3,5);
460 cout << "after resizeKeep(3,5), resizeMem=" << resizeMem << endl;
461 resizeMem.resizeKeep(2,2);
462 cout << "after resizeKeep(2,2), resizeMem=" << resizeMem << endl;
463 resizeMem.resize(3,4);
464 cout << "after resize(3,4), resizeMem=" << resizeMem << endl;
465 cout << "-------- END RESIZE KEEP TEST --------\n\n";
466
467 Mat<3,4,Complex> cm34;
468 for (int i=0; i<3; ++i)
469 for (int j=0; j<4; ++j)
470 cm34(i,j) = mdc[i+j*3];
471 cout << "Mat<3,4,Complex>=" << cm34 << endl;
472 Vec<4,Complex> cv4(&mdc[6]);
473 cout << "Vec<4,Complex>=" << cv4 << endl;
474 cout << "Mat<3,4>*Vec<4>=" << cm34*cv4 << endl;
475 cout << "Mat<3,4>*Mat<4,3>=" << cm34*~cm34;
476
477 complex<float> zzzz(1.,2.);
478 conjugate<float> jjjj(0.3f,0.4f);
479 negator<float> nnnn(7.1);
480 cout << "zzzz=" << zzzz << " jjjj=" << jjjj << " nnnn=" << nnnn << endl;
481 cout << "zzzz*jjjj=" << zzzz*jjjj << endl;
482
483 Matrix_<Complex> cMatrix34(3,4);
484 Vector_<Complex> cVector4(4);
485 for (int i=0; i<3; ++i) for (int j=0; j<4; ++j) cMatrix34(i,j)=cm34(i,j);
486 for (int i=0; i<4; ++i) cVector4[i] = cv4[i];
487 cout << "cMatrix34=" << cMatrix34;
488 cout << "cVector4=" << cVector4 << endl;
489 cout << "cMatrix34*cVector4=" << cMatrix34*cVector4 << endl;
490 cout << "cMatrix34*cMatrix43=" << cMatrix34*~cMatrix34;
491
492 Matrix_<Complex> cMatrix34N = -cMatrix34;
493
494 //TODO: not allowed yet
495 //Matrix_<Complex> cMatrix34H = ~cMatrix34;
496
497 Vector vv(4), ww;
498 vv[0] = 1.; vv[1] = 2.; vv[2] = 3.; vv[3] = 4.;
499 cout << "vv(4)=" << vv << endl;
500 vv *= 9.; cout << "vv*=9:" << vv << endl;
501
502 Vector vvvv;
503 vvvv = vv[2] * vv;
504 cout << "vvvv = vv[2]*vv = " << vvvv << endl;
505
506 Matrix_<Mat<2,2, Mat<2,2,double> > > mmm(2,1);
507 mmm = Mat<2,2, Mat<2,2,double> >(Mat<2,2,double>(1));
508 cout << "*****>>> mmm=" << mmm << endl;
509
510 Matrix mnm(4,2), nn;
511 cout << "mnm(4,2)=" << mnm;
512 mnm(0) = vv;
513 mnm(1) = -0.01 * vv;
514 cout << "mnm(vv,-.01*vv)=" << mnm;
515 cout << "===> mnm.abs()=" << mnm.abs();
516 cout << "mnm(1)=" << mnm(1) << endl;
517 cout << "mnm(1).abs()=" << mnm(1).abs() << endl;
518 cout << "mnm[1]=" << mnm[1] << endl;
519 cout << "mnm[1].abs()=" << mnm[1].abs() << endl;
520
521 ww = vv;
522 ww *= 0.1;
523 vv += ww;
524 cout << "ww=vv, *=0.1:" << ww << endl;
525 cout << "vv+=ww:" << vv << endl;
526
527 const Real rdata[]={1,2,3,
528 9,.1,14,
529 2,6,9};
530 Matrix_<negator<Real> > A(3,3, (negator<Real>*)rdata);
531
532 Matrix AI = A.invert();
533 cout << "A=" << A << "AI=" << AI << " A*AI=" << A*AI;
534
535 cout << "A(1,2).real()=" << A(1,2).real() << endl;
536
537 Matrix AH = ~A;
538 cout << "~A=" << ~A << "inv(~A)=" << (~A).invert() << "~(inv(A))=" << ~AI;
539
540 A.invertInPlace();
541 cout << "after A=inv(A), A=" << A << "norm(A-AI)=" << (A-AI).norm() << endl;
542
543 A.dump("*** A ***");
544 AI.dump("*** AI ***");
545
546 Mat<3,3,negator<Real> > smallNegA((negator<Real>*)rdata);
547 Mat<3,3,negator<Real> >::TInvert smallNegAI(smallNegA.invert());
548
549 cout << "smallNegA=" << smallNegA << " inv(smallNegA)=" << smallNegAI;
550 cout << "smallNegA*inv(smallNegA)=" << smallNegA*smallNegAI << "NORM="
551 << (smallNegA*smallNegAI).norm() << endl;
552 cout << "inverse(smallNegA)=" << inverse(smallNegA) << endl;
553
554 negator<Real> nnn = smallNegA(0,0)-smallNegA(1,1);
555 Real nnnr = smallNegA(0,0)-smallNegA(1,1);
556 cout << "negator nnn=" << nnn << " real nnnr=" << nnnr << endl;
557
558 nnn = smallNegA(0,1)-smallNegA(1,0);
559 nnnr = smallNegA(0,1)-smallNegA(1,0);
560 cout << "negator nnn=" << nnn << " real nnnr=" << nnnr << endl;
561
562 cout << "det(smallNegA)=" << det(smallNegA)
563 << " det(inv(smallNegA))=" << det(smallNegAI) << endl;
564
565 const Real cjdata[]={1,1, 2,2, 3,3, 4,4,
566 9,9, .1,.1, 14,14, 22,22,
567 2,2, 6,6, 9,9, 11,11,
568 .2,.2, .7,.7, 5,5, 10,10};
569
570 // General Lapack inverse.
571 Mat<4,4,conjugate<Real> > smallConjA4((conjugate<Real>*)cjdata);
572 Mat<4,4,conjugate<Real> >::TInvert smallConjAI4(smallConjA4.invert());
573
574
575 cout << "smallConjA4=" << smallConjA4 << " inv(smallConjA4)=" << smallConjAI4;
576 cout << "smallConjA4*inv(smallConjA4)=" << smallConjA4*smallConjAI4;
577 cout << "NORM=" << (smallConjA4*smallConjAI4).norm() << endl;
578 cout << "inverse(smallConjA4)=" << inverse(smallConjA4) << endl;
579 cout << "norm(inverse-lapackInverse4)=" << (inverse(smallConjA4)-lapackInverse(smallConjA4)).norm() << endl;
580
581 cout << "det(smallConjA4)=" << det(smallConjA4)
582 << " det(inv(smallConjA4))=" << det(smallConjAI4) << endl;
583
584 // Specialized inverse.
585 Mat<3,3,conjugate<Real> > smallConjA3((conjugate<Real>*)cjdata);
586 Mat<3,3,conjugate<Real> >::TInvert smallConjAI3(smallConjA3.invert());
587
588 cout << "smallConjA3=" << smallConjA3 << " inv(smallConjA3)=" << smallConjAI3;
589
590 cout << "smallConjA3*inv(smallConjA3)=" << smallConjA3*smallConjAI3 << "NORM="
591 << (smallConjA3*smallConjAI3).norm() << endl;
592
593 cout << "inverse(smallNegA3)=" << inverse(smallConjA3) << endl;
594
595 cout << "norm(inverse-lapackInverse3)=" << (inverse(smallConjA3)-lapackInverse(smallConjA3)).norm() << endl;
596
597 cout << "det(smallConjA3)=" << det(smallConjA3)
598 << " det(inv(smallConjA3))=" << det(smallConjAI3) << endl;
599
600 cout << "Mat<1,1,conj>=" << smallConjA3.getSubMat<1,1>(1,0) << "inv(...)=" <<
601 smallConjA3.getSubMat<1,1>(1,0).invert();
602 cout << "Mat11*inv(Mat11)=" <<
603 smallConjA3.getSubMat<1,1>(1,0)*smallConjA3.getSubMat<1,1>(1,0).invert();
604 cout << "Mat<2,2,conj>=" << smallConjA3.getSubMat<2,2>(0,0) << "inv(...)=" <<
605 smallConjA3.getSubMat<2,2>(0,0).invert();
606 cout << "Mat22*inv(Mat22)=" <<
607 smallConjA3.getSubMat<2,2>(0,0)*smallConjA3.getSubMat<2,2>(0,0).invert();
608
609 try {
610 const double ddd[] = { 11, 12, 13, 14, 15, 16 };
611 const float fddd[] = { 11, 12, 13, 14, 15, 16 };
612 const complex<float> ccc[] = { complex<float>(1.,2.),
613 complex<float>(3.,4.),
614 complex<float>(5.,6.),
615 complex<float>(7.,8.) };
616 Vec< 2,complex<float> > cv2(ccc);
617 cout << "cv2=" << cv2 << endl;
618 cout << "(cv2+cv2)/1000:" << (cv2 + cv2) / complex<float>(1000,0)
619 << endl;
620
621 cout << "cv2: " << cv2 << endl;
622 cout << "cv2T: " << cv2.transpose() << endl;
623 cout << "-cv2: " << -cv2 << endl;
624 cout << "~cv2: " << ~cv2 << endl;
625 cout << "-~cv2: " << -(~cv2) << endl;
626 cout << "~-cv2: " << ~(-cv2) << endl;
627 cout << "~-cv2*10000: " << (~(-cv2))*10000.f << endl;
628
629 (~cv2)[1]=complex<float>(101.1f,202.3f);
630 cout << "after ~cv2[1]=(101.1f,202.3f), cv2= " << cv2 << endl;
631 -(~cv2)[1]=complex<float>(11.1f,22.3f);
632 cout << "after -~cv2[1]=(11.1f,22.3f), cv2= " << cv2 << endl;
633
634 Vec<3,float> dv2(fddd), ddv2(fddd+3);
635 dv2[2] = 1000;
636 cout << 100.* (ddv2 - dv2) / 1000. << endl;
637
638 cout << "dv2=" << dv2 << " dv2.norm()=" << dv2.norm() << endl;
639 cout << "cv2=" << cv2 << " cv2.norm()=" << cv2.norm() << endl;
640
641 dv2 = 100*dv2;
642
643 const Vec<3,float> v3c[] = {Vec<3,float>(fddd),Vec<3,float>(fddd+1)};
644
645 Vector_< Vec<2, Vec<3,float> > > vflt(2);
646 vflt[0] = Vec<2, Vec<3,float> >(v3c);
647 vflt[1] = vflt[0]*100;
648 vflt[1] = 100*vflt[0];
649 cout << "vflt 2xvec3=" << vflt << endl;
650
651
652 cout << "vflt.getNScalarsPerElement()=" << vflt.getNScalarsPerElement() << endl;
653 cout << "vflt.getPackedSizeofElement()=" << vflt.getPackedSizeofElement() << endl;
654 cout << "sizeof(vflt)=" << sizeof(vflt) << " sizeof(vflt[0])=" << sizeof(vflt[0]) << endl;
655 cout << "vflt.hasContiguousData()=" << vflt.hasContiguousData() << endl;
656 cout << "vflt.getContiguousScalarDataLength()=" << vflt.getContiguousScalarDataLength() << endl;
657
658 const float* p = vflt.getContiguousScalarData();
659 cout << "vflt's raw data: " << endl;
660 for (int i=0; i<vflt.getContiguousScalarDataLength(); ++i)
661 cout << " " << p[i];
662 cout << endl;
663
664 float* newData = new float[12];
665 float* oldData;
666 for (int i=0; i<12; ++i) newData[i]=(float)-i;
667 vflt.swapOwnedContiguousScalarData(newData, 12, oldData);
668
669 cout << "after data swap, vflt=" << vflt << endl;
670 cout << "old data =";
671 for (int i=0; i<12; ++i) cout << " " << oldData[i];
672 cout << endl;
673 delete[] oldData;
674
675 }
676 catch(const Exception::Base& b)
677 {
678 cout << b.getMessage() << endl;
679 }
680
681 typedef double P;
682 const int N = 1000;
683 const int LUP = 1;
684 Matrix_<P> big(N,N);
685 for (int j=0; j<N; ++j)
686 for (int i=0; i<N; ++i)
687 big(i,j) = 1+(float)std::rand()/RAND_MAX;
688 cout << "big.norm()=" << big.norm() << endl;
689
690
691 cout << "INVERTING " << LUP << "x" << N << "x" << N
692 << (sizeof(P)==4 ? std::string(" float") : std::string(" double")) << endl << std::flush;
693
694 Matrix_<P> flip(N,N), nrflip(N,N);
695
696 std::clock_t start = std::clock();
697 for (int i=0; i<LUP; ++i)
698 flip = big.invert();
699 cout << "... DONE -- " << (double)(std::clock()-start)/CLOCKS_PER_SEC << " seconds" << endl << std::flush;
700 cout << "initial norm=" << big.norm() << " invert.norm()=" << flip.norm() << endl << std::flush;
701
702 Matrix_<P> nrbig(N,N); nrbig=big;
703 cout << "NOW INVERT WITH NR, initial norm=" << nrbig.norm() << endl << std::flush;
704 start = std::clock();
705 for (int i=0; i<LUP; ++i) {
706 nrbig=big;
707 NR::luinvert(N, &nrbig(0,0), &nrflip(0,0));
708 }
709 cout << "... DONE -- " << (double)(std::clock()-start)/CLOCKS_PER_SEC << " seconds" << endl << std::flush;
710 cout << " nrinverse.norm()=" << nrflip.norm() << endl << std::flush;
711
712 cout << "big.norm()=" << big.norm() << " flip.norm()=" << flip.norm() << endl << std::flush;
713 Matrix_<P> ans(N,N);
714 cout << "Multiplying ..." << endl << std::flush;
715 start = std::clock();
716 for (int i=0; i<LUP; ++i)
717 Lapack::gemm('n','n',N,N,N,P(1),&big(0,0),N,&flip(0,0),N,P(0),&ans(0,0),N);
718 cout << "... DONE -- " << (double)(std::clock()-start)/CLOCKS_PER_SEC << " seconds" << endl << std::flush;
719 cout << "RMS (big*flip).norm() - 1=" << std::sqrt(ans.normSqr()/N)-1. << endl << std::flush;
720
721 cout << "Multiplying with NR result ..." << endl << std::flush;
722 start = std::clock();
723 for (int i=0; i<LUP; ++i)
724 Lapack::gemm('n','n',N,N,N,P(1),&big(0,0),N,&nrflip(0,0),N,P(0),&ans(0,0),N);
725 cout << "... DONE -- " << (double)(std::clock()-start)/CLOCKS_PER_SEC << " seconds" << endl << std::flush;
726 cout << "RMS (big*nrflip).norm() - 1=" << std::sqrt(ans.normSqr()/N)-1. << endl << std::flush;
727
728 cout << "Multiplying by hand ..." << endl << std::flush;
729 start = std::clock();
730 const P* bigp = &big(0,0);
731 const P* flipp = &flip(0,0);
732 P* ansp = &ans(0,0);
733 for (int l=0; l<LUP; ++l) {
734 for (int j=0; j<N; ++j) {
735 for (int i=0; i<N; ++i) {
736 P sum = P(0); const int jN=j*N;
737 for (int k=0; k<N; ++k)
738 sum += bigp[k*N+i] * flipp[jN+k];
739 ansp[jN+i] = sum;
740 }
741 }
742 }
743 cout << "... DONE -- " << (double)(std::clock()-start)/CLOCKS_PER_SEC << " seconds" << endl << std::flush;
744 cout << "RMS (big*flip).norm() - 1=" << std::sqrt(ans.normSqr()/N)-1. << endl << std::flush;
745 return 0;
746 }
747 catch (const std::exception& e) {
748 cout << e.what() << endl;
749 }
750 }
751
752 // Numerical Recipes version 2.11 LU decomp and inversion via backsolve.
753 // This is the C++ version modified to use column-ordered consecutive
754 // storage.
755
756 namespace NR {
757
758 // Return 1d index for column ordered matrix with leading dim N.
759 #define X(i,j) j*N+i
760
761 template <class DP>
luinvert(const int N,DP * a,DP * y)762 void luinvert(const int N, DP* a/*N,N*/, DP* y/*N,N*/) {
763 assert(a && y);
764 int* indx = new int[N];
765 DP d;
766 NR::ludcmp(N, a, indx, d);
767 for (int j=0; j<N; ++j) {
768 DP* col = &y[X(0,j)];
769 for (int i=0; i<N; ++i) col[i]=DP(0);
770 col[j] = DP(1);
771 NR::lubksb(N,a,indx,col); // writes directly into y
772 }
773
774 delete[] indx;
775 }
776
777 template <class DP>
lubksb(const int N,const DP * a,const int * indx,DP * b)778 void lubksb(const int N, const DP* a/*N,N*/, const int* indx/*N*/,
779 DP* b/*N*/)
780 {
781 int i,ii=0,ip,j;
782 DP sum;
783
784 for (i=0;i<N;i++) {
785 ip=indx[i];
786 sum=b[ip];
787 b[ip]=b[i];
788 if (ii != 0)
789 for (j=ii-1;j<i;j++) sum -= a[X(i,j)]*b[j];
790 else if (sum != 0.0)
791 ii=i+1;
792 b[i]=sum;
793 }
794 for (i=N-1;i>=0;i--) {
795 sum=b[i];
796 for (j=i+1;j<N;j++) sum -= a[X(i,j)]*b[j];
797 b[i]=sum/a[X(i,i)];
798 }
799 }
800
801 template <class DP>
ludcmp(const int N,DP * a,int * indx,DP & d)802 void ludcmp(const int N, DP* a/*N,N*/, int* indx/*N*/, DP &d)
803 {
804 const DP TINY=DP(1.0e-20);
805 int i,imax,j,k;
806 DP big,dum,sum,temp;
807
808 DP* vv = new DP[N];
809 d=DP(1);
810 for (i=0;i<N;i++) {
811 big=0.0;
812 for (j=0;j<N;j++)
813 if ((temp=fabs(a[X(i,j)])) > big) big=temp;
814 if (big == 0.0) {
815 std::cerr << "Singular matrix in routine ludcmp" << endl;
816 assert(false);
817 exit(1);
818 }
819 vv[i]=DP(1)/big;
820 }
821 for (j=0;j<N;j++) {
822 for (i=0;i<j;i++) {
823 sum=a[X(i,j)];
824 for (k=0;k<i;k++) sum -= a[X(i,k)]*a[X(k,j)];
825 a[X(i,j)]=sum;
826 }
827 big=0.0;
828 for (i=j;i<N;i++) {
829 sum=a[X(i,j)];
830 for (k=0;k<j;k++) sum -= a[X(i,k)]*a[X(k,j)];
831 a[X(i,j)]=sum;
832 if ((dum=vv[i]*fabs(sum)) >= big) {
833 big=dum;
834 imax=i;
835 }
836 }
837 if (j != imax) {
838 for (k=0;k<N;k++) {
839 dum=a[X(imax,k)];
840 a[X(imax,k)]=a[X(j,k)];
841 a[X(j,k)]=dum;
842 }
843 d = -d;
844 vv[imax]=vv[j];
845 }
846 indx[j]=imax;
847 if (a[X(j,j)] == 0.0) a[X(j,j)]=TINY;
848 if (j != N-1) {
849 dum=DP(1)/(a[X(j,j)]);
850 for (i=j+1;i<N;i++) a[X(i,j)] *= dum;
851 }
852 }
853
854 delete[] vv;
855 }
856
857
858 } // namespace NR
859
860