code/ex15.c

   1
   2 // /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 4    ./ex15 -m 400 -n 400 -ksp_type fgmres
   3
   4
   5
   6 static char help[] = "Solves a linear system in parallel with KSP.  Also\n\
   7 illustrates setting a user-defined shell preconditioner and using the\n\
   8 macro __FUNCT__ to define routine names for use in error handling.\n\
   9 Input parameters include:\n\
  10   -user_defined_pc : Activate a user-defined preconditioner\n\n";
  11
  12 /*T
  13    Concepts: KSP^basic parallel example
  14    Concepts: PC^setting a user-defined shell preconditioner
  15    Concepts: error handling^Using the macro __FUNCT__ to define routine names;
  16    Processors: n
  17 T*/
  18
  19 /*
  20   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
  21   automatically includes:
  22      petscsys.h       - base PETSc routines   petscvec.h - vectors
  23      petscmat.h - matrices
  24      petscis.h     - index sets            petscksp.h - Krylov subspace methods
  25      petscviewer.h - viewers               petscpc.h  - preconditioners
  26 */
  27 #include <petscksp.h>
  28
  29 /* Define context for user-provided preconditioner */
  30 typedef struct {
  31   Vec diag;
  32 } SampleShellPC;
  33
  34
  35
  36 /*
  37    User-defined routines.  Note that immediately before each routine below,
  38    we define the macro __FUNCT__ to be a string containing the routine name.
  39    If defined, this macro is used in the PETSc error handlers to provide a
  40    complete traceback of routine names.  All PETSc library routines use this
  41    macro, and users can optionally employ it as well in their application
  42    codes.  Note that users can get a traceback of PETSc errors regardless of
  43    whether they define __FUNCT__ in application codes; this macro merely
  44    provides the added traceback detail of the application routine names.
  45 */
  46 #undef __FUNCT__
  47 #define __FUNCT__ "main"
  48
  49
  50
  51
  52
  53 int KrylovMinimize(Mat A, Vec b, Vec x) {
  54
  55
  56   //Variables
  57
  58   PetscScalar  gamma, alpha, oldgamma, beta;
  59   PetscReal norm=20, Eprecision=1e-6, cgprec=1e-40;
  60   PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
  61   PetscErrorCode ierr;
  62   PetscScalar T1, T2;
  63   KSP ksp;
  64   PetscInt total=0;
  65   PetscInt size;
  66   PetscInt Istart,Iend;
  67   PetscInt i,its;
  68   Vec       x_old, residu;
  69   Mat S, AS;
  70   PetscScalar *array;
  71   PetscInt *ind_row;
  72   Vec Alpha, p, ss, vect, r, q, Ax;
  73
  74
  75   PetscInt first=1;
  76
  77   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  78   ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
  79   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  80
  81
  82
  83
  84   VecGetSize(b,&size);
  85
  86   ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);
  87
  88   PetscInt aa,bb;
  89   MatGetOwnershipRange(A,&aa,&bb);
  90
  91   //  ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%D %D\n", aa,bb); CHKERRQ(ierr);
  92   //PetscSynchronizedFlush(PETSC_COMM_WORLD);
  93
  94
  95   ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
  96   ierr = MatSetSizes(S, bb-aa,  PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
  97   ierr = MatSetType(S, MATMPIDENSE); CHKERRQ(ierr);
  98   ierr = MatSetUp(S); CHKERRQ(ierr);
  99
 100   ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
 101
 102   ierr = VecCreate(PETSC_COMM_WORLD, &Alpha); CHKERRQ(ierr);
 103   ierr = VecSetSizes(Alpha, PETSC_DECIDE, ColS); CHKERRQ(ierr);
 104   ierr = VecSetFromOptions(Alpha); CHKERRQ(ierr);
 105   ierr = VecDuplicate(Alpha, &vect); CHKERRQ(ierr);
 106   ierr = VecDuplicate(Alpha, &p); CHKERRQ(ierr);
 107   ierr = VecDuplicate(Alpha, &ss); CHKERRQ(ierr);
 108   ierr = VecDuplicate(b, &r); CHKERRQ(ierr);
 109   ierr = VecDuplicate(b, &q); CHKERRQ(ierr);
 110   ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
 111
 112   ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
 113   ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
 114
 115
 116
 117   //indexes of row (these indexes are global)
 118   ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
 119   for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
 120
 121   //Initializations
 122   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
 123   ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
 124   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 125
 126
 127
 128   //GMRES WITH MINIMIZATION
 129   T1 = MPI_Wtime();
 130   ierr = KSPSetUp(ksp); CHKERRQ(ierr);
 131   while(giter<Emaxiter && norm>Eprecision ){
 132     for(col=0; col<ColS  &&  norm>Eprecision; col++){
 133
 134
 135       //Solve
 136       ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
 137
 138       ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
 139       total += its; Iiter ++;
 140
 141
 142
 143       //Build S'
 144       ierr = VecGetArray(x, &array);
 145       ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
 146       VecRestoreArray(x, &array);
 147
 148
 149
 150       KSPGetResidualNorm(ksp,&norm);
 151
 152       /*      //Error
 153       ierr = VecCopy(x, residu); CHKERRQ(ierr);
 154       ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
 155       ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
 156        */
 157
 158
 159       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
 160       ierr = VecCopy(x, x_old); CHKERRQ(ierr);
 161
 162
 163     }
 164
 165
 166     //minimization step
 167     if( norm>Eprecision) {
 168
 169       MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
 170       MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
 171
 172
 173       //Build AS
 174       if(first) {
 175         MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
 176         first=0;
 177       }
 178       else
 179         MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
 180
 181
 182
 183
 184       //Minimization with CGLS
 185       MatMult(AS, Alpha, r);  VecAYPX(r, -1, b); //r_0 = b-AS*x_0
 186
 187
 188       MatMultTranspose(AS, r, p); //p_0 = AS'*r_0
 189
 190
 191
 192
 193       ierr = VecCopy(p, ss); CHKERRQ(ierr); //p_0 = s_0
 194       VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
 195       while(gamma>cgprec && iter<iterations){
 196         MatMult(AS, p, q); //q = AS*p
 197         VecNorm(q, NORM_2, &alpha); alpha = alpha * alpha; //alpha = norm2(q)^2
 198         alpha = gamma / alpha;
 199         VecAXPY(Alpha, alpha, p); //Alpha += alpha*p
 200         VecAXPY(r, -alpha, q); //r -= alpha*q
 201         MatMultTranspose(AS, r, ss); // ss = AS'*r
 202         oldgamma = gamma;
 203         VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
 204         beta = gamma / oldgamma;
 205         VecAYPX(p, beta, ss); //p = s+beta*p;
 206         iter ++;
 207       }
 208
 209       iter = 0; giter ++;
 210       //Minimizer
 211       MatMult(S, Alpha, x); //x = S*Alpha;
 212     }
 213
 214   }
 215   T2 = MPI_Wtime();
 216   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
 217   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
 218
 219   return 0;
 220
 221 }
 222
 223
 224
 225
 226
 227
 228
 229
 230
 231 int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 232
 233
 234   //Variables
 235
 236   PetscScalar  alpha, beta;
 237   PetscReal norm=20, Eprecision=1e-6, tol=1e-40;
 238   PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=20, Iiter=0;
 239   PetscErrorCode ierr;
 240   PetscScalar T1, T2;
 241   KSP ksp;
 242   PetscInt total=0;
 243   PetscInt size;
 244   PetscInt Istart,Iend;
 245   PetscInt i,its;
 246   Vec       x_old, residu;
 247   Mat S, AS;
 248   PetscScalar *array;
 249   PetscInt *ind_row;
 250   Vec  Ax;
 251   PetscScalar norm_ksp;
 252   Vec u,v,d,uu,vt,zero_long,zero_short,x_lsqr;
 253
 254   PetscInt first=1;
 255
 256   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
 257   ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
 258   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
 259
 260
 261
 262
 263   VecGetSize(b,&size);
 264
 265   ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);
 266
 267   PetscInt aa,bb;
 268   MatGetOwnershipRange(A,&aa,&bb);
 269
 270   //  ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%D %D\n", aa,bb); CHKERRQ(ierr);
 271   //PetscSynchronizedFlush(PETSC_COMM_WORLD);
 272
 273
 274   ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
 275   ierr = MatSetSizes(S, bb-aa,  PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
 276   ierr = MatSetType(S, MATMPIDENSE); CHKERRQ(ierr);
 277   ierr = MatSetUp(S); CHKERRQ(ierr);
 278
 279   ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
 280
 281
 282
 283   ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
 284
 285   ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
 286   ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
 287
 288
 289   //long vector
 290   ierr = VecDuplicate(b,&u);CHKERRQ(ierr);
 291
 292
 293   ierr = VecDuplicate(b,&uu);CHKERRQ(ierr);
 294   ierr = VecDuplicate(b,&zero_long);CHKERRQ(ierr);
 295   ierr = VecSet(zero_long,0);CHKERRQ(ierr);
 296
 297   //small vector
 298   ierr = VecCreate(PETSC_COMM_WORLD, &v); CHKERRQ(ierr);
 299   ierr = VecSetSizes(v, PETSC_DECIDE, ColS); CHKERRQ(ierr);
 300   ierr = VecSetFromOptions(v); CHKERRQ(ierr);
 301   ierr = VecDuplicate(v,&zero_short);CHKERRQ(ierr);
 302   ierr = VecSet(zero_short,0);CHKERRQ(ierr);
 303   ierr = VecDuplicate(v,&d);CHKERRQ(ierr);
 304   ierr = VecDuplicate(v,&vt);CHKERRQ(ierr);
 305   ierr = VecDuplicate(v,&x_lsqr);CHKERRQ(ierr);
 306
 307
 308   //indexes of row (these indexes are global)
 309   ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
 310   for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
 311
 312   //Initializations
 313   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
 314   ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
 315   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 316
 317
 318
 319
 320   //GMRES WITH MINIMIZATION
 321   T1 = MPI_Wtime();
 322   ierr = KSPSetUp(ksp); CHKERRQ(ierr);
 323   while(giter<Emaxiter && norm>Eprecision ){
 324     for(col=0; col<ColS  &&  norm>Eprecision; col++){
 325
 326
 327       //Solve
 328       ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
 329
 330       ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
 331       total += its; Iiter ++;
 332
 333
 334
 335       //Build S'
 336       ierr = VecGetArray(x, &array);
 337       ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
 338       VecRestoreArray(x, &array);
 339
 340
 341
 342       KSPGetResidualNorm(ksp,&norm);
 343
 344       /*
 345       //Error
 346       ierr = VecCopy(x, residu); CHKERRQ(ierr);
 347       ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
 348       ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
 349        */
 350
 351
 352       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
 353       ierr = VecCopy(x, x_old); CHKERRQ(ierr);
 354
 355
 356     }
 357
 358
 359     //minimization step
 360     if( norm>Eprecision) {
 361
 362       MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
 363       MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
 364
 365
 366
 367
 368       //Build AS
 369       if(first) {
 370         MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
 371
 372         first=0;
 373       }
 374       else
 375         MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
 376
 377
 378
 379
 380
 381       //LSQR
 382       //LSQR
 383       //LSQR
 384
 385
 386
 387       PetscScalar n2b,tolb,normr,c,s,phibar,normar,norma,thet,rhot,rho,phi;
 388       PetscInt stag;
 389       tolb = tol * n2b;
 390       VecNorm(b, NORM_2, &n2b); //n2b = norm(b);
 391       ierr = VecCopy(b, u); CHKERRQ(ierr);   //u=b
 392       VecNorm(u, NORM_2, &beta); // beta=norm(u)
 393       normr=beta;
 394       if (beta != 0) {
 395         VecAYPX(u,1/beta,zero_long); //  u = u / beta;
 396       }
 397       c=1;
 398       s=0;
 399       phibar=beta;
 400       MatMultTranspose(AS, u, v); //v=A'*u
 401       ierr = VecSet(x_lsqr,0);CHKERRQ(ierr);
 402       VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
 403       if (alpha != 0) {
 404         VecAYPX(v,1/alpha,zero_short); //  v = v / alpha;
 405       }
 406       ierr = VecSet(d,0);CHKERRQ(ierr);
 407       normar = alpha * beta;
 408       norma=0;
 409       //stag=0;
 410       for(i=0;i<iterations;i++) {
 411         MatMult(AS, v, uu); //uu=A*v
 412         VecAYPX(u, -alpha, uu); //u = uu-alpha*u;
 413         VecNorm(u, NORM_2, &beta); // beta=norm(u)
 414         VecAYPX(u,1/beta,zero_long); //  u = u / beta;
 415         norma=sqrt(norma*norma+alpha*alpha+beta*beta); // norma = norm([norma alpha beta]);
 416         thet = - s * alpha;
 417         rhot = c * alpha;
 418         rho = sqrt(rhot*rhot + beta*beta);
 419         c = rhot / rho;
 420         s = - beta / rho;
 421         phi = c * phibar;
 422         if (phi == 0) {             // stagnation of the method
 423           stag = 1;
 424         }
 425         phibar = s * phibar;
 426         VecAYPX(d,-thet,v);       //d = (v - thet * d);
 427         VecAYPX(d,1/rho,zero_short);     //d=d/ rho;
 428
 429         /*
 430         if (normar/(norma*normr) <= tol) { // check for convergence in min{|b-A*x|}
 431           break;
 432         }
 433         if (normr <= tolb) {           // check for convergence in A*x=b
 434           break;
 435         }
 436          */
 437
 438         VecAXPY(x_lsqr,phi,d);     // x_lsqr=x_lsqr+phi*d
 439         normr = abs(s) * normr;
 440         MatMultTranspose(AS, u, vt); //vt=A'*u;
 441         VecAYPX(v,-beta,vt);  // v = vt - beta * v;
 442         VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
 443         VecAYPX(v,1/alpha,zero_short); //  v = v / alpha;
 444         normar = alpha * abs( s * phi);
 445       }
 446
 447
 448
 449       iter = 0; giter ++;
 450       //Minimizer
 451       MatMult(S, x_lsqr, x); //x = S*x_small;
 452     }
 453
 454   }
 455   T2 = MPI_Wtime();
 456   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time LSQR : %g (s)\n", T2-T1); CHKERRQ(ierr);
 457   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations LSQR : %D\n", total); CHKERRQ(ierr);
 458
 459   return 0;
 460
 461 }
 462
 463
 464
 465
 466
 467
 468
 469
 470 int main(int argc,char **args)
 471 {
 472   Vec            x,b,u;   /* approx solution, RHS, exact solution */
 473   Mat            A;         /* linear system matrix */
 474   KSP            ksp;      /* linear solver context */
 475   PC             pc;        /* preconditioner context */
 476   PetscReal      norm;      /* norm of solution error */
 477   PetscScalar    v,one = 1.0;
 478   PetscInt       i,j,Ii,J,Istart,Iend,m = 8,n = 7,its;
 479   PetscErrorCode ierr;
 480
 481   PetscInitialize(&argc,&args,(char*)0,help);
 482   ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr);
 483   ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
 484
 485   PetscMPIInt size;
 486   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 487   PetscPrintf(PETSC_COMM_WORLD,"Number of processors = %d\n",size);
 488
 489   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 490          Compute the matrix and right-hand-side vector that define
 491          the linear system, Ax = b.
 492      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 493   /*
 494      Create parallel matrix, specifying only its global dimensions.
 495      When using MatCreate(), the matrix format can be specified at
 496      runtime. Also, the parallel partioning of the matrix is
 497      determined by PETSc at runtime.
 498   */
 499   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
 500   ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr);
 501   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
 502   ierr = MatSetUp(A);CHKERRQ(ierr);
 503
 504   /*
 505      Currently, all PETSc parallel matrix formats are partitioned by
 506      contiguous chunks of rows across the processors.  Determine which
 507      rows of the matrix are locally owned.
 508   */
 509   ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
 510
 511   /*
 512      Set matrix elements for the 2-D, five-point stencil in parallel.
 513       - Each processor needs to insert only elements that it owns
 514         locally (but any non-local elements will be sent to the
 515         appropriate processor during matrix assembly).
 516       - Always specify global rows and columns of matrix entries.
 517    */
 518   for (Ii=Istart; Ii<Iend; Ii++) {
 519     v = -1.0; i = Ii/n; j = Ii - i*n;
 520     PetscScalar v2=-1.;
 521     if (i>0)   {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v2,INSERT_VALUES);CHKERRQ(ierr);}
 522     if (i<m-1) {J = Ii + n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
 523     if (j>0)   {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
 524     if (j<n-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
 525     v = 4.0; ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);CHKERRQ(ierr);
 526   }
 527
 528   /*
 529      Assemble matrix, using the 2-step process:
 530        MatAssemblyBegin(), MatAssemblyEnd()
 531      Computations can be done while messages are in transition
 532      by placing code between these two statements.
 533   */
 534   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 535   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 536
 537   /*
 538      Create parallel vectors.
 539       - When using VecCreate() VecSetSizes() and VecSetFromOptions(),
 540         we specify only the vector's global
 541         dimension; the parallel partitioning is determined at runtime.
 542       - Note: We form 1 vector from scratch and then duplicate as needed.
 543   */
 544   ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
 545   ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr);
 546   ierr = VecSetFromOptions(u);CHKERRQ(ierr);
 547   ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
 548   ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
 549
 550   /*
 551      Set exact solution; then compute right-hand-side vector.
 552   */
 553   ierr = VecSet(u,one);CHKERRQ(ierr);
 554   ierr = MatMult(A,u,b);CHKERRQ(ierr);
 555
 556   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 557                 Create the linear solver and set various options
 558      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 559
 560   /*
 561      Create linear solver context
 562   */
 563   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
 564
 565   /*
 566      Set operators. Here the matrix that defines the linear system
 567      also serves as the preconditioning matrix.
 568   */
 569   ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
 570
 571   /*
 572      Set linear solver defaults for this problem (optional).
 573      - By extracting the KSP and PC contexts from the KSP context,
 574        we can then directly call any KSP and PC routines
 575        to set various options.
 576   */
 577   ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
 578   ierr = KSPSetTolerances(ksp,1e-7,1e-7,PETSC_DEFAULT,5000000);CHKERRQ(ierr);
 579
 580
 581   //  PC             pc;
 582   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
 583   KSPGetPC(ksp, &pc);
 584   PCType         type;
 585   PCGetType(pc, &type);
 586
 587   PetscPrintf(PETSC_COMM_WORLD, "PC TYPE %s  \n", type);
 588
 589
 590   /*
 591     Set runtime options, e.g.,
 592         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
 593     These options will override those specified above as long as
 594     KSPSetFromOptions() is called _after_ any other customization
 595     routines.
 596   */
 597   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
 598
 599   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 600                       Solve the linear system
 601      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 602   PetscScalar T1,T2;
 603         T1 = MPI_Wtime();
 604         ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 605         T2 = MPI_Wtime();
 606   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 607                       Check solution and clean up
 608      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 609
 610   /*
 611      Check the error
 612   */
 613   Vec sol;
 614   VecDuplicate(b,&sol);
 615   MatMult(A,x,sol);
 616   VecAXPY(sol,-1,b);
 617   VecNorm(sol, NORM_2, &norm);
 618   ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
 619   ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr);
 620   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n\n\n", T2-T1); CHKERRQ(ierr);
 621
 622
 623
 624
 625
 626
 627
 628
 629  {
 630
 631     Vec x2;
 632     Vec sol;
 633     VecDuplicate(b,&x2);
 634     VecDuplicate(b,&sol);
 635
 636     KrylovMinimize(A, b, x2);
 637
 638
 639
 640     MatMult(A,x2,sol);
 641     VecAXPY(sol,-1,b);
 642     VecNorm(sol, NORM_2, &norm);
 643     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Error Krylov Minimization %g\n",norm);
 644
 645   }
 646
 647
 648  {
 649
 650     Vec x2;
 651     Vec sol;
 652     VecDuplicate(b,&x2);
 653     VecDuplicate(b,&sol);
 654
 655     KrylovMinimizeLSQR(A, b, x2);
 656
 657
 658
 659     MatMult(A,x2,sol);
 660     VecAXPY(sol,-1,b);
 661     VecNorm(sol, NORM_2, &norm);
 662     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Error Krylov Minimization LSQR %g\n",norm);
 663
 664   }
 665
 666
 667
 668   /*
 669      Free work space.  All PETSc objects should be destroyed when they
 670      are no longer needed.
 671   */
 672   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
 673   ierr = VecDestroy(&u);CHKERRQ(ierr);  ierr = VecDestroy(&x);CHKERRQ(ierr);
 674   ierr = VecDestroy(&b);CHKERRQ(ierr);  ierr = MatDestroy(&A);CHKERRQ(ierr);
 675
 676   ierr = PetscFinalize();
 677   return 0;
 678
 679 }
 680