1
2 template <class T>
CopyBackM(const tmv::Matrix<std::complex<T>> & m0,tmv::SymMatrixView<std::complex<T>> & m)3 inline void CopyBackM(
4 const tmv::Matrix<std::complex<T> >& m0,
5 tmv::SymMatrixView<std::complex<T> >& m)
6 {
7 if (m.issym()) m = SymMatrixViewOf(m0,m.uplo());
8 else m = HermMatrixViewOf(m0,m.uplo());
9 }
10
11 #define MakeSymList(s,cs,pdc) \
12 std::vector<tmv::SymMatrix<T,tmv::Upper|tmv::RowMajor> > S1; \
13 std::vector<tmv::SymMatrix<T,tmv::Upper|tmv::ColMajor> > S2; \
14 std::vector<tmv::SymMatrix<T,tmv::Lower|tmv::RowMajor> > S3; \
15 std::vector<tmv::SymMatrix<T,tmv::Lower|tmv::ColMajor> > S4; \
16 std::vector<tmv::HermMatrix<T,tmv::Upper|tmv::RowMajor> > S5; \
17 std::vector<tmv::HermMatrix<T,tmv::Upper|tmv::ColMajor> > S6; \
18 std::vector<tmv::HermMatrix<T,tmv::Lower|tmv::RowMajor> > S7; \
19 std::vector<tmv::HermMatrix<T,tmv::Lower|tmv::ColMajor> > S8; \
20 std::vector<tmv::Matrix<T> > S9; \
21 std::vector<tmv::SymMatrix<std::complex<T>,tmv::Upper|tmv::RowMajor> > CS1; \
22 std::vector<tmv::SymMatrix<std::complex<T>,tmv::Upper|tmv::ColMajor> > CS2; \
23 std::vector<tmv::SymMatrix<std::complex<T>,tmv::Lower|tmv::RowMajor> > CS3; \
24 std::vector<tmv::SymMatrix<std::complex<T>,tmv::Lower|tmv::ColMajor> > CS4; \
25 std::vector<tmv::HermMatrix<std::complex<T>,tmv::Upper|tmv::RowMajor> > CS5; \
26 std::vector<tmv::HermMatrix<std::complex<T>,tmv::Upper|tmv::ColMajor> > CS6; \
27 std::vector<tmv::HermMatrix<std::complex<T>,tmv::Lower|tmv::RowMajor> > CS7; \
28 std::vector<tmv::HermMatrix<std::complex<T>,tmv::Lower|tmv::ColMajor> > CS8; \
29 std::vector<tmv::Matrix<std::complex<T> > > CS9; \
30 DoMakeSymList(s,cs,pdc,S1,S2,S3,S4,S5,S6,S7,S8,S9, \
31 CS1,CS2,CS3,CS4,CS5,CS6,CS7,CS8,CS9);
32
33 template <class T>
DoMakeSymList(std::vector<tmv::SymMatrixView<T>> & s,std::vector<tmv::SymMatrixView<std::complex<T>>> & cs,PosDefCode pdc,std::vector<tmv::SymMatrix<T,tmv::Upper|tmv::RowMajor>> & SUR,std::vector<tmv::SymMatrix<T,tmv::Upper|tmv::ColMajor>> & SUC,std::vector<tmv::SymMatrix<T,tmv::Lower|tmv::RowMajor>> & SLR,std::vector<tmv::SymMatrix<T,tmv::Lower|tmv::ColMajor>> & SLC,std::vector<tmv::HermMatrix<T,tmv::Upper|tmv::RowMajor>> & HUR,std::vector<tmv::HermMatrix<T,tmv::Upper|tmv::ColMajor>> & HUC,std::vector<tmv::HermMatrix<T,tmv::Lower|tmv::RowMajor>> & HLR,std::vector<tmv::HermMatrix<T,tmv::Lower|tmv::ColMajor>> & HLC,std::vector<tmv::Matrix<T>> & M,std::vector<tmv::SymMatrix<std::complex<T>,tmv::Upper|tmv::RowMajor>> & CSUR,std::vector<tmv::SymMatrix<std::complex<T>,tmv::Upper|tmv::ColMajor>> & CSUC,std::vector<tmv::SymMatrix<std::complex<T>,tmv::Lower|tmv::RowMajor>> & CSLR,std::vector<tmv::SymMatrix<std::complex<T>,tmv::Lower|tmv::ColMajor>> & CSLC,std::vector<tmv::HermMatrix<std::complex<T>,tmv::Upper|tmv::RowMajor>> & CHUR,std::vector<tmv::HermMatrix<std::complex<T>,tmv::Upper|tmv::ColMajor>> & CHUC,std::vector<tmv::HermMatrix<std::complex<T>,tmv::Lower|tmv::RowMajor>> & CHLR,std::vector<tmv::HermMatrix<std::complex<T>,tmv::Lower|tmv::ColMajor>> & CHLC,std::vector<tmv::Matrix<std::complex<T>>> & CM)34 inline void DoMakeSymList(
35 std::vector<tmv::SymMatrixView<T> >& s,
36 std::vector<tmv::SymMatrixView<std::complex<T> > >& cs,
37 PosDefCode pdc,
38 std::vector<tmv::SymMatrix<T,tmv::Upper|tmv::RowMajor> >& SUR,
39 std::vector<tmv::SymMatrix<T,tmv::Upper|tmv::ColMajor> >& SUC,
40 std::vector<tmv::SymMatrix<T,tmv::Lower|tmv::RowMajor> >& SLR,
41 std::vector<tmv::SymMatrix<T,tmv::Lower|tmv::ColMajor> >& SLC,
42 std::vector<tmv::HermMatrix<T,tmv::Upper|tmv::RowMajor> >& HUR,
43 std::vector<tmv::HermMatrix<T,tmv::Upper|tmv::ColMajor> >& HUC,
44 std::vector<tmv::HermMatrix<T,tmv::Lower|tmv::RowMajor> >& HLR,
45 std::vector<tmv::HermMatrix<T,tmv::Lower|tmv::ColMajor> >& HLC,
46 std::vector<tmv::Matrix<T> >& M,
47 std::vector<tmv::SymMatrix<std::complex<T>,tmv::Upper|tmv::RowMajor> >& CSUR,
48 std::vector<tmv::SymMatrix<std::complex<T>,tmv::Upper|tmv::ColMajor> >& CSUC,
49 std::vector<tmv::SymMatrix<std::complex<T>,tmv::Lower|tmv::RowMajor> >& CSLR,
50 std::vector<tmv::SymMatrix<std::complex<T>,tmv::Lower|tmv::ColMajor> >& CSLC,
51 std::vector<tmv::HermMatrix<std::complex<T>,tmv::Upper|tmv::RowMajor> >& CHUR,
52 std::vector<tmv::HermMatrix<std::complex<T>,tmv::Upper|tmv::ColMajor> >& CHUC,
53 std::vector<tmv::HermMatrix<std::complex<T>,tmv::Lower|tmv::RowMajor> >& CHLR,
54 std::vector<tmv::HermMatrix<std::complex<T>,tmv::Lower|tmv::ColMajor> >& CHLC,
55 std::vector<tmv::Matrix<std::complex<T> > >& CM)
56 {
57 const int N = std::numeric_limits<T>::is_integer ? 6 : 10;
58 const int RESERVE = 20;
59 SUR.reserve(RESERVE);
60 SUC.reserve(RESERVE);
61 SLR.reserve(RESERVE);
62 SLC.reserve(RESERVE);
63 HUR.reserve(RESERVE);
64 HUC.reserve(RESERVE);
65 HLR.reserve(RESERVE);
66 HLC.reserve(RESERVE);
67 M.reserve(RESERVE);
68 CSUR.reserve(RESERVE);
69 CSUC.reserve(RESERVE);
70 CSLR.reserve(RESERVE);
71 CSLC.reserve(RESERVE);
72 CHUR.reserve(RESERVE);
73 CHUC.reserve(RESERVE);
74 CHLR.reserve(RESERVE);
75 CHLC.reserve(RESERVE);
76 CM.reserve(RESERVE);
77
78 tmv::Matrix<T> a1(N,N);
79 tmv::Matrix<std::complex<T> > ca1(N,N);
80 tmv::Matrix<T> a2(2*N,2*N);
81 tmv::Matrix<std::complex<T> > ca2(2*N,2*N);
82
83 if (std::numeric_limits<T>::is_integer) {
84 for (int i=0; i<N; ++i) for (int j=0; j<N; ++j) {
85 a1(i,j) = T(-1+i/2);
86 ca1(i,j) = std::complex<T>(-1+i/2,-2+j/3);
87 }
88 for (int i=0; i<2*N; ++i) for (int j=0; j<2*N; ++j) {
89 a2(i,j) = T(-3+i);
90 ca2(i,j) = std::complex<T>(-3+i,-7+j);
91 }
92 a1.diag().addToAll(3);
93 a2.diag().addToAll(5);
94 ca1.diag().addToAll(3);
95 ca2.diag().addToAll(5);
96 } else {
97 for (int i=0; i<N; ++i) for (int j=0; j<N; ++j) {
98 a1(i,j) = T(3+i-5*j);
99 ca1(i,j) = std::complex<T>(3+i-5*j,2-3*i);
100 }
101 for (int i=0; i<2*N; ++i) for (int j=0; j<2*N; ++j) {
102 a2(i,j) = T(1-3*i+6*j);
103 ca2(i,j) = std::complex<T>(1-3*i+6*j,-4+2*j);
104 }
105
106 if (pdc == PosDef) {
107 a1 /= T(N*N);
108 a2 /= T(N*N);
109 a1.diag().addToAll(T(1));
110 a2.diag().addToAll(T(1));
111 for(int i=0;i<N;++i) a1.diag()(i) += T(i);
112 for(int i=0;i<2*N;++i) a2.diag()(i) += T(i);
113 ca1 /= T(N*N);
114 ca2 /= T(N*N);
115 ca1.diag().addToAll(T(1));
116 ca2.diag().addToAll(T(1));
117 for(int i=0;i<N;++i) ca1.diag()(i) += T(i);
118 for(int i=0;i<2*N;++i) ca2.diag()(i) += T(i);
119 } else if (pdc == InDef) {
120 a1 /= T(N);
121 a2 /= T(N);
122 a1.diag(0,2,4).setZero();
123 a2.diag(0,2,4).setZero();
124 a1(4,0) = a1(0,4) = T(50);
125 a1.diag(-1,6,9).addToAll(T(10));
126 a2.diag(-1,6,9).addToAll(T(10));
127 a1.diag(1,6,9).addToAll(T(10));
128 a2.diag(1,6,9).addToAll(T(10));
129 a1.row(9,0,5).addToAll(T(10));
130 a2.row(9,0,5).addToAll(T(10));
131 a1(3,3) = T(0);
132 a2(3,3) = T(0);
133 if (N > 10) {
134 a1.diag(0,10,N) *= T(0.0001);
135 a2.diag(0,10,N) *= T(0.0001);
136 }
137 ca1 /= T(N);
138 ca2 /= T(N);
139 ca1.diag(0,2,4).setZero();
140 ca2.diag(0,2,4).setZero();
141 ca1(4,0) = ca1(0,4) = T(50);
142 ca1.diag(-1,6,9).addToAll(T(10));
143 ca2.diag(-1,6,9).addToAll(T(10));
144 ca1.diag(1,6,9).addToAll(T(10));
145 ca2.diag(1,6,9).addToAll(T(10));
146 ca1.row(9,0,5).addToAll(T(10));
147 ca2.row(9,0,5).addToAll(T(10));
148 ca1(3,3) = T(0);
149 ca2(3,3) = T(0);
150 if (N > 10) {
151 ca1.diag(0,10,N) *= T(0.0001);
152 ca2.diag(0,10,N) *= T(0.0001);
153 }
154 } else {
155 a1.row(2).setZero();
156 a1.col(2).setZero();
157 a1.row(4) = a1.row(5);
158 a1.col(4) = a1.col(5);
159 a1(4,5) = a1(5,4) = a1(5,5) = a1(4,4);
160 a2.row(2).setZero();
161 a2.col(2).setZero();
162 a2.row(4) = a2.row(5);
163 a2.col(4) = a2.col(5);
164 a2(4,5) = a2(5,4) = a2(5,5) = a2(4,4);
165 ca1.row(2).setZero();
166 ca1.col(2).setZero();
167 ca1.row(4) = ca1.row(5);
168 ca1.col(4) = ca1.col(5);
169 ca1(4,5) = ca1(5,4) = ca1(5,5) = ca1(4,4);
170 ca2.row(2).setZero();
171 ca2.col(2).setZero();
172 ca2.row(4) = ca2.row(5);
173 ca2.col(4) = ca2.col(5);
174 ca2(4,5) = ca2(5,4) = ca2(5,5) = ca2(4,4);
175 }
176 }
177
178 SUR.push_back(tmv::SymMatrix<T,tmv::Upper|tmv::RowMajor>(a1));
179 HUR.push_back(tmv::HermMatrix<T,tmv::Upper|tmv::RowMajor>(a1));
180 CSUR.push_back(tmv::SymMatrix<std::complex<T>,tmv::Upper|tmv::RowMajor>(ca1));
181 CHUR.push_back(tmv::HermMatrix<std::complex<T>,tmv::Upper|tmv::RowMajor>(ca1));
182 s.push_back(SUR.back().view());
183 s.push_back(HUR.back().view());
184 cs.push_back(CSUR.back().view());
185 cs.push_back(CHUR.back().view());
186
187 SUC.push_back(tmv::SymMatrix<T,tmv::Upper|tmv::ColMajor>(a1));
188 HUC.push_back(tmv::HermMatrix<T,tmv::Upper|tmv::ColMajor>(a1));
189 CSUC.push_back(tmv::SymMatrix<std::complex<T>,tmv::Upper|tmv::ColMajor>(ca1));
190 CHUC.push_back(tmv::HermMatrix<std::complex<T>,tmv::Upper|tmv::ColMajor>(ca1));
191 s.push_back(SUC.back().view());
192 s.push_back(HUC.back().view());
193 cs.push_back(CSUC.back().view());
194 cs.push_back(CHUC.back().view());
195
196 #if (XTEST & 2)
197 SLR.push_back(tmv::SymMatrix<T,tmv::Lower|tmv::RowMajor>(a1));
198 HLR.push_back(tmv::HermMatrix<T,tmv::Lower|tmv::RowMajor>(a1));
199 CSLR.push_back(tmv::SymMatrix<std::complex<T>,tmv::Lower|tmv::RowMajor>(ca1));
200 CHLR.push_back(tmv::HermMatrix<std::complex<T>,tmv::Lower|tmv::RowMajor>(ca1));
201 s.push_back(SLR.back().view());
202 s.push_back(HLR.back().view());
203 cs.push_back(CSLR.back().view());
204 cs.push_back(CHLR.back().view());
205
206 SLC.push_back(tmv::SymMatrix<T,tmv::Lower|tmv::ColMajor>(a1));
207 HLC.push_back(tmv::HermMatrix<T,tmv::Lower|tmv::ColMajor>(a1));
208 CSLC.push_back(tmv::SymMatrix<std::complex<T>,tmv::Lower|tmv::ColMajor>(ca1));
209 CHLC.push_back(tmv::HermMatrix<std::complex<T>,tmv::Lower|tmv::ColMajor>(ca1));
210 s.push_back(SLC.back().view());
211 s.push_back(HLC.back().view());
212 cs.push_back(CSLC.back().view());
213 cs.push_back(CHLC.back().view());
214
215 M.push_back(tmv::Matrix<T>(a2));
216 CM.push_back(tmv::Matrix<std::complex<T> >(ca2));
217 s.push_back(SymMatrixViewOf(M.back(),tmv::Upper).subSymMatrix(0,2*N,2));
218 cs.push_back(SymMatrixViewOf(CM.back(),tmv::Upper).subSymMatrix(0,2*N,2));
219
220 M.push_back(tmv::Matrix<T>(a2));
221 CM.push_back(tmv::Matrix<std::complex<T> >(ca2));
222 CM.back().diag().imagPart().setZero();
223 s.push_back(HermMatrixViewOf(M.back(),tmv::Upper).subSymMatrix(0,2*N,2));
224 cs.push_back(HermMatrixViewOf(CM.back(),tmv::Upper).subSymMatrix(0,2*N,2));
225 #endif
226
227 TMVAssert(int(SUR.size()) <= RESERVE);
228 TMVAssert(int(SUC.size()) <= RESERVE);
229 TMVAssert(int(SLR.size()) <= RESERVE);
230 TMVAssert(int(SLC.size()) <= RESERVE);
231 TMVAssert(int(HUR.size()) <= RESERVE);
232 TMVAssert(int(HUC.size()) <= RESERVE);
233 TMVAssert(int(HLR.size()) <= RESERVE);
234 TMVAssert(int(HLC.size()) <= RESERVE);
235 TMVAssert(int(M.size()) <= RESERVE);
236 TMVAssert(int(CSUR.size()) <= RESERVE);
237 TMVAssert(int(CSUC.size()) <= RESERVE);
238 TMVAssert(int(CSLR.size()) <= RESERVE);
239 TMVAssert(int(CSLC.size()) <= RESERVE);
240 TMVAssert(int(CHUR.size()) <= RESERVE);
241 TMVAssert(int(CHUC.size()) <= RESERVE);
242 TMVAssert(int(CHLR.size()) <= RESERVE);
243 TMVAssert(int(CHLC.size()) <= RESERVE);
244 TMVAssert(int(CM.size()) <= RESERVE);
245 }
246
247
248