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