Actual source code: ex2.c

petsc-3.14.0 2020-09-29
Report Typos and Errors

  2: static char help[] = "Tests MatTranspose(), MatNorm(), MatAXPY() and MatAYPX().\n\n";

  4: #include <petscmat.h>

  6: static PetscErrorCode TransposeAXPY(Mat C,PetscScalar alpha,Mat mat,PetscErrorCode (*f)(Mat,Mat*))
  7: {
  8:   Mat            D,E,F,G;

 12:   if (f == MatCreateTranspose) {
 13:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  (C^T)^T = (C^T)^T + alpha * A, C=A, SAME_NONZERO_PATTERN\n");
 14:   } else {
 15:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  (C^H)^H = (C^H)^H + alpha * A, C=A, SAME_NONZERO_PATTERN\n");
 16:   }
 17:   MatDuplicate(mat,MAT_COPY_VALUES,&C);
 18:   f(C,&D);
 19:   f(D,&E);
 20:   MatAXPY(E,alpha,mat,SAME_NONZERO_PATTERN);
 21:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);
 22:   MatDestroy(&E);
 23:   MatDestroy(&D);
 24:   MatDestroy(&C);
 25:   if (f == MatCreateTranspose) {
 26:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  C = C + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n");
 27:   } else {
 28:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  C = C + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n");
 29:   }
 30:   MatDuplicate(mat,MAT_COPY_VALUES,&C);
 31:   /* MATTRANSPOSE should have a MatTranspose_Transpose or MatTranspose_HT implementation */
 32:   if (f == MatCreateTranspose) {
 33:     MatTranspose(mat,MAT_INITIAL_MATRIX,&D);
 34:   } else {
 35:     MatHermitianTranspose(mat,MAT_INITIAL_MATRIX,&D);
 36:   }
 37:   f(D,&E);
 38:   MatAXPY(C,alpha,E,SAME_NONZERO_PATTERN);
 39:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);
 40:   MatDestroy(&E);
 41:   MatDestroy(&D);
 42:   MatDestroy(&C);
 43:   if (f == MatCreateTranspose) {
 44:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  (C^T)^T = (C^T)^T + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n");
 45:   } else {
 46:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  (C^H)^H = (C^H)^H + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n");
 47:   }
 48:   MatDuplicate(mat,MAT_COPY_VALUES,&C);
 49:   f(C,&D);
 50:   f(D,&E);
 51:   f(mat,&F);
 52:   f(F,&G);
 53:   MatAXPY(E,alpha,G,SAME_NONZERO_PATTERN);
 54:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);
 55:   MatDestroy(&G);
 56:   MatDestroy(&F);
 57:   MatDestroy(&E);
 58:   MatDestroy(&D);
 59:   MatDestroy(&C);
 60:   return(0);
 61: }

 63: int main(int argc,char **argv)
 64: {
 65:   Mat            mat,tmat = 0;
 66:   PetscInt       m = 7,n,i,j,rstart,rend,rect = 0;
 68:   PetscMPIInt    size,rank;
 69:   PetscBool      flg;
 70:   PetscScalar    v, alpha;
 71:   PetscReal      normf,normi,norm1;

 73:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 74:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON);
 75:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 76:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 77:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 78:   n    = m;
 79:   PetscOptionsHasName(NULL,NULL,"-rectA",&flg);
 80:   if (flg) {n += 2; rect = 1;}
 81:   PetscOptionsHasName(NULL,NULL,"-rectB",&flg);
 82:   if (flg) {n -= 2; rect = 1;}

 84:   /* ------- Assemble matrix --------- */
 85:   MatCreate(PETSC_COMM_WORLD,&mat);
 86:   MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,m,n);
 87:   MatSetFromOptions(mat);
 88:   MatSetUp(mat);
 89:   MatGetOwnershipRange(mat,&rstart,&rend);
 90:   for (i=rstart; i<rend; i++) {
 91:     for (j=0; j<n; j++) {
 92:       v    = 10.0*i+j;
 93:       MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
 94:     }
 95:   }
 96:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 97:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

 99:   /* ----------------- Test MatNorm()  ----------------- */
