ergo
|
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