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