@@ -173,18 +173,19 @@ void print_matrix(const std::vector<std::vector<double> > & Q){
173
173
}
174
174
std::cout<<" \n " ;
175
175
}
176
+ std::cout<<" --------------------------------------------------------" <<" \n " ;
176
177
177
178
std::cout<<" ------------------------------------------------" <<" \n " ;
178
179
179
180
}
180
181
181
182
182
183
// std::pair<std::vector<double>, std::vector<std::vector<double>> >
183
- void QR_algorithm (std::vector<double > diag, std::vector<double > off_diag, const double toll=1e-8 , const unsigned int max_iter=1000 ){
184
+ void QR_algorithm (std::vector<double > diag, std::vector<double > off_diag, const double toll=1e-8 , const unsigned int max_iter=10000 ){
184
185
185
186
const unsigned int n = diag.size ();
186
187
187
- std::vector<std::vector<double >> Q (n, std::vector<double >(n, 0 )), Q_posterior (n, std::vector< double >(n, 0 )), R (n);
188
+ std::vector<std::vector<double >> Q (n, std::vector<double >(n, 0 )), R (n);
188
189
189
190
for (unsigned int i = 0 ; i < n; i++){
190
191
Q[i][i] = 1 ;
@@ -315,15 +316,16 @@ void QR_algorithm(std::vector<double> diag, std::vector<double> off_diag, cons
315
316
316
317
unsigned j;
317
318
// Uncomment to compute the eigenvalue
318
- # pragma omp parallel for // collapse(2)
319
+ // collapse(2)
319
320
for (unsigned int i=0 ; i<n-1 ; i++){
321
+ // std::vector<std::vector<double>> Q_posterior(Q);
320
322
c=Matrix_trigonometric[i][0 ];
321
323
s=Matrix_trigonometric[i][1 ];
322
-
324
+ # pragma omp parallel for
323
325
for (j=0 ; j<n;j=j+5 ){
324
326
tmp=Q[i][j];
325
327
Q[i][j]=tmp*c-Q[i+1 ][j]*s;
326
- Q_posterior [i+1 ][j]=tmp*s+Q[i+1 ][j]*c;
328
+ Q [i+1 ][j]=tmp*s+Q[i+1 ][j]*c;
327
329
tmp=Q[i][j+1 ];
328
330
Q[i][j+1 ]=tmp*c-Q[i+1 ][j+1 ]*s;
329
331
Q[i+1 ][j+1 ]=tmp*s+Q[i+1 ][j+1 ]*c;
@@ -332,7 +334,7 @@ void QR_algorithm(std::vector<double> diag, std::vector<double> off_diag, cons
332
334
Q[i+1 ][j+2 ]=tmp*s+Q[i+1 ][j+2 ]*c;
333
335
tmp=Q[i][j+3 ];
334
336
Q[i][j+3 ]=tmp*c-Q[i+1 ][j+3 ]*s;
335
- Q[i+1 ][j+3 ]=tmp*s+Q[i+1 ][j+3 ]*c;
337
+ Q[i+1 ][j+3 ]=tmp*s+Q[i+1 ][j+3 ]*c;
336
338
tmp=Q[i][j+4 ];
337
339
Q[i][j+4 ]=tmp*c-Q[i+1 ][j+4 ]*s;
338
340
Q[i+1 ][j+4 ]=tmp*s+Q[i+1 ][j+4 ]*c;
@@ -343,6 +345,7 @@ void QR_algorithm(std::vector<double> diag, std::vector<double> off_diag, cons
343
345
Q[i][j]=tmp*c-Q[i+1 ][j]*s;
344
346
Q[i+1 ][j]=tmp*s+Q[i+1 ][j]*c;
345
347
}
348
+
346
349
347
350
}
348
351
@@ -353,12 +356,10 @@ void QR_algorithm(std::vector<double> diag, std::vector<double> off_diag, cons
353
356
}
354
357
}
355
358
356
- print_matrix (Q);
357
-
358
- for (const auto t: diag){
359
- std::cout<<t<<" \t " ;
359
+ if (iter==max_iter){
360
+ std::cout<<" Converges failed" <<std::endl;
360
361
}
361
- std::cout<<" \n " ;
362
+ std::cout<<" Iteration: " <<iter<<std::endl ;
362
363
363
364
364
365
@@ -368,13 +369,12 @@ void QR_algorithm(std::vector<double> diag, std::vector<double> off_diag, cons
368
369
369
370
int main (){
370
371
371
- std::vector<double > diag (10 , 5 ), offdiag (9 , 20 );
372
+ std::vector<double > diag (2000 , 5 ), offdiag (1999 , 20 );
372
373
373
- QR_algorithm (diag, offdiag);
374
374
375
375
auto start = std::chrono::high_resolution_clock::now ();
376
376
377
- QR_algorithm (diag, offdiag);
377
+ QR_algorithm (diag, offdiag, 1e-8 , 50000 );
378
378
379
379
// Capture the end time
380
380
auto end = std::chrono::high_resolution_clock::now ();
0 commit comments