Actual source code: ex33.c

petsc-3.12.2 2019-11-22
Report Typos and Errors
  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*/