+
+///home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/crashbasis/crashbasis.bin -ksp_type gmres -pc_type none
+//16s
+//temps simi pour lui et nous
+
+
+// /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/parabolic_fem/parabolic_fem.bin -ksp_type gmres -pc_type ilu
+//converge avec nous mais lentement
+//2152s pour lui
+//724s pour nous
+
+
+
+// /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/epb3/epb3.bin -ksp_type fgmres -pc_type sor
+//fonctionne avec ilu asm
+//9s
+//temps simi pour lui et nous
+
+// /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/atmosmodj/atmosmodj.bin -ksp_type fgmres -pc_type none
+//aussi sor
+
+// /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/bfwa398/bfwa398.bin -ksp_type gmres -pc_type none
+
+
+// /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/torso3/torso3.bin -ksp_type gmres -pc_type none
+
+
+static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
+This version first preloads and solves a small system, then loads \n\
+another (larger) system and solves it as well. This example illustrates\n\
+preloading of instructions with the smaller system so that more accurate\n\
+performance monitoring can be done with the larger one (that actually\n\
+is the system of interest). See the 'Performance Hints' chapter of the\n\
+users manual for a discussion of preloading. Input parameters include\n\
+ -f0 <input_file> : first file to load (small system)\n\
+ -f1 <input_file> : second file to load (larger system)\n\n\
+ -nearnulldim <0> : number of vectors in the near-null space immediately following matrix\n\n\
+ -trans : solve transpose system instead\n\n";
+/*
+ This code can be used to test PETSc interface to other packages.\n\
+ Examples of command line options: \n\
+ ./ex10 -f0 <datafile> -ksp_type preonly \n\
+ -help -ksp_view \n\
+ -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
+ -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package superlu or superlu_dist or mumps \n\
+ -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_package mumps \n\
+ mpiexec -n <np> ./ex10 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij
+ \n\n";
+*/
+/*T
+ Concepts: KSP^solving a linear system
+ Processors: n
+T*/
+
+/*
+ Include "petscksp.h" so that we can use KSP solvers. Note that this file
+ automatically includes:
+ petscsys.h - base PETSc routines petscvec.h - vectors
+ petscmat.h - matrices
+ petscis.h - index sets petscksp.h - Krylov subspace methods
+ petscviewer.h - viewers petscpc.h - preconditioners
+*/
+#include <petscksp.h>
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+
+
+
+
+
+
+
+
+
+int KrylovMinimize(Mat A, Vec b, Vec x) {
+
+
+ //Variables
+
+ PetscScalar gamma, alpha, oldgamma, beta;
+ PetscReal norm=20, Eprecision=1e-7, cgprec=1e-40;
+ PetscInt giter=0, ColS=8, col=0, Emaxiter=50000000, iter=0, iterations=20, Iiter=0;
+ PetscErrorCode ierr;
+ PetscScalar T1, T2;
+ KSP ksp;
+ PetscInt total=0;
+ PetscInt size;
+ PetscInt Istart,Iend;
+ PetscInt i,its;
+ Vec x_old, residu;
+ Mat S, AS;
+ PetscScalar *array;
+ PetscInt *ind_row;
+ Vec Alpha, p, ss, vect, r, q, Ax;
+
+
+ PetscInt first=1;
+
+ ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
+ ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
+
+
+
+ ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+
+
+
+
+ VecGetSize(b,&size);
+
+ ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);
+
+ PetscInt aa,bb;
+ MatGetOwnershipRange(A,&aa,&bb);
+
+ // ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%D %D\n", aa,bb); CHKERRQ(ierr);
+ //PetscSynchronizedFlush(PETSC_COMM_WORLD);
+
+
+ ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
+ ierr = MatSetSizes(S, bb-aa, PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
+ ierr = MatSetType(S, MATDENSE); CHKERRQ(ierr);
+ ierr = MatSetUp(S); CHKERRQ(ierr);
+
+ ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
+
+ ierr = VecCreate(PETSC_COMM_WORLD, &Alpha); CHKERRQ(ierr);
+ ierr = VecSetSizes(Alpha, PETSC_DECIDE, ColS); CHKERRQ(ierr);
+ ierr = VecSetFromOptions(Alpha); CHKERRQ(ierr);
+ ierr = VecDuplicate(Alpha, &vect); CHKERRQ(ierr);
+ ierr = VecDuplicate(Alpha, &p); CHKERRQ(ierr);
+ ierr = VecDuplicate(Alpha, &ss); CHKERRQ(ierr);
+ ierr = VecDuplicate(b, &r); CHKERRQ(ierr);
+ ierr = VecDuplicate(b, &q); CHKERRQ(ierr);
+ ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
+
+ ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
+ ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
+
+
+
+ //indexes of row (these indexes are global)
+ ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
+ for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
+
+ //Initializations
+ // ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
+ ierr = KSPSetTolerances(ksp, 1e-10, 1e-10, 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++){
+
+
+ //Solve
+ ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
+
+ ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
+ total += its; Iiter ++;
+
+
+
+ //Build S'
+ ierr = VecGetArray(x, &array);
+ ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
+ VecRestoreArray(x, &array);
+
+
+
+ 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);
+ //ierr = VecCopy(x, x_old); CHKERRQ(ierr);
+
+
+ }
+
+
+ //minimization step
+ if( norm>Eprecision) {
+
+ MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
+
+
+ //Build AS
+ if(first) {
+ MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
+ first=0;
+ }
+ else
+ MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
+
+
+
+
+ //Minimization with CGLS
+ MatMult(AS, Alpha, r); VecAYPX(r, -1, b); //r_0 = b-AS*x_0
+
+
+ MatMultTranspose(AS, r, p); //p_0 = AS'*r_0
+
+
+
+
+ ierr = VecCopy(p, ss); CHKERRQ(ierr); //p_0 = s_0
+ VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
+ while(gamma>cgprec && iter<iterations){
+ MatMult(AS, p, q); //q = AS*p
+ VecNorm(q, NORM_2, &alpha); alpha = alpha * alpha; //alpha = norm2(q)^2
+ alpha = gamma / alpha;
+ VecAXPY(Alpha, alpha, p); //Alpha += alpha*p
+ VecAXPY(r, -alpha, q); //r -= alpha*q
+ MatMultTranspose(AS, r, ss); // ss = AS'*r
+ oldgamma = gamma;
+ VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
+ beta = gamma / oldgamma;
+ VecAYPX(p, beta, ss); //p = s+beta*p;
+ iter ++;
+ }
+
+ iter = 0; giter ++;
+ //Minimizer
+ MatMult(S, Alpha, x); //x = S*Alpha;
+ }
+
+ }
+ T2 = MPI_Wtime();
+ ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
+ ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
+
+ return 0;
+
+}
+
+
+
+
+
+
+
+
+
+int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
+
+
+ //Variables
+
+ PetscScalar alpha, beta;
+ PetscReal norm=20, Eprecision=1e-7, tol=1e-40;
+ PetscInt giter=0, ColS=8, col=0, Emaxiter=50000000, iter=0, iterations=20, Iiter=0;
+ PetscErrorCode ierr;
+ PetscScalar T1, T2;
+ KSP ksp;
+ PetscInt total=0;
+ PetscInt size;
+ PetscInt Istart,Iend;
+ PetscInt i,its;
+ Vec x_old, residu;
+ Mat S, AS;
+ PetscScalar *array;
+ PetscInt *ind_row;
+ Vec Ax;
+ PetscScalar norm_ksp;
+ Vec u,v,d,uu,vt,zero_long,zero_short,x_lsqr;
+
+ PetscInt first=1;
+
+ ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
+ ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
+ ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+
+
+
+
+ VecGetSize(b,&size);
+
+ ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);
+
+ PetscInt aa,bb;
+ MatGetOwnershipRange(A,&aa,&bb);
+
+ // ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%D %D\n", aa,bb); CHKERRQ(ierr);
+ //PetscSynchronizedFlush(PETSC_COMM_WORLD);
+
+
+ ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
+ ierr = MatSetSizes(S, bb-aa, PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
+ ierr = MatSetType(S, MATDENSE); CHKERRQ(ierr);
+ ierr = MatSetUp(S); CHKERRQ(ierr);
+
+ ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
+
+
+
+ ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
+
+ ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
+ ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
+
+
+ //long vector
+ ierr = VecDuplicate(b,&u);CHKERRQ(ierr);
+
+
+ ierr = VecDuplicate(b,&uu);CHKERRQ(ierr);
+ ierr = VecDuplicate(b,&zero_long);CHKERRQ(ierr);
+ ierr = VecSet(zero_long,0);CHKERRQ(ierr);
+
+ //small vector
+ ierr = VecCreate(PETSC_COMM_WORLD, &v); CHKERRQ(ierr);
+ ierr = VecSetSizes(v, PETSC_DECIDE, ColS); CHKERRQ(ierr);
+ ierr = VecSetFromOptions(v); CHKERRQ(ierr);
+ ierr = VecDuplicate(v,&zero_short);CHKERRQ(ierr);
+ ierr = VecSet(zero_short,0);CHKERRQ(ierr);
+ ierr = VecDuplicate(v,&d);CHKERRQ(ierr);
+ ierr = VecDuplicate(v,&vt);CHKERRQ(ierr);
+ ierr = VecDuplicate(v,&x_lsqr);CHKERRQ(ierr);
+
+
+ //indexes of row (these indexes are global)
+ ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
+ for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
+
+ //Initializations
+ // ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
+ ierr = KSPSetTolerances(ksp, 1e-10, 1e-10, 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++){
+
+
+ //Solve
+ ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
+
+ ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
+ total += its; Iiter ++;
+
+
+
+ //Build S'
+ ierr = VecGetArray(x, &array);
+ ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
+ VecRestoreArray(x, &array);
+
+
+
+ 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);
+ //ierr = VecCopy(x, x_old); CHKERRQ(ierr);
+
+
+ }
+
+
+ //minimization step
+ if( norm>Eprecision) {
+
+ MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
+
+
+
+
+ //Build AS
+ if(first) {
+ MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
+
+ first=0;
+ }
+ else
+ MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
+
+
+
+
+
+ //LSQR
+ //LSQR
+ //LSQR
+
+
+
+ PetscScalar n2b,tolb,normr,c,s,phibar,normar,norma,thet,rhot,rho,phi;
+ PetscInt stag;
+ tolb = tol * n2b;
+ VecNorm(b, NORM_2, &n2b); //n2b = norm(b);
+ ierr = VecCopy(b, u); CHKERRQ(ierr); //u=b
+ VecNorm(u, NORM_2, &beta); // beta=norm(u)
+ normr=beta;
+ if (beta != 0) {
+ VecAYPX(u,1/beta,zero_long); // u = u / beta;
+ }
+ c=1;
+ s=0;
+ phibar=beta;
+ MatMultTranspose(AS, u, v); //v=A'*u
+ ierr = VecSet(x_lsqr,0);CHKERRQ(ierr);
+ VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
+ if (alpha != 0) {
+ VecAYPX(v,1/alpha,zero_short); // v = v / alpha;
+ }
+ ierr = VecSet(d,0);CHKERRQ(ierr);
+ normar = alpha * beta;
+ norma=0;
+ //stag=0;
+ for(i=0;i<iterations;i++) {
+ MatMult(AS, v, uu); //uu=A*v
+ VecAYPX(u, -alpha, uu); //u = uu-alpha*u;
+ VecNorm(u, NORM_2, &beta); // beta=norm(u)
+ VecAYPX(u,1/beta,zero_long); // u = u / beta;
+ norma=sqrt(norma*norma+alpha*alpha+beta*beta); // norma = norm([norma alpha beta]);
+ thet = - s * alpha;
+ rhot = c * alpha;
+ rho = sqrt(rhot*rhot + beta*beta);
+ c = rhot / rho;
+ s = - beta / rho;
+ phi = c * phibar;
+ if (phi == 0) { // stagnation of the method
+ stag = 1;
+ }
+ phibar = s * phibar;
+ 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;
+ MatMultTranspose(AS, u, vt); //vt=A'*u;
+ VecAYPX(v,-beta,vt); // v = vt - beta * v;
+ VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
+ VecAYPX(v,1/alpha,zero_short); // v = v / alpha;
+ normar = alpha * abs( s * phi);
+ }
+
+
+
+ iter = 0; giter ++;
+ //Minimizer
+ MatMult(S, x_lsqr, x); //x = S*x_small;
+ }
+
+ }
+ T2 = MPI_Wtime();
+ ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time LSQR : %g (s)\n", T2-T1); CHKERRQ(ierr);
+ ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations LSQR : %D\n", total); CHKERRQ(ierr);
+
+ return 0;
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+int main(int argc,char **args)
+{
+ KSP ksp; /* linear solver context */
+ Mat A,A2; /* matrix */
+ Vec x,b,u; /* approx solution, RHS, exact solution */
+ PetscViewer fd; /* viewer */
+ char file[4][PETSC_MAX_PATH_LEN]; /* input file name */
+ PetscBool table =PETSC_FALSE,flg,trans=PETSC_FALSE,initialguess = PETSC_FALSE;
+ PetscBool outputSoln=PETSC_FALSE;
+ PetscErrorCode ierr;
+ PetscInt its,num_numfac,m,n,M,nearnulldim = 0;
+ PetscReal norm;
+ PetscBool preload=PETSC_TRUE,isSymmetric,cknorm=PETSC_FALSE,initialguessfile = PETSC_FALSE;
+ PetscMPIInt rank;
+ char initialguessfilename[PETSC_MAX_PATH_LEN];
+
+ PetscInitialize(&argc,&args,(char*)0,help);
+ ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
+
+ ierr = PetscOptionsGetString(NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
+ if (flg) {
+ ierr = PetscStrcpy(file[1],file[0]);CHKERRQ(ierr);
+ preload = PETSC_FALSE;
+ } else {
+ ierr = PetscOptionsGetString(NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
+ if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 or -f option");
+ ierr = PetscOptionsGetString(NULL,"-f1",file[1],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
+ if (!flg) preload = PETSC_FALSE; /* don't bother with second system */
+ }
+
+
+ PetscPreLoadBegin(preload,"Load system");
+
+
+
+ ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);CHKERRQ(ierr);
+
+ ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
+ ierr = MatSetFromOptions(A);CHKERRQ(ierr);
+ ierr = MatLoad(A,fd);CHKERRQ(ierr);
+
+ /* char ordering[256] = MATORDERINGRCM;
+ IS rowperm = NULL,colperm = NULL;
+ ierr = MatGetOrdering(A2,ordering,&rowperm,&colperm);CHKERRQ(ierr);
+ ierr = MatPermute(A2,rowperm,colperm,&A);CHKERRQ(ierr);
+ */
+
+ ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
+
+ ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
+
+ ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr);
+ ierr = VecSetSizes(b,m,PETSC_DECIDE);CHKERRQ(ierr);
+ ierr = VecSetFromOptions(b);CHKERRQ(ierr);
+ ierr = VecSet(b,1);CHKERRQ(ierr);
+
+
+
+
+ ierr = MatGetVecs(A,&x,NULL);CHKERRQ(ierr);
+ ierr = VecDuplicate(b,&u);CHKERRQ(ierr);
+
+ ierr = VecSet(x,1.0);CHKERRQ(ierr);
+
+
+ ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
+ ierr = KSPSetInitialGuessNonzero(ksp,initialguess);CHKERRQ(ierr);
+ num_numfac = 1;
+
+
+
+
+
+ ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
+
+ ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+
+ ierr = KSPSetTolerances(ksp,1e-10,1e-10,PETSC_DEFAULT,5000000);CHKERRQ(ierr);
+ PC pc;
+ KSPGetPC(ksp, &pc);
+ PCType type;
+ PCGetType(pc, &type);
+
+ PetscPrintf(PETSC_COMM_WORLD, "PC TYPE %s \n", type);
+
+
+ PetscScalar T1,T2;
+ T1 = MPI_Wtime();
+ ierr = KSPSetUp(ksp);CHKERRQ(ierr);
+ ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
+ T2 = MPI_Wtime();
+
+
+ ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
+ ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
+
+
+
+ ierr = MatMult(A,x,u);CHKERRQ(ierr);
+ ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr);
+ ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr);
+ ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr);
+
+ {
+
+ Vec x2;
+ Vec sol;
+ VecDuplicate(b,&x2);
+ VecDuplicate(b,&sol);
+
+ KrylovMinimize(A, b, x2);
+
+
+
+ MatMult(A,x2,sol);
+ VecAXPY(sol,-1,b);
+ VecNorm(sol, NORM_2, &norm);
+ ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Error Krylov Minimization %g\n",norm);
+
+ }
+
+{
+
+ Vec x2;
+ Vec sol;
+ VecDuplicate(b,&x2);
+ VecDuplicate(b,&sol);
+
+ KrylovMinimizeLSQR(A, b, x2);
+
+
+
+ MatMult(A,x2,sol);
+ VecAXPY(sol,-1,b);
+ VecNorm(sol, NORM_2, &norm);
+ ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Error KrylovLSQR Minimization %g\n",norm);
+
+ }
+
+
+ /*
+ Free work space. All PETSc objects should be destroyed when they
+ are no longer needed.
+ */
+ ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = VecDestroy(&b);CHKERRQ(ierr);
+ ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr);
+ ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+ PetscPreLoadEnd();
+ /* -----------------------------------------------------------
+ End of linear solver loop
+ ----------------------------------------------------------- */
+
+ ierr = PetscFinalize();
+ return 0;
+}
+
