Actual source code: ex123.c
1: static char help[] = "Test MatSetPreallocationCOO and MatSetValuesCOO\n\n";
3: #include <petscmat.h>
4: #define MyMatView(a,b) PetscPrintf(PetscObjectComm((PetscObject)(a)),"LINE %d\n",__LINE__),MatView(a,b);
5: #define MyVecView(a,b) PetscPrintf(PetscObjectComm((PetscObject)(a)),"LINE %d\n",__LINE__),VecView(a,b);
6: int main(int argc,char **args)
7: {
8: Mat A,At,AAt;
9: Vec x,y,z;
10: PetscLayout rmap,cmap;
11: PetscInt n1 = 11, n2 = 9;
12: PetscInt i1[] = { 7, 6, 2, 0, 4, 1, 1, 0, 2, 2, 1 };
13: PetscInt j1[] = { 1, 4, 3, 5, 3, 3, 4, 5, 0, 3, 1 };
14: PetscInt i2[] = { 7, 6, 2, 0, 4, 1, 1, 2, 1 };
15: PetscInt j2[] = { 1, 4, 3, 5, 3, 3, 4, 0, 1 };
16: PetscScalar v1[] = { -1., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.};
17: PetscScalar v2[] = { 1.,-1.,-2.,-3.,-4.,-5.,-6.,-7.,-8.,-9.,-10.};
18: PetscInt N = 6, m = 8, rstart, cstart, i;
19: PetscMPIInt size;
20: PetscBool loc = PETSC_FALSE;
21: PetscBool locdiag = PETSC_TRUE, ismpiaij;
24: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
25: PetscOptionsGetBool(NULL,NULL,"-loc",&loc,NULL);
26: PetscOptionsGetBool(NULL,NULL,"-locdiag",&locdiag,NULL);
28: MatCreate(PETSC_COMM_WORLD,&A);
29: if (loc) {
30: if (locdiag) {
31: MatSetSizes(A,m,N,PETSC_DECIDE,PETSC_DECIDE);
32: } else {
33: MatSetSizes(A,m,m+N,PETSC_DECIDE,PETSC_DECIDE);
34: }
35: } else {
36: MatSetSizes(A,m,PETSC_DECIDE,PETSC_DECIDE,N);
37: }
38: MatSetFromOptions(A);
39: MatGetLayouts(A,&rmap,&cmap);
40: PetscLayoutSetUp(rmap);
41: PetscLayoutSetUp(cmap);
42: MatCreateVecs(A,&x,&y);
43: MatCreateVecs(A,NULL,&z);
44: VecSet(x,1.);
45: VecSet(z,2.);
46: PetscLayoutGetRange(rmap,&rstart,NULL);
47: PetscLayoutGetRange(cmap,&cstart,NULL);
48: for (i = 0; i < n1; i++) i1[i] += rstart;
49: for (i = 0; i < n2; i++) i2[i] += rstart;
50: if (loc) {
51: if (locdiag) {
52: for (i = 0; i < n1; i++) j1[i] += cstart;
53: for (i = 0; i < n2; i++) j2[i] += cstart;
54: } else {
55: for (i = 0; i < n1; i++) j1[i] += cstart + m;
56: for (i = 0; i < n2; i++) j2[i] += cstart + m;
57: }
58: }
60: /* test with repeated entries */
61: MatSetPreallocationCOO(A,n1,i1,j1);
62: MatSetValuesCOO(A,v1,ADD_VALUES);
63: MyMatView(A,NULL);
64: MatMult(A,x,y);
65: MyVecView(y,NULL);
66: MatSetValuesCOO(A,v2,ADD_VALUES);
67: MyMatView(A,NULL);
68: MatMultAdd(A,x,y,y);
69: MyVecView(y,NULL);
70: MatTranspose(A,MAT_INITIAL_MATRIX,&At);
71: MatMatMult(A,At,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);
72: MyMatView(AAt,NULL);
73: MatDestroy(&AAt);
74: MatMatMult(At,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);
75: MyMatView(AAt,NULL);
76: MatDestroy(&AAt);
77: MatDestroy(&At);
78: /* INSERT_VALUES will overwrite matrix entries but
79: still perform the sum of the repeated entries */
80: MatSetValuesCOO(A,v2,INSERT_VALUES);
81: MyMatView(A,NULL);
83: /* test with unique entries */
84: MatSetPreallocationCOO(A,n2,i2,j2);
85: MatSetValuesCOO(A,v1,ADD_VALUES);
86: MyMatView(A,NULL);
87: MatMult(A,x,y);
88: MyVecView(y,NULL);
89: MatSetValuesCOO(A,v2,ADD_VALUES);
90: MyMatView(A,NULL);
91: MatMultAdd(A,x,y,z);
92: MyVecView(z,NULL);
93: MatSetPreallocationCOO(A,n2,i2,j2);
94: MatSetValuesCOO(A,v1,INSERT_VALUES);
95: MyMatView(A,NULL);
96: MatMult(A,x,y);
97: MyVecView(y,NULL);
98: MatSetValuesCOO(A,v2,INSERT_VALUES);
99: MyMatView(A,NULL);
100: MatMultAdd(A,x,y,z);
101: MyVecView(z,NULL);
102: MatTranspose(A,MAT_INITIAL_MATRIX,&At);
103: MatMatMult(A,At,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);
104: MyMatView(AAt,NULL);
105: MatDestroy(&AAt);
106: MatMatMult(At,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);
107: MyMatView(AAt,NULL);
108: MatDestroy(&AAt);
109: MatDestroy(&At);
111: /* test providing diagonal first, the offdiagonal */
112: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
113: PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);
114: if (ismpiaij && size > 1) {
115: Mat lA,lB;
116: const PetscInt *garray,*iA,*jA,*iB,*jB;
117: const PetscScalar *vA,*vB;
118: PetscScalar *coo_v;
119: PetscInt *coo_i,*coo_j;
120: PetscInt i,j,nA,nB,nnz;
121: PetscBool flg;
123: MatMPIAIJGetSeqAIJ(A,&lA,&lB,&garray);
124: MatSeqAIJGetArrayRead(lA,&vA);
125: MatSeqAIJGetArrayRead(lB,&vB);
126: MatGetRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nA,&iA,&jA,&flg);
127: MatGetRowIJ(lB,0,PETSC_FALSE,PETSC_FALSE,&nB,&iB,&jB,&flg);
128: nnz = iA[nA] + iB[nB];
129: PetscMalloc3(nnz,&coo_i,nnz,&coo_j,nnz,&coo_v);
130: nnz = 0;
131: for (i=0;i<nA;i++) {
132: for (j=iA[i];j<iA[i+1];j++,nnz++) {
133: coo_i[nnz] = i+rstart;
134: coo_j[nnz] = jA[j]+cstart;
135: coo_v[nnz] = vA[j];
136: }
137: }
138: for (i=0;i<nB;i++) {
139: for (j=iB[i];j<iB[i+1];j++,nnz++) {
140: coo_i[nnz] = i+rstart;
141: coo_j[nnz] = garray[jB[j]];
142: coo_v[nnz] = vB[j];
143: }
144: }
145: MatRestoreRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nA,&iA,&jA,&flg);
146: MatRestoreRowIJ(lB,0,PETSC_FALSE,PETSC_FALSE,&nB,&iB,&jB,&flg);
147: MatSeqAIJRestoreArrayRead(lA,&vA);
148: MatSeqAIJRestoreArrayRead(lB,&vB);
150: MatSetPreallocationCOO(A,nnz,coo_i,coo_j);
151: MatSetValuesCOO(A,coo_v,ADD_VALUES);
152: MyMatView(A,NULL);
153: MatMult(A,x,y);
154: MyVecView(y,NULL);
155: MatSetValuesCOO(A,coo_v,INSERT_VALUES);
156: MyMatView(A,NULL);
157: MatMult(A,x,y);
158: MyVecView(y,NULL);
159: MatTranspose(A,MAT_INITIAL_MATRIX,&At);
160: MatMatMult(A,At,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);
161: MyMatView(AAt,NULL);
162: MatDestroy(&AAt);
163: MatMatMult(At,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);
164: MyMatView(AAt,NULL);
165: MatDestroy(&AAt);
166: MatDestroy(&At);
168: PetscFree3(coo_i,coo_j,coo_v);
169: }
170: VecDestroy(&z);
171: VecDestroy(&x);
172: VecDestroy(&y);
173: MatDestroy(&A);
174: PetscFinalize();
175: return ierr;
176: }
178: /*TEST
180: test:
181: suffix: 1
182: filter: grep -v type
183: diff_args: -j
184: args: -mat_type {{seqaij mpiaij}}
186: test:
187: requires: cuda
188: suffix: 1_cuda
189: filter: grep -v type
190: diff_args: -j
191: args: -mat_type {{seqaijcusparse mpiaijcusparse}}
192: output_file: output/ex123_1.out
194: test:
195: suffix: 2
196: nsize: 7
197: filter: grep -v type
198: diff_args: -j
199: args: -mat_type mpiaij
201: test:
202: requires: cuda
203: suffix: 2_cuda
204: nsize: 7
205: filter: grep -v type
206: diff_args: -j
207: args: -mat_type mpiaijcusparse
208: output_file: output/ex123_2.out
210: test:
211: suffix: 3
212: nsize: 3
213: filter: grep -v type
214: diff_args: -j
215: args: -mat_type mpiaij -loc
217: test:
218: requires: cuda
219: suffix: 3_cuda
220: nsize: 3
221: filter: grep -v type
222: diff_args: -j
223: args: -mat_type mpiaijcusparse -loc
224: output_file: output/ex123_3.out
226: test:
227: suffix: 4
228: nsize: 4
229: filter: grep -v type
230: diff_args: -j
231: args: -mat_type mpiaij -loc -locdiag 0
233: test:
234: requires: cuda
235: suffix: 4_cuda
236: nsize: 4
237: filter: grep -v type
238: diff_args: -j
239: args: -mat_type mpiaijcusparse -loc -locdiag 0
240: output_file: output/ex123_4.out
242: TEST*/