Actual source code: ex117.c

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

  6: #include <petscmat.h>

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

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

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

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

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

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

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

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

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

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

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