code/ex15.c

   1
   2 static char help[] = "Solves a linear system in parallel with KSP.  Also\n\
   3 illustrates setting a user-defined shell preconditioner and using the\n\
   4 macro __FUNCT__ to define routine names for use in error handling.\n\
   5 Input parameters include:\n\
   6   -user_defined_pc : Activate a user-defined preconditioner\n\n";
   7
   8 /*T
   9    Concepts: KSP^basic parallel example
  10    Concepts: PC^setting a user-defined shell preconditioner
  11    Concepts: error handling^Using the macro __FUNCT__ to define routine names;
  12    Processors: n
  13 T*/
  14
  15 /*
  16   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
  17   automatically includes:
  18      petscsys.h       - base PETSc routines   petscvec.h - vectors
  19      petscmat.h - matrices
  20      petscis.h     - index sets            petscksp.h - Krylov subspace methods
  21      petscviewer.h - viewers               petscpc.h  - preconditioners
  22 */
  23 #include <petscksp.h>
  24
  25 /* Define context for user-provided preconditioner */
  26 typedef struct {
  27   Vec diag;
  28 } SampleShellPC;
  29
  30
  31
  32 /*
  33    User-defined routines.  Note that immediately before each routine below,
  34    we define the macro __FUNCT__ to be a string containing the routine name.
  35    If defined, this macro is used in the PETSc error handlers to provide a
  36    complete traceback of routine names.  All PETSc library routines use this
  37    macro, and users can optionally employ it as well in their application
  38    codes.  Note that users can get a traceback of PETSc errors regardless of
  39    whether they define __FUNCT__ in application codes; this macro merely
  40    provides the added traceback detail of the application routine names.
  41 */
  42 #undef __FUNCT__
  43 #define __FUNCT__ "main"
  44
  45
  46
  47
  48
  49 int KrylovMinimize(Mat A, Vec b, Vec x) {
  50
  51
  52   //Variables
  53
  54   PetscScalar  gamma, alpha, oldgamma, beta, t2;
  55   PetscReal norm=20, Eprecision=1e-8, cgprec=1e-40;
  56   PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
  57   PetscErrorCode ierr;
  58   PetscScalar T1, T2, t1;
  59   KSP ksp;
  60   PetscInt total=0;
  61   PetscInt size;
  62   PetscInt Istart,Iend;
  63   PetscInt i,its;
  64   Vec       x_old, residu;
  65   Mat S, AS;
  66   PetscScalar *array;
  67   PetscInt *ind_row;
  68   Vec Alpha, p, ss, vect, r, q, Ax;
  69   PetscScalar norm_ksp;
  70
  71
  72   PetscInt first=1;
  73
  74   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  75   ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
  76   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  77
  78
  79
  80
  81   VecGetSize(b,&size);
  82
  83   ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);
  84
  85   PetscInt aa,bb;
  86   MatGetOwnershipRange(A,&aa,&bb);
  87
  88   //  ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%D %D\n", aa,bb); CHKERRQ(ierr);
  89   //PetscSynchronizedFlush(PETSC_COMM_WORLD);
  90
  91
  92   ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
  93   ierr = MatSetSizes(S, bb-aa,  PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
  94   ierr = MatSetType(S, MATMPIDENSE); CHKERRQ(ierr);
  95   ierr = MatSetUp(S); CHKERRQ(ierr);
  96
  97   ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
  98
  99   ierr = VecCreate(PETSC_COMM_WORLD, &Alpha); CHKERRQ(ierr);
 100   ierr = VecSetSizes(Alpha, PETSC_DECIDE, ColS); CHKERRQ(ierr);
 101   ierr = VecSetFromOptions(Alpha); CHKERRQ(ierr);
 102   ierr = VecDuplicate(Alpha, &vect); CHKERRQ(ierr);
 103   ierr = VecDuplicate(Alpha, &p); CHKERRQ(ierr);
 104   ierr = VecDuplicate(Alpha, &ss); CHKERRQ(ierr);
 105   ierr = VecDuplicate(b, &r); CHKERRQ(ierr);
 106   ierr = VecDuplicate(b, &q); CHKERRQ(ierr);
 107   ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
 108
 109   ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
 110   ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
 111
 112
 113
 114   //indexes of row (these indexes are global)
 115   ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
 116   for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
 117
 118   //Initializations
 119   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
 120   ierr = KSPSetTolerances(ksp, 1e-10, 1e-10, PETSC_DEFAULT, 16); CHKERRQ(ierr);
 121   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 122
 123
 124
 125   //GMRES WITH MINIMIZATION
 126   T1 = MPI_Wtime();
 127   while(giter<Emaxiter && norm>Eprecision ){
 128     for(col=0; col<ColS  &&  norm>Eprecision; col++){
 129
 130
 131       //Solve
 132       ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
 133
 134       ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
 135       total += its; Iiter ++;
 136
 137
 138
 139       //Build S'
 140       ierr = VecGetArray(x, &array);
 141       ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
 142       VecRestoreArray(x, &array);
 143
 144
 145
 146       //Error
 147       ierr = VecCopy(x, residu); CHKERRQ(ierr);
 148       ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
 149       ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
 150
 151
 152
 153       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
 154       ierr = VecCopy(x, x_old); CHKERRQ(ierr);
 155
 156
 157     }
 158
 159
 160     //minimization step
 161     if( norm>Eprecision) {
 162
 163       MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
 164       MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
 165
 166
 167       //Build AS
 168       if(first) {
 169         MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
 170         first=0;
 171       }
 172       else
 173         MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
 174
 175
 176
 177
 178       //Minimization with CGLS
 179       MatMult(AS, Alpha, r);  VecAYPX(r, -1, b); //r_0 = b-AS*x_0
 180
 181
 182       MatMultTranspose(AS, r, p); //p_0 = AS'*r_0
 183
 184
 185
 186
 187       ierr = VecCopy(p, ss); CHKERRQ(ierr); //p_0 = s_0
 188       VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
 189       while(gamma>cgprec && iter<iterations){
 190         MatMult(AS, p, q); //q = AS*p
 191         VecNorm(q, NORM_2, &alpha); alpha = alpha * alpha; //alpha = norm2(q)^2
 192         alpha = gamma / alpha;
 193         VecAXPY(Alpha, alpha, p); //Alpha += alpha*p
 194         VecAXPY(r, -alpha, q); //r -= alpha*q
 195         MatMultTranspose(AS, r, ss); // ss = AS'*r
 196         oldgamma = gamma;
 197         VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
 198         beta = gamma / oldgamma;
 199         VecAYPX(p, beta, ss); //p = s+beta*p;
 200         iter ++;
 201       }
 202
 203       iter = 0; giter ++;
 204       //Minimizer
 205       MatMult(S, Alpha, x); //x = S*Alpha;
 206     }
 207
 208   }
 209   T2 = MPI_Wtime();
 210   t1 = T2 - T1;
 211   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
 212   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
 213
 214   return 0;
 215
 216 }
 217
 218
 219
 220
 221
 222
 223
 224
 225 int main(int argc,char **args)
 226 {
 227   Vec            x,b,u;   /* approx solution, RHS, exact solution */
 228   Mat            A;         /* linear system matrix */
 229   KSP            ksp;      /* linear solver context */
 230   PC             pc;        /* preconditioner context */
 231   PetscReal      norm;      /* norm of solution error */
 232   SampleShellPC  *shell;    /* user-defined preconditioner context */
 233   PetscScalar    v,one = 1.0,none = -1.0;
 234   PetscInt       i,j,Ii,J,Istart,Iend,m = 8,n = 7,its;
 235   PetscErrorCode ierr;
 236   PetscBool      user_defined_pc = PETSC_FALSE;
 237
 238   PetscInitialize(&argc,&args,(char*)0,help);
 239   ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr);
 240   ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
 241
 242   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 243          Compute the matrix and right-hand-side vector that define
 244          the linear system, Ax = b.
 245      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 246   /*
 247      Create parallel matrix, specifying only its global dimensions.
 248      When using MatCreate(), the matrix format can be specified at
 249      runtime. Also, the parallel partioning of the matrix is
 250      determined by PETSc at runtime.
 251   */
 252   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
 253   ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr);
 254   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
 255   ierr = MatSetUp(A);CHKERRQ(ierr);
 256
 257   /*
 258      Currently, all PETSc parallel matrix formats are partitioned by
 259      contiguous chunks of rows across the processors.  Determine which
 260      rows of the matrix are locally owned.
 261   */
 262   ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
 263
 264   /*
 265      Set matrix elements for the 2-D, five-point stencil in parallel.
 266       - Each processor needs to insert only elements that it owns
 267         locally (but any non-local elements will be sent to the
 268         appropriate processor during matrix assembly).
 269       - Always specify global rows and columns of matrix entries.
 270    */
 271   for (Ii=Istart; Ii<Iend; Ii++) {
 272     v = -1.0; i = Ii/n; j = Ii - i*n;
 273     if (i>0)   {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
 274     if (i<m-1) {J = Ii + n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
 275     if (j>0)   {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
 276     if (j<n-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
 277     v = 4.0; ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);CHKERRQ(ierr);
 278   }
 279
 280   /*
 281      Assemble matrix, using the 2-step process:
 282        MatAssemblyBegin(), MatAssemblyEnd()
 283      Computations can be done while messages are in transition
 284      by placing code between these two statements.
 285   */
 286   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 287   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 288
 289   /*
 290      Create parallel vectors.
 291       - When using VecCreate() VecSetSizes() and VecSetFromOptions(),
 292         we specify only the vector's global
 293         dimension; the parallel partitioning is determined at runtime.
 294       - Note: We form 1 vector from scratch and then duplicate as needed.
 295   */
 296   ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
 297   ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr);
 298   ierr = VecSetFromOptions(u);CHKERRQ(ierr);
 299   ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
 300   ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
 301
 302   /*
 303      Set exact solution; then compute right-hand-side vector.
 304   */
 305   ierr = VecSet(u,one);CHKERRQ(ierr);
 306   ierr = MatMult(A,u,b);CHKERRQ(ierr);
 307
 308   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 309                 Create the linear solver and set various options
 310      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 311
 312   /*
 313      Create linear solver context
 314   */
 315   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
 316
 317   /*
 318      Set operators. Here the matrix that defines the linear system
 319      also serves as the preconditioning matrix.
 320   */
 321   ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
 322
 323   /*
 324      Set linear solver defaults for this problem (optional).
 325      - By extracting the KSP and PC contexts from the KSP context,
 326        we can then directly call any KSP and PC routines
 327        to set various options.
 328   */
 329   ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
 330   ierr = KSPSetTolerances(ksp,1e-9,1e-9,PETSC_DEFAULT,5000000);CHKERRQ(ierr);
 331
 332   /*
 333     Set runtime options, e.g.,
 334         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
 335     These options will override those specified above as long as
 336     KSPSetFromOptions() is called _after_ any other customization
 337     routines.
 338   */
 339   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
 340
 341   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 342                       Solve the linear system
 343      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 344   PetscScalar T1,T2;
 345         T1 = MPI_Wtime();
 346   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 347         T2 = MPI_Wtime();
 348   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 349                       Check solution and clean up
 350      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 351
 352   /*
 353      Check the error
 354   */
 355   Vec sol;
 356   VecDuplicate(b,&sol);
 357   MatMult(A,x,sol);
 358   VecAXPY(sol,-1,b);
 359   VecNorm(sol, NORM_2, &norm);
 360   ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
 361   ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr);
 362   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n\n\n", T2-T1); CHKERRQ(ierr);
 363
 364
 365
 366
 367
 368
 369
 370
 371  {
 372
 373     Vec x2;
 374     Vec sol;
 375     VecDuplicate(b,&x2);
 376     VecDuplicate(b,&sol);
 377
 378     KrylovMinimize(A, b, x2);
 379
 380
 381
 382     MatMult(A,x2,sol);
 383     VecAXPY(sol,-1,b);
 384     VecNorm(sol, NORM_2, &norm);
 385     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Error Krylov Minimization %g\n",norm);
 386
 387   }
 388
 389
 390
 391
 392
 393   /*
 394      Free work space.  All PETSc objects should be destroyed when they
 395      are no longer needed.
 396   */
 397   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
 398   ierr = VecDestroy(&u);CHKERRQ(ierr);  ierr = VecDestroy(&x);CHKERRQ(ierr);
 399   ierr = VecDestroy(&b);CHKERRQ(ierr);  ierr = MatDestroy(&A);CHKERRQ(ierr);
 400
 401   ierr = PetscFinalize();
 402   return 0;
 403
 404 }
 405