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: }