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*/