Actual source code: ex61.c

  1: static char help[] = " * Example code testing SeqDense matrices with an LDA (leading dimension of the user-allocated array) larger than M.\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **argv)
  6: {
  7:   KSP          solver;
  8:   PC           pc;
  9:   Mat          A, B;
 10:   Vec          X, Y, Z;
 11:   MatScalar   *a;
 12:   PetscScalar *b, *x, *y, *z;
 13:   PetscReal    nrm;
 14:   PetscInt     size = 8, lda = 10, i, j;

 16:   PetscFunctionBeginUser;
 17:   PetscCall(PetscInitialize(&argc, &argv, 0, help));
 18:   /* Create matrix and three vectors: these are all normal */
 19:   PetscCall(PetscMalloc1(lda * size, &b));
 20:   for (i = 0; i < size; i++) {
 21:     for (j = 0; j < size; j++) b[i + j * lda] = rand();
 22:   }
 23:   PetscCall(MatCreate(MPI_COMM_SELF, &A));
 24:   PetscCall(MatSetSizes(A, size, size, size, size));
 25:   PetscCall(MatSetType(A, MATSEQDENSE));
 26:   PetscCall(MatSeqDenseSetPreallocation(A, NULL));

 28:   PetscCall(MatDenseGetArray(A, &a));
 29:   for (i = 0; i < size; i++) {
 30:     for (j = 0; j < size; j++) a[i + j * size] = b[i + j * lda];
 31:   }
 32:   PetscCall(MatDenseRestoreArray(A, &a));
 33:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 34:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 36:   PetscCall(MatCreate(MPI_COMM_SELF, &B));
 37:   PetscCall(MatSetSizes(B, size, size, size, size));
 38:   PetscCall(MatSetType(B, MATSEQDENSE));
 39:   PetscCall(MatSeqDenseSetPreallocation(B, b));
 40:   PetscCall(MatDenseSetLDA(B, lda));
 41:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
 42:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

 44:   PetscCall(PetscMalloc1(size, &x));
 45:   for (i = 0; i < size; i++) x[i] = 1.0;
 46:   PetscCall(VecCreateSeqWithArray(MPI_COMM_SELF, 1, size, x, &X));
 47:   PetscCall(VecAssemblyBegin(X));
 48:   PetscCall(VecAssemblyEnd(X));

 50:   PetscCall(PetscMalloc1(size, &y));
 51:   PetscCall(VecCreateSeqWithArray(MPI_COMM_SELF, 1, size, y, &Y));
 52:   PetscCall(VecAssemblyBegin(Y));
 53:   PetscCall(VecAssemblyEnd(Y));

 55:   PetscCall(PetscMalloc1(size, &z));
 56:   PetscCall(VecCreateSeqWithArray(MPI_COMM_SELF, 1, size, z, &Z));
 57:   PetscCall(VecAssemblyBegin(Z));
 58:   PetscCall(VecAssemblyEnd(Z));

 60:   /*
 61:    * Solve with A and B
 62:    */
 63:   PetscCall(KSPCreate(MPI_COMM_SELF, &solver));
 64:   PetscCall(KSPSetType(solver, KSPPREONLY));
 65:   PetscCall(KSPGetPC(solver, &pc));
 66:   PetscCall(PCSetType(pc, PCLU));
 67:   PetscCall(KSPSetOperators(solver, A, A));
 68:   PetscCall(KSPSolve(solver, X, Y));
 69:   PetscCall(KSPSetOperators(solver, B, B));
 70:   PetscCall(KSPSolve(solver, X, Z));
 71:   PetscCall(VecAXPY(Z, -1.0, Y));
 72:   PetscCall(VecNorm(Z, NORM_2, &nrm));
 73:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "Test1; error norm=%e\n", (double)nrm));

 75:   /* Free spaces */
 76:   PetscCall(PetscFree(b));
 77:   PetscCall(PetscFree(x));
 78:   PetscCall(PetscFree(y));
 79:   PetscCall(PetscFree(z));
 80:   PetscCall(VecDestroy(&X));
 81:   PetscCall(VecDestroy(&Y));
 82:   PetscCall(VecDestroy(&Z));
 83:   PetscCall(MatDestroy(&A));
 84:   PetscCall(MatDestroy(&B));
 85:   PetscCall(KSPDestroy(&solver));

 87:   PetscCall(PetscFinalize());
 88:   return 0;
 89: }

 91: /*TEST

 93:    test:

 95: TEST*/