Actual source code: ex2.c
petsc-3.12.4 2020-02-04
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, test MatValid() --------- */
     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:   PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
    176:   /* Free data structures */
    177:   MatDestroy(&mat);
    178:   MatDestroy(&tmat);
    179:   PetscFinalize();
    180:   return ierr;
    181: }

    187: /*TEST

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

    194:    test:
    195:       suffix: 12_A
    196:       args: -mat_type seqdense -rectA
    197:       filter: grep -v "Mat Object"

    199:    test:
    200:       suffix: 11_B
    201:       args: -mat_type seqaij -rectB
    202:       filter: grep -v "Mat Object"

    204:    test:
    205:       suffix: 12_B
    206:       args: -mat_type seqdense -rectB
    207:       filter: grep -v "Mat Object"

    209:    test:
    210:       suffix: 21
    211:       args: -mat_type mpiaij

    213:    test:
    214:       suffix: 22
    215:       args: -mat_type mpidense

    217:    test:
    218:       suffix: 23
    219:       nsize: 3
    220:       args: -mat_type mpiaij
    221:       filter: grep -v type | grep -v "MPI processes"

    223:    test:
    224:       suffix: 24
    225:       nsize: 3
    226:       args: -mat_type mpidense

    228:    test:
    229:       suffix: 2_aijcusparse_1
    230:       args: -mat_type mpiaijcusparse
    231:       output_file: output/ex2_23.out
    232:       requires: cuda
    233:       filter: grep -v type | grep -v "MPI processes"

    235:    test:
    236:       suffix: 2_aijcusparse_2
    237:       nsize: 3
    238:       args: -mat_type mpiaijcusparse
    239:       output_file: output/ex2_23.out
    240:       requires: cuda
    241:       filter: grep -v type | grep -v "MPI processes"

    243:    test:
    244:       suffix: 3
    245:       nsize: 2
    246:       args: -mat_type mpiaij -rectA

    248:    test:
    249:       suffix: 3_aijcusparse
    250:       nsize: 2
    251:       args: -mat_type mpiaijcusparse -rectA
    252:       requires: cuda

    254:    test:
    255:       suffix: 4
    256:       nsize: 2
    257:       args: -mat_type mpidense -rectA

    259:    test:
    260:       suffix: aijcusparse_1
    261:       args: -mat_type seqaijcusparse -rectA
    262:       filter: grep -v "Mat Object"
    263:       output_file: output/ex2_11_A_aijcusparse.out
    264:       requires: cuda

    266:    test:
    267:       suffix: aijcusparse_2
    268:       args: -mat_type seqaijcusparse -rectB
    269:       filter: grep -v "Mat Object"
    270:       output_file: output/ex2_11_B_aijcusparse.out
    271:       requires: cuda

    273: TEST*/