36 #ifndef MAT_PURISTEPINFODEBUG
37 #define MAT_PURISTEPINFODEBUG
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);
Interval< Treal > XmX2EuclNormExact
Definition: PuriStepInfoDebug.h:83
Definition: PuriStepInfoDebug.h:43
PuriStepInfoDebug()
Definition: PuriStepInfoDebug.h:74
Interval< Treal > lminExact
Definition: PuriStepInfoDebug.h:82
Interval< Treal > lmaxExact
Definition: PuriStepInfoDebug.h:81
#define ASSERTALWAYS(x)
Definition: DebugPolicies.h:83
Definition: allocate.cc:30
Treal template_blas_fabs(Treal x)
Definition: Interval.h:44
Interval< Treal > homoExact
Definition: PuriStepInfoDebug.h:79
static void syev(const char *jobz, const char *uplo, const int *n, T *a, const int *lda, T *w, T *work, const int *lwork, int *info)
Definition: gblas.h:380
PuriStepInfoDebug()
Definition: PuriStepInfoDebug.h:55
#define ASSERTDEBUG(x)
Definition: DebugPolicies.h:85
void checkIntervals(Interval< Treal > const &eigInterval, Interval< Treal > const &homo, Interval< Treal > const &lumo, Interval< Treal > const &XmX2EuclNorm, const char *descriptionString) const
Definition: PuriStepInfoDebug.h:46
Definition: DebugPolicies.h:89
Copyright(c) Emanuel Rubensson 2006.
void computeExactValues(Tmatrix const &X, Tmatrix const &X2, int const n, int const nocc)
Definition: PuriStepInfoDebug.h:52
Interval< Treal > lumoExact
Definition: PuriStepInfoDebug.h:80