36 #ifndef MAT_MATRIXTRIDIAGSYMMETRIC
37 #define MAT_MATRIXTRIDIAGSYMMETRIC
44 template<
typename Treal>
50 void increase(Treal
const & alpha, Treal
const & beta) {
67 Treal
const abstol = 0)
const;
73 Treal
const abstol = 0)
const;
86 template<
typename Treal>
94 Treal
const absTol)
const {
95 Treal* eigArray =
new Treal[size];
96 Treal* alphaCopy =
new Treal[size];
97 Treal* betaCopy =
new Treal[size];
98 Treal* work =
new Treal[5 * size];
99 int* iwork =
new int[5 * size];
100 int* ifail =
new int[size];
101 for (
int ind = 0; ind < size; ind++){
102 alphaCopy[ind] = alphaVec[ind];
103 betaCopy[ind] = betaVec[ind];
109 mat::stevx(
"V",
"V", &size, alphaCopy, betaCopy,
110 &lowBound, &uppBound, &dummy, &dummy,
112 &nEigsFound, eigArray, eigVectors, &size,
116 for (
int ind = 0; ind < nEigsFound; ind++) {
117 eigVals[ind] = eigArray[ind];
118 acc[ind] = betaCopy[size - 1] *
129 template<
typename Treal>
136 Treal
const abstol)
const {
137 Treal* eigArray =
new Treal[size];
138 Treal* alphaCopy =
new Treal[size];
139 Treal* betaCopy =
new Treal[size];
140 for (
int ind = 0; ind < size; ind++){
141 alphaCopy[ind] = alphaVec[ind];
142 betaCopy[ind] = betaVec[ind];
153 int const lowIndNew(lowInd + 1);
154 int const uppIndNew(uppInd + 1);
155 int nEigsWanted = uppInd - lowInd + 1;
157 int* isuppz =
new int[2*nEigsWanted];
167 mat::stevr(
"V",
"I", &size, alphaCopy, betaCopy,
168 &dummy, &dummy, &lowIndNew, &uppIndNew,
170 &nEigsFound, eigArray, eigVectors, &size,
172 &work_query, &lwork, &iwork_query, &liwork, &info);
173 lwork = int(work_query);
174 liwork = iwork_query;
175 work =
new Treal[lwork];
176 iwork =
new int[liwork];
177 mat::stevr(
"V",
"I", &size, alphaCopy, betaCopy,
178 &dummy, &dummy, &lowIndNew, &uppIndNew,
180 &nEigsFound, eigArray, eigVectors, &size,
182 work, &lwork, iwork, &liwork, &info);
184 std::cout <<
"info = " << info <<std::endl;
186 assert(nEigsFound == nEigsWanted);
187 for (
int ind = 0; ind < nEigsFound; ind++) {
188 eigVals[ind] = eigArray[ind];
189 acc[ind] = betaCopy[size - 1] *
200 Treal* work =
new Treal[5 * size];
201 int* iwork =
new int[5 * size];
202 int* ifail =
new int[size];
218 int const lowIndNew(lowInd + 1);
219 int const uppIndNew(uppInd + 1);
220 mat::stevx(
"V",
"I", &size, alphaCopy, betaCopy,
221 &dummy, &dummy, &lowIndNew, &uppIndNew,
223 &nEigsFound, eigArray, eigVectors, &size,
227 assert(nEigsFound == uppInd - lowInd + 1);
228 for (
int ind = 0; ind < nEigsFound; ind++) {
229 eigVals[ind] = eigArray[ind];
230 acc[ind] = betaCopy[size - 1] *
236 Treal* eigVectorsTmp =
new Treal[size*size];
237 int const lowIndNew(1);
238 int const uppIndNew(size);
239 mat::stevx(
"V",
"A", &size, alphaCopy, betaCopy,
240 &dummy, &dummy, &lowIndNew, &uppIndNew,
242 &nEigsFound, eigArray, eigVectorsTmp, &size,
246 assert(nEigsFound == size);
247 int nEigsWanted = uppInd - lowInd + 1;
249 for(
int i = 0; i < nEigsWanted; i++)
250 for(
int j = 0; j < size; j++)
251 eigVectors[i*size+j] = eigVectorsTmp[(lowInd+i)*size+j];
252 delete [] eigVectorsTmp;
253 for (
int ind = 0; ind < nEigsWanted; ind++) {
254 eigVals[ind] = eigArray[lowInd+ind];
255 acc[ind] = betaCopy[size - 1] *
270 template<
typename Treal>
273 capacity = newCapacity;
274 Treal* alphaNew =
new Treal[capacity];
275 Treal* betaNew =
new Treal[capacity];
276 for (
int ind = 0; ind < size; ind++){
277 alphaNew[ind] = alphaVec[ind];
278 betaNew[ind] = betaVec[ind];