1 #include "TMV_Test.h"
2 #include "TMV_Test_1.h"
3 #include "TMV.h"
4 #include <fstream>
5 #include <cstdio>
6
7 #define CT std::complex<T>
8
9 #ifdef TMV_MEM_DEBUG
10 // There seems to be something weird about this throw statement. Something gets
11 // allocated when the TriMatrixReadError is constructed that is allocated with
12 // new[], but then at the end of the try block, it is deallocated with plain delete.
13 // The memory manager flags this up as an error, but since it is in not in TMV
14 // code (the module is marked with ??(0000)::??, which means it is in some library
15 // that doesn't include the mmgr.h file), I haven't been able to debug it. I can't
16 // even figure out what object in memory it corresponds to. Maybe something in the
17 // standard library's exception class, or runtime_error? But that's just a guess, really.
18 #define NOTHROW
19 #endif
20
21 template <class T, tmv::DiagType D, tmv::StorageType S>
TestBasicUpperTriMatrix_1()22 static void TestBasicUpperTriMatrix_1()
23 {
24 const int N = 10;
25
26 if (showstartdone) {
27 std::cout<<"Start TestBasicUpperTriMatrix_1\n";
28 std::cout<<"T = "<<tmv::TMV_Text(T())<<std::endl;
29 std::cout<<"D = "<<tmv::TMV_Text(D)<<std::endl;
30 std::cout<<"S = "<<tmv::TMV_Text(S)<<std::endl;
31 std::cout<<"N = "<<N<<std::endl;
32 }
33
34 tmv::UpperTriMatrix<T,D|S> u(N);
35
36 Assert(u.colsize() == N && u.rowsize() == N,
37 "Creating UpperTriMatrix(N)");
38
39 for (int i=0,k=1; i<N; ++i) for (int j=0; j<N; ++j, ++k) {
40 if (i < j || (D==tmv::NonUnitDiag && i==j)) {
41 u(i,j) = T(k);
42 }
43 }
44
45 tmv::UpperTriMatrixView<T> uv = u.view();
46 tmv::ConstUpperTriMatrixView<T> ucv = u.view();
47
48 const tmv::UpperTriMatrix<T,D|S>& ux = u;
49
50 for(int i=0,k=1;i<N;++i) for(int j=0;j<N;++j,++k) {
51 if (i < j) {
52 Assert(u(i,j) == T(k),"Read/Write TriMatrix");
53 Assert(ux(i,j) == T(k),"Access const TriMatrix");
54 Assert(ucv(i,j) == T(k),"Access TriMatrix CV");
55 Assert(uv(i,j) == T(k),"Access TriMatrix V");
56 Assert(u.row(i,i+1,N)(j-i-1) == T(k),"TriMatrix.row1");
57 Assert(ux.row(i,i+1,N)(j-i-1) == T(k),"const TriMatrix.row1");
58 Assert(ucv.row(i,i+1,N)(j-i-1) == T(k),"TriMatrix.row1 CV");
59 Assert(uv.row(i,i+1,N)(j-i-1) == T(k),"TriMatrix.row1 V");
60 Assert(u.row(i,j,N)(0) == T(k),"TriMatrix.row2");
61 Assert(ux.row(i,j,N)(0) == T(k),"const TriMatrix.row2");
62 Assert(ucv.row(i,j,N)(0) == T(k),"TriMatrix.row2 CV");
63 Assert(uv.row(i,j,N)(0) == T(k),"TriMatrix.row2 V");
64 Assert(u.col(j,0,j)(i) == T(k),"TriMatrix.col1");
65 Assert(ux.col(j,0,j)(i) == T(k),"const TriMatrix.col1");
66 Assert(ucv.col(j,0,j)(i) == T(k),"TriMatrix.col1 CV");
67 Assert(uv.col(j,0,j)(i) == T(k),"TriMatrix.col1 V");
68 Assert(u.col(j,i,j)(0) == T(k),"TriMatrix.col2");
69 Assert(ux.col(j,i,j)(0) == T(k),"const TriMatrix.col2");
70 Assert(ucv.col(j,i,j)(0) == T(k),"TriMatrix.col2 CV");
71 Assert(uv.col(j,i,j)(0) == T(k),"TriMatrix.col2 V");
72 Assert(u.diag(j-i)(i) == T(k),"TriMatrix.diag1");
73 Assert(ux.diag(j-i)(i) == T(k),"const TriMatrix.diag1");
74 Assert(ucv.diag(j-i)(i) == T(k),"TriMatrix.diag1 CV ");
75 Assert(uv.diag(j-i)(i) == T(k),"TriMatrix.diag1 V");
76 Assert(u.diag(j-i,i,N-j+i)(0) == T(k),"TriMatrix.diag2");
77 Assert(ux.diag(j-i,i,N-j+i)(0) == T(k),"const TriMatrix.diag2");
78 Assert(ucv.diag(j-i,i,N-j+i)(0) == T(k),"TriMatrix.diag2 CV ");
79 Assert(uv.diag(j-i,i,N-j+i)(0) == T(k),"TriMatrix.diag2 V");
80 } else if (i==j) {
81 if (D == tmv::UnitDiag) {
82 Assert(ux(i,i) == T(1),"Access const TriMatrix");
83 Assert(ucv(i,i) == T(1),"Access TriMatrix CV");
84 } else {
85 Assert(u(i,i) == T(k),"Read/Write TriMatrix");
86 Assert(ux(i,i) == T(k),"Access const TriMatrix");
87 Assert(ucv(i,i) == T(k),"Access TriMatrix CV");
88 Assert(uv(i,i) == T(k),"Access TriMatrix V");
89 Assert(u.row(i,i,N)(0) == T(k),"TriMatrix.row1");
90 Assert(ux.row(i,i,N)(0) == T(k),"const TriMatrix.row1");
91 Assert(ucv.row(i,i,N)(0) == T(k),"TriMatrix.row1 CV");
92 Assert(uv.row(i,i,N)(0) == T(k),"TriMatrix.row1 V");
93 Assert(u.col(i,0,i+1)(i) == T(k),"TriMatrix.col1");
94 Assert(ux.col(i,0,i+1)(i) == T(k),"const TriMatrix.col1");
95 Assert(ucv.col(i,0,i+1)(i) == T(k),"TriMatrix.col1 CV");
96 Assert(uv.col(i,0,i+1)(i) == T(k),"TriMatrix.col1 V");
97 Assert(u.col(i,i,i+1)(0) == T(k),"TriMatrix.col2");
98 Assert(ux.col(i,i,i+1)(0) == T(k),"const TriMatrix.col2");
99 Assert(ucv.col(i,i,i+1)(0) == T(k),"TriMatrix.col2 CV");
100 Assert(uv.col(i,i,i+1)(0) == T(k),"TriMatrix.col2 V");
101 Assert(u.diag()(i) == T(k),"TriMatrix.diag");
102 Assert(ux.diag()(i) == T(k),"const TriMatrix.diag");
103 Assert(ucv.diag()(i) == T(k),"TriMatrix.diag CV ");
104 Assert(uv.diag()(i) == T(k),"TriMatrix.diag V");
105 Assert(u.diag(0)(i) == T(k),"TriMatrix.diag1");
106 Assert(ux.diag(0)(i) == T(k),"const TriMatrix.diag1");
107 Assert(ucv.diag(0)(i) == T(k),"TriMatrix.diag1 CV ");
108 Assert(uv.diag(0)(i) == T(k),"TriMatrix.diag1 V");
109 Assert(u.diag(0,i,N-i+i)(0) == T(k),"TriMatrix.diag2");
110 Assert(ux.diag(0,i,N-i+i)(0) == T(k),"const TriMatrix.diag2");
111 Assert(ucv.diag(0,i,N-i+i)(0) == T(k),"TriMatrix.diag2 CV ");
112 Assert(uv.diag(0,i,N-i+i)(0) == T(k),"TriMatrix.diag2 V");
113 }
114 } else {
115 Assert(ux(i,j) == T(0),"Access const TriMatrix");
116 Assert(ucv(i,j) == T(0),"Access TriMatrix CV");
117 }
118 }
119
120 #if (XTEST & 32)
121 tmv::UpperTriMatrix<T,D|S|tmv::FortranStyle> uf(N);
122 for (int i=0,k=1; i<N; ++i) for (int j=0; j<N; ++j, ++k) {
123 if (i < j || (D==tmv::NonUnitDiag && i==j)) {
124 uf(i+1,j+1) = T(k);
125 }
126 }
127
128 tmv::UpperTriMatrixView<T,tmv::FortranStyle> ufv = uf.view();
129 tmv::ConstUpperTriMatrixView<T,tmv::FortranStyle> ufcv = uf.view();
130
131 const tmv::UpperTriMatrix<T,D|S|tmv::FortranStyle>& ufx = uf;
132
133 for(int i=0,k=1;i<N;++i) for(int j=0;j<N;++j,++k) {
134 if (i < j) {
135 Assert(uf(i+1,j+1) == T(k),"Read/Write TriMatrixF");
136 Assert(ufx(i+1,j+1) == T(k),"Access const TriMatrixF");
137 Assert(ufcv(i+1,j+1) == T(k),"Access TriMatrixF CV");
138 Assert(ufv(i+1,j+1) == T(k),"Access TriMatrixF V");
139 Assert(uf.row(i+1,i+2,N)(j-i) == T(k),"TriMatrixF.row1");
140 Assert(ufx.row(i+1,i+2,N)(j-i) == T(k),"const TriMatrixF.row1");
141 Assert(ufcv.row(i+1,i+2,N)(j-i) == T(k),"TriMatrixF.row1 CV");
142 Assert(ufv.row(i+1,i+2,N)(j-i) == T(k),"TriMatrixF.row1 V");
143 Assert(uf.row(i+1,j+1,N)(1) == T(k),"TriMatrixF.row2");
144 Assert(ufx.row(i+1,j+1,N)(1) == T(k),"const TriMatrixF.row2");
145 Assert(ufcv.row(i+1,j+1,N)(1) == T(k),"TriMatrixF.row2 CV");
146 Assert(ufv.row(i+1,j+1,N)(1) == T(k),"TriMatrixF.row2 V");
147 Assert(uf.col(j+1,1,j)(i+1) == T(k),"TriMatrixF.col1");
148 Assert(ufx.col(j+1,1,j)(i+1) == T(k),"const TriMatrixF.col1");
149 Assert(ufcv.col(j+1,1,j)(i+1) == T(k),"TriMatrixF.col1 CV");
150 Assert(ufv.col(j+1,1,j)(i+1) == T(k),"TriMatrixF.col1 V");
151 Assert(uf.col(j+1,i+1,j)(1) == T(k),"TriMatrixF.col2");
152 Assert(ufx.col(j+1,i+1,j)(1) == T(k),"const TriMatrixF.col2");
153 Assert(ufcv.col(j+1,i+1,j)(1) == T(k),"TriMatrixF.col2 CV");
154 Assert(ufv.col(j+1,i+1,j)(1) == T(k),"TriMatrixF.col2 V");
155 Assert(uf.diag(j-i)(i+1) == T(k),"TriMatrixF.diag1");
156 Assert(ufx.diag(j-i)(i+1) == T(k),"const TriMatrixF.diag1");
157 Assert(ufcv.diag(j-i)(i+1) == T(k),"TriMatrixF.diag1 CV ");
158 Assert(ufv.diag(j-i)(i+1) == T(k),"TriMatrixF.diag1 V");
159 Assert(uf.diag(j-i,i+1,N-j+i)(1) == T(k),"TriMatrixF.diag2");
160 Assert(ufx.diag(j-i,i+1,N-j+i)(1) == T(k),"const TriMatrixF.diag2");
161 Assert(ufcv.diag(j-i,i+1,N-j+i)(1) == T(k),"TriMatrixF.diag2 CV ");
162 Assert(ufv.diag(j-i,i+1,N-j+i)(1) == T(k),"TriMatrixF.diag2 V");
163 } else if (i==j) {
164 if (D == tmv::UnitDiag) {
165 Assert(ufx(i+1,i+1) == T(1),"Access const TriMatrixF");
166 Assert(ufcv(i+1,i+1) == T(1),"Access TriMatrixF CV");
167 } else {
168 Assert(uf(i+1,i+1) == T(k),"Read/Write TriMatrixF");
169 Assert(ufx(i+1,i+1) == T(k),"Access const TriMatrixF");
170 Assert(ufcv(i+1,i+1) == T(k),"Access TriMatrixF CV");
171 Assert(ufv(i+1,i+1) == T(k),"Access TriMatrixF V");
172 Assert(uf.row(i+1,i+1,N)(1) == T(k),"TriMatrixF.row1");
173 Assert(ufx.row(i+1,i+1,N)(1) == T(k),"const TriMatrixF.row1");
174 Assert(ufcv.row(i+1,i+1,N)(1) == T(k),"TriMatrixF.row1 CV");
175 Assert(ufv.row(i+1,i+1,N)(1) == T(k),"TriMatrixF.row1 V");
176 Assert(uf.col(i+1,1,i+1)(i+1) == T(k),"TriMatrixF.col1");
177 Assert(ufx.col(i+1,1,i+1)(i+1) == T(k),"const TriMatrixF.col1");
178 Assert(ufcv.col(i+1,1,i+1)(i+1) == T(k),"TriMatrixF.col1 CV");
179 Assert(ufv.col(i+1,1,i+1)(i+1) == T(k),"TriMatrixF.col1 V");
180 Assert(uf.col(i+1,i+1,i+1)(1) == T(k),"TriMatrixF.col2");
181 Assert(ufx.col(i+1,i+1,i+1)(1) == T(k),"const TriMatrixF.col2");
182 Assert(ufcv.col(i+1,i+1,i+1)(1) == T(k),"TriMatrixF.col2 CV");
183 Assert(ufv.col(i+1,i+1,i+1)(1) == T(k),"TriMatrixF.col2 V");
184 Assert(uf.diag()(i+1) == T(k),"TriMatrixF.diag");
185 Assert(ufx.diag()(i+1) == T(k),"const TriMatrixF.diag");
186 Assert(ufcv.diag()(i+1) == T(k),"TriMatrixF.diag CV ");
187 Assert(ufv.diag()(i+1) == T(k),"TriMatrixF.diag V");
188 Assert(uf.diag(0)(i+1) == T(k),"TriMatrixF.diag1");
189 Assert(ufx.diag(0)(i+1) == T(k),"const TriMatrixF.diag1");
190 Assert(ufcv.diag(0)(i+1) == T(k),"TriMatrixF.diag1 CV ");
191 Assert(ufv.diag(0)(i+1) == T(k),"TriMatrixF.diag1 V");
192 Assert(uf.diag(0,i+1,N-i+i)(1) == T(k),"TriMatrixF.diag2");
193 Assert(ufx.diag(0,i+1,N-i+i)(1) == T(k),"const TriMatrixF.diag2");
194 Assert(ufcv.diag(0,i+1,N-i+i)(1) == T(k),"TriMatrixF.diag2 CV ");
195 Assert(ufv.diag(0,i+1,N-i+i)(1) == T(k),"TriMatrixF.diag2 V");
196 }
197 } else {
198 Assert(ufx(i+1,j+1) == T(0),"Access const TriMatrixF");
199 Assert(ufcv(i+1,j+1) == T(0),"Access TriMatrixF CV");
200 }
201 }
202 #endif
203
204 u.resize(2);
205 Assert(u.colsize() == 2 && u.rowsize() == 2,"u.resize(2)");
206 for (int i=0, k=0; i<2; ++i) for (int j=0; j<2; ++j, ++k)
207 if (i < j || (D==tmv::NonUnitDiag && i==j)) u(i,j) = T(k);
208 for (int i=0, k=0; i<2; ++i) for (int j=0; j<2; ++j, ++k)
209 if (i < j || (D==tmv::NonUnitDiag && i==j))
210 Assert(u(i,j) == k,"Read/Write resized UpperTriMatrix");
211
212 u.resize(2*N);
213 Assert(u.colsize() == 2*N && u.rowsize() == 2*N,"m.resize(2*N)");
214 for (int i=0, k=0; i<2*N; ++i) for (int j=0; j<2*N; ++j, ++k)
215 if (i < j || (D==tmv::NonUnitDiag && i==j)) u(i,j) = T(k);
216 for (int i=0, k=0; i<2*N; ++i) for (int j=0; j<2*N; ++j, ++k)
217 if (i < j || (D==tmv::NonUnitDiag && i==j))
218 Assert(u(i,j) == k,"Read/Write resized UpperTriMatrix");
219
220 }
221
222 template <class T, tmv::DiagType D, tmv::StorageType S>
TestBasicLowerTriMatrix_1()223 static void TestBasicLowerTriMatrix_1()
224 {
225 const int N = 10;
226
227 if (showstartdone) {
228 std::cout<<"Start TestBasicLowerTriMatrix_1\n";
229 std::cout<<"T = "<<tmv::TMV_Text(T())<<std::endl;
230 std::cout<<"D = "<<tmv::TMV_Text(D)<<std::endl;
231 std::cout<<"S = "<<tmv::TMV_Text(S)<<std::endl;
232 std::cout<<"N = "<<N<<std::endl;
233 }
234
235 tmv::LowerTriMatrix<T,D|S> l(N);
236
237 Assert(l.colsize() == N && l.rowsize() == N,
238 "Creating LowerTriMatrix(N)");
239
240 for (int i=0,k=1; i<N; ++i) for (int j=0; j<N; ++j, ++k) {
241 if (j < i || (D==tmv::NonUnitDiag && i==j)) {
242 l(i,j) = T(k);
243 }
244 }
245
246 tmv::LowerTriMatrixView<T> lv = l.view();
247 tmv::ConstLowerTriMatrixView<T> lcv = l.view();
248
249 const tmv::LowerTriMatrix<T,D|S>& lx = l;
250
251 for(int i=0,k=1;i<N;++i) for(int j=0;j<N;++j,++k) {
252 if (i < j) {
253 Assert(lx(i,j) == T(0),"Access const TriMatrix");
254 Assert(lcv(i,j) == T(0),"Access TriMatrix CV");
255 } else if (i==j) {
256 if (D == tmv::UnitDiag) {
257 Assert(lx(i,i) == T(1),"Access const TriMatrix");
258 Assert(lcv(i,i) == T(1),"Access TriMatrix CV");
259 } else {
260 Assert(l(i,i) == T(k),"Read/Write TriMatrix");
261 Assert(lx(i,i) == T(k),"Access const TriMatrix");
262 Assert(lcv(i,i) == T(k),"Access TriMatrix CV");
263 Assert(lv(i,i) == T(k),"Access TriMatrix V");
264 Assert(l.col(i,i,N)(0) == T(k),"TriMatrix.col2");
265 Assert(lx.col(i,i,N)(0) == T(k),"const TriMatrix.col2");
266 Assert(lcv.col(i,i,N)(0) == T(k),"TriMatrix.col2 CV");
267 Assert(lv.col(i,i,N)(0) == T(k),"TriMatrix.col2 V");
268 Assert(l.row(i,0,i+1)(i) == T(k),"TriMatrix.row1");
269 Assert(lx.row(i,0,i+1)(i) == T(k),"const TriMatrix.row1");
270 Assert(lcv.row(i,0,i+1)(i) == T(k),"TriMatrix.row1 CV");
271 Assert(lv.row(i,0,i+1)(i) == T(k),"TriMatrix.row1 V");
272 Assert(l.row(i,i,i+1)(0) == T(k),"TriMatrix.row2");
273 Assert(lx.row(i,i,i+1)(0) == T(k),"const TriMatrix.row2");
274 Assert(lcv.row(i,i,i+1)(0) == T(k),"TriMatrix.row2 CV");
275 Assert(lv.row(i,i,i+1)(0) == T(k),"TriMatrix.row2 V");
276 Assert(l.diag()(i) == T(k),"TriMatrix.diag");
277 Assert(lx.diag()(i) == T(k),"const TriMatrix.diag");
278 Assert(lcv.diag()(i) == T(k),"TriMatrix.diag CV ");
279 Assert(lv.diag()(i) == T(k),"TriMatrix.diag V");
280 Assert(l.diag(0)(i) == T(k),"TriMatrix.diag1");
281 Assert(lx.diag(0)(i) == T(k),"const TriMatrix.diag1");
282 Assert(lcv.diag(0)(i) == T(k),"TriMatrix.diag1 CV ");
283 Assert(lv.diag(0)(i) == T(k),"TriMatrix.diag1 V");
284 Assert(l.diag(0,i,N-i+i)(0) == T(k),"TriMatrix.diag2");
285 Assert(lx.diag(0,i,N-i+i)(0) == T(k),"const TriMatrix.diag2");
286 Assert(lcv.diag(0,i,N-i+i)(0) == T(k),"TriMatrix.diag2 CV ");
287 Assert(lv.diag(0,i,N-i+i)(0) == T(k),"TriMatrix.diag2 V");
288 }
289 } else {
290 Assert(l(i,j) == T(k),"Read/Write TriMatrix");
291 Assert(lx(i,j) == T(k),"Access const TriMatrix");
292 Assert(lcv(i,j) == T(k),"Access TriMatrix CV");
293 Assert(lv(i,j) == T(k),"Access TriMatrix V");
294 Assert(l.col(j,j+1,N)(i-j-1) == T(k),"TriMatrix.col1");
295 Assert(lx.col(j,j+1,N)(i-j-1) == T(k),"const TriMatrix.col1");
296 Assert(lcv.col(j,j+1,N)(i-j-1) == T(k),"TriMatrix.col1 CV");
297 Assert(lv.col(j,j+1,N)(i-j-1) == T(k),"TriMatrix.col1 V");
298 Assert(l.col(j,i,N)(0) == T(k),"TriMatrix.col2");
299 Assert(lx.col(j,i,N)(0) == T(k),"const TriMatrix.col2");
300 Assert(lcv.col(j,i,N)(0) == T(k),"TriMatrix.col2 CV");
301 Assert(lv.col(j,i,N)(0) == T(k),"TriMatrix.col2 V");
302 Assert(l.row(i,0,i)(j) == T(k),"TriMatrix.row1");
303 Assert(lx.row(i,0,i)(j) == T(k),"const TriMatrix.row1");
304 Assert(lcv.row(i,0,i)(j) == T(k),"TriMatrix.row1 CV");
305 Assert(lv.row(i,0,i)(j) == T(k),"TriMatrix.row1 V");
306 Assert(l.row(i,j,i)(0) == T(k),"TriMatrix.row2");
307 Assert(lx.row(i,j,i)(0) == T(k),"const TriMatrix.row2");
308 Assert(lcv.row(i,j,i)(0) == T(k),"TriMatrix.row2 CV");
309 Assert(lv.row(i,j,i)(0) == T(k),"TriMatrix.row2 V");
310 Assert(l.diag(-int(i-j))(j) == T(k),"TriMatrix.diag1");
311 Assert(lx.diag(-int(i-j))(j) == T(k),"const TriMatrix.diag1");
312 Assert(lcv.diag(-int(i-j))(j) == T(k),"TriMatrix.diag1 CV ");
313 Assert(lv.diag(-int(i-j))(j) == T(k),"TriMatrix.diag1 V");
314 Assert(l.diag(-int(i-j),j,N+j-i)(0) == T(k),"TriMatrix.diag2");
315 Assert(lx.diag(-int(i-j),j,N+j-i)(0) == T(k),"const TriMatrix.diag2");
316 Assert(lcv.diag(-int(i-j),j,N+j-i)(0) == T(k),"TriMatrix.diag2 CV ");
317 Assert(lv.diag(-int(i-j),j,N+j-i)(0) == T(k),"TriMatrix.diag2 V");
318 }
319 }
320
321 #if (XTEST & 32)
322 tmv::LowerTriMatrix<T,D|S|tmv::FortranStyle> lf(N);
323
324 for (int i=0,k=1; i<N; ++i) for (int j=0; j<N; ++j, ++k) {
325 if (j < i || (D==tmv::NonUnitDiag && i==j)) {
326 lf(i+1,j+1) = T(k);
327 }
328 }
329
330 tmv::LowerTriMatrixView<T,tmv::FortranStyle> lfv = lf.view();
331 tmv::ConstLowerTriMatrixView<T,tmv::FortranStyle> lfcv = lf.view();
332
333 const tmv::LowerTriMatrix<T,D|S|tmv::FortranStyle>& lfx = lf;
334
335 for(int i=0,k=1;i<N;++i) for(int j=0;j<N;++j,++k) {
336 if (i < j) {
337 Assert(lfx(i+1,j+1) == T(0),"Access const TriMatrixF");
338 Assert(lfcv(i+1,j+1) == T(0),"Access TriMatrixF CV");
339 } else if (i==j) {
340 if (D == tmv::UnitDiag) {
341 Assert(lfx(i+1,i+1) == T(1),"Access const TriMatrixF");
342 Assert(lfcv(i+1,i+1) == T(1),"Access TriMatrixF CV");
343 } else {
344 Assert(lf(i+1,i+1) == T(k),"Read/Write TriMatrixF");
345 Assert(lfx(i+1,i+1) == T(k),"Access const TriMatrixF");
346 Assert(lfcv(i+1,i+1) == T(k),"Access TriMatrixF CV");
347 Assert(lfv(i+1,i+1) == T(k),"Access TriMatrixF V");
348 Assert(lf.col(i+1,i+1,N)(1) == T(k),"TriMatrixF.col2");
349 Assert(lfx.col(i+1,i+1,N)(1) == T(k),"const TriMatrixF.col2");
350 Assert(lfcv.col(i+1,i+1,N)(1) == T(k),"TriMatrixF.col2 CV");
351 Assert(lfv.col(i+1,i+1,N)(1) == T(k),"TriMatrixF.col2 V");
352 Assert(lf.row(i+1,1,i+1)(i+1) == T(k),"TriMatrixF.row1");
353 Assert(lfx.row(i+1,1,i+1)(i+1) == T(k),"const TriMatrixF.row1");
354 Assert(lfcv.row(i+1,1,i+1)(i+1) == T(k),"TriMatrixF.row1 CV");
355 Assert(lfv.row(i+1,1,i+1)(i+1) == T(k),"TriMatrixF.row1 V");
356 Assert(lf.row(i+1,i+1,i+1)(1) == T(k),"TriMatrixF.row2");
357 Assert(lfx.row(i+1,i+1,i+1)(1) == T(k),"const TriMatrixF.row2");
358 Assert(lfcv.row(i+1,i+1,i+1)(1) == T(k),"TriMatrixF.row2 CV");
359 Assert(lfv.row(i+1,i+1,i+1)(1) == T(k),"TriMatrixF.row2 V");
360 Assert(lf.diag()(i+1) == T(k),"TriMatrixF.diag");
361 Assert(lfx.diag()(i+1) == T(k),"const TriMatrixF.diag");
362 Assert(lfcv.diag()(i+1) == T(k),"TriMatrixF.diag CV ");
363 Assert(lfv.diag()(i+1) == T(k),"TriMatrixF.diag V");
364 Assert(lf.diag(0)(i+1) == T(k),"TriMatrixF.diag1");
365 Assert(lfx.diag(0)(i+1) == T(k),"const TriMatrixF.diag1");
366 Assert(lfcv.diag(0)(i+1) == T(k),"TriMatrixF.diag1 CV ");
367 Assert(lfv.diag(0)(i+1) == T(k),"TriMatrixF.diag1 V");
368 Assert(lf.diag(0,i+1,N-i+i)(1) == T(k),"TriMatrixF.diag2");
369 Assert(lfx.diag(0,i+1,N-i+i)(1) == T(k),"const TriMatrixF.diag2");
370 Assert(lfcv.diag(0,i+1,N-i+i)(1) == T(k),"TriMatrixF.diag2 CV ");
371 Assert(lfv.diag(0,i+1,N-i+i)(1) == T(k),"TriMatrixF.diag2 V");
372 }
373 } else {
374 Assert(lf(i+1,j+1) == T(k),"Read/Write TriMatrixF");
375 Assert(lfx(i+1,j+1) == T(k),"Access const TriMatrixF");
376 Assert(lfcv(i+1,j+1) == T(k),"Access TriMatrixF CV");
377 Assert(lfv(i+1,j+1) == T(k),"Access TriMatrixF V");
378 Assert(lf.col(j+1,j+2,N)(i-j) == T(k),"TriMatrixF.col1");
379 Assert(lfx.col(j+1,j+2,N)(i-j) == T(k),"const TriMatrixF.col1");
380 Assert(lfcv.col(j+1,j+2,N)(i-j) == T(k),"TriMatrixF.col1 CV");
381 Assert(lfv.col(j+1,j+2,N)(i-j) == T(k),"TriMatrixF.col1 V");
382 Assert(lf.col(j+1,i+1,N)(1) == T(k),"TriMatrixF.col2");
383 Assert(lfx.col(j+1,i+1,N)(1) == T(k),"const TriMatrixF.col2");
384 Assert(lfcv.col(j+1,i+1,N)(1) == T(k),"TriMatrixF.col2 CV");
385 Assert(lfv.col(j+1,i+1,N)(1) == T(k),"TriMatrixF.col2 V");
386 Assert(lf.row(i+1,1,i)(j+1) == T(k),"TriMatrixF.row1");
387 Assert(lfx.row(i+1,1,i)(j+1) == T(k),"const TriMatrixF.row1");
388 Assert(lfcv.row(i+1,1,i)(j+1) == T(k),"TriMatrixF.row1 CV");
389 Assert(lfv.row(i+1,1,i)(j+1) == T(k),"TriMatrixF.row1 V");
390 Assert(lf.row(i+1,j+1,i)(1) == T(k),"TriMatrixF.row2");
391 Assert(lfx.row(i+1,j+1,i)(1) == T(k),"const TriMatrixF.row2");
392 Assert(lfcv.row(i+1,j+1,i)(1) == T(k),"TriMatrixF.row2 CV");
393 Assert(lfv.row(i+1,j+1,i)(1) == T(k),"TriMatrixF.row2 V");
394 Assert(lf.diag(-int(i-j))(j+1) == T(k),"TriMatrixF.diag1");
395 Assert(lfx.diag(-int(i-j))(j+1) == T(k),"const TriMatrixF.diag1");
396 Assert(lfcv.diag(-int(i-j))(j+1) == T(k),"TriMatrixF.diag1 CV ");
397 Assert(lfv.diag(-int(i-j))(j+1) == T(k),"TriMatrixF.diag1 V");
398 Assert(lf.diag(-int(i-j),j+1,N+j-i)(1) == T(k),"TriMatrixF.diag2");
399 Assert(lfx.diag(-int(i-j),j+1,N+j-i)(1) == T(k),"const TriMatrixF.diag2");
400 Assert(lfcv.diag(-int(i-j),j+1,N+j-i)(1) == T(k),"TriMatrixF.diag2 CV ");
401 Assert(lfv.diag(-int(i-j),j+1,N+j-i)(1) == T(k),"TriMatrixF.diag2 V");
402 }
403 }
404 #endif
405
406 l.resize(2);
407 Assert(l.colsize() == 2 && l.rowsize() == 2,"l.resize(2)");
408 for (int i=0, k=0; i<2; ++i) for (int j=0; j<2; ++j, ++k)
409 if (j < i || (D==tmv::NonUnitDiag && i==j)) l(i,j) = T(k);
410 for (int i=0, k=0; i<2; ++i) for (int j=0; j<2; ++j, ++k)
411 if (j < i || (D==tmv::NonUnitDiag && i==j))
412 Assert(l(i,j) == k,"Read/Write resized UpperTriMatrix");
413
414 l.resize(2*N);
415 Assert(l.colsize() == 2*N && l.rowsize() == 2*N,"m.resize(2*N)");
416 for (int i=0, k=0; i<2*N; ++i) for (int j=0; j<2*N; ++j, ++k)
417 if (j < i || (D==tmv::NonUnitDiag && i==j)) l(i,j) = T(k);
418 for (int i=0, k=0; i<2*N; ++i) for (int j=0; j<2*N; ++j, ++k)
419 if (j < i || (D==tmv::NonUnitDiag && i==j))
420 Assert(l(i,j) == k,"Read/Write resized UpperTriMatrix");
421
422 }
423
424 template <class T, tmv::DiagType D, tmv::StorageType S>
TestBasicTriMatrix_2()425 static void TestBasicTriMatrix_2()
426 {
427 const int N = 10;
428
429 if (showstartdone) {
430 std::cout<<"Start TestBasicTriMatrix_2\n";
431 std::cout<<"T = "<<tmv::TMV_Text(T())<<std::endl;
432 std::cout<<"D = "<<tmv::TMV_Text(D)<<std::endl;
433 std::cout<<"S = "<<tmv::TMV_Text(S)<<std::endl;
434 std::cout<<"N = "<<N<<std::endl;
435 }
436
437 tmv::UpperTriMatrix<T,D|S> u(N);
438 tmv::LowerTriMatrix<T,D|S> l(N);
439 const tmv::UpperTriMatrix<T,D|S>& ux = u;
440 const tmv::LowerTriMatrix<T,D|S>& lx = l;
441
442 Assert(u.colsize() == N && u.rowsize() == N,
443 "Creating UpperTriMatrix(N)");
444 Assert(l.colsize() == N && l.rowsize() == N,
445 "Creating LowerTriMatrix(N)");
446
447 for (int i=0,k=1; i<N; ++i) for (int j=0; j<N; ++j, ++k) {
448 if (i < j || (D==tmv::NonUnitDiag && i==j)) {
449 u(i,j) = T(k);
450 }
451 if (j < i || (D==tmv::NonUnitDiag && i==j)) {
452 l(i,j) = T(k);
453 }
454 }
455
456 tmv::Matrix<T> mu = u;
457 tmv::Matrix<T> ml = l;
458
459 for(int i=0,k=0;i<N;++i) for(int j=0;j<N;++j,++k) {
460 Assert(mu(i,j) == ux(i,j),"TriMatrix -> Matrix");
461 Assert(ml(i,j) == lx(i,j),"TriMatrix -> Matrix");
462 }
463 Assert(u == mu.upperTri(),"TriMatrix == ");
464 tmv::UpperTriMatrix<T,D|S> umu(mu);
465 Assert(u == umu,"TriMatrix == ");
466 Assert(l == ml.lowerTri(),"TriMatrix == ");
467 tmv::LowerTriMatrix<T,D|S> lml(ml);
468 Assert(l == lml,"TriMatrix == ");
469
470 // Test assignments and constructors from arrays
471 const T quarrmnonunit[] = {
472 T(0), T(-1), T(-2),
473 T(1), T(0),
474 T(2)
475 };
476 const T quarrmunit[] = {
477 T(1), T(-1), T(-2),
478 T(1), T(0),
479 T(1)
480 };
481 const T* quarrm = (D == tmv::UnitDiag) ? quarrmunit : quarrmnonunit;
482
483 const T qlarrmnonunit[] = {
484 T(0),
485 T(2), T(1),
486 T(4), T(3), T(2)
487 };
488 const T qlarrmunit[] = {
489 T(1),
490 T(2), T(1),
491 T(4), T(3), T(1)
492 };
493 const T* qlarrm = (D == tmv::UnitDiag) ? qlarrmunit : qlarrmnonunit;
494
495 const T quarcmnonunit[] = {
496 T(0),
497 T(-1), T(1),
498 T(-2), T(0), T(2)
499 };
500 const T quarcmunit[] = {
501 T(1),
502 T(-1), T(1),
503 T(-2), T(0), T(1)
504 };
505 const T* quarcm = (D == tmv::UnitDiag) ? quarcmunit : quarcmnonunit;
506
507 const T qlarcmnonunit[] = {
508 T(0), T(2), T(4),
509 T(1), T(3),
510 T(2)
511 };
512 const T qlarcmunit[] = {
513 T(1), T(2), T(4),
514 T(1), T(3),
515 T(1)
516 };
517 const T* qlarcm = (D == tmv::UnitDiag) ? qlarcmunit : qlarcmnonunit;
518
519 std::vector<T> quvecrm(6);
520 for(int i=0;i<6;i++) quvecrm[i] = quarrm[i];
521 std::vector<T> qlvecrm(6);
522 for(int i=0;i<6;i++) qlvecrm[i] = qlarrm[i];
523 std::vector<T> quveccm(6);
524 for(int i=0;i<6;i++) quveccm[i] = quarcm[i];
525 std::vector<T> qlveccm(6);
526 for(int i=0;i<6;i++) qlveccm[i] = qlarcm[i];
527
528 tmv::UpperTriMatrix<T,D|S> qu1(3);
529 std::copy(quarrm, quarrm+6, qu1.rowmajor_begin());
530 tmv::UpperTriMatrix<T,D|S> qu2(3);
531 std::copy(quarcm, quarcm+6, qu2.colmajor_begin());
532
533 tmv::LowerTriMatrix<T,D|S> ql1(3);
534 std::copy(qlarrm, qlarrm+6, ql1.rowmajor_begin());
535 tmv::LowerTriMatrix<T,D|S> ql2(3);
536 std::copy(qlarcm, qlarcm+6, ql2.colmajor_begin());
537
538 tmv::UpperTriMatrix<T,D|S> qu3(3);
539 std::copy(quvecrm.begin(), quvecrm.end(), qu3.rowmajor_begin());
540 tmv::UpperTriMatrix<T,D|S> qu4(3);
541 std::copy(quveccm.begin(), quveccm.end(), qu4.colmajor_begin());
542
543 tmv::LowerTriMatrix<T,D|S> ql3(3);
544 std::copy(qlvecrm.begin(), qlvecrm.end(), ql3.rowmajor_begin());
545 tmv::LowerTriMatrix<T,D|S> ql4(3);
546 std::copy(qlveccm.begin(), qlveccm.end(), ql4.colmajor_begin());
547
548 tmv::UpperTriMatrix<T,D|S> qu5x(30);
549 tmv::UpperTriMatrixView<T> qu5 = qu5x.subTriMatrix(3,18,5);
550 std::copy(quvecrm.begin(), quvecrm.end(), qu5.rowmajor_begin());
551 tmv::UpperTriMatrix<T,D|S> qu6x(30);
552 tmv::UpperTriMatrixView<T> qu6 = qu6x.subTriMatrix(3,18,5);
553 std::copy(quveccm.begin(), quveccm.end(), qu6.colmajor_begin());
554
555 tmv::LowerTriMatrix<T,D|S> ql5x(30);
556 tmv::LowerTriMatrixView<T> ql5 = ql5x.subTriMatrix(3,18,5);
557 std::copy(qlvecrm.begin(), qlvecrm.end(), ql5.rowmajor_begin());
558 tmv::LowerTriMatrix<T,D|S> ql6x(30);
559 tmv::LowerTriMatrixView<T> ql6 = ql6x.subTriMatrix(3,18,5);
560 std::copy(qlveccm.begin(), qlveccm.end(), ql6.colmajor_begin());
561
562 // Assignment using op<< is always in rowmajor order.
563 tmv::UpperTriMatrix<T,D|S> qu7(3);
564 tmv::LowerTriMatrix<T,D|S> qu8t(3);
565 tmv::UpperTriMatrixView<T> qu8 = qu8t.transpose();
566
567 tmv::LowerTriMatrix<T,D|S> ql7(3);
568 tmv::UpperTriMatrix<T,D|S> ql8t(3);
569 tmv::LowerTriMatrixView<T> ql8 = ql8t.transpose();
570
571 if (D == tmv::UnitDiag) {
572 qu7 <<
573 1, -1, -2,
574 1, 0,
575 1;
576 qu8 <<
577 1, -1, -2,
578 1, 0,
579 1;
580
581 ql7 <<
582 1,
583 2, 1,
584 4, 3, 1;
585 ql8 <<
586 1,
587 2, 1,
588 4, 3, 1;
589 } else {
590 qu7 <<
591 0, -1, -2,
592 1, 0,
593 2;
594 qu8 <<
595 0, -1, -2,
596 1, 0,
597 2;
598
599 ql7 <<
600 0,
601 2, 1,
602 4, 3, 2;
603 ql8 <<
604 0,
605 2, 1,
606 4, 3, 2;
607 }
608
609 // Can also view memory directly
610 // (Diagonal elements in memory do not have to be 1 if D == UnitDiag.)
611 T qarrmfull[] = {
612 T(0), T(-1), T(-2),
613 T(2), T(1), T(0),
614 T(4), T(3), T(2),
615 };
616 T qarcmfull[] = {
617 T(0), T(2), T(4),
618 T(-1), T(1), T(3),
619 T(-2), T(0), T(2),
620 };
621 T* qarfull = (S == tmv::RowMajor) ? qarrmfull : qarcmfull;
622 const int Si = (S == tmv::RowMajor) ? 3 : 1;
623 const int Sj = (S == tmv::RowMajor) ? 1 : 3;
624 const tmv::ConstUpperTriMatrixView<T> qu9 =
625 tmv::UpperTriMatrixViewOf(qarfull,3,S,D);
626 const tmv::ConstUpperTriMatrixView<T> qu10 =
627 tmv::UpperTriMatrixViewOf(qarfull,3,Si,Sj,D);
628 const tmv::ConstLowerTriMatrixView<T> ql9 =
629 tmv::LowerTriMatrixViewOf(qarfull,3,S,D);
630 const tmv::ConstLowerTriMatrixView<T> ql10 =
631 tmv::LowerTriMatrixViewOf(qarfull,3,Si,Sj,D);
632
633 if (showacc) {
634 std::cout<<"qu1 = "<<qu1<<std::endl;
635 std::cout<<"qu2 = "<<qu2<<std::endl;
636 std::cout<<"qu3 = "<<qu3<<std::endl;
637 std::cout<<"qu4 = "<<qu4<<std::endl;
638 std::cout<<"qu5 = "<<qu5<<std::endl;
639 std::cout<<"qu6 = "<<qu6<<std::endl;
640 std::cout<<"qu7 = "<<qu7<<std::endl;
641 std::cout<<"qu8 = "<<qu8<<std::endl;
642 std::cout<<"qu9 = "<<qu9<<std::endl;
643 std::cout<<"qu10 = "<<qu10<<std::endl;
644
645 std::cout<<"ql1 = "<<ql1<<std::endl;
646 std::cout<<"ql2 = "<<ql2<<std::endl;
647 std::cout<<"ql3 = "<<ql3<<std::endl;
648 std::cout<<"ql4 = "<<ql4<<std::endl;
649 std::cout<<"ql5 = "<<ql5<<std::endl;
650 std::cout<<"ql6 = "<<ql6<<std::endl;
651 std::cout<<"ql7 = "<<ql7<<std::endl;
652 std::cout<<"ql8 = "<<ql8<<std::endl;
653 std::cout<<"ql9 = "<<ql9<<std::endl;
654 std::cout<<"ql10 = "<<ql10<<std::endl;
655 }
656
657 for(int i=0;i<3;i++) for(int j=0;j<3;j++) {
658 T val = (D == tmv::UnitDiag && i==j) ? T(1) : T(2*i-j);
659 if (j>=i) {
660 Assert(qu1(i,j) == val,"Create UpperTriMatrix from T* rm");
661 Assert(qu2(i,j) == val,"Create UpperTriMatrix from T* cm");
662 Assert(qu3(i,j) == val,"Create UpperTriMatrix from vector rm");
663 Assert(qu4(i,j) == val,"Create UpperTriMatrix from vector cm");
664 Assert(qu5(i,j) == val,"Create UpperTriMatrixView from vector rm");
665 Assert(qu6(i,j) == val,"Create UpperTriMatrixView from vector cm");
666 Assert(qu7(i,j) == val,"Create UpperTriMatrix from << list");
667 Assert(qu8(i,j) == val,"Create UpperTriMatrixView from << list");
668 Assert(qu9(i,j) == val,"Create UpperTriMatrixView of T* (S)");
669 Assert(qu10(i,j) == val,"Create UpperTriMatrixView of T* (Si,Sj)");
670 }
671 if (i>=j) {
672 Assert(ql1(i,j) == val,"Create UpperTriMatrix from T* rm");
673 Assert(ql2(i,j) == val,"Create UpperTriMatrix from T* cm");
674 Assert(ql3(i,j) == val,"Create UpperTriMatrix from vector rm");
675 Assert(ql4(i,j) == val,"Create UpperTriMatrix from vector cm");
676 Assert(ql5(i,j) == val,"Create UpperTriMatrixView from vector rm");
677 Assert(ql6(i,j) == val,"Create UpperTriMatrixView from vector cm");
678 Assert(ql7(i,j) == val,"Create UpperTriMatrix from << list");
679 Assert(ql8(i,j) == val,"Create UpperTriMatrixView from << list");
680 Assert(ql9(i,j) == val,"Create UpperTriMatrixView of T* (S)");
681 Assert(ql10(i,j) == val,"Create UpperTriMatrixView of T* (Si,Sj)");
682 }
683 }
684
685 // Test the span of the iteration (i.e. the validity of begin(), end())
686 const tmv::UpperTriMatrix<T,D|S>& qu1_const = qu1;
687 tmv::UpperTriMatrixView<T> qu1_view = qu1.view();
688 tmv::ConstUpperTriMatrixView<T> qu1_constview = qu1_const.view();
689 tmv::ConstUpperTriMatrixView<T> qu5_const = qu5;
690
691 typename tmv::UpperTriMatrix<T,D|S>::rowmajor_iterator rmitu1 =
692 qu1.rowmajor_begin();
693 typename tmv::UpperTriMatrix<T,D|S>::const_rowmajor_iterator rmitu2 =
694 qu1_const.rowmajor_begin();
695 typename tmv::UpperTriMatrixView<T>::rowmajor_iterator rmitu3 =
696 qu1_view.rowmajor_begin();
697 typename tmv::ConstUpperTriMatrixView<T>::const_rowmajor_iterator rmitu4 =
698 qu1_constview.rowmajor_begin();
699 typename tmv::UpperTriMatrixView<T>::rowmajor_iterator rmitu5 =
700 qu5.rowmajor_begin();
701 typename tmv::ConstUpperTriMatrixView<T>::const_rowmajor_iterator rmitu6 =
702 qu5_const.rowmajor_begin();
703 int ii = 0;
704 while (rmitu1 != qu1.rowmajor_end()) {
705 Assert(*rmitu1++ == quarrm[ii], "RowMajor iteration 1");
706 Assert(*rmitu2++ == quarrm[ii], "RowMajor iteration 2");
707 Assert(*rmitu3++ == quarrm[ii], "RowMajor iteration 3");
708 Assert(*rmitu4++ == quarrm[ii], "RowMajor iteration 4");
709 Assert(*rmitu5++ == quarrm[ii], "RowMajor iteration 5");
710 Assert(*rmitu6++ == quarrm[ii], "RowMajor iteration 6");
711 ++ii;
712 }
713 Assert(ii == 6, "RowMajor iteration number of elements");
714 Assert(rmitu2 == qu1_const.rowmajor_end(), "rmitu2 reaching end");
715 Assert(rmitu3 == qu1_view.rowmajor_end(), "rmitu3 reaching end");
716 Assert(rmitu4 == qu1_constview.rowmajor_end(), "rmitu4 reaching end");
717 Assert(rmitu5 == qu5.rowmajor_end(), "rmitu5 reaching end");
718 Assert(rmitu6 == qu5_const.rowmajor_end(), "rmitu6 reaching end");
719
720 const tmv::LowerTriMatrix<T,D|S>& ql1_const = ql1;
721 tmv::LowerTriMatrixView<T> ql1_view = ql1.view();
722 tmv::ConstLowerTriMatrixView<T> ql1_constview = ql1_const.view();
723 tmv::ConstLowerTriMatrixView<T> ql5_const = ql5;
724
725 typename tmv::LowerTriMatrix<T,D|S>::rowmajor_iterator rmitl1 =
726 ql1.rowmajor_begin();
727 typename tmv::LowerTriMatrix<T,D|S>::const_rowmajor_iterator rmitl2 =
728 ql1_const.rowmajor_begin();
729 typename tmv::LowerTriMatrixView<T>::rowmajor_iterator rmitl3 =
730 ql1_view.rowmajor_begin();
731 typename tmv::ConstLowerTriMatrixView<T>::const_rowmajor_iterator rmitl4 =
732 ql1_constview.rowmajor_begin();
733 typename tmv::LowerTriMatrixView<T>::rowmajor_iterator rmitl5 =
734 ql5.rowmajor_begin();
735 typename tmv::ConstLowerTriMatrixView<T>::const_rowmajor_iterator rmitl6 =
736 ql5_const.rowmajor_begin();
737 ii = 0;
738 while (rmitl1 != ql1.rowmajor_end()) {
739 Assert(*rmitl1++ == qlarrm[ii], "RowMajor iteration 1");
740 Assert(*rmitl2++ == qlarrm[ii], "RowMajor iteration 2");
741 Assert(*rmitl3++ == qlarrm[ii], "RowMajor iteration 3");
742 Assert(*rmitl4++ == qlarrm[ii], "RowMajor iteration 4");
743 Assert(*rmitl5++ == qlarrm[ii], "RowMajor iteration 5");
744 Assert(*rmitl6++ == qlarrm[ii], "RowMajor iteration 6");
745 ++ii;
746 }
747 Assert(ii == 6, "RowMajor iteration number of elements");
748 Assert(rmitl2 == ql1_const.rowmajor_end(), "rmitl2 reaching end");
749 Assert(rmitl3 == ql1_view.rowmajor_end(), "rmitl3 reaching end");
750 Assert(rmitl4 == ql1_constview.rowmajor_end(), "rmitl4 reaching end");
751 Assert(rmitl5 == ql5.rowmajor_end(), "rmitl5 reaching end");
752 Assert(rmitl6 == ql5_const.rowmajor_end(), "rmitl6 reaching end");
753
754 typename tmv::UpperTriMatrix<T,D|S>::colmajor_iterator cmitu1 =
755 qu1.colmajor_begin();
756 typename tmv::UpperTriMatrix<T,D|S>::const_colmajor_iterator cmitu2 =
757 qu1_const.colmajor_begin();
758 typename tmv::UpperTriMatrixView<T>::colmajor_iterator cmitu3 =
759 qu1_view.colmajor_begin();
760 typename tmv::ConstUpperTriMatrixView<T>::const_colmajor_iterator cmitu4 =
761 qu1_constview.colmajor_begin();
762 typename tmv::UpperTriMatrixView<T>::colmajor_iterator cmitu5 =
763 qu5.colmajor_begin();
764 typename tmv::ConstUpperTriMatrixView<T>::const_colmajor_iterator cmitu6 =
765 qu5_const.colmajor_begin();
766 ii = 0;
767 while (cmitu1 != qu1.colmajor_end()) {
768 Assert(*cmitu1++ == quarcm[ii], "ColMajor iteration 1");
769 Assert(*cmitu2++ == quarcm[ii], "ColMajor iteration 2");
770 Assert(*cmitu3++ == quarcm[ii], "ColMajor iteration 3");
771 Assert(*cmitu4++ == quarcm[ii], "ColMajor iteration 4");
772 Assert(*cmitu5++ == quarcm[ii], "ColMajor iteration 5");
773 Assert(*cmitu6++ == quarcm[ii], "ColMajor iteration 6");
774 ++ii;
775 }
776 Assert(ii == 6, "ColMajor iteration number of elements");
777 Assert(cmitu2 == qu1_const.colmajor_end(), "cmitu2 reaching end");
778 Assert(cmitu3 == qu1_view.colmajor_end(), "cmitu3 reaching end");
779 Assert(cmitu4 == qu1_constview.colmajor_end(), "cmitu4 reaching end");
780 Assert(cmitu5 == qu5.colmajor_end(), "cmitu5 reaching end");
781 Assert(cmitu6 == qu5_const.colmajor_end(), "cmitu6 reaching end");
782
783 typename tmv::LowerTriMatrix<T,D|S>::colmajor_iterator cmitl1 =
784 ql1.colmajor_begin();
785 typename tmv::LowerTriMatrix<T,D|S>::const_colmajor_iterator cmitl2 =
786 ql1_const.colmajor_begin();
787 typename tmv::LowerTriMatrixView<T>::colmajor_iterator cmitl3 =
788 ql1_view.colmajor_begin();
789 typename tmv::ConstLowerTriMatrixView<T>::const_colmajor_iterator cmitl4 =
790 ql1_constview.colmajor_begin();
791 typename tmv::LowerTriMatrixView<T>::colmajor_iterator cmitl5 =
792 ql5.colmajor_begin();
793 typename tmv::ConstLowerTriMatrixView<T>::const_colmajor_iterator cmitl6 =
794 ql5_const.colmajor_begin();
795 ii = 0;
796 while (cmitl1 != ql1.colmajor_end()) {
797 Assert(*cmitl1++ == qlarcm[ii], "ColMajor iteration 1");
798 Assert(*cmitl2++ == qlarcm[ii], "ColMajor iteration 2");
799 Assert(*cmitl3++ == qlarcm[ii], "ColMajor iteration 3");
800 Assert(*cmitl4++ == qlarcm[ii], "ColMajor iteration 4");
801 Assert(*cmitl5++ == qlarcm[ii], "ColMajor iteration 5");
802 Assert(*cmitl6++ == qlarcm[ii], "ColMajor iteration 6");
803 ++ii;
804 }
805 Assert(ii == 6, "ColMajor iteration number of elements");
806 Assert(cmitl2 == ql1_const.colmajor_end(), "cmitl2 reaching end");
807 Assert(cmitl3 == ql1_view.colmajor_end(), "cmitl3 reaching end");
808 Assert(cmitl4 == ql1_constview.colmajor_end(), "cmitl4 reaching end");
809 Assert(cmitl5 == ql5.colmajor_end(), "cmitl5 reaching end");
810 Assert(cmitl6 == ql5_const.colmajor_end(), "cmitl6 reaching end");
811
812
813 // Test Basic Arithmetic
814 tmv::Matrix<T> a(N,N);
815 for(int i=0;i<N;++i) for(int j=0;j<N;++j) a(i,j) = T(12+3*i-5*j);
816
817 tmv::UpperTriMatrix<T,D|S> u1(a);
818 tmv::UpperTriMatrix<T,D|S> u2(a.transpose());
819 tmv::LowerTriMatrix<T,D|S> l1(a);
820 tmv::LowerTriMatrix<T,D|S> l2(a.transpose());
821 const tmv::UpperTriMatrix<T,D|S>& u1x = u1;
822 const tmv::UpperTriMatrix<T,D|S>& u2x = u2;
823 const tmv::LowerTriMatrix<T,D|S>& l1x = l1;
824 const tmv::LowerTriMatrix<T,D|S>& l2x = l2;
825
826 tmv::UpperTriMatrix<T> u4 = u1+u1;
827 tmv::UpperTriMatrix<T> u5 = u1+u2;
828 tmv::UpperTriMatrix<T> u6 = u2+u2;
829 tmv::LowerTriMatrix<T> l4 = l1+l1;
830 tmv::LowerTriMatrix<T> l5 = l1+l2;
831 tmv::LowerTriMatrix<T> l6 = l2+l2;
832 tmv::Matrix<T> m1 = l1+u1;
833 tmv::Matrix<T> m2 = l1+u2;
834 tmv::Matrix<T> m3 = l2+u1;
835 tmv::Matrix<T> m4 = l2+u2;
836
837 const tmv::UpperTriMatrix<T>& u4x = u4;
838 const tmv::UpperTriMatrix<T>& u5x = u5;
839 const tmv::UpperTriMatrix<T>& u6x = u6;
840 const tmv::LowerTriMatrix<T>& l4x = l4;
841 const tmv::LowerTriMatrix<T>& l5x = l5;
842 const tmv::LowerTriMatrix<T>& l6x = l6;
843
844 for(int i=0;i<N;i++) for(int j=0;j<N;j++) {
845 Assert(u4x(i,j) == 2*u1x(i,j),"Add triMatrices");
846 Assert(u5x(i,j) == u1x(i,j) + u2x(i,j),"Add triMatrices");
847 Assert(u6x(i,j) == 2*u2x(i,j),"Add triMatrices");
848 Assert(l4x(i,j) == 2*l1x(i,j),"Add triMatrices");
849 Assert(l5x(i,j) == l1x(i,j) + l2x(i,j),"Add triMatrices");
850 Assert(l6x(i,j) == 2*l2x(i,j),"Add triMatrices");
851 Assert(m1(i,j) == l1x(i,j) + u1x(i,j),"Add triMatrices");
852 Assert(m2(i,j) == l1x(i,j) + u2x(i,j),"Add triMatrices");
853 Assert(m3(i,j) == l2x(i,j) + u1x(i,j),"Add triMatrices");
854 Assert(m4(i,j) == l2x(i,j) + u2x(i,j),"Add triMatrices");
855 }
856
857 u4 = u1-u1;
858 u5 = u1-u2;
859 u6 = u2-u2;
860 l4 = l1-l1;
861 l5 = l1-l2;
862 l6 = l2-l2;
863 m1 = l1-u1;
864 m2 = l1-u2;
865 m3 = l2-u1;
866 m4 = l2-u2;
867
868 for(int i=0;i<N;i++) for(int j=0;j<N;j++) {
869 Assert(u4x(i,j) == T(0),"Subtract TriMatrices");
870 Assert(u5x(i,j) == u1x(i,j) - u2x(i,j),"Subtract TriMatrices");
871 Assert(u6x(i,j) == T(0),"Subtract TriMatrices");
872 Assert(l4x(i,j) == T(0),"Subtract TriMatrices");
873 Assert(l5x(i,j) == l1x(i,j) - l2x(i,j),"Subtract TriMatrices");
874 Assert(l6x(i,j) == T(0),"Subtract TriMatrices");
875 Assert(m1(i,j) == l1x(i,j) - u1x(i,j),"Subtract TriMatrices");
876 Assert(m2(i,j) == l1x(i,j) - u2x(i,j),"Subtract TriMatrices");
877 Assert(m3(i,j) == l2x(i,j) - u1x(i,j),"Subtract TriMatrices");
878 Assert(m4(i,j) == l2x(i,j) - u2x(i,j),"Subtract TriMatrices");
879 }
880
881 tmv::UpperTriMatrix<std::complex<T>,tmv::NonUnitDiag|tmv::RowMajor> cunr = u;
882 if (D == tmv::UnitDiag)
883 cunr.offDiag() *= std::complex<T>(1,2);
884 else
885 cunr *= std::complex<T>(1,2);
886 tmv::UpperTriMatrix<std::complex<T>,D|S> cu = cunr;
887 Assert(cunr == cu,"TriMatrix == TriMatrix<N,R>");
888 Assert(cunr.view() == cu.view(),"TriMatrix View");
889 Assert(cunr.transpose() == cu.transpose(),"TriMatrix transpose");
890 Assert(cunr.conjugate() == cu.conjugate(),"TriMatrix conjugate");
891 Assert(cunr.adjoint() == cu.adjoint(),"TriMatrix adjoint");
892 Assert(cunr.offDiag() == cu.offDiag(),"TriMatrix offDiag");
893 Assert(cunr.realPart() == cu.realPart(),"TriMatrix realPart");
894 if (D == tmv::NonUnitDiag)
895 Assert(cunr.imagPart() == cu.imagPart(),"TriMatrix imagPart");
896 Assert(cunr.subMatrix(0,N/2,N/2,N) == cu.subMatrix(0,N/2,N/2,N),"TriMatrix subMatrix");
897 Assert(cunr.subTriMatrix(0,N/2) == cu.subTriMatrix(0,N/2),"TriMatrix subTriMatrix");
898 Assert(cunr.subTriMatrix(N/2,N) == cu.subTriMatrix(N/2,N),"TriMatrix subTriMatrix");
899
900 tmv::LowerTriMatrix<std::complex<T>,tmv::NonUnitDiag|tmv::RowMajor> clnr = l;
901 if (D == tmv::UnitDiag)
902 clnr.offDiag() *= std::complex<T>(1,2);
903 else
904 clnr *= std::complex<T>(1,2);
905 tmv::LowerTriMatrix<std::complex<T>,D|S> cl = clnr;
906 Assert(clnr == cl,"TriMatrix == TriMatrix<N,R>");
907 Assert(clnr.view() == cl.view(),"TriMatrix View");
908 Assert(clnr.transpose() == cl.transpose(),"TriMatrix transpose");
909 Assert(clnr.conjugate() == cl.conjugate(),"TriMatrix conjugate");
910 Assert(clnr.adjoint() == cl.adjoint(),"TriMatrix adjoint");
911 Assert(clnr.offDiag() == cl.offDiag(),"TriMatrix offDiag");
912 Assert(clnr.realPart() == cl.realPart(),"TriMatrix realPart");
913 if (D == tmv::NonUnitDiag)
914 Assert(clnr.imagPart() == cl.imagPart(),"TriMatrix imagPart");
915 Assert(clnr.subMatrix(N/2,N,0,N/2) == cl.subMatrix(N/2,N,0,N/2),"TriMatrix subMatrix");
916 Assert(clnr.subTriMatrix(0,N/2) == cl.subTriMatrix(0,N/2),"TriMatrix subTriMatrix");
917 Assert(clnr.subTriMatrix(N/2,N) == cl.subTriMatrix(N/2,N),"TriMatrix subTriMatrix");
918 }
919
920 template <class T, tmv::DiagType D, tmv::StorageType S>
TestBasicTriMatrix_IO()921 static void TestBasicTriMatrix_IO()
922 {
923 const int N = 10;
924
925 if (showstartdone) {
926 std::cout<<"Start TestBasicTriMatrix_IO\n";
927 std::cout<<"T = "<<tmv::TMV_Text(T())<<std::endl;
928 std::cout<<"D = "<<tmv::TMV_Text(D)<<std::endl;
929 std::cout<<"S = "<<tmv::TMV_Text(S)<<std::endl;
930 std::cout<<"N = "<<N<<std::endl;
931 }
932
933 tmv::UpperTriMatrix<T,D|S> u(N);
934 tmv::UpperTriMatrix<CT,D|S> cu(N);
935 tmv::LowerTriMatrix<T,D|S> l(N);
936 tmv::LowerTriMatrix<CT,D|S> cl(N);
937
938 for (int i=0,k=1; i<N; ++i) for (int j=0; j<N; ++j, ++k) {
939 if (i < j || (D==tmv::NonUnitDiag && i==j)) {
940 u(i,j) = T(k);
941 cu(i,j) = CT(k,k+1000);
942 }
943 if (j < i || (D==tmv::NonUnitDiag && i==j)) {
944 l(i,j) = T(k);
945 cl(i,j) = CT(k,k+1000);
946 }
947 }
948 u(1,3) = l(3,1) = T(1.e-30);
949 cu(1,3) = cl(3,1) = CT(T(1.e-30),T(1.e-30));
950 u(5,6) = l(6,5) = T(9.e-3);
951 cu(5,6) = cl(6,5) = CT(T(9.e-3),T(9.e-3));
952 cu(5,7) = cl(7,5) = CT(T(9),T(9.e-3));
953 u(4,7) = l(7,4) = T(0.123456789);
954 cu(4,7) = cl(7,4) = CT(T(3.123456789),T(6.987654321));
955
956 // First check clipping function...
957 tmv::UpperTriMatrix<T> u2 = u;
958 tmv::UpperTriMatrix<CT> cu2 = cu;
959 tmv::LowerTriMatrix<T> l2 = l;
960 tmv::LowerTriMatrix<CT> cl2 = cl;
961 if (!std::numeric_limits<T>::is_integer) {
962 u2.clip(T(1.e-2));
963 cu2.clip(T(1.e-2));
964 l2.clip(T(1.e-2));
965 cl2.clip(T(1.e-2));
966 }
967 tmv::UpperTriMatrix<T> u3 = u;
968 tmv::UpperTriMatrix<CT> cu3 = cu;
969 tmv::LowerTriMatrix<T> l3 = l;
970 tmv::LowerTriMatrix<CT> cl3 = cl;
971 u3(1,3) = l3(3,1) = T(0);
972 cu3(1,3) = cl3(3,1) = T(0);
973 u3(5,6) = l3(6,5) = T(0); // Others shouldn't get clipped.
974 Assert(u2 == u3,"UpperTriMatrix clip");
975 Assert(cu2 == cu3,"Complex UpperTriMatrix clip");
976 Assert(l2 == l3,"LowerTriMatrix clip");
977 Assert(cl2 == cl3,"Complex LowerTriMatrix clip");
978
979 // However, ThreshIO for complex works slightly differently than clip.
980 // It clips _either_ the real or imag component, so now cu2(5,6) and
981 // cu2(5,7) need to be modified.
982 cu2(5,6) = cu3(5,6) = cl2(6,5) = cl3(6,5) = T(0);
983 cu2(5,7) = cu3(5,7) = cl2(7,5) = cl3(7,5) = T(9);
984
985 // Write matrices with 4 different styles
986 std::ofstream fout("tmvtest_trimatrix_io.dat");
987 Assert(bool(fout),"Couldn't open tmvtest_trimatrix_io.dat for output");
988 fout << u << std::endl;
989 fout << l << std::endl;
990 fout << cu << std::endl;
991 fout << cl << std::endl;
992 fout << tmv::CompactIO() << u << std::endl;
993 fout << tmv::CompactIO() << l << std::endl;
994 fout << tmv::CompactIO() << cu << std::endl;
995 fout << tmv::CompactIO() << cl << std::endl;
996 fout << tmv::ThreshIO(1.e-2).setPrecision(12) << u << std::endl;
997 fout << tmv::ThreshIO(1.e-2).setPrecision(12) << l << std::endl;
998 fout << tmv::ThreshIO(1.e-2).setPrecision(12) << cu << std::endl;
999 fout << tmv::ThreshIO(1.e-2).setPrecision(12) << cl << std::endl;
1000 tmv::IOStyle myStyle =
1001 tmv::CompactIO().setThresh(1.e-2).setPrecision(4).
1002 markup("Start","[",",","]","---","Done");
1003 fout << myStyle << u << std::endl;
1004 fout << myStyle << l << std::endl;
1005 fout << myStyle << cu << std::endl;
1006 fout << myStyle << cl << std::endl;
1007 fout.close();
1008
1009 // When using (the default) prec(6), these will be the values read in.
1010 u(4,7) = l(7,4) = T(0.123457);
1011 cu(4,7) = cl(7,4) = CT(T(3.12346),T(6.98765));
1012
1013 // When using prec(12), the full correct values will be read in.
1014
1015 // When using prec(4), these will be the values read in.
1016 u3(4,7) = l3(7,4) = T(0.1235);
1017 if (std::numeric_limits<T>::is_integer) cu3(4,7) = cl3(7,4) = CT(3,6);
1018 else cu3(4,7) = cl3(7,4) = CT(T(3.123),T(6.988));
1019
1020 // Read them back in
1021 tmv::UpperTriMatrix<T,D|tmv::RowMajor> xu1(N);
1022 tmv::LowerTriMatrix<T,D|tmv::RowMajor> xl1(N);
1023 tmv::UpperTriMatrix<CT,D|tmv::RowMajor> xcu1(N);
1024 tmv::LowerTriMatrix<CT,D|tmv::RowMajor> xcl1(N);
1025 std::ifstream fin("tmvtest_trimatrix_io.dat");
1026 Assert(bool(fin),"Couldn't open tmvtest_trimatrix_io.dat for input");
1027 fin >> xu1 >> xl1 >> xcu1 >> xcl1;
1028 Assert(EqualIO(u,xu1,EPS),"UpperTriMatrix I/O check normal");
1029 Assert(EqualIO(l,xl1,EPS),"LowerTriMatrix I/O check normal");
1030 Assert(EqualIO(cu,xcu1,EPS),"CUpperTriMatrix I/O check normal");
1031 Assert(EqualIO(cl,xcl1,EPS),"CLowerTriMatrix I/O check normal");
1032 fin >> tmv::CompactIO() >> xu1 >> tmv::CompactIO() >> xl1;
1033 fin >> tmv::CompactIO() >> xcu1 >> tmv::CompactIO() >> xcl1;
1034 Assert(EqualIO(u,xu1,EPS),"UpperTriMatrix I/O check compact");
1035 Assert(EqualIO(l,xl1,EPS),"LowerTriMatrix I/O check compact");
1036 Assert(EqualIO(cu,xcu1,EPS),"CUpperTriMatrix I/O check compact");
1037 Assert(EqualIO(cl,xcl1,EPS),"CLowerTriMatrix I/O check compact");
1038 fin >> xu1.view() >> xl1.view() >> xcu1.view() >> xcl1.view();
1039 Assert(EqualIO(u2,xu1,EPS),"UpperTriMatrix I/O check thresh");
1040 Assert(EqualIO(l2,xl1,EPS),"LowerTriMatrix I/O check thresh");
1041 Assert(EqualIO(cu2,xcu1,EPS),"CUpperTriMatrix I/O check thresh");
1042 Assert(EqualIO(cl2,xcl1,EPS),"CLowerTriMatrix I/O check thresh");
1043 fin >> myStyle >> xu1.view() >> myStyle >> xl1.view();
1044 fin >> myStyle >> xcu1.view() >> myStyle >> xcl1.view();
1045 Assert(EqualIO(u3,xu1,EPS),"UpperTriMatrix I/O check compact thresh & prec(4)");
1046 Assert(EqualIO(l3,xl1,EPS),"LowerTriMatrix I/O check compact thresh & prec(4)");
1047 Assert(EqualIO(cu3,xcu1,EPS),"CUpperTriMatrix I/O check compact thresh & prec(4)");
1048 Assert(EqualIO(cl3,xcl1,EPS),"CLowerTriMatrix I/O check compact thresh & prec(4)");
1049 fin.close();
1050
1051 // Repeat for column major
1052 tmv::UpperTriMatrix<T,D|tmv::ColMajor> xu2(N);
1053 tmv::LowerTriMatrix<T,D|tmv::ColMajor> xl2(N);
1054 tmv::UpperTriMatrix<CT,D|tmv::ColMajor> xcu2(N);
1055 tmv::LowerTriMatrix<CT,D|tmv::ColMajor> xcl2(N);
1056 fin.open("tmvtest_trimatrix_io.dat");
1057 Assert(bool(fin),"Couldn't open tmvtest_trimatrix_io.dat for input");
1058 fin >> xu2.view() >> xl2.view() >> xcu2.view() >> xcl2.view();
1059 Assert(EqualIO(u,xu2,EPS),"UpperTriMatrix I/O check normal");
1060 Assert(EqualIO(l,xl2,EPS),"LowerTriMatrix I/O check normal");
1061 Assert(EqualIO(cu,xcu2,EPS),"CUpperTriMatrix I/O check normal");
1062 Assert(EqualIO(cl,xcl2,EPS),"CLowerTriMatrix I/O check normal");
1063 fin >> tmv::CompactIO() >> xu2.view() >> tmv::CompactIO() >> xl2.view();
1064 fin >> tmv::CompactIO() >> xcu2.view() >> tmv::CompactIO() >> xcl2.view();
1065 Assert(EqualIO(u,xu2,EPS),"UpperTriMatrix I/O check compact");
1066 Assert(EqualIO(l,xl2,EPS),"LowerTriMatrix I/O check compact");
1067 Assert(EqualIO(cu,xcu2,EPS),"CUpperTriMatrix I/O check compact");
1068 Assert(EqualIO(cl,xcl2,EPS),"CLowerTriMatrix I/O check compact");
1069 fin >> xu2 >> xl2 >> xcu2 >> xcl2;
1070 Assert(EqualIO(u2,xu2,EPS),"UpperTriMatrix I/O check thresh");
1071 Assert(EqualIO(l2,xl2,EPS),"LowerTriMatrix I/O check thresh");
1072 Assert(EqualIO(cu2,xcu2,EPS),"CUpperTriMatrix I/O check thresh");
1073 Assert(EqualIO(cl2,xcl2,EPS),"CLowerTriMatrix I/O check thresh");
1074 fin >> myStyle >> xu2 >> myStyle >> xl2;
1075 fin >> myStyle >> xcu2 >> myStyle >> xcl2;
1076 Assert(EqualIO(u3,xu2,EPS),"UpperTriMatrix I/O check compact thresh & prec(4)");
1077 Assert(EqualIO(l3,xl2,EPS),"LowerTriMatrix I/O check compact thresh & prec(4)");
1078 Assert(EqualIO(cu3,xcu2,EPS),"CUpperTriMatrix I/O check compact thresh & prec(4)");
1079 Assert(EqualIO(cl3,xcl2,EPS),"CLowerTriMatrix I/O check compact thresh & prec(4)");
1080 fin.close();
1081
1082 // And repeat for matrices that need to be resized.
1083 // Also check switching the default IOStyle.
1084 tmv::CompactIO().makeDefault();
1085 tmv::UpperTriMatrix<T,D> zu1,zu2,zu3,zu4;
1086 tmv::LowerTriMatrix<T,D> zl1,zl2,zl3,zl4;
1087 tmv::UpperTriMatrix<CT,D> zcu1,zcu2,zcu3,zcu4;
1088 tmv::LowerTriMatrix<CT,D> zcl1,zcl2,zcl3,zcl4;
1089 fin.open("tmvtest_trimatrix_io.dat");
1090 Assert(bool(fin),"Couldn't open tmvtest_trimatrix_io.dat for input");
1091 fin >> tmv::NormalIO() >> zu1 >> tmv::NormalIO() >> zl1;
1092 fin >> tmv::NormalIO() >> zcu1 >> tmv::NormalIO() >> zcl1;
1093 Assert(EqualIO(u,zu1,EPS),"UpperTriMatrix I/O check normal");
1094 Assert(EqualIO(l,zl1,EPS),"LowerTriMatrix I/O check normal");
1095 Assert(EqualIO(cu,zcu1,EPS),"CUpperTriMatrix I/O check normal");
1096 Assert(EqualIO(cl,zcl1,EPS),"CLowerTriMatrix I/O check normal");
1097 fin >> zu2 >> zl2 >> zcu2 >> zcl2;
1098 Assert(EqualIO(u,zu2,EPS),"UpperTriMatrix I/O check compact");
1099 Assert(EqualIO(l,zl2,EPS),"LowerTriMatrix I/O check compact");
1100 Assert(EqualIO(cu,zcu2,EPS),"CUpperTriMatrix I/O check compact");
1101 Assert(EqualIO(cl,zcl2,EPS),"CLowerTriMatrix I/O check compact");
1102 fin >> tmv::NormalIO() >> zu3 >> tmv::NormalIO() >> zl3;
1103 fin >> tmv::NormalIO() >> zcu3 >> tmv::NormalIO() >> zcl3;
1104 Assert(EqualIO(u2,zu3,EPS),"UpperTriMatrix I/O check thresh");
1105 Assert(EqualIO(l2,zl3,EPS),"LowerTriMatrix I/O check thresh");
1106 Assert(EqualIO(cu2,zcu3,EPS),"CUpperTriMatrix I/O check thresh");
1107 Assert(EqualIO(cl2,zcl3,EPS),"CLowerTriMatrix I/O check thresh");
1108 fin >> myStyle >> zu4 >> myStyle >> zl4;
1109 fin >> myStyle >> zcu4 >> myStyle >> zcl4;
1110 Assert(EqualIO(u3,zu4,EPS),"UpperTriMatrix I/O check compact thresh & prec(4)");
1111 Assert(EqualIO(l3,zl4,EPS),"LowerTriMatrix I/O check compact thresh & prec(4)");
1112 Assert(EqualIO(cu3,zcu4,EPS),"CUpperTriMatrix I/O check compact thresh & prec(4)");
1113 Assert(EqualIO(cl3,zcl4,EPS),"CLowerTriMatrix I/O check compact thresh & prec(4)");
1114 fin.close();
1115
1116 // Also try reading with the opposite DiagType.
1117 // This will succeed if doing Unit -> NonUnit,
1118 // but it will throw an exception if doing NonUnit -> Unit.
1119 fin.open("tmvtest_trimatrix_io.dat");
1120 Assert(bool(fin),"Couldn't open tmvtest_trimatrix_io.dat for input");
1121 if (D==tmv::UnitDiag) {
1122 tmv::UpperTriMatrix<T,tmv::NonUnitDiag> zu5,zu6,zu7,zu8;
1123 tmv::LowerTriMatrix<T,tmv::NonUnitDiag> zl5,zl6,zl7,zl8;
1124 tmv::UpperTriMatrix<CT,tmv::NonUnitDiag> zcu5,zcu6,zcu7,zcu8;
1125 tmv::LowerTriMatrix<CT,tmv::NonUnitDiag> zcl5,zcl6,zcl7,zcl8;
1126 fin >> tmv::NormalIO() >> zu5 >> tmv::NormalIO() >> zl5;
1127 fin >> tmv::NormalIO() >> zcu5 >> tmv::NormalIO() >> zcl5;
1128 Assert(EqualIO(u,zu5,EPS),"UpperTriMatrix I/O check normal Unit->NonUnit");
1129 Assert(EqualIO(l,zl5,EPS),"LowerTriMatrix I/O check normal Unit->NonUnit");
1130 Assert(EqualIO(cu,zcu5,EPS),"CUpperTriMatrix I/O check normal Unit->NonUnit");
1131 Assert(EqualIO(cl,zcl5,EPS),"CLowerTriMatrix I/O check normal Unit->NonUnit");
1132 fin >> zu6 >> zl6 >> zcu6 >> zcl6;
1133 Assert(EqualIO(u,zu6,EPS),"UpperTriMatrix I/O check compact Unit->NonUnit");
1134 Assert(EqualIO(l,zl6,EPS),"LowerTriMatrix I/O check compact Unit->NonUnit");
1135 Assert(EqualIO(cu,zcu6,EPS),"CUpperTriMatrix I/O check compact Unit->NonUnit");
1136 Assert(EqualIO(cl,zcl6,EPS),"CLowerTriMatrix I/O check compact Unit->NonUnit");
1137 fin >> tmv::NormalIO() >> zu7 >> tmv::NormalIO() >> zl7;
1138 fin >> tmv::NormalIO() >> zcu7 >> tmv::NormalIO() >> zcl7;
1139 Assert(EqualIO(u2,zu7,EPS),"UpperTriMatrix I/O check thresh Unit->NonUnit");
1140 Assert(EqualIO(l2,zl7,EPS),"LowerTriMatrix I/O check thresh Unit->NonUnit");
1141 Assert(EqualIO(cu2,zcu7,EPS),"CUpperTriMatrix I/O check thresh Unit->NonUnit");
1142 Assert(EqualIO(cl2,zcl7,EPS),"CLowerTriMatrix I/O check thresh Unit->NonUnit");
1143 fin >> myStyle >> zu8 >> myStyle >> zl8;
1144 fin >> myStyle >> zcu8 >> myStyle >> zcl8;
1145 Assert(EqualIO(u3,zu8,EPS),"UpperTriMatrix I/O check compact thresh Unit->NonUnit");
1146 Assert(EqualIO(l3,zl8,EPS),"LowerTriMatrix I/O check compact thresh Unit->NonUnit");
1147 Assert(EqualIO(cu3,zcu8,EPS),"CUpperTriMatrix I/O check compact thresh Unit->NonUnit");
1148 Assert(EqualIO(cl3,zcl8,EPS),"CLowerTriMatrix I/O check compact thresh Unit->NonUnit");
1149 } else {
1150 #ifndef NOTHROW
1151 tmv::UpperTriMatrix<T,tmv::UnitDiag> zu5;
1152 try {
1153 fin >> tmv::NormalIO() >> zu5;
1154 Assert(false,"Throw ReadError for UnitDiag read of NonUnitDiag");
1155 } catch (tmv::ReadError&) {
1156 Assert(true,"Catch ReadError for UnitDiag read of NonUnitDiag");
1157 }
1158 // stream is already corrupted at this point, so no point
1159 // in continuing on with other matrices.
1160 #endif
1161 }
1162 fin.close();
1163 tmv::IOStyle::revertDefault();
1164
1165 // Finally, check that the NormalIO can be read in as a regular matrix.
1166 tmv::Matrix<T> zm1,zm2;
1167 tmv::Matrix<CT> zcm1,zcm2;
1168 fin.open("tmvtest_trimatrix_io.dat");
1169 Assert(bool(fin),"Couldn't open tmvtest_trimatrix_io.dat for input");
1170 fin >> zm1 >> zm2 >> zcm1 >> zcm2;
1171 Assert(EqualIO(u,zm1,EPS),"UpperTriMatrix -> Matrix I/O check");
1172 Assert(EqualIO(l,zm2,EPS),"LowerTriMatrix -> Matrix I/O check");
1173 Assert(EqualIO(cu,zcm1,EPS),"CUpperTriMatrix -> CMatrix I/O check");
1174 Assert(EqualIO(cl,zcm2,EPS),"CLowerTriMatrix -> CMatrix I/O check");
1175 fin.close();
1176
1177 #if XTEST == 0
1178 std::remove("tmvtest_trimatrix_io.dat");
1179 #endif
1180 }
1181
1182 template <class T, tmv::DiagType D, tmv::StorageType S>
TestBasicTriMatrix()1183 static void TestBasicTriMatrix()
1184 {
1185 TestBasicUpperTriMatrix_1<T,D,S>();
1186 TestBasicLowerTriMatrix_1<T,D,S>();
1187 TestBasicTriMatrix_2<T,D,S>();
1188 TestBasicTriMatrix_IO<T,D,S>();
1189 }
1190
TestTriMatrix()1191 template <class T> void TestTriMatrix()
1192 {
1193 #if 1
1194 TestBasicTriMatrix<T,tmv::UnitDiag,tmv::RowMajor>();
1195 TestBasicTriMatrix<T,tmv::NonUnitDiag,tmv::RowMajor>();
1196 TestBasicTriMatrix<T,tmv::UnitDiag,tmv::ColMajor>();
1197 TestBasicTriMatrix<T,tmv::NonUnitDiag,tmv::ColMajor>();
1198 std::cout<<"TriMatrix<"<<tmv::TMV_Text(T())<<"> passed all basic tests\n";
1199 #endif
1200
1201 #if 1
1202 TestAllAliasMultUL<T>();
1203 std::cout<<"TriMatrix<"<<tmv::TMV_Text(T())<<"> passed all aliased multiplication tests\n";
1204 #endif
1205
1206 #if 1
1207 TestTriMatrixArith_A1a<T>();
1208 TestTriMatrixArith_A1b<T>();
1209 TestTriMatrixArith_A2<T>();
1210 TestTriMatrixArith_A3<T>();
1211 TestTriMatrixArith_A4a<T>();
1212 TestTriMatrixArith_A4b<T>();
1213 TestTriMatrixArith_A5a<T>();
1214 TestTriMatrixArith_A5b<T>();
1215 TestTriMatrixArith_A6a<T>();
1216 TestTriMatrixArith_A6b<T>();
1217 TestTriMatrixArith_A6c<T>();
1218 std::cout<<"TriMatrix<"<<tmv::TMV_Text(T())<<"> (Tri/Tri) Arithmetic passed all tests\n";
1219 #endif
1220 #if 1
1221 TestTriMatrixArith_B4a<T>();
1222 TestTriMatrixArith_B4b<T>();
1223 TestTriMatrixArith_B5a<T>();
1224 TestTriMatrixArith_B5b<T>();
1225 TestTriMatrixArith_B6a<T>();
1226 TestTriMatrixArith_B6b<T>();
1227 std::cout<<"TriMatrix<"<<tmv::TMV_Text(T())<<"> (Matrix/Tri) Arithmetic passed all tests\n";
1228 #endif
1229 #if 1
1230 TestTriMatrixArith_C4a<T>();
1231 TestTriMatrixArith_C4b<T>();
1232 TestTriMatrixArith_C5a<T>();
1233 TestTriMatrixArith_C5b<T>();
1234 TestTriMatrixArith_C6a<T>();
1235 TestTriMatrixArith_C6b<T>();
1236 std::cout<<"TriMatrix<"<<tmv::TMV_Text(T())<<"> (Diag/Tri) Arithmetic passed all tests\n";
1237 #endif
1238 }
1239
1240 #ifdef TEST_DOUBLE
1241 template void TestTriMatrix<double>();
1242 #endif
1243 #ifdef TEST_FLOAT
1244 template void TestTriMatrix<float>();
1245 #endif
1246 #ifdef TEST_LONGDOUBLE
1247 template void TestTriMatrix<long double>();
1248 #endif
1249 #ifdef TEST_INT
1250 template void TestTriMatrix<int>();
1251 #endif
1252