Actual source code: ex58.c

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

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

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

 21:   PetscFunctionBeginUser;
 22:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 23:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 24:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
 25:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 26:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-nonzero_guess", &nonzeroguess, NULL));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

190: /*TEST

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

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

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

206: TEST*/