Actual source code: ex33.c

petsc-3.5.4 2015-05-23
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_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(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(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(NULL,"-n",&n,NULL);
 51:     PetscOptionsGetInt(NULL,"-m",&m,&flag);
 52:     if (!flag) 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:     MatSetUp(A);

 60:     MatSetOption(A,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);
 61:     MatGetOwnershipRange(A,&Istart,&Iend);
 62:     for (II=Istart; II<Iend; II++) {
 63:       v = -1.0; i = II/n; j = II-i*n;
 64:       if (i>0) { J=II-n; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
 65:       if (i<m-1) { J=II+n; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
 66:       if (j>0) { J=II-1; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
 67:       if (j<n-1) { J=II+1; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
 68:       v=4.0; MatSetValues(A,1,&II,1,&II,&v,INSERT_VALUES);

 70:     }
 71:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 72:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 73:   }
 74:   /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */

 76:   if (!loadB) {
 77:     MatGetLocalSize(A,&m,&n);
 78:     MatCreate(PETSC_COMM_WORLD,&B);
 79:     MatSetSizes(B,m,n,PETSC_DECIDE,PETSC_DECIDE);
 80:     MatSetType(B,MATSBAIJ);
 81:     MatSetFromOptions(B);
 82:     MatSetUp(B);
 83:     MatSetOption(B,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);
 84:     MatGetOwnershipRange(A,&Istart,&Iend);

 86:     for (II=Istart; II<Iend; II++) {
 87:       /* v=4.0; MatSetValues(B,1,&II,1,&II,&v,INSERT_VALUES); */
 88:       v=1.0; MatSetValues(B,1,&II,1,&II,&v,INSERT_VALUES);
 89:     }
 90:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 91:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 92:   }
 93:   /* MatView(B,PETSC_VIEWER_STDOUT_WORLD); */

 95:   /* Set a shift: A = A - sigma*B */
 96:   PetscOptionsGetScalar(NULL,"-sigma",&sigma,&flag);
 97:   if (flag) {
 98:     sigma = -1.0 * sigma;
 99:     MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN); /* A <- A - sigma*B */
100:     /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
101:   }

103:   /* Test MatGetInertia() */
104:   KSPCreate(PETSC_COMM_WORLD,&ksp);
105:   KSPSetType(ksp,KSPPREONLY);
106:   KSPSetOperators(ksp,A,A);

108:   KSPGetPC(ksp,&pc);
109:   PCSetType(pc,PCCHOLESKY);
110:   PCSetFromOptions(pc);

112:   PCSetUp(pc);
113:   PCFactorGetMatrix(pc,&F);
114:   MatGetInertia(F,&nneg,&nzero,&npos);
115:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
116:   if (!rank) {
117:     PetscPrintf(PETSC_COMM_SELF," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);
118:   }

120:   /* Destroy */
121:   KSPDestroy(&ksp);
122:   MatDestroy(&A);
123:   MatDestroy(&B);
124:   PetscFinalize();
125:   return 0;
126: }