Actual source code: ex2.c
petsc-3.9.4 2018-09-11
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 MatInsertValues(), MatMult(),
6: * and MatMultTranspose().
7: */
8: /*T
9: requires: x
10: T*/
12: #include <petscmat.h>
14: int main(int argc,char **argv)
15: {
16: Mat A,A11,A12,A21,A22;
17: Vec X,X1,X2,Y,Z,Z1,Z2;
18: PetscScalar *a,*b,*x,*y,*z,v,one=1;
19: PetscReal nrm;
21: PetscInt size=8,size1=6,size2=2, i,j;
23: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
25: /*
26: * Create matrix and three vectors: these are all normal
27: */
28: PetscMalloc1(size*size,&a);
29: PetscMalloc1(size*size,&b);
30: for (i=0; i<size; i++) {
31: for (j=0; j<size; j++) {
32: a[i+j*size] = rand(); 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: PetscPrintf(PETSC_COMM_WORLD,"Test1; error norm=%g\n",(double)nrm);
112: PetscPrintf(PETSC_COMM_WORLD,"MatMult the usual way:\n");
113: VecView(Y,0);
114: PetscPrintf(PETSC_COMM_WORLD,"MatMult by subblock:\n");
115: VecView(Z,0);
117: /*
118: * Next test: change both matrices
119: */
120: v = rand(); i=1; 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: v = rand(); i=j=size1+1;
125: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
126: i =j=1;
127: MatSetValues(A22,1,&i,1,&j,&v,INSERT_VALUES);
128: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
129: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
130: MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
131: MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);
132: MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
133: MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);
135: MatMult(A,X,Y);
136: MatMult(A11,X1,Z1);
137: MatMultAdd(A12,X2,Z1,Z1);
138: MatMult(A22,X2,Z2);
139: MatMultAdd(A21,X1,Z2,Z2);
140: VecAXPY(Z,-1.0,Y);
141: VecNorm(Z,NORM_2,&nrm);
142: PetscPrintf(PETSC_COMM_WORLD,"Test2; error norm=%g\n",(double)nrm);
144: /*
145: * Transpose product
146: */
147: MatMultTranspose(A,X,Y);
148: MatMultTranspose(A11,X1,Z1);
149: MatMultTransposeAdd(A21,X2,Z1,Z1);
150: MatMultTranspose(A22,X2,Z2);
151: MatMultTransposeAdd(A12,X1,Z2,Z2);
152: VecAXPY(Z,-1.0,Y);
153: VecNorm(Z,NORM_2,&nrm);
154: PetscPrintf(PETSC_COMM_WORLD,"Test3; error norm=%g\n",(double)nrm);
156: PetscFree(a);
157: PetscFree(b);
158: PetscFree(x);
159: PetscFree(y);
160: PetscFree(z);
161: MatDestroy(&A);
162: MatDestroy(&A11);
163: MatDestroy(&A12);
164: MatDestroy(&A21);
165: MatDestroy(&A22);
167: VecDestroy(&X);
168: VecDestroy(&Y);
169: VecDestroy(&Z);
171: VecDestroy(&X1);
172: VecDestroy(&X2);
173: VecDestroy(&Z1);
174: VecDestroy(&Z2);
176: PetscFinalize();
177: return ierr;
178: }
181: /*TEST
183: test:
184: requires: x
185: TODO: Need to implement test
187: TEST*/