Actual source code: ex2.c
1: /*
2: * Example code testing SeqDense matrices with an LDA (leading dimension
3: * of the user-allocated arrray) larger than M.
4: * This example tests the functionality of MatInsertValues, MatMult,
5: * and MatMultTranspose.
6: */
7: #include <stdlib.h>
8: #include petscmat.h
12: int main(int argc,char **argv)
13: {
14: Mat A,A11,A12,A21,A22;
15: Vec X,X1,X2,Y,Z,Z1,Z2;
16: PetscScalar *a,*b,*x,*y,*z,v,one=1;
17: PetscReal nrm;
19: PetscInt size=8,size1=6,size2=2, i,j;
21: PetscInitialize(&argc,&argv,0,0);
23: /*
24: * Create matrix and three vectors: these are all normal
25: */
26: PetscMalloc(size*size*sizeof(PetscScalar),&a);
27: PetscMalloc(size*size*sizeof(PetscScalar),&b);
28: for (i=0; i<size; i++) {
29: for (j=0; j<size; j++) {
30: a[i+j*size] = rand(); b[i+j*size] = a[i+j*size];
31: }
32: }
33: MatCreate(MPI_COMM_SELF,&A);
34: MatSetSizes(A,size,size,size,size);
35: MatSetType(A,MATSEQDENSE);
36: MatSeqDenseSetPreallocation(A,a);
37: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
38: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
40: PetscMalloc(size*sizeof(PetscScalar),&x);
41: for (i=0; i<size; i++) {
42: x[i] = one;
43: }
44: VecCreateSeqWithArray(MPI_COMM_SELF,size,x,&X);
45: VecAssemblyBegin(X);
46: VecAssemblyEnd(X);
48: PetscMalloc(size*sizeof(PetscScalar),&y);
49: VecCreateSeqWithArray(MPI_COMM_SELF,size,y,&Y);
50: VecAssemblyBegin(Y);
51: VecAssemblyEnd(Y);
53: PetscMalloc(size*sizeof(PetscScalar),&z);
54: VecCreateSeqWithArray(MPI_COMM_SELF,size,z,&Z);
55: VecAssemblyBegin(Z);
56: VecAssemblyEnd(Z);
58: /*
59: * Now create submatrices and subvectors
60: */
61: MatCreate(MPI_COMM_SELF,&A11);
62: MatSetSizes(A11,size1,size1,size1,size1);
63: MatSetType(A11,MATSEQDENSE);
64: MatSeqDenseSetPreallocation(A11,b);
65: MatSeqDenseSetLDA(A11,size);
66: MatAssemblyBegin(A11,MAT_FINAL_ASSEMBLY);
67: MatAssemblyEnd(A11,MAT_FINAL_ASSEMBLY);
69: MatCreate(MPI_COMM_SELF,&A12);
70: MatSetSizes(A12,size1,size2,size1,size2);
71: MatSetType(A12,MATSEQDENSE);
72: MatSeqDenseSetPreallocation(A12,b+size1*size);
73: MatSeqDenseSetLDA(A12,size);
74: MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
75: MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);
77: MatCreate(MPI_COMM_SELF,&A21);
78: MatSetSizes(A21,size2,size1,size2,size1);
79: MatSetType(A21,MATSEQDENSE);
80: MatSeqDenseSetPreallocation(A21,b+size1);
81: MatSeqDenseSetLDA(A21,size);
82: MatAssemblyBegin(A21,MAT_FINAL_ASSEMBLY);
83: MatAssemblyEnd(A21,MAT_FINAL_ASSEMBLY);
85: MatCreate(MPI_COMM_SELF,&A22);
86: MatSetSizes(A22,size2,size2,size2,size2);
87: MatSetType(A22,MATSEQDENSE);
88: MatSeqDenseSetPreallocation(A22,b+size1*size+size1);
89: MatSeqDenseSetLDA(A22,size);
90: MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
91: MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);
93: VecCreateSeqWithArray(MPI_COMM_SELF,size1,x,&X1);
94: VecCreateSeqWithArray(MPI_COMM_SELF,size2,x+size1,&X2);
95: VecCreateSeqWithArray(MPI_COMM_SELF,size1,z,&Z1);
96: VecCreateSeqWithArray(MPI_COMM_SELF,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: printf("Test1; error norm=%e\n",nrm);
110: /*
111: printf("MatMult the usual way:\n"); VecView(Y,0);
112: printf("MatMult by subblock:\n"); VecView(Z,0);
113: */
115: /*
116: * Next test: change both matrices
117: */
118: v = rand(); i=1; j=size-2;
119: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
120: j -= size1;
121: MatSetValues(A12,1,&i,1,&j,&v,INSERT_VALUES);
122: v = rand(); i=j=size1+1;
123: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
124: i=j=1;
125: MatSetValues(A22,1,&i,1,&j,&v,INSERT_VALUES);
126: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
127: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
128: MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
129: MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);
130: MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
131: MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);
133: MatMult(A,X,Y);
134: MatMult(A11,X1,Z1);
135: MatMultAdd(A12,X2,Z1,Z1);
136: MatMult(A22,X2,Z2);
137: MatMultAdd(A21,X1,Z2,Z2);
138: VecAXPY(Z,-1.0,Y);
139: VecNorm(Z,NORM_2,&nrm);
140: printf("Test2; error norm=%e\n",nrm);
142: /*
143: * Transpose product
144: */
145: MatMultTranspose(A,X,Y);
146: MatMultTranspose(A11,X1,Z1);
147: MatMultTransposeAdd(A21,X2,Z1,Z1);
148: MatMultTranspose(A22,X2,Z2);
149: MatMultTransposeAdd(A12,X1,Z2,Z2);
150: VecAXPY(Z,-1.0,Y);
151: VecNorm(Z,NORM_2,&nrm);
152: printf("Test3; error norm=%e\n",nrm);
154: PetscFree(a);
155: PetscFree(b);
156: PetscFree(x);
157: PetscFree(y);
158: PetscFree(z);
159: MatDestroy(A);
160: MatDestroy(A11);
161: MatDestroy(A12);
162: MatDestroy(A21);
163: MatDestroy(A22);
165: VecDestroy(X);
166: VecDestroy(Y);
167: VecDestroy(Z);
169: VecDestroy(X1);
170: VecDestroy(X2);
171: VecDestroy(Z1);
172: VecDestroy(Z2);
175: PetscFinalize();
176: return 0;
177: }