100:   MatNorm(mat,NORM_FROBENIUS,&normf);
101:   MatNorm(mat,NORM_1,&norm1);
102:   MatNorm(mat,NORM_INFINITY,&normi);
103:   PetscPrintf(PETSC_COMM_WORLD,"original A: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",(double)normf,(double)norm1,(double)normi);
104:   MatView(mat,PETSC_VIEWER_STDOUT_WORLD);

106:   /* --------------- Test MatTranspose()  -------------- */
107:   PetscOptionsHasName(NULL,NULL,"-in_place",&flg);
108:   if (!rect && flg) {
109:     MatTranspose(mat,MAT_REUSE_MATRIX,&mat);   /* in-place transpose */
110:     tmat = mat; mat = 0;
111:   } else {      /* out-of-place transpose */
112:     MatTranspose(mat,MAT_INITIAL_MATRIX,&tmat);
113:   }

115:   /* ----------------- Test MatNorm()  ----------------- */
116:   /* Print info about transpose matrix */
117:   MatNorm(tmat,NORM_FROBENIUS,&normf);
118:   MatNorm(tmat,NORM_1,&norm1);
119:   MatNorm(tmat,NORM_INFINITY,&normi);
120:   PetscPrintf(PETSC_COMM_WORLD,"B = A^T: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",(double)normf,(double)norm1,(double)normi);
121:   MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);

123:   /* ----------------- Test MatAXPY(), MatAYPX()  ----------------- */
124:   if (mat && !rect) {
125:     alpha = 1.0;
126:     PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL);
127:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  B = B + alpha * A\n");
128:     MatAXPY(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN);
129:     MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);

131:     PetscPrintf(PETSC_COMM_WORLD,"MatAYPX:  B = alpha*B + A\n");
132:     MatAYPX(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN);
133:     MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
134:   }

136:   {
137:     Mat C;
138:     alpha = 1.0;
139:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n");
140:     MatDuplicate(mat,MAT_COPY_VALUES,&C);
141:     MatAXPY(C,alpha,mat,SAME_NONZERO_PATTERN);
142:     MatView(C,PETSC_VIEWER_STDOUT_WORLD);
143:     MatDestroy(&C);
144:     TransposeAXPY(C,alpha,mat,MatCreateTranspose);
145:     TransposeAXPY(C,alpha,mat,MatCreateHermitianTranspose);
146:   }

148:   {
149:     Mat matB;
150:     /* get matB that has nonzeros of mat in all even numbers of row and col */
151:     MatCreate(PETSC_COMM_WORLD,&matB);
152:     MatSetSizes(matB,PETSC_DECIDE,PETSC_DECIDE,m,n);
153:     MatSetFromOptions(matB);
154:     MatSetUp(matB);
155:     MatGetOwnershipRange(matB,&rstart,&rend);
156:     if (rstart % 2 != 0) rstart++;
157:     for (i=rstart; i<rend; i += 2) {
158:       for (j=0; j<n; j += 2) {
159:         v    = 10.0*i+j;
160:         MatSetValues(matB,1,&i,1,&j,&v,INSERT_VALUES);
161:       }
162:     }
163:     MatAssemblyBegin(matB,MAT_FINAL_ASSEMBLY);
164:     MatAssemblyEnd(matB,MAT_FINAL_ASSEMBLY);
165:     PetscPrintf(PETSC_COMM_WORLD," A: original matrix:\n");
166:     MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
167:     PetscPrintf(PETSC_COMM_WORLD," B(a subset of A):\n");
168:     MatView(matB,PETSC_VIEWER_STDOUT_WORLD);
169:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  B = B + alpha * A, SUBSET_NONZERO_PATTERN\n");
170:     MatAXPY(mat,alpha,matB,SUBSET_NONZERO_PATTERN);
171:     MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
172:     MatDestroy(&matB);
173:   }

175:   /* Test MatZeroRows */
176:   j = rstart - 1;
177:   if (j < 0) j = m-1;
178:   MatZeroRows(mat,1,&j,0.0,NULL,NULL);
179:   MatView(mat,PETSC_VIEWER_STDOUT_WORLD);

