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;

 17:   PetscInitialize(&argc, &argv, 0, help);
 18:   /* Create matrix and three vectors: these are all normal */
 19:   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:   MatCreate(MPI_COMM_SELF, &A);
 24:   MatSetSizes(A, size, size, size, size);
 25:   MatSetType(A, MATSEQDENSE);
 26:   MatSeqDenseSetPreallocation(A, NULL);

 28:   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:   MatDenseRestoreArray(A, &a);
 33:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 34:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

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

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

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

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

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

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

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

 91: /*TEST

 93:    test:

 95: TEST*/