36 #ifndef MAT_PURISTEPINFODEBUG
37 #define MAT_PURISTEPINFODEBUG
40 #define __ID__ "$Id: PuriStepInfoDebug.h 4437 2012-07-05 09:01:18Z elias $"
42 template<
typename Treal,
typename TdebugPolicy>
50 const char* descriptionString)
const {}
51 template<
typename Tmatrix>
53 int const n,
int const nocc) {}
62 template<
typename Treal>
69 const char* descriptionString)
const;
70 template<
typename Tmatrix>
72 int const n,
int const nocc);
75 :homoExact(), lumoExact(),
76 lmaxExact(), lminExact(),
87 template<
typename Treal>
93 const char* descriptionString)
const {
98 std::cout <<
"PuriStepInfoDebug::checkIntervals failed, descriptionString = '"
99 << descriptionString <<
"'" << std::endl;
103 std::cout<<std::endl<<
" eigInterval "<<
104 std::setprecision(25)<<eigInterval<<std::endl
105 <<
" lminExact "<<std::setprecision(25)
106 <<lminExact<<std::endl;
109 std::cout<<std::endl<<
" eigInterval "<<
110 std::setprecision(25)<<eigInterval<<std::endl
111 <<
" lmaxExact "<<lmaxExact<<std::endl;
118 std::cout<<std::endl<<
" homo "<<
119 std::setprecision(25)<<homo<<std::endl
120 <<
" homoExact "<<homoExact<<std::endl;
122 std::cout <<
"PuriStepInfoDebug::checkIntervals failed due to intersect(homo, homoExact) , descriptionString = '"
123 << descriptionString <<
"'" << std::endl;
128 std::cout<<std::endl<<
" lumo "<<
129 std::setprecision(25)<<lumo<<std::endl
130 <<
" lumoExact "<<lumoExact<<std::endl;
134 intersect(XmX2EuclNorm, XmX2EuclNormExact).empty()) {
135 std::cout<<std::endl<<
" XmX2EuclNorm "<<
136 std::setprecision(25)<<XmX2EuclNorm<<std::endl
137 <<
" XmX2EuclNormExact "<<XmX2EuclNormExact<<std::endl;
138 std::cout <<
"PuriStepInfoDebug::checkIntervals failed, descriptionString = '"
139 << descriptionString <<
"'" << std::endl;
142 intersect(XmX2EuclNorm, XmX2EuclNormExact).empty());
145 template<
typename Treal>
146 template<
typename Tmatrix>
149 int const n,
int const nocc) {
151 std::vector<Treal> full(n*n);
153 Treal* eigvals =
new Treal[n];
155 Treal* work =
new Treal[lwork];
157 Treal euclNormMatrix;
161 syev(
"N",
"U", &n, &full[0], &n, eigvals, work, &lwork, &info);
165 precSyev = euclNormMatrix * getRelPrecision<Treal>();
167 eigvals[n-nocc-1] + precSyev);
169 eigvals[n-nocc] + precSyev);
171 eigvals[0] + precSyev);
173 eigvals[n-1] + precSyev);
177 XmX2.fullMatrix(full);
178 syev(
"N",
"U", &n, &full[0], &n, eigvals, work, &lwork, &info);
182 precSyev = euclNormMatrix * getRelPrecision<Treal>();
185 std::cout<<
"EIGS XmX2: "<<std::endl;
186 for(
int ind = 0; ind < n; ++ind)
187 std::cout<<std::setprecision(20)<<eigvals[ind]<<std::endl;
188 std::cout<<
"end EIGS XmX2: "<<std::endl<<std::endl;
191 Treal lowBound = euclNormMatrix - precSyev > 0 ?
192 euclNormMatrix - precSyev : 0;
193 XmX2EuclNormExact =
Interval<Treal>(lowBound, euclNormMatrix + precSyev);