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