Actual source code: ex58.c


  2: static char help[] = "Solves a tridiagonal linear system with KSP.\n\n";

  4: /*
  5:   Modified from ex1.c for testing matrix operations when matrix structure is changed.
  6:   Contributed by Jose E. Roman, Feb. 2012.
  7: */
  8: #include <petscksp.h>

 10: int main(int argc, char **args)
 11: {
 12:   Vec         x, b, u; /* approx solution, RHS, exact solution */
 13:   Mat         A, B, C; /* linear system matrix */
 14:   KSP         ksp;     /* linear solver context */
 15:   PC          pc;      /* preconditioner context */
 16:   PetscReal   norm;    /* norm of solution error */
 17:   PetscInt    i, n = 20, col[3], its;
 18:   PetscMPIInt size;
 19:   PetscScalar one          = 1.0, value[3];
 20:   PetscBool   nonzeroguess = PETSC_FALSE;

 23:   PetscInitialize(&argc, &args, (char *)0, help);
 24:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 26:   PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
 27:   PetscOptionsGetBool(NULL, NULL, "-nonzero_guess", &nonzeroguess, NULL);

 29:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 30:          Compute the matrix and right-hand-side vector that define
 31:          the linear system, Ax = b.
 32:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 34:   /*
 35:      Create vectors.  Note that we form 1 vector from scratch and
 36:      then duplicate as needed.
 37:   */
 38:   VecCreate(PETSC_COMM_WORLD, &x);
 39:   PetscObjectSetName((PetscObject)x, "Solution");
 40:   VecSetSizes(x, PETSC_DECIDE, n);
 41:   VecSetFromOptions(x);
 42:   VecDuplicate(x, &b);
 43:   VecDuplicate(x, &u);

 45:   /*
 46:      Create matrix.  When using MatCreate(), the matrix format can
 47:      be specified at runtime.

 49:      Performance tuning note:  For problems of substantial size,
 50:      preallocation of matrix memory is crucial for attaining good
 51:      performance. See the matrix chapter of the users manual for details.
 52:   */
 53:   MatCreate(PETSC_COMM_WORLD, &A);
 54:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
 55:   MatSetFromOptions(A);
 56:   MatSetUp(A);

 58:   /*
 59:      Assemble matrix
 60:   */
 61:   value[0] = -1.0;
 62:   value[1] = 2.0;
 63:   value[2] = -1.0;
 64:   for (i = 1; i < n - 1; i++) {
 65:     col[0] = i - 1;
 66:     col[1] = i;
 67:     col[2] = i + 1;
 68:     MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
 69:   }
 70:   i      = n - 1;
 71:   col[0] = n - 2;
 72:   col[1] = n - 1;
 73:   MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
 74:   i        = 0;
 75:   col[0]   = 0;
 76:   col[1]   = 1;
 77:   value[0] = 2.0;
 78:   value[1] = -1.0;
 79:   MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
 80:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 81:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

 83:   /*
 84:      jroman: added matrices
 85:   */
 86:   MatCreate(PETSC_COMM_WORLD, &B);
 87:   MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, n, n);
 88:   MatSetFromOptions(B);
 89:   MatSetUp(B);
 90:   for (i = 0; i < n; i++) {
 91:     MatSetValue(B, i, i, value[1], INSERT_VALUES);
 92:     if (n - i + n / 3 < n) {
 93:       MatSetValue(B, n - i + n / 3, i, value[0], INSERT_VALUES);
 94:       MatSetValue(B, i, n - i + n / 3, value[0], INSERT_VALUES);
 95:     }
 96:   }
 97:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
 98:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
 99:   MatDuplicate(A, MAT_COPY_VALUES, &C);
100:   MatAXPY(C, 2.0, B, DIFFERENT_NONZERO_PATTERN);

102:   /*
103:      Set exact solution; then compute right-hand-side vector.
104:   */
105:   VecSet(u, one);
106:   MatMult(C, u, b);

108:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109:                 Create the linear solver and set various options
110:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111:   /*
112:      Create linear solver context
113:   */
114:   KSPCreate(PETSC_COMM_WORLD, &ksp);

116:   /*
117:      Set operators. Here the matrix that defines the linear system
118:      also serves as the preconditioning matrix.
119:   */
120:   KSPSetOperators(ksp, C, C);

122:   /*
123:      Set linear solver defaults for this problem (optional).
124:      - By extracting the KSP and PC contexts from the KSP context,
125:        we can then directly call any KSP and PC routines to set
126:        various options.
127:      - The following four statements are optional; all of these
128:        parameters could alternatively be specified at runtime via
129:        KSPSetFromOptions();
130:   */
131:   KSPGetPC(ksp, &pc);
132:   PCSetType(pc, PCJACOBI);
133:   KSPSetTolerances(ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);

135:   /*
136:     Set runtime options, e.g.,
137:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
138:     These options will override those specified above as long as
139:     KSPSetFromOptions() is called _after_ any other customization
140:     routines.
141:   */
142:   KSPSetFromOptions(ksp);

144:   if (nonzeroguess) {
145:     PetscScalar p = .5;
146:     VecSet(x, p);
147:     KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
148:   }

150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151:                       Solve the linear system
152:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153:   /*
154:      Solve linear system
155:   */
156:   KSPSolve(ksp, b, x);

158:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159:                       Check solution and clean up
160:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161:   /*
162:      Check the error
163:   */
164:   VecAXPY(x, -1.0, u);
165:   VecNorm(x, NORM_2, &norm);
166:   KSPGetIterationNumber(ksp, &its);
167:   PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);

169:   /*
170:      Free work space.  All PETSc objects should be destroyed when they
171:      are no longer needed.
172:   */
173:   VecDestroy(&x);
174:   VecDestroy(&u);
175:   VecDestroy(&b);
176:   MatDestroy(&A);
177:   MatDestroy(&B);
178:   MatDestroy(&C);
179:   KSPDestroy(&ksp);

181:   /*
182:      Always call PetscFinalize() before exiting a program.  This routine
183:        - finalizes the PETSc libraries as well as MPI
184:        - provides summary and diagnostic information if certain runtime
185:          options are chosen (e.g., -log_view).
186:   */
187:   PetscFinalize();
188:   return 0;
189: }

191: /*TEST

193:    test:
194:       args: -mat_type aij
195:       output_file: output/ex58.out

197:    test:
198:       suffix: baij
199:       args: -mat_type baij
200:       output_file: output/ex58.out

202:    test:
203:       suffix: sbaij
204:       args: -mat_type sbaij
205:       output_file: output/ex58.out

207: TEST*/