00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00036 #ifndef MAT_PURISTEPINFODEBUG
00037 #define MAT_PURISTEPINFODEBUG
00038 #include "matInclude.h"
00039 #include "Interval.h"
00040 #define __ID__ "$Id: PuriStepInfoDebug.h 4437 2012-07-05 09:01:18Z elias $"
00041 namespace mat {
00042 template<typename Treal, typename TdebugPolicy>
00043 class PuriStepInfoDebug : public TdebugPolicy {
00044 public:
00045
00046 inline void checkIntervals(Interval<Treal> const & eigInterval,
00047 Interval<Treal> const & homo,
00048 Interval<Treal> const & lumo,
00049 Interval<Treal> const & XmX2EuclNorm,
00050 const char* descriptionString) const {}
00051 template<typename Tmatrix>
00052 inline void computeExactValues(Tmatrix const & X, Tmatrix const & X2,
00053 int const n, int const nocc) {}
00054 protected:
00055 PuriStepInfoDebug(){}
00056 private:
00057 };
00058
00059
00060
00061
00062 template<typename Treal>
00063 class PuriStepInfoDebug<Treal, DebugLevelHigh> : public DebugLevelHigh {
00064 public:
00065 void checkIntervals(Interval<Treal> const & eigInterval,
00066 Interval<Treal> const & homo,
00067 Interval<Treal> const & lumo,
00068 Interval<Treal> const & XmX2EuclNorm,
00069 const char* descriptionString) const;
00070 template<typename Tmatrix>
00071 void computeExactValues(Tmatrix const & X, Tmatrix const & X2,
00072 int const n, int const nocc);
00073 protected:
00074 PuriStepInfoDebug()
00075 :homoExact(), lumoExact(),
00076 lmaxExact(), lminExact(),
00077 XmX2EuclNormExact()
00078 {}
00079 Interval<Treal> homoExact;
00080 Interval<Treal> lumoExact;
00081 Interval<Treal> lmaxExact;
00082 Interval<Treal> lminExact;
00083 Interval<Treal> XmX2EuclNormExact;
00084 private:
00085 };
00086
00087 template<typename Treal>
00088 void PuriStepInfoDebug<Treal, DebugLevelHigh>::
00089 checkIntervals(Interval<Treal> const & eigInterval,
00090 Interval<Treal> const & homo,
00091 Interval<Treal> const & lumo,
00092 Interval<Treal> const & XmX2EuclNorm,
00093 const char* descriptionString) const {
00094
00095
00096
00097 if(lminExact.empty())
00098 std::cout << "PuriStepInfoDebug::checkIntervals failed, descriptionString = '"
00099 << descriptionString << "'" << std::endl;
00100 ASSERTALWAYS(!lminExact.empty());
00101 Interval<Treal> zeroOneInt(0.0,1.0);
00102 if (Interval<Treal>::intersect(eigInterval, lminExact).empty())
00103 std::cout<<std::endl<<" eigInterval "<<
00104 std::setprecision(25)<<eigInterval<<std::endl
00105 <<" lminExact "<<std::setprecision(25)
00106 <<lminExact<<std::endl;
00107 ASSERTDEBUG(!Interval<Treal>::intersect(eigInterval, lminExact).empty());
00108 if (Interval<Treal>::intersect(eigInterval, lmaxExact).empty())
00109 std::cout<<std::endl<<" eigInterval "<<
00110 std::setprecision(25)<<eigInterval<<std::endl
00111 <<" lmaxExact "<<lmaxExact<<std::endl;
00112 ASSERTDEBUG(!Interval<Treal>::intersect(eigInterval, lmaxExact).empty());
00113
00114
00115
00116 if (!Interval<Treal>::intersect(zeroOneInt, homoExact).empty()) {
00117 if (Interval<Treal>::intersect(homo, homoExact).empty())
00118 std::cout<<std::endl<<" homo "<<
00119 std::setprecision(25)<<homo<<std::endl
00120 <<" homoExact "<<homoExact<<std::endl;
00121 if(Interval<Treal>::intersect(homo, homoExact).empty())
00122 std::cout << "PuriStepInfoDebug::checkIntervals failed due to intersect(homo, homoExact) , descriptionString = '"
00123 << descriptionString << "'" << std::endl;
00124 ASSERTDEBUG(!Interval<Treal>::intersect(homo, homoExact).empty());
00125 }
00126 if (!Interval<Treal>::intersect(zeroOneInt, lumoExact).empty()) {
00127 if (Interval<Treal>::intersect(lumo, lumoExact).empty())
00128 std::cout<<std::endl<<" lumo "<<
00129 std::setprecision(25)<<lumo<<std::endl
00130 <<" lumoExact "<<lumoExact<<std::endl;
00131 ASSERTDEBUG(!Interval<Treal>::intersect(lumo, lumoExact).empty());
00132 }
00133 if (Interval<Treal>::
00134 intersect(XmX2EuclNorm, XmX2EuclNormExact).empty()) {
00135 std::cout<<std::endl<<" XmX2EuclNorm "<<
00136 std::setprecision(25)<<XmX2EuclNorm<<std::endl
00137 <<" XmX2EuclNormExact "<<XmX2EuclNormExact<<std::endl;
00138 std::cout << "PuriStepInfoDebug::checkIntervals failed, descriptionString = '"
00139 << descriptionString << "'" << std::endl;
00140 }
00141 ASSERTDEBUG(!Interval<Treal>::
00142 intersect(XmX2EuclNorm, XmX2EuclNormExact).empty());
00143 }
00144
00145 template<typename Treal>
00146 template<typename Tmatrix>
00147 void PuriStepInfoDebug<Treal, DebugLevelHigh>::
00148 computeExactValues(Tmatrix const & X, Tmatrix const & X2,
00149 int const n, int const nocc) {
00150
00151 std::vector<Treal> full(n*n);
00152
00153 Treal* eigvals = new Treal[n];
00154 int lwork = 3*n-1;
00155 Treal* work = new Treal[lwork];
00156 int info = 0;
00157 Treal euclNormMatrix;
00158 Treal precSyev;
00159
00160 X.fullMatrix(full);
00161 syev("N", "U", &n, &full[0], &n, eigvals, work, &lwork, &info);
00162 ASSERTALWAYS(!info);
00163 euclNormMatrix = template_blas_fabs(eigvals[0]) > template_blas_fabs(eigvals[n-1]) ?
00164 template_blas_fabs(eigvals[0]) : template_blas_fabs(eigvals[n-1]);
00165 precSyev = euclNormMatrix * getRelPrecision<Treal>();
00166 lumoExact = Interval<Treal>(eigvals[n-nocc-1] - precSyev,
00167 eigvals[n-nocc-1] + precSyev);
00168 homoExact = Interval<Treal>(eigvals[n-nocc] - precSyev,
00169 eigvals[n-nocc] + precSyev);
00170 lminExact = Interval<Treal>(eigvals[0] - precSyev,
00171 eigvals[0] + precSyev);
00172 lmaxExact = Interval<Treal>(eigvals[n-1] - precSyev,
00173 eigvals[n-1] + precSyev);
00174
00175 Tmatrix XmX2(X);
00176 XmX2 += -1.0 * X2;
00177 XmX2.fullMatrix(full);
00178 syev("N", "U", &n, &full[0], &n, eigvals, work, &lwork, &info);
00179 ASSERTALWAYS(!info);
00180 euclNormMatrix = template_blas_fabs(eigvals[0]) > template_blas_fabs(eigvals[n-1]) ?
00181 template_blas_fabs(eigvals[0]) : template_blas_fabs(eigvals[n-1]);
00182 precSyev = euclNormMatrix * getRelPrecision<Treal>();
00183
00184 #if 0
00185 std::cout<<"EIGS XmX2: "<<std::endl;
00186 for(int ind = 0; ind < n; ++ind)
00187 std::cout<<std::setprecision(20)<<eigvals[ind]<<std::endl;
00188 std::cout<<"end EIGS XmX2: "<<std::endl<<std::endl;
00189 #endif
00190
00191 Treal lowBound = euclNormMatrix - precSyev > 0 ?
00192 euclNormMatrix - precSyev : 0;
00193 XmX2EuclNormExact = Interval<Treal>(lowBound, euclNormMatrix + precSyev);
00194
00195 delete[] eigvals;
00196 delete[] work;
00197 }
00198
00199
00200 }
00201 #undef __ID__
00202 #endif