+// /home/couturie/work/petsc-3.5.1_old/arch-linux2-c-debug/bin/mpirun -np 4 -machinefile archi ./ex15 -m 1000 -n 1000 -ksp_type fgmres -pc_type mg
-// /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 4 ./ex15 -m 400 -n 400
+
+// /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 4 ./ex15 -m 400 -n 400 -ksp_type fgmres
//Variables
PetscScalar gamma, alpha, oldgamma, beta;
- PetscReal norm=20, Eprecision=1e-8, cgprec=1e-40;
+ PetscReal norm=20, Eprecision=1e-3, cgprec=1e-40;
PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
PetscErrorCode ierr;
PetscScalar T1, T2;
//Initializations
// ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
- ierr = KSPSetTolerances(ksp, 1e-10, 1e-10, PETSC_DEFAULT, 16); CHKERRQ(ierr);
+ ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
//GMRES WITH MINIMIZATION
T1 = MPI_Wtime();
+ ierr = KSPSetUp(ksp); CHKERRQ(ierr);
while(giter<Emaxiter && norm>Eprecision ){
for(col=0; col<ColS && norm>Eprecision; col++){
- //Error
+ KSPGetResidualNorm(ksp,&norm);
+
+ /* //Error
ierr = VecCopy(x, residu); CHKERRQ(ierr);
ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
-
+ */
ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
//Variables
PetscScalar alpha, beta;
- PetscReal norm=20, Eprecision=1e-8, tol=1e-40;
- PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
+ PetscReal norm=20, Eprecision=1e-3, tol=1e-40;
+ PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=20, Iiter=0;
PetscErrorCode ierr;
PetscScalar T1, T2;
KSP ksp;
//Initializations
// ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
- ierr = KSPSetTolerances(ksp, 1e-10, 1e-10, PETSC_DEFAULT, 16); CHKERRQ(ierr);
+ ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
//GMRES WITH MINIMIZATION
T1 = MPI_Wtime();
+ ierr = KSPSetUp(ksp); CHKERRQ(ierr);
while(giter<Emaxiter && norm>Eprecision ){
for(col=0; col<ColS && norm>Eprecision; col++){
+ KSPGetResidualNorm(ksp,&norm);
+
+ /*
//Error
ierr = VecCopy(x, residu); CHKERRQ(ierr);
ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
-
+ */
ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
VecAYPX(d,-thet,v); //d = (v - thet * d);
VecAYPX(d,1/rho,zero_short); //d=d/ rho;
-
+ /*
if (normar/(norma*normr) <= tol) { // check for convergence in min{|b-A*x|}
break;
}
if (normr <= tolb) { // check for convergence in A*x=b
break;
}
-
+ */
VecAXPY(x_lsqr,phi,d); // x_lsqr=x_lsqr+phi*d
normr = abs(s) * normr;
ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
+ PetscMPIInt size;
+ MPI_Comm_size(PETSC_COMM_WORLD,&size);
+ PetscPrintf(PETSC_COMM_WORLD,"Number of processors = %d\n",size);
+
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Compute the matrix and right-hand-side vector that define
the linear system, Ax = b.
*/
for (Ii=Istart; Ii<Iend; Ii++) {
v = -1.0; i = Ii/n; j = Ii - i*n;
- if (i>0) {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
+ PetscScalar v2=-1.;
+ if (i>0) {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v2,INSERT_VALUES);CHKERRQ(ierr);}
if (i<m-1) {J = Ii + n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
if (j>0) {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
if (j<n-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
to set various options.
*/
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
- ierr = KSPSetTolerances(ksp,1e-9,1e-9,PETSC_DEFAULT,5000000);CHKERRQ(ierr);
+ ierr = KSPSetTolerances(ksp,4e-6,4e-6,PETSC_DEFAULT,500000);CHKERRQ(ierr);
+
+
+ // PC pc;
+ ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+ KSPGetPC(ksp, &pc);
+ PCType type;
+ PCGetType(pc, &type);
+
+ PetscPrintf(PETSC_COMM_WORLD, "PC TYPE %s \n", type);
+
/*
Set runtime options, e.g.,
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve the linear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
- PetscScalar T1,T2;
- T1 = MPI_Wtime();
- ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
- T2 = MPI_Wtime();
- /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- Check solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
- /*
- Check the error
- */
- Vec sol;
- VecDuplicate(b,&sol);
- MatMult(A,x,sol);
- VecAXPY(sol,-1,b);
- VecNorm(sol, NORM_2, &norm);
- ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
- ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr);
- ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n\n\n", T2-T1); CHKERRQ(ierr);
-
-
+ /* PetscScalar T1,T2; */
+ /* T1 = MPI_Wtime(); */
+ /* ierr = KSPSetUp(ksp); CHKERRQ(ierr); */
+ /* ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); */
+ /* T2 = MPI_Wtime(); */
+ /* /\* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+ /* Check solution and clean up */
+ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *\/ */
+
+ /* /\* */
+ /* Check the error */
+ /* *\/ */
+ /* Vec sol; */
+ /* VecDuplicate(b,&sol); */
+ /* MatMult(A,x,sol); */
+ /* VecAXPY(sol,-1,b); */
+ /* VecNorm(sol, NORM_2, &norm); */
+ /* ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); */
+ /* ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr); */
+ /* ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n\n\n", T2-T1); CHKERRQ(ierr); */
+
+
+
+
+ //version to control the error
+ {
+ Vec x2;
+ Vec sol;
+ VecDuplicate(b,&x2);
+ VecDuplicate(b,&sol);
+
+ PetscScalar norm=100;
+ PetscScalar T1,T2;
+ PetscInt total,its;
+ ierr = KSPSetTolerances(ksp,1e-10,1e-10,PETSC_DEFAULT,30);CHKERRQ(ierr);
+ ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
+ T1 = MPI_Wtime();
+ while(norm>1e-3) {
+ ierr = KSPSolve(ksp,b,x2);CHKERRQ(ierr);
+ KSPGetResidualNorm(ksp,&norm);
+ ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
+ total += its;
+ ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g\n", norm); CHKERRQ(ierr);
+ }
+ T2 = MPI_Wtime();
+ MatMult(A,x2,sol);
+ VecAXPY(sol,-1,b);
+ VecNorm(sol, NORM_2, &norm);
+ ierr = PetscPrintf(PETSC_COMM_WORLD,"Computed norm of error %g iterations %D\n",(double)norm,total);CHKERRQ(ierr);
+ ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time NORMAL GMRES : %g (s)\n\n\n", T2-T1); CHKERRQ(ierr);
+ }
+
{
Vec x2;
VecDuplicate(b,&x2);
VecDuplicate(b,&sol);
- KrylovMinimize(A, b, x2);
+ KrylovMinimizeLSQR(A, b, x2);
ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Error Krylov Minimization LSQR %g\n",norm);
}
-
-
-
+
+
/*
Free work space. All PETSc objects should be destroyed when they
are no longer needed.
