Actual source code: ex33.c
petsc-3.4.2 2013-07-02

      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,DIFFERENT_NONZERO_PATTERN);

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