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

      1: static char help[] = "Tests the use of MatTranspose_Nest\n";

      3:  #include <petscmat.h>

      5: PetscErrorCode TestInitialMatrix(void)
      6: {
      7:   const PetscInt  nr = 2,nc = 3;
      8:   const PetscInt  arow[2*3] = { 2,2,2,3,3,3 };
      9:   const PetscInt  acol[2*3] = { 3,2,4,3,2,4 };
     10:   Mat             A,Atranspose;
     11:   Mat             subs[2*3], **block;
     12:   Vec             x,y,Ax,ATy;
     13:   PetscInt        i,j;
     14:   PetscScalar     dot1,dot2,zero = 0.0, one = 1.0;
     15:   PetscRandom     rctx;
     16:   PetscErrorCode  ierr;

     19:   PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
     20:   /* Force the random numbers to have imaginary part 0 so printed results are the same for --with-scalar-type=real or --with-scalar-type=complex */
     21:   PetscRandomSetInterval(rctx,zero,one);
     22:   PetscRandomSetFromOptions(rctx);
     23:   for (i=0; i<(nr * nc); i++) {
     24:     MatCreateSeqDense(PETSC_COMM_WORLD,arow[i],acol[i],NULL,&subs[i]);
     25:   }
     26:   MatCreateNest(PETSC_COMM_WORLD,nr,NULL,nc,NULL,subs,&A);
     27:   MatCreateVecs(A, &x, NULL);
     28:   MatCreateVecs(A, NULL, &y);
     29:   VecDuplicate(x, &ATy);
     30:   VecDuplicate(y, &Ax);
     31:   MatSetRandom(A,rctx);
     32:   MatTranspose(A, MAT_INITIAL_MATRIX, &Atranspose);

     34:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);
     35:   MatNestGetSubMats(A, NULL, NULL, &block);
     36:   for (i=0; i<nr; i++) {
     37:     for (j=0; j<nc; j++) {
     38:       MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
     39:     }
     40:   }

     42:   MatView(Atranspose, PETSC_VIEWER_STDOUT_WORLD);
     43:   MatNestGetSubMats(Atranspose, NULL, NULL, &block);
     44:   for (i=0; i<nc; i++) {
     45:     for (j=0; j<nr; j++) {
     46:       MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
     47:     }
     48:   }

     50:   /* Check <Ax, y> = <x, A^Ty> */
     51:   for (i=0; i<10; i++) {
     52:     VecSetRandom(x,rctx);
     53:     VecSetRandom(y,rctx);

     55:     MatMult(A, x, Ax);
     56:     VecDot(Ax, y, &dot1);
     57:     MatMult(Atranspose, y, ATy);
     58:     VecDot(ATy, x, &dot2);

     60:     PetscPrintf(PETSC_COMM_WORLD, "<Ax, y> = %g\n", (double)PetscRealPart(dot1));
     61:     PetscPrintf(PETSC_COMM_WORLD, "<x, A^Ty> = %g\n",(double)PetscRealPart(dot2));
     62:   }

     64:   for (i=0; i<(nr * nc); i++) {
     65:     MatDestroy(&subs[i]);
     66:   }
     67:   MatDestroy(&A);
     68:   MatDestroy(&Atranspose);
     69:   VecDestroy(&x);
     70:   VecDestroy(&y);
     71:   VecDestroy(&Ax);
     72:   VecDestroy(&ATy);
     73:   PetscRandomDestroy(&rctx);
     74:   return(0);
     75: }

     77: PetscErrorCode TestReuseMatrix(void)
     78: {
     79:   const PetscInt  n = 2;
     80:   Mat             A;
     81:   Mat             subs[2*2],**block;
     82:   PetscInt        i,j;
     83:   PetscRandom     rctx;
     84:   PetscErrorCode  ierr;
     85:   PetscScalar     zero = 0.0, one = 1.0;

     88:   PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
     89:   PetscRandomSetInterval(rctx,zero,one);
     90:   PetscRandomSetFromOptions(rctx);
     91:   for (i=0; i<(n * n); i++) {
     92:     MatCreateSeqDense(PETSC_COMM_WORLD,n,n,NULL,&subs[i]);
     93:   }
     94:   MatCreateNest(PETSC_COMM_WORLD,n,NULL,n,NULL,subs,&A);
     95:   MatSetRandom(A,rctx);

     97:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);
     98:   MatNestGetSubMats(A, NULL, NULL, &block);
     99:   for (i=0; i<n; i++) {
    100:     for (j=0; j<n; j++) {
    101:       MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
    102:     }
    103:   }
    104:   MatTranspose(A,MAT_INPLACE_MATRIX,&A);
    105:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);
    106:   MatNestGetSubMats(A, NULL, NULL, &block);
    107:   for (i=0; i<n; i++) {
    108:     for (j=0; j<n; j++) {
    109:       MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
    110:     }
    111:   }

    113:   for (i=0; i<(n * n); i++) {
    114:     MatDestroy(&subs[i]);
    115:   }
    116:   MatDestroy(&A);
    117:   PetscRandomDestroy(&rctx);
    118:   return(0);
    119: }

    121: int main(int argc,char **args)
    122: {
    123:   PetscErrorCode      ierr;

    125:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
    126:   TestInitialMatrix();
    127:   TestReuseMatrix();
    128:   PetscFinalize();
    129:   return ierr;
    130: }

    132: /*TEST

    134:    test:
    135:       args: -malloc_dump

    137: TEST*/