1 static char help[] = " * Example code testing SeqDense matrices with an LDA (leading dimension of the user-allocated array) larger than M.\n";
2
3 #include <petscksp.h>
4
main(int argc,char ** argv)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;
15
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));
27
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));
35
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));
43
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));
49
50 PetscCall(PetscMalloc1(size, &y));
51 PetscCall(VecCreateSeqWithArray(MPI_COMM_SELF, 1, size, y, &Y));
52 PetscCall(VecAssemblyBegin(Y));
53 PetscCall(VecAssemblyEnd(Y));
54
55 PetscCall(PetscMalloc1(size, &z));
56 PetscCall(VecCreateSeqWithArray(MPI_COMM_SELF, 1, size, z, &Z));
57 PetscCall(VecAssemblyBegin(Z));
58 PetscCall(VecAssemblyEnd(Z));
59
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));
74
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));
86
87 PetscCall(PetscFinalize());
88 return 0;
89 }
90
91 /*TEST
92
93 test:
94
95 TEST*/
96