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*/