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;
21: PetscInitialize(&argc, &argv, 0, help);
22: PetscRandomCreate(PETSC_COMM_SELF, &rnd);
24: /*
25: * Create matrix and three vectors: these are all normal
26: */
27: PetscMalloc1(size * size, &a);
28: PetscMalloc1(size * size, &b);
29: for (i = 0; i < size; i++) {
30: for (j = 0; j < size; j++) {
31: PetscRandomGetValue(rnd, &a[i + j * size]);
32: b[i + j * size] = a[i + j * size];
33: }
34: }
35: MatCreate(PETSC_COMM_SELF, &A);
36: MatSetSizes(A, size, size, size, size);
37: MatSetType(A, MATSEQDENSE);
38: MatSeqDenseSetPreallocation(A, a);
39: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
40: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
42: PetscMalloc1(size, &x);
43: for (i = 0; i < size; i++) x[i] = one;
44: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size, x, &X);
45: VecAssemblyBegin(X);
46: VecAssemblyEnd(X);
48: PetscMalloc1(size, &y);
49: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size, y, &Y);
50: VecAssemblyBegin(Y);
51: VecAssemblyEnd(Y);
53: PetscMalloc1(size, &z);
54: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size, z, &Z);
55: VecAssemblyBegin(Z);
56: VecAssemblyEnd(Z);
58: /*
59: * Now create submatrices and subvectors
60: */
61: MatCreate(PETSC_COMM_SELF, &A11);
62: MatSetSizes(A11, size1, size1, size1, size1);
63: MatSetType(A11, MATSEQDENSE);
64: MatSeqDenseSetPreallocation(A11, b);
65: MatDenseSetLDA(A11, size);
66: MatAssemblyBegin(A11, MAT_FINAL_ASSEMBLY);
67: MatAssemblyEnd(A11, MAT_FINAL_ASSEMBLY);
69: MatCreate(PETSC_COMM_SELF, &A12);
70: MatSetSizes(A12, size1, size2, size1, size2);
71: MatSetType(A12, MATSEQDENSE);
72: MatSeqDenseSetPreallocation(A12, b + size1 * size);
73: MatDenseSetLDA(A12, size);
74: MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY);
75: MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY);
77: MatCreate(PETSC_COMM_SELF, &A21);
78: MatSetSizes(A21, size2, size1, size2, size1);
79: MatSetType(A21, MATSEQDENSE);
80: MatSeqDenseSetPreallocation(A21, b + size1);
81: MatDenseSetLDA(A21, size);
82: MatAssemblyBegin(A21, MAT_FINAL_ASSEMBLY);
83: MatAssemblyEnd(A21, MAT_FINAL_ASSEMBLY);
85: MatCreate(PETSC_COMM_SELF, &A22);
86: MatSetSizes(A22, size2, size2, size2, size2);
87: MatSetType(A22, MATSEQDENSE);
88: MatSeqDenseSetPreallocation(A22, b + size1 * size + size1);
89: MatDenseSetLDA(A22, size);
90: MatAssemblyBegin(A22, MAT_FINAL_ASSEMBLY);
91: MatAssemblyEnd(A22, MAT_FINAL_ASSEMBLY);
93: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size1, x, &X1);
94: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size2, x + size1, &X2);
95: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size1, z, &Z1);
96: 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: MatMult(A, X, Y);
103: MatMult(A11, X1, Z1);
104: MatMultAdd(A12, X2, Z1, Z1);
105: MatMult(A22, X2, Z2);
106: MatMultAdd(A21, X1, Z2, Z2);
107: VecAXPY(Z, -1.0, Y);
108: VecNorm(Z, NORM_2, &nrm);
109: if (nrm > 100.0 * PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD, "Test1; error norm=%g\n", (double)nrm);
111: /*
112: * Next test: change both matrices
113: */
114: PetscRandomGetValue(rnd, &v);
115: i = 1;
116: j = size - 2;
117: MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES);
118: j -= size1;
119: MatSetValues(A12, 1, &i, 1, &j, &v, INSERT_VALUES);
120: PetscRandomGetValue(rnd, &v);
121: i = j = size1 + 1;
122: MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES);
123: i = j = 1;
124: MatSetValues(A22, 1, &i, 1, &j, &v, INSERT_VALUES);
125: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
126: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
127: MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY);
128: MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY);
129: MatAssemblyBegin(A22, MAT_FINAL_ASSEMBLY);
130: MatAssemblyEnd(A22, MAT_FINAL_ASSEMBLY);
132: MatMult(A, X, Y);
133: MatMult(A11, X1, Z1);
134: MatMultAdd(A12, X2, Z1, Z1);
135: MatMult(A22, X2, Z2);
136: MatMultAdd(A21, X1, Z2, Z2);
137: VecAXPY(Z, -1.0, Y);
138: VecNorm(Z, NORM_2, &nrm);
139: if (nrm > 100.0 * PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD, "Test2; error norm=%g\n", (double)nrm);
141: /*
142: * Transpose product
143: */
144: MatMultTranspose(A, X, Y);
145: MatMultTranspose(A11, X1, Z1);
146: MatMultTransposeAdd(A21, X2, Z1, Z1);
147: MatMultTranspose(A22, X2, Z2);
148: MatMultTransposeAdd(A12, X1, Z2, Z2);
149: VecAXPY(Z, -1.0, Y);
150: VecNorm(Z, NORM_2, &nrm);
151: if (nrm > 100.0 * PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD, "Test3; error norm=%g\n", (double)nrm);
152: PetscFree(a);
153: PetscFree(b);
154: PetscFree(x);
155: PetscFree(y);
156: PetscFree(z);
157: MatDestroy(&A);
158: MatDestroy(&A11);
159: MatDestroy(&A12);
160: MatDestroy(&A21);
161: MatDestroy(&A22);
163: VecDestroy(&X);
164: VecDestroy(&Y);
165: VecDestroy(&Z);
167: VecDestroy(&X1);
168: VecDestroy(&X2);
169: VecDestroy(&Z1);
170: VecDestroy(&Z2);
171: PetscRandomDestroy(&rnd);
172: PetscFinalize();
173: return 0;
174: }
176: /*TEST
178: test:
180: TEST*/