Actual source code: ex93.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      1: static char help[] = "Test MatMatMult() and MatPtAP() for AIJ matrices.\n\n";

      3:  #include <petscmat.h>

      5: extern PetscErrorCode testPTAPRectangular(void);

      7: int main(int argc,char **argv)
      8: {
      9:   Mat            A,B,C,D;
     10:   PetscScalar    a[] ={1.,1.,0.,0.,1.,1.,0.,0.,1.};
     11:   PetscInt       ij[]={0,1,2};
     12:   PetscScalar    none=-1.;
     14:   PetscReal      fill=4;
     15:   PetscReal      norm;
     16:   PetscMPIInt    size,rank;
     17:   PetscBool      test_hypre=PETSC_FALSE;

     19:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     20: #if defined(PETSC_HAVE_HYPRE)
     21:   PetscOptionsGetBool(NULL,NULL,"-test_hypre",&test_hypre,NULL);
     22: #endif
     23:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     24:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

     26:   MatCreate(PETSC_COMM_WORLD,&A);
     27:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);
     28:   MatSetType(A,MATAIJ);
     29:   MatSetFromOptions(A);
     30:   MatSetUp(A);
     31:   MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);

     33:   if (!rank) {
     34:     MatSetValues(A,3,ij,3,ij,a,ADD_VALUES);
     35:   }
     36:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
     37:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
     38:   MatSetOptionsPrefix(A,"A_");
     39:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);
     40:   PetscPrintf(PETSC_COMM_WORLD,"\n");

     42:   /* Test MatMatMult() */
     43:   MatTranspose(A,MAT_INITIAL_MATRIX,&B);      /* B = A^T */
     44:   MatSetOptionsPrefix(B,"B_");
     45:   MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A */
     46:   MatMatMultNumeric(B,A,C);                   /* recompute C=B*A */
     47:   MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C);   /* recompute C=B*A */
     48:   MatSetOptionsPrefix(C,"C_");
     49:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);
     50:   if (!rank) {PetscPrintf(PETSC_COMM_SELF,"\n");}

     52:   MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&D);
     53:   MatMatMultNumeric(C,A,D);  /* D = C*A = (A^T*A)*A */
     54:   MatSetOptionsPrefix(D,"D_");
     55:   MatView(D,PETSC_VIEWER_STDOUT_WORLD);
     56:   if (!rank) {PetscPrintf(PETSC_COMM_SELF,"\n");}

     58:   /* Repeat the numeric product to test reuse of the previous symbolic product */
     59:   MatMatMultNumeric(C,A,D);
     60:   MatView(D,PETSC_VIEWER_STDOUT_WORLD);
     61:   if (!rank) {PetscPrintf(PETSC_COMM_SELF,"\n");}

     63:   MatDestroy(&B);
     64:   MatDestroy(&C);

     66:   /* Test PtAP routine. */
     67:   MatDuplicate(A,MAT_COPY_VALUES,&B);      /* B = A */
     68:   MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C); /* C = B^T*A*B */
     69:   if (!test_hypre) {
     70:     MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
     71:     MatNorm(D,NORM_FROBENIUS,&norm);
     72:     if (norm > 1.e-15 && !rank) {
     73:       PetscPrintf(PETSC_COMM_SELF,"Error in MatPtAP: %g\n",norm);
     74:     }
     75:   }
     76:   MatDestroy(&C);
     77:   MatDestroy(&D);

     79:   /* Repeat PtAP to test symbolic/numeric separation for reuse of the symbolic product */
     80:   MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C);
     81:   MatSetOptionsPrefix(C,"C=BtAB_");
     82:   if (!test_hypre) {
     83:     MatView(C,PETSC_VIEWER_STDOUT_WORLD);
     84:     PetscPrintf(PETSC_COMM_WORLD,"\n");
     85:   }

     87:   if (!test_hypre) {
     88:     MatPtAPSymbolic(A,B,fill,&D);
     89:     MatPtAPNumeric(A,B,D);
     90:     MatSetOptionsPrefix(D,"D=BtAB_");
     91:     MatView(D,PETSC_VIEWER_STDOUT_WORLD);
     92:     PetscPrintf(PETSC_COMM_WORLD,"\n");

     94:     /* Repeat numeric product to test reuse of the previous symbolic product */
     95:     MatPtAPNumeric(A,B,D);
     96:     MatAXPY(D,none,C,DIFFERENT_NONZERO_PATTERN);
     97:     MatNorm(D,NORM_FROBENIUS,&norm);
     98:     if (norm > 1.e-15) {
     99:       PetscPrintf(PETSC_COMM_WORLD,"Error in symbolic/numeric MatPtAP: %g\n",norm);
    100:     }
    101:     MatDestroy(&D);
    102:   }
    103:   MatDestroy(&C);
    104:   MatDestroy(&B);

    106:   if (size == 1) {
    107:     /* A test contributed by Tobias Neckel <neckel@in.tum.de> */
    108:     testPTAPRectangular();

    110:     /* test MatMatTransposeMult(): A*B^T */
    111:     MatMatTransposeMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A*A^T */
    112:     MatSetOptionsPrefix(D,"D=A*A^T_");
    113:     MatView(D,PETSC_VIEWER_STDOUT_WORLD);
    114:     PetscPrintf(PETSC_COMM_WORLD,"\n");

    116:     MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
    117:     MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C); /* C=A*B */
    118:     MatSetOptionsPrefix(C,"D=A*B=A*A^T_");
    119:     MatView(C,PETSC_VIEWER_STDOUT_WORLD);
    120:     PetscPrintf(PETSC_COMM_WORLD,"\n");
    121:   }

    123:   MatDestroy(&A);
    124:   MatDestroy(&B);
    125:   MatDestroy(&C);
    126:   MatDestroy(&D);
    127:   PetscFinalize();
    128:   return ierr;
    129: }

    131: /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */
    132: PetscErrorCode testPTAPRectangular(void)
    133: {

    135:   const int      rows = 3;
    136:   const int      cols = 5;
    138:   int            i;
    139:   Mat            A;
    140:   Mat            P;
    141:   Mat            C;

    144:   /* set up A  */
    145:   MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows,1, NULL, &A);
    146:   for (i=0; i<rows; i++) {
    147:     MatSetValue(A, i, i, 1.0, INSERT_VALUES);
    148:   }
    149:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    150:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

    152:   /* set up P */
    153:   MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols,5, NULL, &P);
    154:   MatSetValue(P, 0, 0,  1.0, INSERT_VALUES);
    155:   MatSetValue(P, 0, 1,  2.0, INSERT_VALUES);
    156:   MatSetValue(P, 0, 2,  0.0, INSERT_VALUES);

    158:   MatSetValue(P, 0, 3, -1.0, INSERT_VALUES);

    160:   MatSetValue(P, 1, 0,  0.0, INSERT_VALUES);
    161:   MatSetValue(P, 1, 1, -1.0, INSERT_VALUES);
    162:   MatSetValue(P, 1, 2,  1.0, INSERT_VALUES);

    164:   MatSetValue(P, 2, 0,  3.0, INSERT_VALUES);
    165:   MatSetValue(P, 2, 1,  0.0, INSERT_VALUES);
    166:   MatSetValue(P, 2, 2, -3.0, INSERT_VALUES);

    168:   MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
    169:   MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);

    171:   /* compute C */
    172:   MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C);

    174:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
    175:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

    177:   /* compare results */
    178:   /*
    179:   printf("C:\n");
    180:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);

    182:   blitz::Array<double,2> actualC(cols, cols);
    183:   actualC = 0.0;
    184:   for (int i=0; i<cols; i++) {
    185:     for (int j=0; j<cols; j++) {
    186:       MatGetValues(C, 1, &i, 1, &j, &actualC(i,j));
    187:        ;
    188:     }
    189:   }
    190:   blitz::Array<double,2> expectedC(cols, cols);
    191:   expectedC = 0.0;

    193:   expectedC(0,0) = 10.0;
    194:   expectedC(0,1) =  2.0;
    195:   expectedC(0,2) = -9.0;
    196:   expectedC(0,3) = -1.0;
    197:   expectedC(1,0) =  2.0;
    198:   expectedC(1,1) =  5.0;
    199:   expectedC(1,2) = -1.0;
    200:   expectedC(1,3) = -2.0;
    201:   expectedC(2,0) = -9.0;
    202:   expectedC(2,1) = -1.0;
    203:   expectedC(2,2) = 10.0;
    204:   expectedC(2,3) =  0.0;
    205:   expectedC(3,0) = -1.0;
    206:   expectedC(3,1) = -2.0;
    207:   expectedC(3,2) =  0.0;
    208:   expectedC(3,3) =  1.0;

    210:   int check = areBlitzArrays2NumericallyEqual(actualC,expectedC);
    211:   validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check));
    212:   */

    214:   MatDestroy(&A);
    215:   MatDestroy(&P);
    216:   MatDestroy(&C);
    217:   return(0);
    218: }

    220: /*TEST

    222:    test:

    224:    test:
    225:       suffix: 2
    226:       nsize: 2
    227:       args: -B_matmatmult_via nonscalable

    229:    test:
    230:       suffix: 3
    231:       nsize: 2
    232:       output_file: output/ex93_2.out

    234:    test:
    235:       suffix: 4
    236:       nsize: 2
    237:       args: -A_matptap_via scalable
    238:       output_file: output/ex93_2.out

    240:    test:
    241:       suffix: btheap
    242:       args: -B_matmatmult_via btheap
    243:       output_file: output/ex93_1.out

    245:    test:
    246:       suffix: heap
    247:       args: -B_matmatmult_via heap
    248:       output_file: output/ex93_1.out

    250:    #HYPRE PtAP is broken for complex numbers
    251:    test:
    252:       suffix: hypre
    253:       nsize: 3
    254:       requires: hypre !complex
    255:       args: -B_matmatmult_via hypre -A_matptap_via hypre -test_hypre

    257:    test:
    258:       suffix: llcondensed
    259:       args: -B_matmatmult_via llcondensed
    260:       output_file: output/ex93_1.out

    262:    test:
    263:       suffix: scalable
    264:       args: -B_matmatmult_via scalable
    265:       output_file: output/ex93_1.out

    267:    test:
    268:       suffix: scalable_fast
    269:       args: -B_matmatmult_via scalable_fast
    270:       output_file: output/ex93_1.out

    272: TEST*/