181:   PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
182:   /* Free data structures */
183:   MatDestroy(&mat);
184:   MatDestroy(&tmat);
185:   PetscFinalize();
186:   return ierr;
187: }

189: /*TEST

191:    test:
192:       suffix: 11_A
193:       args: -mat_type seqaij -rectA
194:       filter: grep -v "Mat Object"

196:    test:
197:       suffix: 12_A
198:       args: -mat_type seqdense -rectA
199:       filter: grep -v type | grep -v "Mat Object"

201:    test:
202:       requires: cuda
203:       suffix: 12_A_cuda
204:       args: -mat_type seqdensecuda -rectA
205:       output_file: output/ex2_12_A.out
206:       filter: grep -v type | grep -v "Mat Object"

208:    test:
209:       suffix: 11_B
210:       args: -mat_type seqaij -rectB
211:       filter: grep -v "Mat Object"

213:    test:
214:       suffix: 12_B
215:       args: -mat_type seqdense -rectB
216:       filter: grep -v type | grep -v "Mat Object"

218:    test:
219:       requires: cuda
220:       suffix: 12_B_cuda
221:       args: -mat_type seqdensecuda -rectB
222:       output_file: output/ex2_12_B.out
223:       filter: grep -v type | grep -v "Mat Object"

225:    test:
226:       suffix: 21
227:       args: -mat_type mpiaij
228:       filter: grep -v type | grep -v "MPI processes"

230:    test:
231:       suffix: 22
232:       args: -mat_type mpidense
233:       filter: grep -v type | grep -v "Mat Object"

235:    test:
236:       requires: cuda
237:       suffix: 22_cuda
238:       output_file: output/ex2_22.out
239:       args: -mat_type mpidensecuda
240:       filter: grep -v type | grep -v "Mat Object"

242:    test:
243:       suffix: 23
244:       nsize: 3
245:       args: -mat_type mpiaij
246:       filter: grep -v type | grep -v "MPI processes"

248:    test:
249:       suffix: 24
250:       nsize: 3
251:       args: -mat_type mpidense
252:       filter: grep -v type | grep -v "Mat Object"

254:    test:
255:       requires: cuda
256:       suffix: 24_cuda
257:       nsize: 3
258:       output_file: output/ex2_24.out
259:       args: -mat_type mpidensecuda
260:       filter: grep -v type | grep -v "Mat Object"

262:    test:
263:       suffix: 2_aijcusparse_1
264:       args: -mat_type mpiaijcusparse
265:       output_file: output/ex2_21.out
266:       requires: cuda
267:       filter: grep -v type | grep -v "MPI processes"

269:    test:
270:       suffix: 2_aijcusparse_2
271:       nsize: 3
272:       args: -mat_type mpiaijcusparse
273:       output_file: output/ex2_23.out
274:       requires: cuda
275:       filter: grep -v type | grep -v "MPI processes"

277:    test:
278:       suffix: 3
279:       nsize: 2
280:       args: -mat_type mpiaij -rectA

282:    test:
283:       suffix: 3_aijcusparse
284:       nsize: 2
285:       args: -mat_type mpiaijcusparse -rectA
286:       requires: cuda

288:    test:
289:       suffix: 4
290:       nsize: 2
291:       args: -mat_type mpidense -rectA
292:       filter: grep -v type | grep -v "MPI processes"

294:    test:
295:       requires: cuda
296:       suffix: 4_cuda
297:       nsize: 2
298:       output_file: output/ex2_4.out
299:       args: -mat_type mpidensecuda -rectA
300:       filter: grep -v type | grep -v "MPI processes"

302:    test:
303:       suffix: aijcusparse_1
304:       args: -mat_type seqaijcusparse -rectA
305:       filter: grep -v "Mat Object"
306:       output_file: output/ex2_11_A_aijcusparse.out
307:       requires: cuda

309:    test:
310:       suffix: aijcusparse_2
311:       args: -mat_type seqaijcusparse -rectB
312:       filter: grep -v "Mat Object"
313:       output_file: output/ex2_11_B_aijcusparse.out
314:       requires: cuda

316: TEST*/