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 */
8
9 #include <petscmat.h>
10
main(int argc,char ** argv)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;
19
20 PetscFunctionBeginUser;
21 PetscCall(PetscInitialize(&argc, &argv, 0, help));
22 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rnd));
23
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));
41
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));
47
48 PetscCall(PetscMalloc1(size, &y));
49 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size, y, &Y));
50 PetscCall(VecAssemblyBegin(Y));
51 PetscCall(VecAssemblyEnd(Y));
52
53 PetscCall(PetscMalloc1(size, &z));
54 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size, z, &Z));
55 PetscCall(VecAssemblyBegin(Z));
56 PetscCall(VecAssemblyEnd(Z));
57
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));
68
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));
76
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));
84
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));
92
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));
97
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));
110
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));
131
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));
140
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));
162
163 PetscCall(VecDestroy(&X));
164 PetscCall(VecDestroy(&Y));
165 PetscCall(VecDestroy(&Z));
166
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 }
175
176 /*TEST
177
178 test:
179 output_file: output/empty.out
180
181 TEST*/
182