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_MATRIXTRIDIAGSYMMETRIC
00037 #define MAT_MATRIXTRIDIAGSYMMETRIC
00038 #include "gblas.h"
00039 namespace mat {
00040 namespace arn {
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,
00089 Treal* eigVectors,
00090 Treal* acc,
00091 int & nEigsFound,
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
00108
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,
00132 Treal* eigVectors,
00133 Treal* acc,
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
00146
00147
00148
00149
00150
00151
00152
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
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
00208
00209
00210
00211
00212
00213
00214
00215 #if 1
00216
00217
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
00235
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
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 }
00288 }
00289 #endif