MatrixTridiagSymmetric.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_MATRIXTRIDIAGSYMMETRIC
00037 #define MAT_MATRIXTRIDIAGSYMMETRIC
00038 #include "gblas.h"
00039 namespace mat { /* Matrix namespace */
00040   namespace arn { /* Arnoldi type methods namespace */
00044     template<typename Treal>
00045       class MatrixTridiagSymmetric {
00046     public:
00047       explicit MatrixTridiagSymmetric(int k = 100) 
00048         : alphaVec(new Treal[k]), betaVec(new Treal[k]),
00049         size(0), capacity(k) {}
00050         void increase(Treal const & alpha, Treal const & beta) {
00051           if (size + 1 > capacity)
00052             increaseCapacity(capacity * 2);
00053           ++size;
00054           alphaVec[size - 1] = alpha;
00055           betaVec[size - 1] = beta;
00056         }
00057         virtual ~MatrixTridiagSymmetric() {
00058           delete[] alphaVec;
00059           delete[] betaVec;
00060         }
00061         void getEigsByInterval(Treal* eigVals,
00062                                Treal* eigVectors, 
00063                                Treal* acc, 
00064                                int & nEigsFound, 
00065                                Treal const lowBound,
00066                                Treal const uppBound,
00067                                Treal const abstol = 0) const;
00068         void getEigsByIndex(Treal* eigVals, 
00069                             Treal* eigVectors, 
00070                             Treal* acc, 
00071                             int const lowInd,
00072                             int const uppInd,
00073                             Treal const abstol = 0) const;
00074         inline void clear() {
00075           size = 0;
00076         }
00077     protected:
00078         Treal* alphaVec;
00079         Treal* betaVec;
00080         int size;
00081         int capacity;
00082         void increaseCapacity(int const newCapacity);
00083     private:
00084     };
00085 
00086     template<typename Treal>
00087       void MatrixTridiagSymmetric<Treal>::
00088       getEigsByInterval(Treal* eigVals, /* length: >= nEigsFound */
00089                         Treal* eigVectors, /* length: >= size * nEigsFound */
00090                         Treal* acc, /* length: size */
00091                         int & nEigsFound, /* The number of found eigenpairs. */
00092                         Treal const lowBound,
00093                         Treal const uppBound,
00094                         Treal const absTol) const {
00095       Treal* eigArray  = new Treal[size];
00096       Treal* alphaCopy = new Treal[size];
00097       Treal* betaCopy  = new Treal[size];
00098       Treal* work = new Treal[5 * size];
00099       int* iwork = new int[5 * size];
00100       int* ifail = new int[size];
00101       for (int ind = 0; ind < size; ind++){
00102         alphaCopy[ind] = alphaVec[ind];
00103         betaCopy[ind] = betaVec[ind];
00104       }
00105       int dummy = -1;
00106       int info;
00107       /* Find eigenvalues */
00108       /* FIXME: change to stevr */
00109       mat::stevx("V", "V", &size, alphaCopy, betaCopy, 
00110                  &lowBound, &uppBound, &dummy, &dummy,
00111                  &absTol, 
00112                  &nEigsFound, eigArray, eigVectors, &size, 
00113                  work, iwork, ifail, 
00114                  &info);
00115       assert(info == 0);
00116       for (int ind = 0; ind < nEigsFound; ind++) {
00117         eigVals[ind] = eigArray[ind];
00118         acc[ind] = betaCopy[size - 1] * 
00119           template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9;
00120       }
00121       delete[] eigArray;
00122       delete[] alphaCopy;
00123       delete[] betaCopy;
00124       delete[] work;
00125       delete[] iwork;
00126       delete[] ifail;
00127     }
00128 
00129     template<typename Treal>
00130       void MatrixTridiagSymmetric<Treal>::
00131       getEigsByIndex(Treal* eigVals, /* length: uppInd-lowInd+1 */
00132                      Treal* eigVectors, /* length: size*(uppInd-lowInd+1) */
00133                      Treal* acc, /* length: size */
00134                      int const lowInd,
00135                      int const uppInd,
00136                      Treal const abstol) const {
00137       Treal* eigArray  = new Treal[size];
00138       Treal* alphaCopy = new Treal[size];
00139       Treal* betaCopy  = new Treal[size];
00140       for (int ind = 0; ind < size; ind++){
00141         alphaCopy[ind] = alphaVec[ind];
00142         betaCopy[ind] = betaVec[ind];
00143       }
00144 #if 1
00145       // Emanuel note 2010-03-14:
00146       // The following code uses stevr instead of stevx for two reasons: 
00147       // 1) Due to a bug in LAPACK we previously computed all
00148       // eigenvalues (see Elias' note below) which turned out to be
00149       // too time consuming in some cases.
00150       // 2) Contrary to stevx, stevr should never fail to compute the
00151       // desired eigenpairs unless there is a bug in the implementation
00152       // or erroneous input.
00153       int const lowIndNew(lowInd + 1);
00154       int const uppIndNew(uppInd + 1);
00155       int nEigsWanted = uppInd - lowInd + 1;
00156       int nEigsFound = 0;
00157       int* isuppz = new int[2*nEigsWanted];
00158       Treal* work;
00159       int lwork = -1;
00160       int* iwork;
00161       int liwork = -1;
00162       Treal dummy = -1.0;
00163       int info;
00164       // First do a workspace query:
00165       Treal work_query;
00166       int iwork_query;
00167       mat::stevr("V", "I", &size, alphaCopy, betaCopy, 
00168                  &dummy, &dummy, &lowIndNew, &uppIndNew,
00169                  &abstol,
00170                  &nEigsFound, eigArray, eigVectors, &size, 
00171                  isuppz, 
00172                  &work_query, &lwork, &iwork_query, &liwork, &info);
00173       lwork = int(work_query);
00174       liwork = iwork_query;
00175       work = new Treal[lwork];
00176       iwork = new int[liwork];
00177       mat::stevr("V", "I", &size, alphaCopy, betaCopy, 
00178                  &dummy, &dummy, &lowIndNew, &uppIndNew,
00179                  &abstol,
00180                  &nEigsFound, eigArray, eigVectors, &size, 
00181                  isuppz, 
00182                  work, &lwork, iwork, &liwork, &info);
00183       if (info)
00184         std::cout << "info = " << info <<std::endl;
00185       assert(info == 0);
00186       assert(nEigsFound == nEigsWanted);
00187       for (int ind = 0; ind < nEigsFound; ind++) {
00188         eigVals[ind] = eigArray[ind];
00189         acc[ind] = betaCopy[size - 1] * 
00190           template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9;
00191       }
00192       delete[] eigArray;
00193       delete[] alphaCopy;
00194       delete[] betaCopy;
00195       delete[] isuppz;
00196       delete[] work;
00197       delete[] iwork;
00198       
00199 #else
00200       Treal* work = new Treal[5 * size];
00201       int* iwork = new int[5 * size];
00202       int* ifail = new int[size];
00203       Treal dummy = -1.0;
00204       int info;
00205       int nEigsFound = 0;
00206       /* 
00207          Elias note 2007-07-02:
00208          There have been some problems with stevx returning with info=0
00209          at the same time as nEigsFound != uppInd - lowInd + 1.
00210          This is due to a bug in LAPACK which has been reported and confirmed.
00211          To avoid this problem Elias changed the code so that ALL eigenvectors
00212          are computed and then the desired ones are selected from the
00213          complete list. 
00214       */
00215 #if 1
00216       /* Original version of code calling stevx to get only the 
00217          desired eigenvalues/vectors. */
00218       int const lowIndNew(lowInd + 1);
00219       int const uppIndNew(uppInd + 1);
00220       mat::stevx("V", "I", &size, alphaCopy, betaCopy, 
00221                  &dummy, &dummy, &lowIndNew, &uppIndNew,
00222                  &abstol, 
00223                  &nEigsFound, eigArray, eigVectors, &size, 
00224                  work, iwork, ifail, 
00225                  &info);
00226       assert(info == 0);
00227       assert(nEigsFound == uppInd - lowInd + 1);
00228       for (int ind = 0; ind < nEigsFound; ind++) {
00229         eigVals[ind] = eigArray[ind];
00230         acc[ind] = betaCopy[size - 1] * 
00231           template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9;
00232       }
00233 #else
00234       /* Modified version of code calling stevx to get ALL 
00235          eigenvalues/vectors, and then picking out the desired ones. */
00236       Treal* eigVectorsTmp = new Treal[size*size];
00237       int const lowIndNew(1);
00238       int const uppIndNew(size);
00239       mat::stevx("V", "A", &size, alphaCopy, betaCopy, 
00240                  &dummy, &dummy, &lowIndNew, &uppIndNew, 
00241                  &abstol, 
00242                  &nEigsFound, eigArray, eigVectorsTmp, &size, 
00243                  work, iwork, ifail, 
00244                  &info);
00245       assert(info == 0);
00246       assert(nEigsFound == size);
00247       int nEigsWanted = uppInd - lowInd + 1;
00248       /* Copy desired eigenvectors from eigVectorsTmp to eigVectors. */
00249       for(int i = 0; i < nEigsWanted; i++)
00250         for(int j = 0; j < size; j++)
00251           eigVectors[i*size+j] = eigVectorsTmp[(lowInd+i)*size+j];
00252       delete [] eigVectorsTmp;
00253       for (int ind = 0; ind < nEigsWanted; ind++) {
00254         eigVals[ind] = eigArray[lowInd+ind];
00255         acc[ind] = betaCopy[size - 1] * 
00256           template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9;
00257       }
00258 #endif
00259       delete[] eigArray;
00260       delete[] alphaCopy;
00261       delete[] betaCopy;
00262       delete[] work;
00263       delete[] iwork;
00264       delete[] ifail;
00265 #endif
00266     }
00267 
00268 
00269 
00270     template<typename Treal>
00271       void MatrixTridiagSymmetric<Treal>::
00272       increaseCapacity(int const newCapacity) {
00273       capacity = newCapacity;
00274       Treal* alphaNew = new Treal[capacity];
00275       Treal* betaNew  = new Treal[capacity];
00276       for (int ind = 0; ind < size; ind++){
00277         alphaNew[ind] = alphaVec[ind];
00278         betaNew[ind] = betaVec[ind];
00279       }
00280       delete[] alphaVec;
00281       delete[] betaVec;
00282       alphaVec = alphaNew;
00283       betaVec = betaNew;
00284     }
00285 
00286 
00287   } /* end namespace arn */
00288 } /* end namespace mat */
00289 #endif

Generated on Mon Sep 17 14:32:56 2012 for ergo by  doxygen 1.4.7