@@ -173,6 +173,7 @@ 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
}
178
179
@@ -320,31 +321,40 @@ void QR_algorithm(std::vector<double> diag, std::vector<double> off_diag, cons
320
321
321
322
for (j=0 ; j<n;j=j+5 ){
322
323
tmp=Q[i][j];
323
- Q_posterior[i][j]=tmp*c-Q[i+1 ][j]*s;
324
- Q_posterior[i+1 ][j]=tmp*s+Q[i+1 ][j]*c;
324
+ Q[i][j]=tmp*c-Q[i+1 ][j]*s;
325
+ Q[i+1 ][j]=tmp*s+Q[i+1 ][j]*c;
326
+ // print_matrix(Q_posterior);
325
327
tmp=Q[i][j+1 ];
326
- Q_posterior[i][j+1 ]=tmp*c-Q[i+1 ][j+1 ]*s;
327
- Q_posterior[i+1 ][j+1 ]=tmp*s+Q[i+1 ][j+1 ]*c;
328
+ Q[i][j+1 ]=tmp*c-Q[i+1 ][j+1 ]*s;
329
+ Q[i+1 ][j+1 ]=tmp*s+Q[i+1 ][j+1 ]*c;
330
+ // print_matrix(Q_posterior);
328
331
tmp=Q[i][j+2 ];
329
- Q_posterior[i][j+2 ]=tmp*c-Q[i+1 ][j+2 ]*s;
330
- Q_posterior[i+1 ][j+2 ]=tmp*s+Q[i+1 ][j+2 ]*c;
332
+ Q[i][j+2 ]=tmp*c-Q[i+1 ][j+2 ]*s;
333
+ Q[i+1 ][j+2 ]=tmp*s+Q[i+1 ][j+2 ]*c;
334
+ // print_matrix(Q_posterior);
331
335
tmp=Q[i][j+3 ];
332
- Q_posterior[i][j+3 ]=tmp*c-Q[i+1 ][j+3 ]*s;
333
- Q_posterior[i+1 ][j+3 ]=tmp*s+Q[i+1 ][j+3 ]*c;
336
+ Q[i][j+3 ]=tmp*c-Q[i+1 ][j+3 ]*s;
337
+ Q[i+1 ][j+3 ]=tmp*s+Q[i+1 ][j+3 ]*c;
338
+ // print_matrix(Q_posterior);
334
339
tmp=Q[i][j+4 ];
335
- Q_posterior[i][j+4 ]=tmp*c-Q[i+1 ][j+4 ]*s;
336
- Q_posterior[i+1 ][j+4 ]=tmp*s+Q[i+1 ][j+4 ]*c;
340
+ Q[i][j+4 ]=tmp*c-Q[i+1 ][j+4 ]*s;
341
+ Q[i+1 ][j+4 ]=tmp*s+Q[i+1 ][j+4 ]*c;
342
+ // print_matrix(Q_posterior);
337
343
}
338
344
for (; j < n; j++)
339
345
{
340
346
tmp=Q[i][j];
341
- Q_posterior [i][j]=tmp*c-Q[i+1 ][j]*s;
342
- Q_posterior [i+1 ][j]=tmp*s+Q[i+1 ][j]*c;
347
+ Q [i][j]=tmp*c-Q[i+1 ][j]*s;
348
+ Q [i+1 ][j]=tmp*s+Q[i+1 ][j]*c;
343
349
}
344
350
351
+
345
352
}
346
353
347
- std::swap (Q, Q_posterior);
354
+
355
+
356
+ // std::swap(Q, Q_posterior);
357
+ // print_matrix(Q);
348
358
349
359
350
360
// for(unsigned int i=0;i<n_processor; i++){
@@ -406,12 +416,16 @@ void QR_algorithm(std::vector<double> diag, std::vector<double> off_diag, cons
406
416
}
407
417
}
408
418
409
- print_matrix (Q);
419
+ // print_matrix(Q);
410
420
411
- for (const auto t: diag){
412
- std::cout<<t<<" \t " ;
421
+ // for(const auto t: diag){
422
+ // std::cout<<t<<"\t";
423
+ // }
424
+ // std::cout<<"\n";
425
+ if (iter==max_iter){
426
+ std::cout<<" Converges failed" <<std::endl;
413
427
}
414
- std::cout<<" \n " ;
428
+ std::cout<<" Iteration: " <<iter<<std::endl ;
415
429
416
430
417
431
@@ -420,13 +434,13 @@ void QR_algorithm(std::vector<double> diag, std::vector<double> off_diag, cons
420
434
421
435
int main (){
422
436
423
- std::vector<double > diag (10 , 5 ), offdiag (9 , 20 );
437
+ std::vector<double > diag (5000 , 5 ), offdiag (4999 , 20 );
424
438
425
439
QR_algorithm (diag, offdiag);
426
440
427
441
auto start = std::chrono::high_resolution_clock::now ();
428
442
429
- QR_algorithm (diag, offdiag);
443
+ QR_algorithm (diag, offdiag, 1e-8 , 50000 );
430
444
431
445
// Capture the end time
432
446
auto end = std::chrono::high_resolution_clock::now ();
0 commit comments