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,mone=-1;
17: PetscReal nrm;
18: int ierr,size=8,size1=6,size2=2, i,j;
20: PetscInitialize(&argc,&argv,0,0);
22: /*
23: * Create matrix and three vectors: these are all normal
24: */
25: PetscMalloc(size*size*sizeof(PetscScalar),&a);
26: PetscMalloc(size*size*sizeof(PetscScalar),&b);
27: for (i=0; i<size; i++) {
28: for (j=0; j<size; j++) {
29: a[i+j*size] = rand(); b[i+j*size] = a[i+j*size];
30: }
31: }
32: MatCreate(MPI_COMM_SELF,size,size,size,size,&A);
33: MatSetType(A,MATSEQDENSE);
34: MatSeqDenseSetPreallocation(A,a);
35: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
36: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
38: PetscMalloc(size*sizeof(PetscScalar),&x);
39: for (i=0; i<size; i++) {
40: x[i] = one;
41: }
42: VecCreateSeqWithArray(MPI_COMM_SELF,size,x,&X);
43: VecAssemblyBegin(X);
44: VecAssemblyEnd(X);
46: PetscMalloc(size*sizeof(PetscScalar),&y);
47: VecCreateSeqWithArray(MPI_COMM_SELF,size,y,&Y);
48: VecAssemblyBegin(Y);
49: VecAssemblyEnd(Y);
51: PetscMalloc(size*sizeof(PetscScalar),&z);
52: VecCreateSeqWithArray(MPI_COMM_SELF,size,z,&Z);
53: VecAssemblyBegin(Z);
54: VecAssemblyEnd(Z);
56: /*
57: * Now create submatrices and subvectors
58: */
59: MatCreate(MPI_COMM_SELF,size1,size1,size1,size1,&A11);
60: MatSetType(A11,MATSEQDENSE);
61: MatSeqDenseSetPreallocation(A11,b);
62: MatSeqDenseSetLDA(A11,size);
63: MatAssemblyBegin(A11,MAT_FINAL_ASSEMBLY);
64: MatAssemblyEnd(A11,MAT_FINAL_ASSEMBLY);
66: MatCreate(MPI_COMM_SELF,size1,size2,size1,size2,&A12);
67: MatSetType(A12,MATSEQDENSE);
68: MatSeqDenseSetPreallocation(A12,b+size1*size);
69: MatSeqDenseSetLDA(A12,size);
70: MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
71: MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);
73: MatCreate(MPI_COMM_SELF,size2,size1,size2,size1,&A21);
74: MatSetType(A21,MATSEQDENSE);
75: MatSeqDenseSetPreallocation(A21,b+size1);
76: MatSeqDenseSetLDA(A21,size);
77: MatAssemblyBegin(A21,MAT_FINAL_ASSEMBLY);
78: MatAssemblyEnd(A21,MAT_FINAL_ASSEMBLY);
80: MatCreate(MPI_COMM_SELF,size2,size2,size2,size2,&A22);
81: MatSetType(A22,MATSEQDENSE);
82: MatSeqDenseSetPreallocation(A22,b+size1*size+size1);
83: MatSeqDenseSetLDA(A22,size);
84: MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
85: MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);
87: VecCreateSeqWithArray(MPI_COMM_SELF,size1,x,&X1);
88: VecCreateSeqWithArray(MPI_COMM_SELF,size2,x+size1,&X2);
89: VecCreateSeqWithArray(MPI_COMM_SELF,size1,z,&Z1);
90: VecCreateSeqWithArray(MPI_COMM_SELF,size2,z+size1,&Z2);
92: /*
93: * Now multiple matrix times input in two ways;
94: * compare the results
95: */
96: MatMult(A,X,Y);
97: MatMult(A11,X1,Z1);
98: MatMultAdd(A12,X2,Z1,Z1);
99: MatMult(A22,X2,Z2);
100: MatMultAdd(A21,X1,Z2,Z2);
101: VecAXPY(&mone,Y,Z);
102: VecNorm(Z,NORM_2,&nrm);
103: printf("Test1; error norm=%e\n",nrm);
104: /*
105: printf("MatMult the usual way:\n"); VecView(Y,0);
106: printf("MatMult by subblock:\n"); VecView(Z,0);
107: */
109: /*
110: * Next test: change both matrices
111: */
112: v = rand(); i=1; j=size-2;
113: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
114: j -= size1;
115: MatSetValues(A12,1,&i,1,&j,&v,INSERT_VALUES);
116: v = rand(); i=j=size1+1;
117: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
118: i=j=1;
119: MatSetValues(A22,1,&i,1,&j,&v,INSERT_VALUES);
120: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
121: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
122: MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
123: MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);
124: MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
125: MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);
127: MatMult(A,X,Y);
128: MatMult(A11,X1,Z1);
129: MatMultAdd(A12,X2,Z1,Z1);
130: MatMult(A22,X2,Z2);
131: MatMultAdd(A21,X1,Z2,Z2);
132: VecAXPY(&mone,Y,Z);
133: VecNorm(Z,NORM_2,&nrm);
134: printf("Test2; error norm=%e\n",nrm);
136: /*
137: * Transpose product
138: */
139: MatMultTranspose(A,X,Y);
140: MatMultTranspose(A11,X1,Z1);
141: MatMultTransposeAdd(A21,X2,Z1,Z1);
142: MatMultTranspose(A22,X2,Z2);
143: MatMultTransposeAdd(A12,X1,Z2,Z2);
144: VecAXPY(&mone,Y,Z);
145: VecNorm(Z,NORM_2,&nrm);
146: printf("Test3; error norm=%e\n",nrm);
148: PetscFinalize();
149: return 0;
150: }