1
2 #include "TMV_Test.h"
3 #include "TMV_Test_1.h"
4 #include "TMV.h"
5 #include <fstream>
6 #include <cstdio>
7 #include <vector>
8
9 #include "TMV_TestVectorArith.h"
10
11 #define CT std::complex<T>
12
13 template <class T>
TestVectorReal()14 static void TestVectorReal()
15 {
16 const int N = 100;
17
18 if (showstartdone) {
19 std::cout<<"Start Test Real Vector"<<std::endl;
20 std::cout<<"T = "<<tmv::TMV_Text(T())<<std::endl;
21 std::cout<<"N = "<<N<<std::endl;
22 }
23
24 tmv::Vector<T> v(N);
25 Assert(v.size() == N,"Creating Vector");
26 for (int i=0; i<N; ++i) v(i) = T(i);
27 for (int i=0; i<N; ++i) Assert(v(i) == T(i),"Setting Vector");
28
29 tmv::VectorView<T> v2 = v.subVector(0,N,2);
30 for (int i=0; i<N/2; ++i)
31 Assert(v2(i) == T(2*i),"Reading Vector with stride = 2");
32
33 for (int i=0; i<N/2; ++i) v2[i] = T(i + 1000);
34 for (int i=0; i<N/2; ++i)
35 Assert(v(2*i) == T(i+1000),"Writing Vector with stride = 2");
36
37 tmv::Vector<T> v3 = v2;
38 for (int i=0; i<N/2; ++i)
39 Assert(v3[i] == v2[i],"Copying Vector with stride = 2");
40
41 for (int i=0; i<N; ++i) v[i] = T(i);
42 v.swap(2,5);
43 Assert(v(2) == T(5) && v(5) == T(2),"Swapping elements of Vector");
44 v.swap(2,5);
45 Assert(v(2) == T(2) && v(5) == T(5),"Swapping elements of Vector");
46
47 T sum = N*(N-1)/2;
48 if (showacc) {
49 std::cout<<"SumElements = "<<SumElements(v)<<std::endl;
50 std::cout<<"expected "<<sum<<std::endl;
51 std::cout<<"== "<<(SumElements(v)==sum)<<std::endl;
52 std::cout<<"diff = "<<(SumElements(v)-sum)<<std::endl;
53 }
54 Assert(SumElements(v) == sum,"Vector SumElements(v)");
55
56 v.reverseSelf();
57 for (int i=0; i<N; ++i) Assert(v(i) == T(N-i-1),"Reversing Vector");
58
59 for (int i=0; i<N; ++i) v(i) = T(i+10);
60 v(23) = T(10*N);
61 v(42) = T(0.25);
62 v(15) = T(-20*N);
63 int imax,imin;
64 if (showacc) {
65 std::cout<<"v = "<<v<<std::endl;
66 std::cout<<"v.MaxAbs = "<<v.maxAbsElement(&imax)<<std::endl;
67 std::cout<<"imax = "<<imax<<std::endl;
68 std::cout<<"v.MinAbs = "<<v.minAbsElement(&imin)<<std::endl;
69 std::cout<<"imin = "<<imin<<std::endl;
70 }
71 Assert(Equal2(v.maxAbsElement(&imax),T(20*N),EPS),
72 "MaxAbsElement of Vector did not return correct value");
73 Assert(imax == 15,
74 "MaxAbsElement of Vector did not return correct index");
75 Assert(Equal2(v.minAbsElement(&imin),T(0.25),EPS),
76 "MinAbsElement of Vector did not return correct value");
77 Assert(imin == 42,
78 "MinAbsElement of Vector did not return correct index");
79 Assert(Equal2(v.maxElement(&imax),T(10*N),EPS),
80 "MaxElement of Vector did not return correct value");
81 Assert(imax == 23,
82 "MaxElement of Vector did not return correct index");
83 Assert(Equal2(v.minElement(&imin),T(-20*N),EPS),
84 "MinElement of Vector did not return correct value");
85 Assert(imin == 15,
86 "MinElement of Vector did not return correct index");
87 Assert(Equal2(v.maxAbs2Element(&imax),T(20*N),EPS),
88 "MaxAbs2Element of Vector did not return correct value");
89 Assert(imax == 15,
90 "MaxAbs2Element of Vector did not return correct index");
91 Assert(Equal2(v.minAbs2Element(&imin),T(0.25),EPS),
92 "MinAbs2Element of Vector did not return correct value");
93 Assert(imin == 42,
94 "MinAbs2Element of Vector did not return correct index");
95
96 tmv::Vector<T> a(N);
97 tmv::Vector<T> b(N);
98 for (int i=0; i<N; ++i) a(i) = T(3+i);
99
100 b = a;
101 for (int i=0; i<N; ++i) Assert(a(i) == b(i),"Vector1 = Vector2");
102
103 Assert(a == b,"Testing Equality of Vectors");
104
105 b(4) = T(0);
106 Assert(a != b,"Vector = Vector copied address, not values");
107
108 tmv::Vector<T,tmv::FortranStyle> af(N);
109 for (int i=1; i<=N; ++i) {
110 af(i) = T(3+i-1);
111 }
112 for (int i=1; i<=N; ++i)
113 Assert(af(i) == a(i-1),"FortranStyle Vector access");
114 tmv::ConstVectorView<T,tmv::FortranStyle> afcv = af.view();
115 for (int i=1; i<=N; ++i)
116 Assert(afcv(i) == a(i-1),"FortranStyle Vector CV access");
117 tmv::VectorView<T,tmv::FortranStyle> afv = af.view();
118 for (int i=1; i<=N; ++i)
119 Assert(afv(i) == a(i-1),"FortranStyle Vector V access");
120 Assert(a == af,"FortransStyle Vector == CStyle Vector");
121 tmv::ConstVectorView<T> afcv_c = afcv;
122 Assert(afcv_c == a,"CStyle View of FortransStyle Vector == CStyle Vector");
123 Assert(afcv == a,"FortranStyle View of Vector == CStyle Vector");
124
125 // Test assignments and constructors from arrays
126 std::vector<T> qv(6);
127 T qvar[] = { T(8), T(6), T(4), T(2), T(0), T(-2) };
128 T qvar2[] = { T(8), T(7), T(6), T(5), T(4), T(3), T(2), T(1), T(0), T(-1), T(-2) };
129 for(int i=0;i<6;++i) qv[i] = qvar[i];
130
131 // Construct from C array
132 tmv::Vector<T> q1(6);
133 std::copy(qvar,qvar+6,q1.begin());
134 // Construct from vector
135 tmv::Vector<T> q2(6);
136 std::copy(qv.begin(),qv.end(),q2.begin());
137 // Assign VectorView from vector
138 tmv::Vector<T> q3x(50);
139 tmv::VectorView<T> q3 = q3x.subVector(4,34,5);
140 std::copy(qv.begin(),qv.end(),q3.begin());
141 // Use op<< assignment
142 tmv::Vector<T> q4(6);
143 q4 << 8, 6, 4, 2, 0, -2;
144 // Directly view an array
145 tmv::VectorView<T> q5 = tmv::VectorViewOf(qvar,6);
146 // Directly view an array with step = 2
147 tmv::VectorView<T> q6 = tmv::VectorViewOf(qvar2,6,2);
148
149 if (showacc) {
150 std::cout<<"q1 = "<<q1<<std::endl;
151 std::cout<<"q2 = "<<q2<<std::endl;
152 std::cout<<"q3 = "<<q3<<std::endl;
153 std::cout<<"q4 = "<<q4<<std::endl;
154 std::cout<<"q5 = "<<q5<<std::endl;
155 std::cout<<"q6 = "<<q6<<std::endl;
156 }
157
158 for(int i=0;i<6;++i) {
159 Assert(q1(i) == T(8-2*i),"Create Vector from T*");
160 Assert(q2(i) == T(8-2*i),"Create Vector from std::vector");
161 Assert(q3(i) == T(8-2*i),"Create VectorView from std::vector");
162 Assert(q4(i) == T(8-2*i),"Create Vector from << list");
163 Assert(q5(i) == T(8-2*i),"Create VectorView of T* ");
164 Assert(q6(i) == T(8-2*i),"Create VectorView of T* with step");
165 }
166
167 // Test Basic Arithmetic
168 for (int i=0; i<N; ++i) b(i) = T(5+2*i);
169
170 v = a+b;
171 for (int i=0; i<N; ++i) Assert(v(i) == T(8+3*i),"Adding Vectors");
172
173 v = a-b;
174 if (showacc) {
175 std::cout<<"a = "<<a<<std::endl;
176 std::cout<<"b = "<<b<<std::endl;
177 std::cout<<"a-b = "<<v<<std::endl;
178 }
179 for (int i=0; i<N; ++i) Assert(v(i) == T(-2-i),"Subtracting Vectors");
180
181 // b(i) = 5+2i
182 // a(i) = 3+i
183 // a(i)*b(i) = 15+11i+2i^2
184 // Sum = 15N + 11N(N-1)/2 + 2*N*(N-1)*(2N-1)/6
185 T prod = 15*N + 11*N*(N-1)/2 + 2*N*(N-1)*(2*N-1)/6;
186 if (showacc) {
187 std::cout<<"a = "<<a<<std::endl;
188 std::cout<<"b = "<<b<<std::endl;
189 std::cout<<"a*b = "<<a*b<<std::endl;
190 std::cout<<"prod = "<<prod<<std::endl;
191 }
192 Assert(a*b == prod,"Multiplying Vectors");
193
194 tmv::Vector<T> c(5);
195 c = v.subVector(10,70,12);
196 for (int i=0; i<5; ++i) Assert(c(i) == v(10+12*i),"SubVector");
197
198 for(int i=0;i<N;++i) a(i) = T(i+10);
199 for(int i=0;i<N;++i) b(i) = T(-3*i+191);
200
201 prod = 2900;
202 T normsqsum = 1373700;
203 T normsqdiff = 1362100;
204 if (showacc) {
205 std::cout<<"a*b = "<<a*b<<std::endl;
206 std::cout<<"expected prod = "<<prod<<std::endl;
207 std::cout<<"abs(diff) = "<<tmv::TMV_ABS(a*b-prod)<<std::endl;
208 std::cout<<"eps = "<<EPS<<std::endl;
209 std::cout<<"a+b = "<<a+b<<std::endl;
210 std::cout<<"NormSq(a+b) = "<<NormSq(a+b)<<std::endl;
211 std::cout<<"expected normsqsum = "<<normsqsum<<std::endl;
212 std::cout<<"abs(diff) = "<<tmv::TMV_ABS(NormSq(a+b)-normsqsum)<<std::endl;
213 std::cout<<"eps = "<<EPS*tmv::TMV_ABS(Norm1(a)+Norm1(b))<<std::endl;
214 std::cout<<"Norm1(a) = "<<Norm1(a)<<std::endl;
215 std::cout<<"Norm1(b) = "<<Norm1(b)<<std::endl;
216 }
217 T eps = EPS;
218 if (!std::numeric_limits<T>::is_integer) eps *= Norm(a) * Norm(b);
219 Assert(Equal2(a*b,prod,eps),"Inner Product");
220 T eps2 = EPS * tmv::TMV_ABS2(Norm1(a)+Norm1(b));
221 Assert(Equal2(NormSq(a+b),normsqsum,eps2),"Vector Sum");
222 Assert(Equal2(NormSq(a-b),normsqdiff,eps2),"Vector Diff");
223
224 const int NN=20;
225 tmv::Vector<T> w(NN);
226 w << 33,12,54,-12,43,-94,0,-20,40,-115,
227 -120,140,330,10,-93,-39,49,100,-310,1;
228
229 tmv::Vector<T> origw = w;
230 tmv::Vector<T> w2 = w;
231 tmv::Permutation P(NN);
232 if (showacc) std::cout<<"unsorted w = "<<w<<std::endl;
233
234 w.sort(P);
235 if (showacc) std::cout<<"sorted w = "<<w<<std::endl;
236 for(int i=1;i<NN;++i) {
237 Assert(w(i-1) <= w(i),"Sort real Vector");
238 }
239 w2 = P * w2;
240 Assert(w2 == w,"Sort real Vector -- perm");
241 w = origw;
242 w.sort();
243 Assert(w2 == w,"Sort real Vector -- without perm");
244 w = P.inverse() * w;
245 if (showacc) std::cout<<"reverse permute sorted Vector = "<<w<<std::endl;
246 Assert(w == origw,"Reverse permute sorted Vector = orig");
247 w2 = origw;
248
249 w.sort(P,tmv::Ascend,tmv::AbsComp);
250 if (showacc) std::cout<<"sorted w abs = "<<w<<std::endl;
251 for(int i=1;i<NN;++i) {
252 Assert(tmv::TMV_ABS(w(i-1)) <= tmv::TMV_ABS(w(i)),
253 "Sort real Vector abs");
254 }
255 w2 = P * w2;
256 if (showacc) std::cout<<"permuted w2 = "<<w<<std::endl;
257 Assert(w2 == w,"Sort real Vector abs -- perm");
258 w = origw;
259 w.sort(tmv::Ascend,tmv::AbsComp);
260 if (showacc)
261 std::cout<<"sorted w abs (without perm) = "<<w<<std::endl;
262 Assert(w2 == w,"Sort real Vector abs -- without perm");
263 w = w2 = origw;
264
265 w.sort(P,tmv::Descend);
266 if (showacc) std::cout<<"sorted w desc = "<<w<<std::endl;
267 for(int i=1;i<NN;++i) {
268 Assert(w(i-1) >= w(i),"Sort real Vector desc");
269 }
270 w2 = P * w2;
271 Assert(w2 == w,"Sort real Vector desc -- perm");
272 w = origw;
273 w.sort(tmv::Descend);
274 Assert(w2 == w,"Sort real Vector desc -- without perm");
275 w = w2 = origw;
276
277 w.sort(P,tmv::Descend,tmv::AbsComp);
278 if (showacc) std::cout<<"sorted w desc abs = "<<w<<std::endl;
279 for(int i=1;i<NN;++i) {
280 Assert(tmv::TMV_ABS(w(i-1)) >= tmv::TMV_ABS(w(i)),
281 "Sort real Vector desc abs");
282 }
283 w2 = P * w2;
284 Assert(w2 == w,"Sort real Vector desc abs -- perm");
285 w = origw;
286 w.sort(tmv::Descend,tmv::AbsComp);
287 Assert(w2 == w,"Sort real Vector desc abs -- without perm");
288 w = w2 = origw;
289
290 v.resize(2);
291 Assert(v.size() == 2,"v.resize(2)");
292 for (int i=0; i<2; ++i) v(i) = T(i);
293 for (int i=0; i<2; ++i) Assert(v(i) == T(i),"Setting resized Vector");
294
295 v.resize(2*N);
296 Assert(int(v.size()) == 2*N,"v.resize(2*N)");
297 for (int i=0; i<2*N; ++i) v(i) = T(i);
298 for (int i=0; i<2*N; ++i) Assert(v(i) == T(i),"Setting resized Vector");
299
300 if (showstartdone) {
301 std::cout<<"Done Test Real Vector"<<std::endl;
302 }
303 }
304
305 template <class T>
TestVectorComplex()306 static void TestVectorComplex()
307 {
308 const int N = 100;
309
310 if (showstartdone) {
311 std::cout<<"Start Test Complex Vector"<<std::endl;
312 std::cout<<"T = "<<tmv::TMV_Text(T())<<std::endl;
313 std::cout<<"N = "<<N<<std::endl;
314 }
315
316 tmv::Vector<CT > v(N);
317 for (int i=0; i<N; ++i) v(i) = CT(T(i),T(i+1234));
318
319 for (int i=0; i<N; ++i)
320 Assert(v(i).real() == T(i), "CVector set");
321 for (int i=0; i<N; ++i)
322 Assert(v(i).imag() == T(i+1234), "CVector set");
323
324 tmv::VectorView<CT > v1(v.subVector(0,N,2));
325 for (int i=0; i<N/2; ++i)
326 Assert(v1(i) == CT(T(2*i),T(2*i+1234)),
327 "CVector stride=2");
328
329 for (int i=0; i<N/2; ++i) v1[i] = CT(T(i),T(i+1234));
330 for (int i=0; i<N/2; ++i)
331 Assert(v[2*i] == CT(T(i),T(i+1234)),
332 "setting CVector with stride = 2");
333
334 for (int i=0; i<N; ++i) v(i) = CT(T(i),T(i+1234));
335
336 v.swap(2,5);
337 Assert(v[2] == CT(5,5+1234),"Swap in CVector");
338 Assert(v[5] == CT(2,2+1234),"Swap in CVector");
339 v.swap(2,5);
340
341 tmv::Vector<CT > v2 = v.conjugate();
342
343 for (int i=0; i<N; ++i)
344 Assert(v2(i) == CT(T(i),T(-i-1234)), "Conjugate CVector");
345 Assert(v2 == v.conjugate(),"Conjugate == CVector");
346
347 tmv::Vector<CT > v3(N);
348 for (int i=0; i<N; ++i) v3(i) = CT(i+10,2*i);
349 v3(23) = CT(40*N,9*N);
350 v3(42) = CT(0,1);
351 v3(15) = CT(-32*N,24*N);
352 int imax,imin;
353
354 if (!std::numeric_limits<T>::is_integer) {
355 if (showacc) {
356 std::cout<<"v = "<<v3<<std::endl;
357 std::cout<<"v.MaxAbs = "<<v3.maxAbsElement(&imax)<<std::endl;
358 std::cout<<"imax = "<<imax<<std::endl;
359 std::cout<<"v.MinAbs = "<<v3.minAbsElement(&imin)<<std::endl;
360 std::cout<<"imin = "<<imin<<std::endl;
361 }
362 Assert(Equal2(v3.maxAbsElement(&imax),T(41*N),EPS),
363 "MaxAbsElement of Vector did not return correct value");
364 Assert(imax == 23,
365 "MaxAbsElement of Vector did not return correct index");
366 Assert(Equal2(v3.minAbsElement(&imin),T(1),EPS),
367 "MinAbsElement of Vector did not return correct value");
368 Assert(imin == 42,
369 "MinAbsElement of Vector did not return correct index");
370 }
371 Assert(Equal2(v3.maxAbs2Element(&imax),T(56*N),EPS),
372 "MaxAbs2Element of complex Vector did not return correct value");
373 Assert(imax == 15,
374 "MaxAbs2Element of complex Vector did not return correct index");
375 Assert(Equal2(v3.minAbs2Element(&imin),T(1),EPS),
376 "MinAbs2Element of complex Vector did not return correct value");
377 Assert(imin == 42,
378 "MinAbs2Element of complex Vector did not return correct index");
379
380 CT prod_act(0);
381 for (int i=0; i<N; ++i) prod_act += v[i] * v2[i];
382 CT prod = v*v2;
383 Assert(Equal2(prod,prod_act,EPS*tmv::TMV_ABS2(prod_act)),
384 "CVector * CVector");
385 prod = v*v.conjugate();
386 prod_act = T(0);
387 for (int i=0; i<N; ++i) prod_act += v[i] * std::conj(v[i]);
388 Assert(Equal2(prod.imag(),T(0),EPS),"prod is real");
389 Assert(Equal2(prod,prod_act,EPS*tmv::TMV_ABS2(prod_act)),
390 "CVector * conj(CVector)");
391
392 if (!std::numeric_limits<T>::is_integer) {
393 T norm1 = tmv::TMV_SQRT(prod.real());
394 T norm2 = Norm(v);
395 if (showacc) {
396 std::cout<<"v = "<<v<<std::endl;
397 std::cout<<"v2 = "<<v2<<std::endl;
398 std::cout<<"v*v2 = "<<v*v2<<std::endl;
399 std::cout<<"norm1 = "<<norm1<<std::endl;
400 std::cout<<"norm2 = "<<norm2<<std::endl;
401 }
402 Assert(Equal2(norm1,norm2,EPS*norm1),"Norm CVector");
403 }
404
405 CT sum_act(0);
406 for (int i=0; i<N; ++i) sum_act += v[i];
407 CT sumel = v.sumElements();
408 if (showacc) {
409 std::cout<<"sumel = "<<sumel<<std::endl;
410 std::cout<<"sumact = "<<sum_act<<std::endl;
411 std::cout<<"diff = "<<tmv::TMV_ABS(sumel-sum_act)<<std::endl;
412 }
413 Assert(Equal2(sumel,sum_act,EPS*tmv::TMV_ABS2(sum_act)),
414 "CVector SumElements");
415
416 if (!std::numeric_limits<T>::is_integer) {
417 T sumabs_act(0);
418 for (int i=0; i<N; ++i) sumabs_act += tmv::TMV_ABS(v[i]);
419 T sumabsel = v.sumAbsElements();
420 if (showacc) {
421 std::cout<<"sumabsel = "<<sumabsel<<std::endl;
422 std::cout<<"sumabs_act = "<<sumabs_act<<std::endl;
423 std::cout<<"diff = "<<tmv::TMV_ABS(sumabsel-sumabs_act)<<std::endl;
424 }
425 Assert(Equal2(sumabsel,sumabs_act,EPS*tmv::TMV_ABS2(sumabs_act)),
426 "CVector SumAbsElements");
427 }
428 T sumabs2_act(0);
429 for (int i=0; i<N; ++i) sumabs2_act += tmv::TMV_ABS2(v[i]);
430 T sumabs2el = v.sumAbs2Elements();
431 Assert(Equal2(sumabs2el,sumabs2_act,EPS*tmv::TMV_ABS2(sumabs2_act)),
432 "CVector SumAbs2Elements");
433
434 v.conjugateSelf();
435 Assert(v == v2,"ConjugateSelf CVector");
436 v = v.conjugate();
437 Assert(v == v2.conjugate(),"v = v.conjugate() CVector");
438
439 tmv::Vector<T> a(N);
440 for(int i=0;i<N;++i) a(i) = T(i+10);
441 tmv::Vector<T> b(N);
442 for(int i=0;i<N;++i) b(i) = T(-3*i+191);
443
444 tmv::Vector<CT > ca = a;
445 Assert(Equal(ca,a,EPS),"Copy real V -> complex V");
446
447 ca *= CT(3,4);
448 tmv::Vector<CT > cb = b*CT(3,4);
449
450 prod = T(29)*T(25)*CT(-28,96);
451 T normsqsum = 34342500;
452 T normsqdiff = 34052500;
453 T eps = EPS;
454 if (!std::numeric_limits<T>::is_integer) eps *= Norm(ca) * Norm(cb);
455 if (showacc) {
456 std::cout<<"ca*cb = "<<ca*cb<<std::endl;
457 std::cout<<"expected prod = "<<prod<<std::endl;
458 std::cout<<"abs(diff) = "<<tmv::TMV_ABS(ca*cb-prod)<<std::endl;
459 std::cout<<"eps = "<<eps<<std::endl;
460 }
461 Assert(Equal2(ca*cb,prod,eps),"CInner Product");
462 T eps2 = EPS;
463 if (!std::numeric_limits<T>::is_integer) eps2 *= Norm1(ca) * Norm1(cb);
464 Assert(Equal2(NormSq(ca+cb),normsqsum,eps2),"CVector Sum");
465 Assert(Equal2(NormSq(ca-cb),normsqdiff,eps2),"CVector Diff");
466
467 const int NN=20;
468 tmv::Vector<CT > w(NN);
469 w << 33,12,54,-12,43,-94,0,-20,40,-115,
470 -120,140,330,10,-93,-39,49,100,-310,1;
471
472 tmv::Vector<T> iw(NN);
473 iw << 14,98,-2,-86,30,-44,30,90,-19,-114,
474 111,-1400,-230,110,52,-39,48,990,-710,-5;
475 w.imagPart() = iw;
476
477 tmv::Vector<CT > origw = w;
478 tmv::Vector<CT > w2 = w;
479 tmv::Permutation P(NN);
480 if (showacc) std::cout<<"unsorted w = "<<w<<std::endl;
481
482 w.sort(P);
483 if (showacc) std::cout<<"sorted w = "<<w<<std::endl;
484 for(int i=1;i<NN;++i) {
485 Assert(real(w(i-1)) <= real(w(i)),"Sort complex Vector");
486 }
487 w2 = P * w2;
488 Assert(w2 == w,"Sort complex Vector -- perm");
489 w = origw;
490 w.sort();
491 Assert(w2 == w,"Sort complex Vector -- without perm");
492 w = P.inverse() * w;
493 if (showacc) std::cout<<"reverse permute sorted Vector = "<<w<<std::endl;
494 Assert(w == origw,"Reverse permute sorted Vector = orig");
495 w2 = origw;
496
497 if (!std::numeric_limits<T>::is_integer) {
498 w.sort(P,tmv::Ascend,tmv::AbsComp);
499 if (showacc) std::cout<<"sorted w abs = "<<w<<std::endl;
500 for(int i=1;i<NN;++i) {
501 Assert(tmv::TMV_ABS(w(i-1)) <= tmv::TMV_ABS(w(i)),
502 "Sort complex Vector abs");
503 }
504 w2 = P * w2;
505 Assert(w2 == w,"Sort complex Vector abs -- perm");
506 w = origw;
507 w.sort(tmv::Ascend,tmv::AbsComp);
508 Assert(w2 == w,"Sort complex Vector abs -- without perm");
509 w = w2 = origw;
510
511 w.sort(P,tmv::Ascend,tmv::ArgComp);
512 if (showacc) std::cout<<"sorted w arg = "<<w<<std::endl;
513 for(int i=1;i<NN;++i) {
514 Assert(tmv::TMV_ARG(w(i-1)) <= tmv::TMV_ARG(w(i)),
515 "Sort complex Vector arg");
516 }
517 w2 = P * w2;
518 Assert(w2 == w,"Sort complex Vector arg -- perm");
519 w = origw;
520 w.sort(tmv::Ascend,tmv::ArgComp);
521 Assert(w2 == w,"Sort complex Vector arg -- without perm");
522 w = w2 = origw;
523 }
524
525 w.sort(P,tmv::Ascend,tmv::ImagComp);
526 if (showacc) std::cout<<"sorted w imag = "<<w<<std::endl;
527 for(int i=1;i<NN;++i) {
528 Assert(imag(w(i-1)) <= imag(w(i)),"Sort complex Vector imag");
529 }
530 w2 = P * w2;
531 Assert(w2 == w,"Sort complex Vector imag -- perm");
532 w = origw;
533 w.sort(tmv::Ascend,tmv::ImagComp);
534 Assert(w2 == w,"Sort complex Vector imag -- without perm");
535 w = w2 = origw;
536
537 w.sort(P,tmv::Descend);
538 if (showacc) std::cout<<"sorted w desc = "<<w<<std::endl;
539 for(int i=1;i<NN;++i) {
540 Assert(real(w(i-1)) >= real(w(i)),"Sort complex Vector desc");
541 }
542 w2 = P * w2;
543 Assert(w2 == w,"Sort complex Vector desc -- perm");
544 w = origw;
545 w.sort(tmv::Descend);
546 Assert(w2 == w,"Sort complex Vector desc -- without perm");
547 w = w2 = origw;
548
549 if (!std::numeric_limits<T>::is_integer) {
550 w.sort(P,tmv::Descend,tmv::AbsComp);
551 if (showacc) std::cout<<"sorted w desc abs = "<<w<<std::endl;
552 for(int i=1;i<NN;++i) {
553 Assert(tmv::TMV_ABS(w(i-1)) >= tmv::TMV_ABS(w(i)),
554 "Sort complex Vector desc abs");
555 }
556 w2 = P * w2;
557 Assert(w2 == w,"Sort complex Vector desc abs -- perm");
558 w = origw;
559 w.sort(tmv::Descend,tmv::AbsComp);
560 Assert(w2 == w,"Sort complex Vector desc abs -- without perm");
561 w = w2 = origw;
562
563 w.sort(P,tmv::Descend,tmv::ArgComp);
564 if (showacc) std::cout<<"sorted w desc arg = "<<w<<std::endl;
565 for(int i=1;i<NN;++i) {
566 Assert(tmv::TMV_ARG(w(i-1)) >= tmv::TMV_ARG(w(i)),
567 "Sort complex Vector desc arg");
568 }
569 w2 = P * w2;
570 Assert(w2 == w,"Sort complex Vector desc arg -- perm");
571 w = origw;
572 w.sort(tmv::Descend,tmv::ArgComp);
573 Assert(w2 == w,"Sort complex Vector desc arg -- without perm");
574 w = w2 = origw;
575 }
576
577 w.sort(P,tmv::Descend,tmv::ImagComp);
578 if (showacc)
579 std::cout<<"sorted w desc imag = "<<w<<std::endl;
580 for(int i=1;i<NN;++i) {
581 Assert(imag(w(i-1)) >= imag(w(i)),"Sort complex Vector desc imag");
582 }
583 w2 = P * w2;
584 Assert(w2 == w,"Sort complex Vector desc imag -- perm");
585 w = origw;
586 w.sort(tmv::Descend,tmv::ImagComp);
587 Assert(w2 == w,"Sort complex Vector desc imag -- without perm");
588 w = w2 = origw;
589
590 if (showstartdone) {
591 std::cout<<"Done Test Complex Vector"<<std::endl;
592 }
593 }
594
595 template <class T>
TestVectorArith()596 static void TestVectorArith()
597 {
598 typedef tmv::VectorView<T> V;
599 typedef tmv::VectorView<CT > CV;
600
601 const int N = 100;
602
603 if (showstartdone) {
604 std::cout<<"Start Test Vector Arith"<<std::endl;
605 std::cout<<"T = "<<tmv::TMV_Text(T())<<std::endl;
606 std::cout<<"N = "<<N<<std::endl;
607 }
608
609 tmv::Vector<T> a(N);
610 for(int i=0;i<N;++i) a(i) = T(i+10);
611 tmv::Vector<CT > ca = a*CT(2,-1);;
612 V aa = a.view();
613 CV caa = ca.view();
614 TestVectorArith1(aa,caa,"Vector C");
615
616 tmv::Vector<T> b(N);
617 for(int i=0;i<N;++i) b(i) = T(-3*i+2);
618 tmv::Vector<CT > cb = b*CT(-5,1);
619 V bb = b.view();
620 CV cbb = cb.view();
621 TestVectorArith2(aa,caa,b,cbb,"Vector CC");
622
623 tmv::Vector<T> a10(10*N);
624 V as = a10.subVector(0,10*N,10);
625 as = a;
626 tmv::Vector<CT > ca10(10*N);
627 CV cas = ca10.subVector(0,10*N,10);
628 cas = ca;
629 TestVectorArith1(as,cas,"Vector C Step");
630
631 tmv::Vector<T> b10(10*N);
632 V bs = b10.subVector(0,10*N,10);
633 bs = b;
634 tmv::Vector<CT > cb10(10*N);
635 CV cbs = cb10.subVector(0,10*N,10);
636 cbs = cb;
637 TestVectorArith2(as,cas,bb,cbb,"Vector C StepA");
638 TestVectorArith2(aa,caa,bs,cbs,"Vector C StepB");
639 TestVectorArith2(as,cas,bs,cbs,"Vector C StepAB");
640
641 V ar = aa.reverse();
642 CV car = caa.reverse();
643 TestVectorArith1(ar,car,"Vector C Rev");
644
645 V br = bb.reverse();
646 CV cbr = cbb.reverse();
647 TestVectorArith2(ar,car,bb,cbb,"Vector C RevA");
648 TestVectorArith2(aa,caa,br,cbr,"Vector C RevB");
649 TestVectorArith2(ar,car,br,cbr,"Vector C RevAB");
650
651 #if (XTEST & 32)
652 typedef tmv::VectorView<T,tmv::FortranStyle> VF;
653 typedef tmv::VectorView<CT,tmv::FortranStyle> CVF;
654 VF af = a.fView();
655 CVF caf = ca.fView();
656 TestVectorArith1(af,caf,"Vector F");
657
658 VF bf = b.fView();
659 CVF cbf = cb.fView();
660 TestVectorArith2(af,caf,b,cbb,"Vector FC");
661 TestVectorArith2(aa,caa,bf,cbf,"Vector CF");
662 TestVectorArith2(af,caf,bf,cbf,"Vector FF");
663
664 VF asf = as;
665 CVF casf = cas;
666 TestVectorArith1(asf,casf,"Vector F Step");
667
668 VF bsf = bs;
669 CVF cbsf = cbs;
670 TestVectorArith2(asf,casf,bf,cbf,"Vector F StepA");
671 TestVectorArith2(af,caf,bsf,cbsf,"Vector F StepB");
672 TestVectorArith2(asf,casf,bsf,cbsf,"Vector F StepAB");
673
674 VF arf = ar;
675 VF brf = br;
676 CVF carf = car;
677 CVF cbrf = cbr;
678 TestVectorArith1(arf,carf,"Vector F Rev");
679 TestVectorArith2(arf,carf,bf,cbf,"Vector F RevA");
680 TestVectorArith2(af,caf,brf,cbrf,"Vector F RevB");
681 TestVectorArith2(arf,carf,brf,cbrf,"Vector F RevAB");
682 #endif
683
684 if (showstartdone) {
685 std::cout<<"Done Test Vector Arith"<<std::endl;
686 }
687 }
688
689 template <class T>
TestVectorIO()690 static void TestVectorIO()
691 {
692 const int N = 20;
693
694 if (showstartdone) {
695 std::cout<<"Start Test Vector I/O"<<std::endl;
696 std::cout<<"T = "<<tmv::TMV_Text(T())<<std::endl;
697 std::cout<<"N = "<<N<<std::endl;
698 }
699
700 tmv::Vector<T> v(N);
701 tmv::Vector<CT > cv(N);
702
703 for (int i=0; i<N; ++i) {
704 v(i) = T(i+34);
705 cv(i) = CT(T(i),T(N-i));
706 }
707 v(3) = T(1.e-30);
708 cv(3) = CT(T(1.e-30),T(1.e-30));
709 v(8) = T(9.e-3);
710 cv(8) = CT(T(9.e-3),T(9.e-3));
711 cv(9) = CT(T(9),T(9.e-3));
712 v(12) = T(0.123456789);
713 cv(12) = CT(T(3.123456789),T(6.987654321));
714
715 // First check clipping function...
716 tmv::Vector<T> v2 = v;
717 tmv::Vector<CT > cv2 = cv;
718 if (!std::numeric_limits<T>::is_integer) {
719 v2.clip(T(1.e-2));
720 cv2.clip(T(1.e-2));
721 }
722 tmv::Vector<T> v3 = v;
723 tmv::Vector<CT > cv3 = cv;
724 v3(3) = T(0);
725 cv3(3) = T(0);
726 v3(8) = T(0);
727 // Others, esp. cv3(8) and cv3(9) shouldn't get clipped.
728 Assert(v2 == v3,"Vector clip");
729 Assert(cv2 == cv3,"Complex Vector clip");
730
731 // However, ThreshIO for complex works slightly differently than clip.
732 // It clips _either_ the real or imag component, so now cv2(8) and
733 // cv2(9) need to be modified.
734 cv2(8) = cv3(8) = T(0);
735 cv2(9) = cv3(9) = T(9);
736
737 // Write vectors with 4 different styles
738 std::ofstream fout("tmvtest_vector_io.dat");
739 Assert(bool(fout),"Couldn't open tmvtest_vector_io.dat for output");
740 fout << v << std::endl;
741 fout << cv << std::endl;
742 fout << tmv::CompactIO() << v << std::endl;
743 fout << tmv::CompactIO() << cv << std::endl;
744 fout << tmv::ThreshIO(1.e-2).setPrecision(12) << v << std::endl;
745 fout << tmv::ThreshIO(1.e-2).setPrecision(12) << cv << std::endl;
746 // Not a very pretty IO style, but it tests being able to read
747 // a style that has no whitespace and has more than one
748 // character for some of the markup elements.
749 tmv::IOStyle myStyle =
750 tmv::CompactIO().setThresh(1.e-2).setPrecision(4).
751 markup("Start","[",",","]","---","Done");
752 fout << myStyle << v << std::endl;
753 fout << myStyle << cv << std::endl;
754 fout.close();
755
756 // When using (the default) prec(6), these will be the values read in.
757 v(12) = T(0.123457);
758 cv(12) = CT(T(3.12346),T(6.98765));
759
760 // When using prec(12), the full correct values will be read in. (v2,cv2)
761
762 // When using prec(4), these will be the values read in.
763 v3(12) = T(0.1235);
764 if (std::numeric_limits<T>::is_integer) cv3(12) = CT(3,6);
765 else cv3(12) = CT(T(3.123),T(6.988));
766
767 // Read them back in
768 tmv::Vector<T> xv(N);
769 tmv::Vector<CT > xcv(N);
770 std::ifstream fin("tmvtest_vector_io.dat");
771 Assert(bool(fin),"Couldn't open tmvtest_vector_io.dat for input");
772 fin >> xv >> xcv;
773 if (showacc) {
774 std::cout<<"v = "<<v<<std::endl;
775 std::cout<<"xv = "<<xv<<std::endl;
776 std::cout<<"xcv = "<<xcv<<std::endl;
777 std::cout<<"v-xv = "<<(v-xv)<<std::endl;
778 std::cout<<"cv-xcv = "<<(cv-xcv)<<std::endl;
779 }
780 Assert(EqualIO(v,xv,EPS),"Vector I/O check normal");
781 Assert(EqualIO(cv,xcv,EPS),"CVector I/O check normal");
782 fin >> tmv::CompactIO() >> xv >> tmv::CompactIO() >> xcv;
783 Assert(EqualIO(v,xv,EPS),"Vector I/O check compact");
784 Assert(EqualIO(cv,xcv,EPS),"CVector I/O check compact");
785 fin >> xv.view() >> xcv.view();
786 Assert(EqualIO(v2,xv,EPS),"Vector I/O check thresh");
787 Assert(EqualIO(cv2,xcv,EPS),"CVector I/O check thresh");
788 fin >> myStyle >> xv.view() >> myStyle >> xcv.view();
789 Assert(EqualIO(v3,xv,EPS),"Vector I/O check compact thresh & prec(4)");
790 Assert(EqualIO(cv3,xcv,EPS),"CVector I/O check compact thresh & prec(4)");
791 fin.close();
792
793 // Read into vectors that need to be resized.
794 // Also check switching the default IOStyle.
795 tmv::CompactIO().makeDefault();
796 tmv::Vector<T> zv1, zv2, zv3, zv4;
797 tmv::Vector<CT > zcv1, zcv2, zcv3, zcv4;
798 fin.open("tmvtest_vector_io.dat");
799 Assert(bool(fin),"Couldn't open tmvtest_vector_io.dat for input");
800 fin >> tmv::NormalIO() >> zv1 >> tmv::NormalIO() >> zcv1;
801 Assert(EqualIO(v,zv1,EPS),"Vector I/O check normal with resize");
802 Assert(EqualIO(cv,zcv1,EPS),"CVector I/O check normal with resize");
803 fin >> zv2 >> zcv2;
804 Assert(EqualIO(v,zv2,EPS),"Vector I/O check compact with resize");
805 Assert(EqualIO(cv,zcv2,EPS),"CVector I/O check compact with resize");
806 fin >> tmv::NormalIO() >> zv3 >> tmv::NormalIO() >> zcv3;
807 Assert(EqualIO(v2,zv3,EPS),"Vector I/O check thresh with resize");
808 Assert(EqualIO(cv2,zcv3,EPS),"CVector I/O check thresh with resize");
809 fin >> myStyle >> zv4 >> myStyle >> zcv4;
810 Assert(EqualIO(v3,zv4,EPS),"Vector I/O check compact thresh with resize");
811 Assert(EqualIO(cv3,zcv4,EPS),"CVector I/O check compact thresh with resize");
812 fin.close();
813 // Switch it back.
814 tmv::IOStyle::revertDefault();
815
816 #if XTEST == 0
817 std::remove("tmvtest_vector_io.dat");
818 #endif
819
820 if (showstartdone) {
821 std::cout<<"Done Test Vector I/O"<<std::endl;
822 }
823 }
824
TestVector()825 template <class T> void TestVector()
826 {
827 #if 1
828 TestVectorReal<T>();
829 TestVectorComplex<T>();
830 TestVectorIO<T>();
831 #endif
832
833 #if 1
834 TestVectorArith<T>();
835 #endif
836
837 std::cout<<"Vector<"<<tmv::TMV_Text(T())<<"> passed all tests\n";
838 }
839
840 #ifdef TEST_DOUBLE
841 template void TestVector<double>();
842 #endif
843 #ifdef TEST_FLOAT
844 template void TestVector<float>();
845 #endif
846 #ifdef TEST_LONGDOUBLE
847 template void TestVector<long double>();
848 #endif
849 #ifdef TEST_INT
850 template void TestVector<int>();
851 #endif
852