Actual source code: ex49.c

petsc-3.12.2 2019-11-22
Report Typos and Errors

  2: static char help[] = "Tests SeqSBAIJ factorizations 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;
 15:   PetscScalar    v;
 16: 
 17:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 18:   PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
 19:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);

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

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

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

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

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

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

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

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

 76:   KSPSetFromOptions(ksp);

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

 82:   KSPSolve(ksp,b,x);

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

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

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

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

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


126: /*TEST

128:    test:
129:       args: -bs {{1 2 3 4 5 6 7 8 9 10 11 12}} -pc_type cholesky

131: TEST*/