Mentions légales du service

Skip to content
Snippets Groups Projects

cg_creux

Open MELLOUKI Bilal requested to merge patch-3 into master
1 file
+ 120
0
Compare changes
  • Side-by-side
  • Inline
creux_cg.cpp 0 → 100644
+ 120
0
#include <iostream>
#include <Eigen/Sparse>
using namespace Eigen;
using namespace std;
SparseMatrix<double> generer_sparse_sdp(int n){
SparseMatrix<double> mat(n, n);
SparseMatrix<double> sdp_sparse(n, n);
for (int i = 0; i < n-1; i++) {
mat.coeffRef(i, i) += (rand() % 10)/10. -.5;
mat.coeffRef(i, i+1) += (rand() % 10)/10.;
}
sdp_sparse = mat.transpose()*mat;
return sdp_sparse;
}
SparseVector<double> generer_sparse_vector(int n){
SparseVector<double> V(n);
for(int i = 0; i < n/2; i++){
V.coeffRef( rand()%n ) = 1;
}
return V;
}
SparseVector<double> cg(SparseMatrix<double> A, SparseVector<double> b, SparseVector<double> X, float eps, int Nmax){
SparseVector<double> q;
SparseVector<double> r = b - A*X;
SparseVector<double> p = r;
float alpha;
float beta;
float norm_r_old = r.dot(r);
float norm_r_new;
int n = 0;
for (int i = 0; i < Nmax; i++){
q = A * p;
alpha = norm_r_old/p.dot(q);
X += alpha * p;
r -= alpha * q ;
//cout << r.norm() << endl;
if (r.norm() < eps){
break ;
}
norm_r_new = r.dot(r);
beta = norm_r_new/norm_r_old;
norm_r_old = norm_r_new;
p = beta * p + r;
n++;
}
return X;
}
int main(){
int n = 4;
float eps = 0.01;
int Nmax = 10;
SparseMatrix<double> A(3,3);
A.coeffRef(0,0) = 1;
A.coeffRef(1,1) = 0.4;
A.coeffRef(2,2) = 0.25;
A.coeffRef(0,1) = A.coeffRef(0,1) = 0;
A.coeffRef(0,2) = A.coeffRef(2,0) = 0;
A.coeffRef(1,2) = A.coeffRef(2,1) = 0.1;
SparseVector<double> B(3);
B.coeffRef(0) = 0 ;
B.coeffRef(1) = 1 ;
B.coeffRef(2) = 5 ;
SparseVector<double> X(3);
X.coeffRef(0) = 0 ;
X.coeffRef(1) = 0 ;
X.coeffRef(2) = 0 ;
SparseVector<double> solution = cg(A,B,X,eps,Nmax);
//-------->Test_1<----------:
cout << "La matrice creuse est :" << endl << A << endl;
cout << "Le vecteur creux est :" << endl << B << endl;
cout << "La solution du systeme AX = B est :" << endl << solution << endl;
cout << "Verification:" << endl << A * solution << endl;
//-------->Test_avec_generer<----------:
SparseMatrix<double> A2 = generer_sparse_sdp(3);
cout << "La matrice creuse est :" << endl << A2 << endl;
cout << "Le vecteur creux est :" << endl << B << endl;
cout << "La solution du systeme AX = B est :" << endl << solution << endl;
cout << "Verification:" << endl << A * solution << endl;
}
Loading