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_package mumps -mat_mumps_icntl_13 1
6: ./ex33 -sigma <shift> -fA <matrix_file>
7: */
9: #include <petscksp.h>
12: int main(int argc,char **args)
13: {
14: Mat A,B,F;
16: KSP ksp;
17: PC pc;
18: PetscInt N, n=10, m, Istart, Iend, II, J, i,j;
19: PetscInt nneg, nzero, npos;
20: PetscScalar v,sigma;
21: PetscBool flag,loadA=PETSC_FALSE,loadB=PETSC_FALSE;
22: char file[2][PETSC_MAX_PATH_LEN];
23: PetscViewer viewer;
24: PetscMPIInt rank;
26: PetscInitialize(&argc,&args,(char *)0,help);
27: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
28: Compute the matrices that define the eigensystem, Ax=kBx
29: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
30: PetscOptionsGetString(PETSC_NULL,"-fA",file[0],PETSC_MAX_PATH_LEN,&loadA);
31: if (loadA) {
32: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer);
33: MatCreate(PETSC_COMM_WORLD,&A);
34: MatSetType(A,MATSBAIJ);
35: MatLoad(A,viewer);
36: PetscViewerDestroy(&viewer);
38: PetscOptionsGetString(PETSC_NULL,"-fB",file[1],PETSC_MAX_PATH_LEN,&loadB);
39: if (loadB){
40: /* load B to get A = A + sigma*B */
41: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer);
42: MatCreate(PETSC_COMM_WORLD,&B);
43: MatSetType(B,MATSBAIJ);
44: MatLoad(B,viewer);
45: PetscViewerDestroy(&viewer);
46: }
47: }
49: if (!loadA) { /* Matrix A is copied from slepc-3.0.0-p6/src/examples/ex13.c. */
50: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
51: PetscOptionsGetInt(PETSC_NULL,"-m",&m,&flag);
52: if( flag==PETSC_FALSE ) m=n;
53: N = n*m;
54: MatCreate(PETSC_COMM_WORLD,&A);
55: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
56: MatSetType(A,MATSBAIJ);
57: MatSetFromOptions(A);
58: MatSetOption(A,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);
59: MatGetOwnershipRange(A,&Istart,&Iend);
60: for( II=Istart; II<Iend; II++ ) {
61: v = -1.0; i = II/n; j = II-i*n;
62: if(i>0) { J=II-n; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
63: if(i<m-1) { J=II+n; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
64: if(j>0) { J=II-1; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
65: if(j<n-1) { J=II+1; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
66: v=4.0; MatSetValues(A,1,&II,1,&II,&v,INSERT_VALUES);
67:
68: }
69: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
70: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
71: }
72: /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
73:
74: if (!loadB) {
75: MatGetLocalSize(A,&m,&n);
76: MatCreate(PETSC_COMM_WORLD,&B);
77: MatSetSizes(B,m,n,PETSC_DECIDE,PETSC_DECIDE);
78: MatSetType(B,MATSBAIJ);
79: MatSetFromOptions(B);
80: MatSetOption(B,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);
81: MatGetOwnershipRange(A,&Istart,&Iend);
82:
83: for( II=Istart; II<Iend; II++ ) {
84: /* v=4.0; MatSetValues(B,1,&II,1,&II,&v,INSERT_VALUES); */
85: v=1.0; MatSetValues(B,1,&II,1,&II,&v,INSERT_VALUES);
86: }
87: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
88: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
89: }
90: /* MatView(B,PETSC_VIEWER_STDOUT_WORLD); */
92: /* Set a shift: A = A - sigma*B */
93: PetscOptionsGetScalar(PETSC_NULL,"-sigma",&sigma,&flag);
94: if (flag){
95: sigma = -1.0 * sigma;
96: MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN); /* A <- A - sigma*B */
97: /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
98: }
100: /* Test MatGetInertia() */
101: KSPCreate(PETSC_COMM_WORLD,&ksp);
102: KSPSetType(ksp,KSPPREONLY);
103: KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
105: KSPGetPC(ksp,&pc);
106: PCSetType(pc,PCCHOLESKY);
107: PCSetFromOptions(pc);
109: PCSetUp(pc);
110: PCFactorGetMatrix(pc,&F);
111: MatGetInertia(F,&nneg,&nzero,&npos);
112: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
113: if (!rank){
114: PetscPrintf(PETSC_COMM_SELF," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);
115: }
117: /* Destroy */
118: KSPDestroy(&ksp);
119: MatDestroy(&A);
120: MatDestroy(&B);
121: PetscFinalize();
122: return 0;
123: }