Actual source code: ex2.c
1: static char help[] = "testing SeqDense matrices with an LDA (leading dimension of the user-allocated array) larger than M.\n";
2: /*
3: * Example code testing SeqDense matrices with an LDA (leading dimension
4: * of the user-allocated array) larger than M.
5: * This example tests the functionality of MatSetValues(), MatMult(),
6: * and MatMultTranspose().
7: */
9: #include <petscmat.h>
11: int main(int argc, char **argv)
12: {
13: Mat A, A11, A12, A21, A22;
14: Vec X, X1, X2, Y, Z, Z1, Z2;
15: PetscScalar *a, *b, *x, *y, *z, v, one = 1;
16: PetscReal nrm;
17: PetscInt size = 8, size1 = 6, size2 = 2, i, j;
18: PetscRandom rnd;
20: PetscFunctionBeginUser;
21: PetscCall(PetscInitialize(&argc, &argv, 0, help));
22: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rnd));
24: /*
25: * Create matrix and three vectors: these are all normal
26: */
27: PetscCall(PetscMalloc1(size * size, &a));
28: PetscCall(PetscMalloc1(size * size, &b));
29: for (i = 0; i < size; i++) {
30: for (j = 0; j < size; j++) {
31: PetscCall(PetscRandomGetValue(rnd, &a[i + j * size]));
32: b[i + j * size] = a[i + j * size];
33: }
34: }
35: PetscCall(MatCreate(PETSC_COMM_SELF, &A));
36: PetscCall(MatSetSizes(A, size, size, size, size));
37: PetscCall(MatSetType(A, MATSEQDENSE));
38: PetscCall(MatSeqDenseSetPreallocation(A, a));
39: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
40: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
42: PetscCall(PetscMalloc1(size, &x));
43: for (i = 0; i < size; i++) x[i] = one;
44: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size, x, &X));
45: PetscCall(VecAssemblyBegin(X));
46: PetscCall(VecAssemblyEnd(X));
48: PetscCall(PetscMalloc1(size, &y));
49: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size, y, &Y));
50: PetscCall(VecAssemblyBegin(Y));
51: PetscCall(VecAssemblyEnd(Y));
53: PetscCall(PetscMalloc1(size, &z));
54: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size, z, &Z));
55: PetscCall(VecAssemblyBegin(Z));
56: PetscCall(VecAssemblyEnd(Z));
58: /*
59: * Now create submatrices and subvectors
60: */
61: PetscCall(MatCreate(PETSC_COMM_SELF, &A11));
62: PetscCall(MatSetSizes(A11, size1, size1, size1, size1));
63: PetscCall(MatSetType(A11, MATSEQDENSE));
64: PetscCall(MatSeqDenseSetPreallocation(A11, b));
65: PetscCall(MatDenseSetLDA(A11, size));
66: PetscCall(MatAssemblyBegin(A11, MAT_FINAL_ASSEMBLY));
67: PetscCall(MatAssemblyEnd(A11, MAT_FINAL_ASSEMBLY));
69: PetscCall(MatCreate(PETSC_COMM_SELF, &A12));
70: PetscCall(MatSetSizes(A12, size1, size2, size1, size2));
71: PetscCall(MatSetType(A12, MATSEQDENSE));
72: PetscCall(MatSeqDenseSetPreallocation(A12, b + size1 * size));
73: PetscCall(MatDenseSetLDA(A12, size));
74: PetscCall(MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY));
75: PetscCall(MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY));
77: PetscCall(MatCreate(PETSC_COMM_SELF, &A21));
78: PetscCall(MatSetSizes(A21, size2, size1, size2, size1));
79: PetscCall(MatSetType(A21, MATSEQDENSE));
80: PetscCall(MatSeqDenseSetPreallocation(A21, b + size1));
81: PetscCall(MatDenseSetLDA(A21, size));
82: PetscCall(MatAssemblyBegin(A21, MAT_FINAL_ASSEMBLY));
83: PetscCall(MatAssemblyEnd(A21, MAT_FINAL_ASSEMBLY));
85: PetscCall(MatCreate(PETSC_COMM_SELF, &A22));
86: PetscCall(MatSetSizes(A22, size2, size2, size2, size2));
87: PetscCall(MatSetType(A22, MATSEQDENSE));
88: PetscCall(MatSeqDenseSetPreallocation(A22, b + size1 * size + size1));
89: PetscCall(MatDenseSetLDA(A22, size));
90: PetscCall(MatAssemblyBegin(A22, MAT_FINAL_ASSEMBLY));
91: PetscCall(MatAssemblyEnd(A22, MAT_FINAL_ASSEMBLY));
93: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size1, x, &X1));
94: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size2, x + size1, &X2));
95: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size1, z, &Z1));
96: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size2, z + size1, &Z2));
98: /*
99: * Now multiple matrix times input in two ways;
100: * compare the results
101: */
102: PetscCall(MatMult(A, X, Y));
103: PetscCall(MatMult(A11, X1, Z1));
104: PetscCall(MatMultAdd(A12, X2, Z1, Z1));
105: PetscCall(MatMult(A22, X2, Z2));
106: PetscCall(MatMultAdd(A21, X1, Z2, Z2));
107: PetscCall(VecAXPY(Z, -1.0, Y));
108: PetscCall(VecNorm(Z, NORM_2, &nrm));
109: if (nrm > 100.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test1; error norm=%g\n", (double)nrm));
111: /*
112: * Next test: change both matrices
113: */
114: PetscCall(PetscRandomGetValue(rnd, &v));
115: i = 1;
116: j = size - 2;
117: PetscCall(MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES));
118: j -= size1;
119: PetscCall(MatSetValues(A12, 1, &i, 1, &j, &v, INSERT_VALUES));
120: PetscCall(PetscRandomGetValue(rnd, &v));
121: i = j = size1 + 1;
122: PetscCall(MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES));
123: i = j = 1;
124: PetscCall(MatSetValues(A22, 1, &i, 1, &j, &v, INSERT_VALUES));
125: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
126: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
127: PetscCall(MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY));
128: PetscCall(MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY));
129: PetscCall(MatAssemblyBegin(A22, MAT_FINAL_ASSEMBLY));
130: PetscCall(MatAssemblyEnd(A22, MAT_FINAL_ASSEMBLY));
132: PetscCall(MatMult(A, X, Y));
133: PetscCall(MatMult(A11, X1, Z1));
134: PetscCall(MatMultAdd(A12, X2, Z1, Z1));
135: PetscCall(MatMult(A22, X2, Z2));
136: PetscCall(MatMultAdd(A21, X1, Z2, Z2));
137: PetscCall(VecAXPY(Z, -1.0, Y));
138: PetscCall(VecNorm(Z, NORM_2, &nrm));
139: if (nrm > 100.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test2; error norm=%g\n", (double)nrm));
141: /*
142: * Transpose product
143: */
144: PetscCall(MatMultTranspose(A, X, Y));
145: PetscCall(MatMultTranspose(A11, X1, Z1));
146: PetscCall(MatMultTransposeAdd(A21, X2, Z1, Z1));
147: PetscCall(MatMultTranspose(A22, X2, Z2));
148: PetscCall(MatMultTransposeAdd(A12, X1, Z2, Z2));
149: PetscCall(VecAXPY(Z, -1.0, Y));
150: PetscCall(VecNorm(Z, NORM_2, &nrm));
151: if (nrm > 100.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test3; error norm=%g\n", (double)nrm));
152: PetscCall(PetscFree(a));
153: PetscCall(PetscFree(b));
154: PetscCall(PetscFree(x));
155: PetscCall(PetscFree(y));
156: PetscCall(PetscFree(z));
157: PetscCall(MatDestroy(&A));
158: PetscCall(MatDestroy(&A11));
159: PetscCall(MatDestroy(&A12));
160: PetscCall(MatDestroy(&A21));
161: PetscCall(MatDestroy(&A22));
163: PetscCall(VecDestroy(&X));
164: PetscCall(VecDestroy(&Y));
165: PetscCall(VecDestroy(&Z));
167: PetscCall(VecDestroy(&X1));
168: PetscCall(VecDestroy(&X2));
169: PetscCall(VecDestroy(&Z1));
170: PetscCall(VecDestroy(&Z2));
171: PetscCall(PetscRandomDestroy(&rnd));
172: PetscCall(PetscFinalize());
173: return 0;
174: }
176: /*TEST
178: test:
180: TEST*/