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

      1: static char help[] = "Test MatGetInertia() for Hermitian matrix. \n\n";
      2: /*
      3:   Example of usage
      4:     ./ex36 -check_Hermitian -display_mat -display_vec
      5:     mpiexec -n 2 ./ex36

      7:   This example is modified from src/mat/examples/tests/ex127.c
      8: */

     10: #include <petscksp.h>

     14: PetscInt main(PetscInt argc,char **args)
     15: {
     16:   Mat            A,As;
     17:   PetscBool      flg,disp_mat=PETSC_FALSE;
     19:   PetscMPIInt    size,rank;
     20:   PetscInt       i,j;
     21:   PetscScalar    v,sigma2;
     22:   PetscRandom    rctx;
     23:   PetscReal      h2,sigma1=100.0;
     24:   PetscInt       dim,Ii,J,n = 3,use_random,rstart,rend;
     25:   KSP            ksp;
     26:   PC             pc;
     27:   Mat            F;
     28:   PetscInt       nneg, nzero, npos;

     30:   PetscInitialize(&argc,&args,(char*)0,help);
     31: #if !defined(PETSC_USE_COMPLEX)
     32:   SETERRQ(PETSC_COMM_WORLD,1,"This example requires complex numbers");
     33: #endif
     34:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     35:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
     36:   PetscOptionsHasName(NULL, "-display_mat", &disp_mat);

     38:   PetscOptionsGetReal(NULL,"-sigma1",&sigma1,NULL);
     39:   PetscOptionsGetInt(NULL,"-n",&n,NULL);
     40:   dim  = n*n;

     42:   MatCreate(PETSC_COMM_WORLD,&A);
     43:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
     44:   MatSetType(A,MATAIJ);
     45:   MatSetFromOptions(A);

     47:   PetscOptionsHasName(NULL,"-norandom",&flg);
     48:   if (flg) use_random = 0;
     49:   else use_random = 1;
     50:   if (use_random) {
     51:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
     52:     PetscRandomSetFromOptions(rctx);
     53:     PetscRandomSetInterval(rctx,0.0,PETSC_i);
     54:     PetscRandomGetValue(rctx,&sigma2); /* RealPart(sigma2) == 0.0 */
     55:   } else {
     56:     sigma2 = 10.0*PETSC_i;
     57:   }
     58:   h2 = 1.0/((n+1)*(n+1));

     60:   MatGetOwnershipRange(A,&rstart,&rend);
     61:   for (Ii=rstart; Ii<rend; Ii++) {
     62:     v = -1.0; i = Ii/n; j = Ii - i*n;
     63:     if (i>0) {
     64:       J = Ii-n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
     65:     }
     66:     if (i<n-1) {
     67:       J = Ii+n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
     68:     }
     69:     if (j>0) {
     70:       J = Ii-1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
     71:     }
     72:     if (j<n-1) {
     73:       J = Ii+1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
     74:     }
     75:     v    = 4.0 - sigma1*h2;
     76:     MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
     77:   }
     78:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
     79:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

     81:   /* Check whether A is symmetric */
     82:   PetscOptionsHasName(NULL, "-check_symmetric", &flg);
     83:   if (flg) {
     84:     Mat Trans;
     85:     MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
     86:     MatEqual(A, Trans, &flg);
     87:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"A is not symmetric");
     88:     MatDestroy(&Trans);
     89:   }
     90:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);

     92:   /* make A complex Hermitian */
     93:   Ii = 0; J = dim-1;
     94:   if (Ii >= rstart && Ii < rend) {
     95:     v    = sigma2*h2; /* RealPart(v) = 0.0 */
     96:     MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
     97:     v    = -sigma2*h2;
     98:     MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
     99:   }

    101:   Ii = dim-2; J = dim-1;
    102:   if (Ii >= rstart && Ii < rend) {
    103:     v    = sigma2*h2; /* RealPart(v) = 0.0 */
    104:     MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
    105:     v    = -sigma2*h2;
    106:     MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
    107:   }

    109:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    110:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

    112:   /* Check whether A is Hermitian */
    113:   PetscOptionsHasName(NULL, "-check_Hermitian", &flg);
    114:   if (flg) {
    115:     Mat Hermit;
    116:     if (disp_mat) {
    117:       if (!rank) printf(" A:\n");
    118:       MatView(A,PETSC_VIEWER_STDOUT_WORLD);
    119:     }
    120:     MatHermitianTranspose(A,MAT_INITIAL_MATRIX, &Hermit);
    121:     if (disp_mat) {
    122:       if (!rank) printf(" A_Hermitian:\n");
    123:       MatView(Hermit,PETSC_VIEWER_STDOUT_WORLD);
    124:     }
    125:     MatEqual(A, Hermit, &flg);
    126:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"A is not Hermitian");
    127:     MatDestroy(&Hermit);
    128:   }
    129:   MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);

    131:   /* Create a Hermitian matrix As in sbaij format */
    132:   MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&As);
    133:   if (disp_mat) {
    134:     if (!rank) {PetscPrintf(PETSC_COMM_SELF," As:\n");}
    135:     MatView(As,PETSC_VIEWER_STDOUT_WORLD);
    136:   }

    138:   /* Test MatGetInertia() */
    139:   KSPCreate(PETSC_COMM_WORLD,&ksp);
    140:   KSPSetType(ksp,KSPPREONLY);
    141:   KSPSetOperators(ksp,As,As,DIFFERENT_NONZERO_PATTERN);

    143:   KSPGetPC(ksp,&pc);
    144:   PCSetType(pc,PCCHOLESKY);
    145:   PCSetFromOptions(pc);

    147:   PCSetUp(pc);
    148:   PCFactorGetMatrix(pc,&F);
    149:   MatGetInertia(F,&nneg,&nzero,&npos);
    150:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    151:   if (!rank) {
    152:     PetscPrintf(PETSC_COMM_SELF," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);
    153:   }

    155:   /* Free spaces */
    156:   KSPDestroy(&ksp);
    157:   if (use_random) {PetscRandomDestroy(&rctx);}
    158:   MatDestroy(&A);
    159:   MatDestroy(&As);
    160:   PetscFinalize();
    161:   return 0;
    162: }