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: PetscInitialize(&argc,&argv,0,help);
21: PetscRandomCreate(PETSC_COMM_SELF,&rnd);
23: /*
24: * Create matrix and three vectors: these are all normal
25: */
26: PetscMalloc1(size*size,&a);
27: PetscMalloc1(size*size,&b);
28: for (i=0; i<size; i++) {
29: for (j=0; j<size; j++) {
30: PetscRandomGetValue(rnd,&a[i+j*size]);
31: b[i+j*size] = a[i+j*size];
32: }
33: }
34: MatCreate(PETSC_COMM_SELF,&A);
35: MatSetSizes(A,size,size,size,size);
36: MatSetType(A,MATSEQDENSE);
37: MatSeqDenseSetPreallocation(A,a);
38: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
39: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
41: PetscMalloc1(size,&x);
42: for (i=0; i<size; i++) {
43: x[i] = one;
44: }
45: VecCreateSeqWithArray(PETSC_COMM_SELF,1,size,x,&X);
46: VecAssemblyBegin(X);
47: VecAssemblyEnd(X);
49: PetscMalloc1(size,&y);
50: VecCreateSeqWithArray(PETSC_COMM_SELF,1,size,y,&Y);
51: VecAssemblyBegin(Y);
52: VecAssemblyEnd(Y);
54: PetscMalloc1(size,&z);
55: VecCreateSeqWithArray(PETSC_COMM_SELF,1,size,z,&Z);
56: VecAssemblyBegin(Z);
57: VecAssemblyEnd(Z);
59: /*
60: * Now create submatrices and subvectors
61: */
62: MatCreate(PETSC_COMM_SELF,&A11);
63: MatSetSizes(A11,size1,size1,size1,size1);
64: MatSetType(A11,MATSEQDENSE);
65: MatSeqDenseSetPreallocation(A11,b);
66: MatDenseSetLDA(A11,size);
67: MatAssemblyBegin(A11,MAT_FINAL_ASSEMBLY);
68: MatAssemblyEnd(A11,MAT_FINAL_ASSEMBLY);
70: MatCreate(PETSC_COMM_SELF,&A12);
71: MatSetSizes(A12,size1,size2,size1,size2);
72: MatSetType(A12,MATSEQDENSE);
73: MatSeqDenseSetPreallocation(A12,b+size1*size);
74: MatDenseSetLDA(A12,size);
75: MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
76: MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);
78: MatCreate(PETSC_COMM_SELF,&A21);
79: MatSetSizes(A21,size2,size1,size2,size1);
80: MatSetType(A21,MATSEQDENSE);
81: MatSeqDenseSetPreallocation(A21,b+size1);
82: MatDenseSetLDA(A21,size);
83: MatAssemblyBegin(A21,MAT_FINAL_ASSEMBLY);
84: MatAssemblyEnd(A21,MAT_FINAL_ASSEMBLY);
86: MatCreate(PETSC_COMM_SELF,&A22);
87: MatSetSizes(A22,size2,size2,size2,size2);
88: MatSetType(A22,MATSEQDENSE);
89: MatSeqDenseSetPreallocation(A22,b+size1*size+size1);
90: MatDenseSetLDA(A22,size);
91: MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
92: MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);
94: VecCreateSeqWithArray(PETSC_COMM_SELF,1,size1,x,&X1);
95: VecCreateSeqWithArray(PETSC_COMM_SELF,1,size2,x+size1,&X2);
96: VecCreateSeqWithArray(PETSC_COMM_SELF,1,size1,z,&Z1);
97: VecCreateSeqWithArray(PETSC_COMM_SELF,1,size2,z+size1,&Z2);
99: /*
100: * Now multiple matrix times input in two ways;
101: * compare the results
102: */
103: MatMult(A,X,Y);
104: MatMult(A11,X1,Z1);
105: MatMultAdd(A12,X2,Z1,Z1);
106: MatMult(A22,X2,Z2);
107: MatMultAdd(A21,X1,Z2,Z2);
108: VecAXPY(Z,-1.0,Y);
109: VecNorm(Z,NORM_2,&nrm);
110: if (nrm > 100.0*PETSC_MACHINE_EPSILON) {
111: PetscPrintf(PETSC_COMM_WORLD,"Test1; error norm=%g\n",(double)nrm);
112: }
114: /*
115: * Next test: change both matrices
116: */
117: PetscRandomGetValue(rnd,&v);
118: i = 1;
119: j = size-2;
120: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
121: j -= size1;
122: MatSetValues(A12,1,&i,1,&j,&v,INSERT_VALUES);
123: PetscRandomGetValue(rnd,&v);
124: 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: if (nrm > 100.0*PETSC_MACHINE_EPSILON) {
143: PetscPrintf(PETSC_COMM_WORLD,"Test2; error norm=%g\n",(double)nrm);
144: }
146: /*
147: * Transpose product
148: */
149: MatMultTranspose(A,X,Y);
150: MatMultTranspose(A11,X1,Z1);
151: MatMultTransposeAdd(A21,X2,Z1,Z1);
152: MatMultTranspose(A22,X2,Z2);
153: MatMultTransposeAdd(A12,X1,Z2,Z2);
154: VecAXPY(Z,-1.0,Y);
155: VecNorm(Z,NORM_2,&nrm);
156: if (nrm > 100.0*PETSC_MACHINE_EPSILON) {
157: PetscPrintf(PETSC_COMM_WORLD,"Test3; error norm=%g\n",(double)nrm);
158: }
159: PetscFree(a);
160: PetscFree(b);
161: PetscFree(x);
162: PetscFree(y);
163: PetscFree(z);
164: MatDestroy(&A);
165: MatDestroy(&A11);
166: MatDestroy(&A12);
167: MatDestroy(&A21);
168: MatDestroy(&A22);
170: VecDestroy(&X);
171: VecDestroy(&Y);
172: VecDestroy(&Z);
174: VecDestroy(&X1);
175: VecDestroy(&X2);
176: VecDestroy(&Z1);
177: VecDestroy(&Z2);
178: PetscRandomDestroy(&rnd);
179: PetscFinalize();
180: return 0;
181: }
183: /*TEST
185: test:
187: TEST*/