Actual source code: ex8.c
2: static char help[] = "Solves a linear system in parallel with KSP. \n\
3: Contributed by Jose E. Roman, SLEPc developer, for testing repeated call of KSPSetOperators(), 2014 \n\n";
5: #include <petscksp.h>
6: int main(int argc,char **args)
7: {
8: Vec x,b,u; /* approx solution, RHS, exact solution */
9: Mat A; /* linear system matrix */
10: KSP ksp; /* linear solver context */
11: PetscRandom rctx; /* random number generator context */
12: PetscInt i,j,Ii,J,Istart,Iend,m = 8,n = 7;
14: PetscBool flg = PETSC_FALSE;
15: PetscScalar v;
16: PC pc;
17: PetscInt in;
18: Mat F,B;
19: PetscBool solve=PETSC_FALSE,sameA=PETSC_FALSE,setfromoptions_first=PETSC_FALSE;
20: #if defined(PETSC_USE_LOG)
21: PetscLogStage stage;
22: #endif
23: #if !defined(PETSC_HAVE_MUMPS)
24: PetscMPIInt size;
25: #endif
27: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
28: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
29: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
30: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
31: Compute the matrix and right-hand-side vector that define
32: the linear system, Ax = b.
33: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
34: MatCreate(PETSC_COMM_WORLD,&A);
35: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
36: MatSetFromOptions(A);
37: MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
38: MatSeqAIJSetPreallocation(A,5,NULL);
39: MatSetUp(A);
41: MatGetOwnershipRange(A,&Istart,&Iend);
43: PetscLogStageRegister("Assembly", &stage);
44: PetscLogStagePush(stage);
45: for (Ii=Istart; Ii<Iend; Ii++) {
46: v = -1.0; i = Ii/n; j = Ii - i*n;
47: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
48: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
49: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
50: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
51: v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
52: }
53: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
54: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
55: PetscLogStagePop();
57: /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
58: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
60: /* Create parallel vectors. */
61: VecCreate(PETSC_COMM_WORLD,&u);
62: VecSetSizes(u,PETSC_DECIDE,m*n);
63: VecSetFromOptions(u);
64: VecDuplicate(u,&b);
65: VecDuplicate(b,&x);
67: /*
68: Set exact solution; then compute right-hand-side vector.
69: By default we use an exact solution of a vector with all
70: elements of 1.0; Alternatively, using the runtime option
71: -random_sol forms a solution vector with random components.
72: */
73: PetscOptionsGetBool(NULL,NULL,"-random_exact_sol",&flg,NULL);
74: if (flg) {
75: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
76: PetscRandomSetFromOptions(rctx);
77: VecSetRandom(u,rctx);
78: PetscRandomDestroy(&rctx);
79: } else {
80: VecSet(u,1.0);
81: }
82: MatMult(A,u,b);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Create the linear solver and set various options
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: /* Create linear solver context */
88: KSPCreate(PETSC_COMM_WORLD,&ksp);
90: /* Set operators. */
91: KSPSetOperators(ksp,A,A);
93: KSPSetTolerances(ksp,1.e-2/((m+1)*(n+1)),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
95: PetscOptionsGetBool(NULL,NULL,"-setfromoptions_first",&setfromoptions_first,NULL);
96: if (setfromoptions_first) {
97: /* code path for changing from KSPLSQR to KSPREONLY */
98: KSPSetFromOptions(ksp);
99: }
100: KSPSetType(ksp,KSPPREONLY);
101: KSPGetPC(ksp, &pc);
102: PCSetType(pc,PCCHOLESKY);
103: #if defined(PETSC_HAVE_MUMPS)
104: #if defined(PETSC_USE_COMPLEX)
105: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Spectrum slicing with MUMPS is not available for complex scalars");
106: #endif
107: PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
108: /*
109: must use runtime option '-mat_mumps_icntl_13 1' (turn off ScaLAPACK for
110: matrix inertia), currently there is no better way of setting this in program
111: */
112: PetscOptionsInsertString(NULL,"-mat_mumps_icntl_13 1");
113: #else
114: MPI_Comm_size(PETSC_COMM_WORLD,&size);
115: if (size>1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Configure with MUMPS if you want to run this example in parallel");
116: #endif
118: if (!setfromoptions_first) {
119: /* when -setfromoptions_first is true, do not call KSPSetFromOptions() again and stick to KSPPREONLY */
120: KSPSetFromOptions(ksp);
121: }
123: /* get inertia */
124: PetscOptionsGetBool(NULL,NULL,"-solve",&solve,NULL);
125: PetscOptionsGetBool(NULL,NULL,"-sameA",&sameA,NULL);
126: KSPSetUp(ksp);
127: PCFactorGetMatrix(pc,&F);
128: MatGetInertia(F,&in,NULL,NULL);
129: PetscPrintf(PETSC_COMM_WORLD,"INERTIA=%D\n",in);
130: if (solve) {
131: PetscPrintf(PETSC_COMM_WORLD,"Solving the intermediate KSP\n");
132: KSPSolve(ksp,b,x);
133: } else {PetscPrintf(PETSC_COMM_WORLD,"NOT Solving the intermediate KSP\n");}
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Solve the linear system
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: MatDuplicate(A,MAT_COPY_VALUES,&B);
139: if (sameA) {
140: PetscPrintf(PETSC_COMM_WORLD,"Setting A\n");
141: MatAXPY(A,1.1,B,DIFFERENT_NONZERO_PATTERN);
142: KSPSetOperators(ksp,A,A);
143: } else {
144: PetscPrintf(PETSC_COMM_WORLD,"Setting B\n");
145: MatAXPY(B,1.1,A,DIFFERENT_NONZERO_PATTERN);
146: KSPSetOperators(ksp,B,B);
147: }
148: KSPSetUp(ksp);
149: PCFactorGetMatrix(pc,&F);
150: MatGetInertia(F,&in,NULL,NULL);
151: PetscPrintf(PETSC_COMM_WORLD,"INERTIA=%D\n",in);
152: KSPSolve(ksp,b,x);
153: MatDestroy(&B);
155: /* Free work space.*/
156: KSPDestroy(&ksp);
157: VecDestroy(&u); VecDestroy(&x);
158: VecDestroy(&b); MatDestroy(&A);
160: PetscFinalize();
161: return ierr;
162: }
164: /*TEST
166: build:
167: requires: !complex
169: test:
170: args:
172: test:
173: suffix: 2
174: args: -sameA
176: test:
177: suffix: 3
178: args: -ksp_lsqr_monitor -ksp_type lsqr -setfromoptions_first {{0 1}separate output}
180: TEST*/