1
2 #include "TMV_Test.h"
3 #include "TMV_Test_1.h"
4 #include "TMV.h"
5
6 template <class T>
TestMatrixDet()7 void TestMatrixDet()
8 {
9 tmv::Matrix<T> m0(0,0);
10
11 tmv::Matrix<T> m1(1,1);
12 m1 <<
13 6;
14
15 tmv::Matrix<T> m2(2,2);
16 m2 <<
17 6, 1,
18 3, 2;
19
20 tmv::Matrix<T> m3(3,3);
21 m3 <<
22 6, 1, 2,
23 3, 2, 1,
24 5, 6, 7;
25
26 tmv::Matrix<T> m4(4,4);
27 m4 <<
28 6, 1, 2, 4,
29 3, 2, 1, 8,
30 5, 6, 7, 1,
31 2, 9, 8, 3;
32
33 tmv::Matrix<T> m5(5,5);
34 m5 <<
35 6, 1, 2, 4, 1,
36 3, 2, 1, 8, 4,
37 5, 6, 7, 1, 2,
38 2, 9, 8, 3, 7,
39 1, 4, 5, 9, 3;
40
41 tmv::Matrix<T> m6(5,5);
42 m6 <<
43 6, 5, 2, 4, 2,
44 3, 0, 9, 3, 6,
45 5, 5, 7, 1, 8,
46 2, 10, 8, 3, 4,
47 1, 15, 5, 9, 4;
48
49 tmv::Matrix<T> m7(5,5);
50 m7 <<
51 6, 1, 0, 4, 1,
52 3, 0, 0, 8, 0,
53 5, 6, 7, 1, 2,
54 2, 9, 0, 3, 0,
55 0, 0, 0, 9, 0;
56
57 tmv::Matrix<T> m8(5,5);
58 m8 <<
59 6, 1, 2, 4, 1,
60 3, 2, 8, 8, 4,
61 5, 6, 4, 1, 2,
62 2, 9, 8, 3, 4,
63 1, 4, 6, 9, 3;
64
65 tmv::Matrix<T> m9 = m5.row(0) ^ m5.row(0);
66
67 tmv::Matrix<T> m10 = m5.rowRange(0,4).transpose() * m5.rowRange(0,4);
68
69 tmv::Matrix<T> m11(10,10);
70 m11 <<
71 0, 2, 3, 0, 0, 2, 6, 0, 1, 0,
72 0, 0, 2, 6, 0, 1, 0, 2, 1, 0,
73 1, 4, 0, 0, 3, 0, 3, 0, 0, 3,
74 7, 0, 0, 0, 8, 0, 4, 0, 0, 5,
75 0, 0, 0, 3, 0, 1, 0, 3, 1, 0,
76 0, 2, 1, 0, 4, 0, 9, 0, 0, 6,
77 0, 4, 7, 6, 0, 1, 0, 0, 4, 0,
78 0, 0, 0, 1, 0, 5, 0, 1, 0, 0,
79 4, 0, 0, 0, 6, 0, 0, 0, 0, 4,
80 0, 1, 8, 0, 0, 0, 3, 0, 9, 2;
81
82 tmv::Matrix<T> m12(10,10);
83 m12 <<
84 9, 2, 3, 5, 3, 2, 6, 1, 1, 8,
85 5, 6, 2, 6, 3, 1, 2, 2, 1, 8,
86 1, 4, 4, 7, 3, 3, 3, 1, 9, 3,
87 7, 8, 4, 5, 8, 8, 4, 6, 3, 5,
88 2, 1, 2, 3, 8, 1, 4, 3, 1, 6,
89 4, 2, 1, 3, 4, 3, 9, 1, 4, 6,
90 7, 4, 7, 6, 8, 1, 3, 8, 4, 4,
91 8, 8, 3, 1, 1, 5, 8, 1, 3, 4,
92 4, 1, 4, 2, 6, 1, 5, 8, 3, 4,
93 3, 1, 8, 8, 3, 2, 3, 8, 9, 2;
94
95 tmv::Matrix<std::complex<T> > cm1(5,5);
96 cm1.realPart() = m5;
97 cm1.imagPart() = m7.transpose();
98
99 tmv::Matrix<std::complex<T> > cm2(10,10);
100 for(int i=0;i<10;++i) for(int j=0;j<10;++j)
101 if (m11(i,j) == T(0))
102 cm2(i,j) = T(0);
103 else
104 cm2(i,j) = std::complex<T>(m11(i,j)-5, 3-m11(i,j));
105
106 tmv::Matrix<std::complex<T> > cm3(10,10);
107 cm3.realPart() = m12; cm3.realPart().addToAll(-4);
108 cm3.imagPart() = T(2)*m12.transpose();
109 for(int i=0;i<10;++i) for(int j=0;j<10;++j) {
110 // The intent here is to get the determinant as close to MAX_INT
111 // for a 32 bit integer as possible (in either component) without
112 // going over. i.e. close to 2e9.
113 // This is to prove that the algorithm works for any matrix that
114 // has its answer expressible as an int (or complex<int> in this case).
115 if (cm3.realPart()(i,j) >= T(5)) cm3.realPart()(i,j) = T(0);
116 if (cm3.imagPart()(i,j) > T(5)) cm3.imagPart()(i,j) -= T(10);
117 }
118
119 if (showacc) {
120 std::cout<<"m0: "<<m0.det()<<std::endl;
121 std::cout<<"diff = "<<tmv::TMV_ABS2(m0.det())<<" "<<EPS<<std::endl;
122 std::cout<<"m1: "<<m1.det()<<std::endl;
123 std::cout<<"diff = "<<tmv::TMV_ABS2(m1.det()-6)<<" "<<
124 6*NormSq(m1)*EPS<<std::endl;
125 std::cout<<"m2: "<<m2.det()<<std::endl;
126 std::cout<<"diff = "<<tmv::TMV_ABS2(m2.det()-9)<<" "<<
127 9*NormSq(m2)*EPS<<std::endl;
128 std::cout<<"m3: "<<m3.det()<<std::endl;
129 std::cout<<"diff = "<<tmv::TMV_ABS2(m3.det()-48)<<" "<<
130 48*NormSq(m3)*EPS<<std::endl;
131 std::cout<<"m4: "<<m4.det()<<std::endl;
132 std::cout<<"diff = "<<tmv::TMV_ABS2(m4.det()+66)<<" "<<
133 66*NormSq(m4)*EPS<<std::endl;
134 std::cout<<"m5: "<<m5.det()<<std::endl;
135 std::cout<<"diff = "<<tmv::TMV_ABS2(m5.det()+1118)<<" "<<
136 1000*NormSq(m5)*EPS<<std::endl;
137 std::cout<<"m6: "<<m6.det()<<std::endl;
138 std::cout<<"diff = "<<tmv::TMV_ABS2(m6.det()+15360)<<" "<<
139 15000*NormSq(m6)*EPS<<std::endl;
140 std::cout<<"m7: "<<m7.det()<<std::endl;
141 std::cout<<"diff = "<<tmv::TMV_ABS2(m7.det()+1701)<<" "<<
142 1700*NormSq(m7)*EPS<<std::endl;
143 std::cout<<"m8: "<<m8.det()<<std::endl;
144 std::cout<<"diff = "<<tmv::TMV_ABS2(m8.det())<<" "<<
145 1000*NormSq(m8)*EPS<<std::endl;
146 std::cout<<"m9: "<<m9.det()<<std::endl;
147 std::cout<<"diff = "<<tmv::TMV_ABS2(m9.det())<<" "<<
148 1000*NormSq(m9)*EPS<<std::endl;
149 std::cout<<"m10: "<<m10.det()<<std::endl;
150 std::cout<<"diff = "<<tmv::TMV_ABS2(m10.det())<<" "<<
151 1000*NormSq(m10)*EPS<<std::endl;
152 std::cout<<"m11: "<<m11.det()<<std::endl;
153 std::cout<<"diff = "<<tmv::TMV_ABS2(m11.det()-251656)<<" "<<
154 250000*NormSq(m11)*EPS<<std::endl;
155 std::cout<<"m12: "<<m12.det()<<std::endl;
156 std::cout<<"diff = "<<tmv::TMV_ABS2(m12.det()-20762581)<<" "<<
157 20000000*NormSq(m12)*EPS<<std::endl;
158 std::cout<<"cm1: "<<cm1.det()<<std::endl;
159 std::cout<<"diff = "<<
160 tmv::TMV_ABS2(cm1.det()-std::complex<T>(-19796,-2165))<<" "<<
161 20000*NormSq(cm1)*EPS<<std::endl;
162 std::cout<<"cm2: "<<cm2.det()<<std::endl;
163 std::cout<<"diff = "<<
164 tmv::TMV_ABS2(cm2.det()-std::complex<T>(2386176,-1084032))<<" "<<
165 3000000*NormSq(cm2)*EPS<<std::endl;
166 std::cout<<"cm3: "<<cm3.det()<<std::endl;
167 std::cout<<"diff = "<<
168 tmv::TMV_ABS2(cm3.det()-std::complex<T>(-1923395548,811584412))<<" "<<
169 2000000000*NormSq(cm3)*EPS<<std::endl;
170 }
171 Assert(Equal2(m0.det(),1,EPS),"0x0 determinant");
172 Assert(Equal2(m1.det(),6,6*NormSq(m1)*EPS),"1x1 determinant");
173 Assert(Equal2(m2.det(),9,9*NormSq(m2)*EPS),"2x2 determinant");
174 Assert(Equal2(m3.det(),48,48*NormSq(m3)*EPS),"3x3 determinant");
175 Assert(Equal2(m4.det(),-66,66*NormSq(m4)*EPS),"4x4 determinant");
176 Assert(Equal2(m5.det(),-1118,1000*NormSq(m5)*EPS),"5x5 determinant");
177 Assert(Equal2(m6.det(),-15360,15000*NormSq(m6)*EPS),
178 "5x5 determinant with GCD factors");
179 Assert(Equal2(m7.det(),-1701,1700*NormSq(m7)*EPS),
180 "5x5 determinant permuted triangle");
181 Assert(Equal2(m8.det(),0,1000*NormSq(m8)*EPS),"5x5 determinant singular");
182 Assert(Equal2(m9.det(),0,1000*NormSq(m9)*EPS),
183 "5x5 determinant rank 1 outer product");
184 Assert(Equal2(m10.det(),0,1000*NormSq(m10)*EPS),
185 "5x5 determinant rank 4 outer product");
186 Assert(Equal2(m11.det(),251656,250000*NormSq(m11)*EPS),
187 "10x10 determinant permuted band");
188 Assert(Equal2(m12.det(),20762581,20000000*NormSq(m11)*EPS),
189 "10x10 determinant full");
190 Assert(Equal2(cm1.det(),std::complex<T>(-19796,-2165),
191 20000*NormSq(cm1)*EPS),
192 "5x5 determinant complex");
193 Assert(Equal2(cm2.det(),std::complex<T>(2386176,-1084032),
194 4000000*NormSq(cm2)*EPS),
195 "10x10 determinant complex permuted band");
196 Assert(Equal2(cm3.det(),std::complex<T>(-1923395548,811584412),
197 2000000000*NormSq(cm3)*EPS),
198 "10x10 determinant complex full");
199
200 std::cout<<"Matrix<"<<tmv::TMV_Text(T())<<"> passed all determinant tests\n";
201 }
202
203 #ifdef TEST_DOUBLE
204 template void TestMatrixDet<double>();
205 #endif
206 #ifdef TEST_FLOAT
207 template void TestMatrixDet<float>();
208 #endif
209 #ifdef TEST_LONGDOUBLE
210 template void TestMatrixDet<long double>();
211 #endif
212 #ifdef TEST_INT
213 template void TestMatrixDet<int>();
214 #endif
215