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
00028
00029
00038 #ifndef MAT_MATRIXTRIDIAGSYMMETRIC
00039 #define MAT_MATRIXTRIDIAGSYMMETRIC
00040 #include "mat_gblas.h"
00041 namespace mat {
00042 namespace arn {
00046 template<typename Treal>
00047 class MatrixTridiagSymmetric {
00048 public:
00049 explicit MatrixTridiagSymmetric(int k = 100)
00050 : alphaVec(new Treal[k]), betaVec(new Treal[k]),
00051 size(0), capacity(k) {}
00052 void increase(Treal const & alpha, Treal const & beta) {
00053 if (size + 1 > capacity)
00054 increaseCapacity(capacity * 2);
00055 ++size;
00056 alphaVec[size - 1] = alpha;
00057 betaVec[size - 1] = beta;
00058 }
00059 virtual ~MatrixTridiagSymmetric() {
00060 delete[] alphaVec;
00061 delete[] betaVec;
00062 }
00063
00064
00065 void update_beta(Treal const & beta)
00066 {
00067 betaVec[size-1] = beta;
00068 }
00069
00070 void getEigsByInterval(Treal* eigVals,
00071 Treal* eigVectors,
00072 Treal* acc,
00073 int & nEigsFound,
00074 Treal const lowBound,
00075 Treal const uppBound,
00076 Treal const abstol = 0) const;
00077 void getEigsByIndex(Treal* eigVals,
00078 Treal* eigVectors,
00079 Treal* acc,
00080 int const lowInd,
00081 int const uppInd,
00082 Treal const abstol = 0) const;
00083 inline void clear() {
00084 size = 0;
00085 }
00086 protected:
00087 Treal* alphaVec;
00088 Treal* betaVec;
00089 int size;
00090 int capacity;
00091 void increaseCapacity(int const newCapacity);
00092 private:
00093 };
00094
00095 template<typename Treal>
00096 void MatrixTridiagSymmetric<Treal>::
00097 getEigsByInterval(Treal* eigVals,
00098 Treal* eigVectors,
00099 Treal* acc,
00100 int & nEigsFound,
00101 Treal const lowBound,
00102 Treal const uppBound,
00103 Treal const absTol) const {
00104 Treal* eigArray = new Treal[size];
00105 Treal* alphaCopy = new Treal[size];
00106 Treal* betaCopy = new Treal[size];
00107 Treal* work = new Treal[5 * size];
00108 int* iwork = new int[5 * size];
00109 int* ifail = new int[size];
00110 for (int ind = 0; ind < size; ind++){
00111 alphaCopy[ind] = alphaVec[ind];
00112 betaCopy[ind] = betaVec[ind];
00113 }
00114 int dummy = -1;
00115 int info;
00116
00117
00118 mat::stevx("V", "V", &size, alphaCopy, betaCopy,
00119 &lowBound, &uppBound, &dummy, &dummy,
00120 &absTol,
00121 &nEigsFound, eigArray, eigVectors, &size,
00122 work, iwork, ifail,
00123 &info);
00124 assert(info == 0);
00125 for (int ind = 0; ind < nEigsFound; ind++) {
00126 eigVals[ind] = eigArray[ind];
00127 acc[ind] = betaCopy[size - 1] *
00128 template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9;
00129 }
00130 delete[] eigArray;
00131 delete[] alphaCopy;
00132 delete[] betaCopy;
00133 delete[] work;
00134 delete[] iwork;
00135 delete[] ifail;
00136 }
00137
00138 template<typename Treal>
00139 void MatrixTridiagSymmetric<Treal>::
00140 getEigsByIndex(Treal* eigVals,
00141 Treal* eigVectors,
00142 Treal* acc,
00143 int const lowInd,
00144 int const uppInd,
00145 Treal const abstol) const {
00146 Treal* eigArray = new Treal[size];
00147 Treal* alphaCopy = new Treal[size];
00148 Treal* betaCopy = new Treal[size];
00149 for (int ind = 0; ind < size; ind++){
00150 alphaCopy[ind] = alphaVec[ind];
00151 betaCopy[ind] = betaVec[ind];
00152 }
00153 #if 1
00154
00155
00156
00157
00158
00159
00160
00161
00162 int const lowIndNew(lowInd + 1);
00163 int const uppIndNew(uppInd + 1);
00164 int nEigsWanted = uppInd - lowInd + 1;
00165 int nEigsFound = 0;
00166 int* isuppz = new int[2*nEigsWanted];
00167 Treal* work;
00168 int lwork = -1;
00169 int* iwork;
00170 int liwork = -1;
00171 Treal dummy = -1.0;
00172 int info;
00173
00174 Treal work_query;
00175 int iwork_query;
00176 mat::stevr("V", "I", &size, alphaCopy, betaCopy,
00177 &dummy, &dummy, &lowIndNew, &uppIndNew,
00178 &abstol,
00179 &nEigsFound, eigArray, eigVectors, &size,
00180 isuppz,
00181 &work_query, &lwork, &iwork_query, &liwork, &info);
00182 lwork = int(work_query);
00183 liwork = iwork_query;
00184 work = new Treal[lwork];
00185 iwork = new int[liwork];
00186 mat::stevr("V", "I", &size, alphaCopy, betaCopy,
00187 &dummy, &dummy, &lowIndNew, &uppIndNew,
00188 &abstol,
00189 &nEigsFound, eigArray, eigVectors, &size,
00190 isuppz,
00191 work, &lwork, iwork, &liwork, &info);
00192 if (info)
00193 std::cout << "info = " << info <<std::endl;
00194 assert(info == 0);
00195 assert(nEigsFound == nEigsWanted);
00196 for (int ind = 0; ind < nEigsFound; ind++) {
00197 eigVals[ind] = eigArray[ind];
00198 acc[ind] = betaCopy[size - 1] *
00199 template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9;
00200 }
00201 delete[] eigArray;
00202 delete[] alphaCopy;
00203 delete[] betaCopy;
00204 delete[] isuppz;
00205 delete[] work;
00206 delete[] iwork;
00207
00208 #else
00209 Treal* work = new Treal[5 * size];
00210 int* iwork = new int[5 * size];
00211 int* ifail = new int[size];
00212 Treal dummy = -1.0;
00213 int info;
00214 int nEigsFound = 0;
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224 #if 1
00225
00226
00227 int const lowIndNew(lowInd + 1);
00228 int const uppIndNew(uppInd + 1);
00229 mat::stevx("V", "I", &size, alphaCopy, betaCopy,
00230 &dummy, &dummy, &lowIndNew, &uppIndNew,
00231 &abstol,
00232 &nEigsFound, eigArray, eigVectors, &size,
00233 work, iwork, ifail,
00234 &info);
00235 assert(info == 0);
00236 assert(nEigsFound == uppInd - lowInd + 1);
00237 for (int ind = 0; ind < nEigsFound; ind++) {
00238 eigVals[ind] = eigArray[ind];
00239 acc[ind] = betaCopy[size - 1] *
00240 template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9;
00241 }
00242 #else
00243
00244
00245 Treal* eigVectorsTmp = new Treal[size*size];
00246 int const lowIndNew(1);
00247 int const uppIndNew(size);
00248 mat::stevx("V", "A", &size, alphaCopy, betaCopy,
00249 &dummy, &dummy, &lowIndNew, &uppIndNew,
00250 &abstol,
00251 &nEigsFound, eigArray, eigVectorsTmp, &size,
00252 work, iwork, ifail,
00253 &info);
00254 assert(info == 0);
00255 assert(nEigsFound == size);
00256 int nEigsWanted = uppInd - lowInd + 1;
00257
00258 for(int i = 0; i < nEigsWanted; i++)
00259 for(int j = 0; j < size; j++)
00260 eigVectors[i*size+j] = eigVectorsTmp[(lowInd+i)*size+j];
00261 delete [] eigVectorsTmp;
00262 for (int ind = 0; ind < nEigsWanted; ind++) {
00263 eigVals[ind] = eigArray[lowInd+ind];
00264 acc[ind] = betaCopy[size - 1] *
00265 template_blas_fabs(eigVectors[(ind * size) + size - 1]) / 0.9;
00266 }
00267 #endif
00268 delete[] eigArray;
00269 delete[] alphaCopy;
00270 delete[] betaCopy;
00271 delete[] work;
00272 delete[] iwork;
00273 delete[] ifail;
00274 #endif
00275 }
00276
00277
00278
00279 template<typename Treal>
00280 void MatrixTridiagSymmetric<Treal>::
00281 increaseCapacity(int const newCapacity) {
00282 capacity = newCapacity;
00283 Treal* alphaNew = new Treal[capacity];
00284 Treal* betaNew = new Treal[capacity];
00285 for (int ind = 0; ind < size; ind++){
00286 alphaNew[ind] = alphaVec[ind];
00287 betaNew[ind] = betaVec[ind];
00288 }
00289 delete[] alphaVec;
00290 delete[] betaVec;
00291 alphaVec = alphaNew;
00292 betaVec = betaNew;
00293 }
00294
00295
00296 }
00297 }
00298 #endif