Actual source code: ex50.c


  2: static char help[] = "Tests SeqBAIJ point block Jacobi for different block sizes\n\n";

  4: #include <petscksp.h>

  6: int main(int argc,char **args)
  7: {
  8:   Vec            x,b,u;
  9:   Mat            A;        /* linear system matrix */
 10:   KSP            ksp;     /* linear solver context */
 11:   PetscRandom    rctx;     /* random number generator context */
 12:   PetscReal      norm;     /* norm of solution error */
 13:   PetscInt       i,j,k,l,n = 27,its,bs = 2,Ii,J;
 14:   PetscScalar    v;

 16:   PetscInitialize(&argc,&args,(char*)0,help);
 17:   PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
 18:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);

 20:   MatCreate(PETSC_COMM_WORLD,&A);
 21:   MatSetSizes(A,n*bs,n*bs,PETSC_DETERMINE,PETSC_DETERMINE);
 22:   MatSetBlockSize(A,bs);
 23:   MatSetType(A,MATSEQBAIJ);
 24:   MatSetFromOptions(A);
 25:   MatSetUp(A);

 27:   /*
 28:      Don't bother to preallocate matrix
 29:   */
 30:   PetscRandomCreate(PETSC_COMM_SELF,&rctx);
 31:   for (i=0; i<n; i++) {
 32:     for (j=0; j<n; j++) {
 33:       PetscRandomGetValue(rctx,&v);
 34:       if (PetscRealPart(v) < .25 || i == j) {
 35:         for (k=0; k<bs; k++) {
 36:           for (l=0; l<bs; l++) {
 37:             PetscRandomGetValue(rctx,&v);
 38:             Ii = i*bs + k;
 39:             J = j*bs + l;
 40:             if (Ii == J) v += 10.;
 41:             MatSetValue(A,Ii,J,v,INSERT_VALUES);
 42:           }
 43:         }
 44:       }
 45:     }
 46:   }

 48:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 49:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 51:   VecCreate(PETSC_COMM_WORLD,&u);
 52:   VecSetSizes(u,PETSC_DECIDE,n*bs);
 53:   VecSetFromOptions(u);
 54:   VecDuplicate(u,&b);
 55:   VecDuplicate(b,&x);

 57:   VecSet(u,1.0);
 58:   MatMult(A,u,b);

 60:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 61:                 Create the linear solver and set various options
 62:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 64:   /*
 65:      Create linear solver context
 66:   */
 67:   KSPCreate(PETSC_COMM_WORLD,&ksp);

 69:   /*
 70:      Set operators. Here the matrix that defines the linear system
 71:      also serves as the preconditioning matrix.
 72:   */
 73:   KSPSetOperators(ksp,A,A);

 75:   KSPSetFromOptions(ksp);

 77:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 78:                       Solve the linear system
 79:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 81:   KSPSolve(ksp,b,x);

 83:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 84:                       Check solution and clean up
 85:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 87:   /*
 88:      Check the error
 89:   */
 90:   VecAXPY(x,-1.0,u);
 91:   VecNorm(x,NORM_2,&norm);
 92:   KSPGetIterationNumber(ksp,&its);

 94:   /*
 95:      Print convergence information.  PetscPrintf() produces a single
 96:      print statement from all processes that share a communicator.
 97:      An alternative is PetscFPrintf(), which prints to a file.
 98:   */
 99:   if (norm > .1) {
100:     PetscPrintf(PETSC_COMM_WORLD,"Norm of residual %g iterations %D bs %D\n",(double)norm,its,bs);
101:   }

103:   /*
104:      Free work space.  All PETSc objects should be destroyed when they
105:      are no longer needed.
106:   */
107:   KSPDestroy(&ksp);
108:   VecDestroy(&u);
109:   VecDestroy(&x);
110:   VecDestroy(&b);
111:   MatDestroy(&A);
112:   PetscRandomDestroy(&rctx);

114:   /*
115:      Always call PetscFinalize() before exiting a program.  This routine
116:        - finalizes the PETSc libraries as well as MPI
117:        - provides summary and diagnostic information if certain runtime
118:          options are chosen (e.g., -log_view).
119:   */
120:   PetscFinalize();
121:   return 0;
122: }

124: /*TEST

126:    test:
127:       args: -bs {{1 2 3 4 5 6 7}} -pc_type pbjacobi

129:    test:
130:       suffix: 2
131:       args: -bs {{8 9 10 11 12 13 14 15}} -pc_type ilu

133: TEST*/