PuriStepInfoDebug.h

Go to the documentation of this file.
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure
00002  * calculations.
00003  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
00004  * 
00005  * This program is free software: you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation, either version 3 of the License, or
00008  * (at your option) any later version.
00009  * 
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  * 
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00017  * 
00018  * Primary academic reference:
00019  * Kohn−Sham Density Functional Theory Electronic Structure Calculations 
00020  * with Linearly Scaling Computational Time and Memory Usage,
00021  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
00022  * J. Chem. Theory Comput. 7, 340 (2011),
00023  * <http://dx.doi.org/10.1021/ct100611z>
00024  * 
00025  * For further information about Ergo, see <http://www.ergoscf.org>.
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 /* Specialization for high debug level */
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     /*  Failure here probably means that checkIntervals was called before
00095      *  the exact values were calculated.
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     /* The homo and lumo intervals only provides information if 
00114      * the exact values lie in [0, 1] otherwise no check is done.
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     /* Set up work space */
00151     std::vector<Treal> full(n*n);
00152     //    Treal* full = new Treal[n*n];
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     /* Compute lumoExact, homoExact, lminExact, and lmaxExact.  */
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     /* Compute largest magnitude eigenvalue of X-X^2 */
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     /* Free memory */
00195     delete[] eigvals;
00196     delete[] work;
00197   }
00198 
00199 
00200 } /* end namespace mat */
00201 #undef __ID__
00202 #endif

Generated on Wed Nov 21 09:32:11 2012 for ergo by  doxygen 1.4.7