Actual source code: ex33.c
petsc-3.12.0 2019-09-29
1: static char help[] = "Test MatGetInertia().\n\n";
3: /*
4: Examples of command line options:
5: ./ex33 -sigma 2.0 -pc_factor_mat_solver_type mumps -mat_mumps_icntl_13 1 -mat_mumps_icntl_24 1
6: ./ex33 -sigma <shift> -fA <matrix_file>
7: */
9: #include <petscksp.h>
10: int main(int argc,char **args)
11: {
12: Mat A,B,F;
14: KSP ksp;
15: PC pc;
16: PetscInt N, n=10, m, Istart, Iend, II, J, i,j;
17: PetscInt nneg, nzero, npos;
18: PetscScalar v,sigma;
19: PetscBool flag,loadA=PETSC_FALSE,loadB=PETSC_FALSE;
20: char file[2][PETSC_MAX_PATH_LEN];
21: PetscViewer viewer;
22: PetscMPIInt rank;
24: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
25: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
26: Compute the matrices that define the eigensystem, Ax=kBx
27: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
28: PetscOptionsGetString(NULL,NULL,"-fA",file[0],PETSC_MAX_PATH_LEN,&loadA);
29: if (loadA) {
30: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer);
31: MatCreate(PETSC_COMM_WORLD,&A);
32: MatLoad(A,viewer);
33: PetscViewerDestroy(&viewer);
35: PetscOptionsGetString(NULL,NULL,"-fB",file[1],PETSC_MAX_PATH_LEN,&loadB);
36: if (loadB) {
37: /* load B to get A = A + sigma*B */
38: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer);
39: MatCreate(PETSC_COMM_WORLD,&B);
40: MatLoad(B,viewer);
41: PetscViewerDestroy(&viewer);
42: }
43: }
45: if (!loadA) { /* Matrix A is copied from slepc-3.0.0-p6/src/examples/ex13.c. */
46: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
47: PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag);
48: if (!flag) m=n;
49: N = n*m;
50: MatCreate(PETSC_COMM_WORLD,&A);
51: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
52: MatSetFromOptions(A);
53: MatSetUp(A);
55: MatGetOwnershipRange(A,&Istart,&Iend);
56: for (II=Istart; II<Iend; II++) {
57: v = -1.0; i = II/n; j = II-i*n;
58: if (i>0) { J=II-n; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
59: if (i<m-1) { J=II+n; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
60: if (j>0) { J=II-1; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
61: if (j<n-1) { J=II+1; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
62: v=4.0; MatSetValues(A,1,&II,1,&II,&v,INSERT_VALUES);
64: }
65: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
66: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
67: }
68: /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
70: if (!loadB) {
71: MatGetLocalSize(A,&m,&n);
72: MatCreate(PETSC_COMM_WORLD,&B);
73: MatSetSizes(B,m,n,PETSC_DECIDE,PETSC_DECIDE);
74: MatSetFromOptions(B);
75: MatSetUp(B);
76: MatGetOwnershipRange(A,&Istart,&Iend);
78: for (II=Istart; II<Iend; II++) {
79: /* v=4.0; MatSetValues(B,1,&II,1,&II,&v,INSERT_VALUES); */
80: v=1.0; MatSetValues(B,1,&II,1,&II,&v,INSERT_VALUES);
81: }
82: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
83: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
84: }
85: /* MatView(B,PETSC_VIEWER_STDOUT_WORLD); */
87: /* Set a shift: A = A - sigma*B */
88: PetscOptionsGetScalar(NULL,NULL,"-sigma",&sigma,&flag);
89: if (flag) {
90: sigma = -1.0 * sigma;
91: MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN); /* A <- A - sigma*B */
92: /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
93: }
95: /* Test MatGetInertia() */
96: /* if A is symmetric, set its flag -- required by MatGetInertia() */
97: MatIsSymmetric(A,0.0,&flag);
99: KSPCreate(PETSC_COMM_WORLD,&ksp);
100: KSPSetType(ksp,KSPPREONLY);
101: KSPSetOperators(ksp,A,A);
103: KSPGetPC(ksp,&pc);
104: PCSetType(pc,PCCHOLESKY);
105: PCSetFromOptions(pc);
107: PCSetUp(pc);
108: PCFactorGetMatrix(pc,&F);
109: MatGetInertia(F,&nneg,&nzero,&npos);
111: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
112: if (!rank) {
113: PetscPrintf(PETSC_COMM_SELF," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);
114: }
116: /* Destroy */
117: KSPDestroy(&ksp);
118: MatDestroy(&A);
119: MatDestroy(&B);
120: PetscFinalize();
121: return ierr;
122: }
124: /*TEST
126: test:
127: args: -sigma 2.0
128: requires: !complex
129: output_file: output/ex33.out
131: test:
132: suffix: mumps
133: args: -sigma 2.0 -pc_factor_mat_solver_type mumps -mat_mumps_icntl_13 1 -mat_mumps_icntl_24 1
134: requires: mumps !complex
135: output_file: output/ex33.out
137: test:
138: suffix: mumps_2
139: args: -sigma 2.0 -pc_factor_mat_solver_type mumps -mat_mumps_icntl_13 1 -mat_mumps_icntl_24 1
140: requires: mumps !complex
141: nsize: 3
142: output_file: output/ex33.out
144: test:
145: suffix: mkl_pardiso
146: args: -sigma 2.0 -pc_factor_mat_solver_type mkl_pardiso -mat_type sbaij
147: requires: mkl_pardiso
148: output_file: output/ex33.out
150: test:
151: suffix: superlu_dist
152: args: -sigma 2.0 -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
153: requires: superlu_dist !complex
154: output_file: output/ex33.out
156: test:
157: suffix: superlu_dist_2
158: args: -sigma 2.0 -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
159: requires: superlu_dist !complex
160: nsize: 3
161: output_file: output/ex33.out
163: TEST*/