Actual source code: ex33.c
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;
13: KSP ksp;
14: PC pc;
15: PetscInt N, n=10, m, Istart, Iend, II, J, i,j;
16: PetscInt nneg, nzero, npos;
17: PetscScalar v,sigma;
18: PetscBool flag,loadA=PETSC_FALSE,loadB=PETSC_FALSE;
19: char file[2][PETSC_MAX_PATH_LEN];
20: PetscViewer viewer;
21: PetscMPIInt rank;
23: PetscInitialize(&argc,&args,(char*)0,help);
24: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
25: Compute the matrices that define the eigensystem, Ax=kBx
26: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
27: PetscOptionsGetString(NULL,NULL,"-fA",file[0],sizeof(file[0]),&loadA);
28: if (loadA) {
29: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer);
30: MatCreate(PETSC_COMM_WORLD,&A);
31: MatLoad(A,viewer);
32: PetscViewerDestroy(&viewer);
34: PetscOptionsGetString(NULL,NULL,"-fB",file[1],sizeof(file[1]),&loadB);
35: if (loadB) {
36: /* load B to get A = A + sigma*B */
37: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer);
38: MatCreate(PETSC_COMM_WORLD,&B);
39: MatLoad(B,viewer);
40: PetscViewerDestroy(&viewer);
41: }
42: }
44: if (!loadA) { /* Matrix A is copied from slepc-3.0.0-p6/src/examples/ex13.c. */
45: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
46: PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag);
47: if (!flag) m=n;
48: N = n*m;
49: MatCreate(PETSC_COMM_WORLD,&A);
50: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
51: MatSetFromOptions(A);
52: MatSetUp(A);
54: MatGetOwnershipRange(A,&Istart,&Iend);
55: for (II=Istart; II<Iend; II++) {
56: v = -1.0; i = II/n; j = II-i*n;
57: if (i>0) { J=II-n; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
58: if (i<m-1) { J=II+n; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
59: if (j>0) { J=II-1; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
60: if (j<n-1) { J=II+1; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
61: v=4.0; MatSetValues(A,1,&II,1,&II,&v,INSERT_VALUES);
63: }
64: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
65: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
66: }
67: /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
69: if (!loadB) {
70: MatGetLocalSize(A,&m,&n);
71: MatCreate(PETSC_COMM_WORLD,&B);
72: MatSetSizes(B,m,n,PETSC_DECIDE,PETSC_DECIDE);
73: MatSetFromOptions(B);
74: MatSetUp(B);
75: MatGetOwnershipRange(A,&Istart,&Iend);
77: for (II=Istart; II<Iend; II++) {
78: v=1.0; MatSetValues(B,1,&II,1,&II,&v,INSERT_VALUES);
79: }
80: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
81: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
82: }
83: /* MatView(B,PETSC_VIEWER_STDOUT_WORLD); */
85: /* Set a shift: A = A - sigma*B */
86: PetscOptionsGetScalar(NULL,NULL,"-sigma",&sigma,&flag);
87: if (flag) {
88: sigma = -1.0 * sigma;
89: MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN); /* A <- A - sigma*B */
90: /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
91: }
93: /* Test MatGetInertia() */
94: /* if A is symmetric, set its flag -- required by MatGetInertia() */
95: MatIsSymmetric(A,0.0,&flag);
97: KSPCreate(PETSC_COMM_WORLD,&ksp);
98: KSPSetType(ksp,KSPPREONLY);
99: KSPSetOperators(ksp,A,A);
101: KSPGetPC(ksp,&pc);
102: PCSetType(pc,PCCHOLESKY);
103: PCSetFromOptions(pc);
105: PCSetUp(pc);
106: PCFactorGetMatrix(pc,&F);
107: MatGetInertia(F,&nneg,&nzero,&npos);
109: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
110: if (rank == 0) {
111: PetscPrintf(PETSC_COMM_SELF," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);
112: }
114: /* Destroy */
115: KSPDestroy(&ksp);
116: MatDestroy(&A);
117: MatDestroy(&B);
118: PetscFinalize();
119: return 0;
120: }
122: /*TEST
124: test:
125: args: -sigma 2.0
126: requires: !complex
127: output_file: output/ex33.out
129: test:
130: suffix: mumps
131: args: -sigma 2.0 -pc_factor_mat_solver_type mumps -mat_mumps_icntl_13 1 -mat_mumps_icntl_24 1
132: requires: mumps !complex
133: output_file: output/ex33.out
135: test:
136: suffix: mumps_2
137: args: -sigma 2.0 -pc_factor_mat_solver_type mumps -mat_mumps_icntl_13 1 -mat_mumps_icntl_24 1
138: requires: mumps !complex
139: nsize: 3
140: output_file: output/ex33.out
142: test:
143: suffix: mkl_pardiso
144: args: -sigma 2.0 -pc_factor_mat_solver_type mkl_pardiso -mat_type sbaij
145: requires: mkl_pardiso !complex
146: output_file: output/ex33.out
148: test:
149: suffix: superlu_dist
150: args: -sigma 2.0 -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
151: requires: superlu_dist !complex
152: output_file: output/ex33.out
154: test:
155: suffix: superlu_dist_2
156: args: -sigma 2.0 -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
157: requires: superlu_dist !complex
158: nsize: 3
159: output_file: output/ex33.out
161: TEST*/