Actual source code: ex117.c


  2: static char help[] = "Tests Cholesky factorization for a SBAIJ matrix, (bs=2).\n";
  3: /*
  4:   This code is modified from the code contributed by JUNWANG@uwm.edu on Apr 13, 2007
  5: */

  7: #include <petscmat.h>

  9: int main(int argc, char **args)
 10: {
 11:   Mat           mat, fact, B;
 12:   PetscInt      ind1[2], ind2[2];
 13:   PetscScalar   temp[4];
 14:   PetscInt      nnz[3];
 15:   IS            perm, colp;
 16:   MatFactorInfo info;
 17:   PetscMPIInt   size;

 19:   PetscFunctionBeginUser;
 20:   PetscCall(PetscInitialize(&argc, &args, 0, help));
 21:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 22:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

 24:   nnz[0] = 2;
 25:   nnz[1] = 1;
 26:   nnz[2] = 1;
 27:   PetscCall(MatCreateSeqSBAIJ(PETSC_COMM_SELF, 2, 6, 6, 0, nnz, &mat));

 29:   ind1[0] = 0;
 30:   ind1[1] = 1;
 31:   temp[0] = 3;
 32:   temp[1] = 2;
 33:   temp[2] = 0;
 34:   temp[3] = 3;
 35:   PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES));
 36:   ind2[0] = 4;
 37:   ind2[1] = 5;
 38:   temp[0] = 1;
 39:   temp[1] = 1;
 40:   temp[2] = 2;
 41:   temp[3] = 1;
 42:   PetscCall(MatSetValues(mat, 2, ind1, 2, ind2, temp, INSERT_VALUES));
 43:   ind1[0] = 2;
 44:   ind1[1] = 3;
 45:   temp[0] = 4;
 46:   temp[1] = 1;
 47:   temp[2] = 1;
 48:   temp[3] = 5;
 49:   PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES));
 50:   ind1[0] = 4;
 51:   ind1[1] = 5;
 52:   temp[0] = 5;
 53:   temp[1] = 1;
 54:   temp[2] = 1;
 55:   temp[3] = 6;
 56:   PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES));

 58:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
 59:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));

 61:   PetscCall(MatDuplicate(mat, MAT_SHARE_NONZERO_PATTERN, &B));
 62:   ind1[0] = 0;
 63:   ind1[1] = 1;
 64:   temp[0] = 3;
 65:   temp[1] = 2;
 66:   temp[2] = 0;
 67:   temp[3] = 3;
 68:   PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES));
 69:   ind2[0] = 4;
 70:   ind2[1] = 5;
 71:   temp[0] = 1;
 72:   temp[1] = 1;
 73:   temp[2] = 2;
 74:   temp[3] = 1;
 75:   PetscCall(MatSetValues(mat, 2, ind1, 2, ind2, temp, INSERT_VALUES));
 76:   ind1[0] = 2;
 77:   ind1[1] = 3;
 78:   temp[0] = 4;
 79:   temp[1] = 1;
 80:   temp[2] = 1;
 81:   temp[3] = 5;
 82:   PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES));
 83:   ind1[0] = 4;
 84:   ind1[1] = 5;
 85:   temp[0] = 5;
 86:   temp[1] = 1;
 87:   temp[2] = 1;
 88:   temp[3] = 6;
 89:   PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES));

 91:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
 92:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));

 94:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mat: \n"));
 95:   PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_SELF));

 97:   /* begin cholesky factorization */
 98:   PetscCall(MatGetOrdering(mat, MATORDERINGNATURAL, &perm, &colp));
 99:   PetscCall(ISDestroy(&colp));

101:   info.fill = 1.0;
102:   PetscCall(MatGetFactor(mat, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &fact));
103:   PetscCall(MatCholeskyFactorSymbolic(fact, mat, perm, &info));
104:   PetscCall(MatCholeskyFactorNumeric(fact, mat, &info));

106:   PetscCall(ISDestroy(&perm));
107:   PetscCall(MatDestroy(&mat));
108:   PetscCall(MatDestroy(&fact));
109:   PetscCall(MatDestroy(&B));
110:   PetscCall(PetscFinalize());
111:   return 0;
112: }