38 #ifndef MAT_MATRIXTRIDIAGSYMMETRIC
39 #define MAT_MATRIXTRIDIAGSYMMETRIC
46 template<
typename Treal>
52 void increase(Treal
const & alpha, Treal
const & beta) {
76 Treal
const abstol = 0)
const;
82 Treal
const abstol = 0)
const;
95 template<
typename Treal>
101 Treal
const lowBound,
102 Treal
const uppBound,
103 Treal
const absTol)
const {
104 Treal* eigArray =
new Treal[size];
105 Treal* alphaCopy =
new Treal[size];
106 Treal* betaCopy =
new Treal[size];
107 Treal* work =
new Treal[5 * size];
108 int* iwork =
new int[5 * size];
109 int* ifail =
new int[size];
110 for (
int ind = 0; ind < size; ind++){
111 alphaCopy[ind] = alphaVec[ind];
112 betaCopy[ind] = betaVec[ind];
118 mat::stevx(
"V",
"V", &size, alphaCopy, betaCopy,
119 &lowBound, &uppBound, &dummy, &dummy,
121 &nEigsFound, eigArray, eigVectors, &size,
125 for (
int ind = 0; ind < nEigsFound; ind++) {
126 eigVals[ind] = eigArray[ind];
127 acc[ind] = betaCopy[size - 1] *
138 template<
typename Treal>
145 Treal
const abstol)
const {
146 Treal* eigArray =
new Treal[size];
147 Treal* alphaCopy =
new Treal[size];
148 Treal* betaCopy =
new Treal[size];
149 for (
int ind = 0; ind < size; ind++){
150 alphaCopy[ind] = alphaVec[ind];
151 betaCopy[ind] = betaVec[ind];
162 int const lowIndNew(lowInd + 1);
163 int const uppIndNew(uppInd + 1);
164 int nEigsWanted = uppInd - lowInd + 1;
166 int* isuppz =
new int[2*nEigsWanted];
176 mat::stevr(
"V",
"I", &size, alphaCopy, betaCopy,
177 &dummy, &dummy, &lowIndNew, &uppIndNew,
179 &nEigsFound, eigArray, eigVectors, &size,
181 &work_query, &lwork, &iwork_query, &liwork, &info);
182 lwork = int(work_query);
183 liwork = iwork_query;
184 work =
new Treal[lwork];
185 iwork =
new int[liwork];
186 mat::stevr(
"V",
"I", &size, alphaCopy, betaCopy,
187 &dummy, &dummy, &lowIndNew, &uppIndNew,
189 &nEigsFound, eigArray, eigVectors, &size,
191 work, &lwork, iwork, &liwork, &info);
193 std::cout <<
"info = " << info <<std::endl;
195 assert(nEigsFound == nEigsWanted);
196 for (
int ind = 0; ind < nEigsFound; ind++) {
197 eigVals[ind] = eigArray[ind];
198 acc[ind] = betaCopy[size - 1] *
209 Treal* work =
new Treal[5 * size];
210 int* iwork =
new int[5 * size];
211 int* ifail =
new int[size];
227 int const lowIndNew(lowInd + 1);
228 int const uppIndNew(uppInd + 1);
229 mat::stevx(
"V",
"I", &size, alphaCopy, betaCopy,
230 &dummy, &dummy, &lowIndNew, &uppIndNew,
232 &nEigsFound, eigArray, eigVectors, &size,
236 assert(nEigsFound == uppInd - lowInd + 1);
237 for (
int ind = 0; ind < nEigsFound; ind++) {
238 eigVals[ind] = eigArray[ind];
239 acc[ind] = betaCopy[size - 1] *
245 Treal* eigVectorsTmp =
new Treal[size*size];
246 int const lowIndNew(1);
247 int const uppIndNew(size);
248 mat::stevx(
"V",
"A", &size, alphaCopy, betaCopy,
249 &dummy, &dummy, &lowIndNew, &uppIndNew,
251 &nEigsFound, eigArray, eigVectorsTmp, &size,
255 assert(nEigsFound == size);
256 int nEigsWanted = uppInd - lowInd + 1;
258 for(
int i = 0; i < nEigsWanted; i++)
259 for(
int j = 0; j < size; j++)
260 eigVectors[i*size+j] = eigVectorsTmp[(lowInd+i)*size+j];
261 delete [] eigVectorsTmp;
262 for (
int ind = 0; ind < nEigsWanted; ind++) {
263 eigVals[ind] = eigArray[lowInd+ind];
264 acc[ind] = betaCopy[size - 1] *
279 template<
typename Treal>
282 capacity = newCapacity;
283 Treal* alphaNew =
new Treal[capacity];
284 Treal* betaNew =
new Treal[capacity];
285 for (
int ind = 0; ind < size; ind++){
286 alphaNew[ind] = alphaVec[ind];
287 betaNew[ind] = betaVec[ind];