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: }