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*/