Actual source code: ex2.c
petsc-3.12.5 2020-03-29
1: static char help[] = "testing SeqDense matrices with an LDA (leading dimension of the user-allocated arrray) larger than M.\n";
2: /*
3: * Example code testing SeqDense matrices with an LDA (leading dimension
4: * of the user-allocated arrray) 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;
18: PetscInt size=8,size1=6,size2=2, i,j;
19: PetscRandom rnd;
21: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
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(MPI_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++) {
44: x[i] = one;
45: }
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: * Now create submatrices and subvectors
62: */
63: MatCreate(MPI_COMM_SELF,&A11);
64: MatSetSizes(A11,size1,size1,size1,size1);
65: MatSetType(A11,MATSEQDENSE);
66: MatSeqDenseSetPreallocation(A11,b);
67: MatSeqDenseSetLDA(A11,size);
68: MatAssemblyBegin(A11,MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(A11,MAT_FINAL_ASSEMBLY);
71: MatCreate(MPI_COMM_SELF,&A12);
72: MatSetSizes(A12,size1,size2,size1,size2);
73: MatSetType(A12,MATSEQDENSE);
74: MatSeqDenseSetPreallocation(A12,b+size1*size);
75: MatSeqDenseSetLDA(A12,size);
76: MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
77: MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);
79: MatCreate(MPI_COMM_SELF,&A21);
80: MatSetSizes(A21,size2,size1,size2,size1);
81: MatSetType(A21,MATSEQDENSE);
82: MatSeqDenseSetPreallocation(A21,b+size1);
83: MatSeqDenseSetLDA(A21,size);
84: MatAssemblyBegin(A21,MAT_FINAL_ASSEMBLY);
85: MatAssemblyEnd(A21,MAT_FINAL_ASSEMBLY);
87: MatCreate(MPI_COMM_SELF,&A22);
88: MatSetSizes(A22,size2,size2,size2,size2);
89: MatSetType(A22,MATSEQDENSE);
90: MatSeqDenseSetPreallocation(A22,b+size1*size+size1);
91: MatSeqDenseSetLDA(A22,size);
92: MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
93: MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);
95: VecCreateSeqWithArray(MPI_COMM_SELF,1,size1,x,&X1);
96: VecCreateSeqWithArray(MPI_COMM_SELF,1,size2,x+size1,&X2);
97: VecCreateSeqWithArray(MPI_COMM_SELF,1,size1,z,&Z1);
98: VecCreateSeqWithArray(MPI_COMM_SELF,1,size2,z+size1,&Z2);
100: /*
101: * Now multiple matrix times input in two ways;
102: * compare the results
103: */
104: MatMult(A,X,Y);
105: MatMult(A11,X1,Z1);
106: MatMultAdd(A12,X2,Z1,Z1);
107: MatMult(A22,X2,Z2);
108: MatMultAdd(A21,X1,Z2,Z2);
109: VecAXPY(Z,-1.0,Y);
110: VecNorm(Z,NORM_2,&nrm);
111: if (nrm > 100.0*PETSC_MACHINE_EPSILON) {
112: PetscPrintf(PETSC_COMM_WORLD,"Test1; error norm=%g\n",(double)nrm);
113: }
115: /*
116: * Next test: change both matrices
117: */
118: PetscRandomGetValue(rnd,&v);
119: i = 1;
120: j = size-2;
121: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
122: j -= size1;
123: MatSetValues(A12,1,&i,1,&j,&v,INSERT_VALUES);
124: PetscRandomGetValue(rnd,&v);
125: i = j = size1+1;
126: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
127: i =j = 1;
128: MatSetValues(A22,1,&i,1,&j,&v,INSERT_VALUES);
129: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
130: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
131: MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
132: MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);
133: MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
134: MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);
136: MatMult(A,X,Y);
137: MatMult(A11,X1,Z1);
138: MatMultAdd(A12,X2,Z1,Z1);
139: MatMult(A22,X2,Z2);
140: MatMultAdd(A21,X1,Z2,Z2);
141: VecAXPY(Z,-1.0,Y);
142: VecNorm(Z,NORM_2,&nrm);
143: if (nrm > 100.0*PETSC_MACHINE_EPSILON) {
144: PetscPrintf(PETSC_COMM_WORLD,"Test2; error norm=%g\n",(double)nrm);
145: }
147: /*
148: * Transpose product
149: */
150: MatMultTranspose(A,X,Y);
151: MatMultTranspose(A11,X1,Z1);
152: MatMultTransposeAdd(A21,X2,Z1,Z1);
153: MatMultTranspose(A22,X2,Z2);
154: MatMultTransposeAdd(A12,X1,Z2,Z2);
155: VecAXPY(Z,-1.0,Y);
156: VecNorm(Z,NORM_2,&nrm);
157: if (nrm > 100.0*PETSC_MACHINE_EPSILON) {
158: PetscPrintf(PETSC_COMM_WORLD,"Test3; error norm=%g\n",(double)nrm);
159: }
160: PetscFree(a);
161: PetscFree(b);
162: PetscFree(x);
163: PetscFree(y);
164: PetscFree(z);
165: MatDestroy(&A);
166: MatDestroy(&A11);
167: MatDestroy(&A12);
168: MatDestroy(&A21);
169: MatDestroy(&A22);
171: VecDestroy(&X);
172: VecDestroy(&Y);
173: VecDestroy(&Z);
175: VecDestroy(&X1);
176: VecDestroy(&X2);
177: VecDestroy(&Z1);
178: VecDestroy(&Z2);
179: PetscRandomDestroy(&rnd);
180: PetscFinalize();
181: return ierr;
182: }
185: /*TEST
187: test:
189: TEST*/