Actual source code: ex202.c

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

  3: #include <petscmat.h>

  5: PetscErrorCode TestInitialMatrix(void)
  6: {
  7:   const PetscInt  nr = 2,nc = 3,nk = 10;
  8:   PetscInt        n,N,m,M;
  9:   const PetscInt  arow[2*3] = { 2,2,2,3,3,3 };
 10:   const PetscInt  acol[2*3] = { 3,2,4,3,2,4 };
 11:   Mat             A,Atranspose,B,C;
 12:   Mat             subs[2*3],**block;
 13:   Vec             x,y,Ax,ATy;
 14:   PetscInt        i,j;
 15:   PetscScalar     dot1,dot2,zero = 0.0,one = 1.0,*valsB,*valsC;
 16:   PetscReal       norm;
 17:   PetscRandom     rctx;
 18:   PetscBool       equal;

 20:   PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
 21:   /* 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 */
 22:   PetscRandomSetInterval(rctx,zero,one);
 23:   PetscRandomSetFromOptions(rctx);
 24:   for (i=0; i<(nr * nc); i++) {
 25:     MatCreateSeqDense(PETSC_COMM_WORLD,arow[i],acol[i],NULL,&subs[i]);
 26:   }
 27:   MatCreateNest(PETSC_COMM_WORLD,nr,NULL,nc,NULL,subs,&A);
 28:   MatCreateVecs(A, &x, NULL);
 29:   MatCreateVecs(A, NULL, &y);
 30:   VecDuplicate(x, &ATy);
 31:   VecDuplicate(y, &Ax);
 32:   MatSetRandom(A,rctx);
 33:   MatTranspose(A, MAT_INITIAL_MATRIX, &Atranspose);

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

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

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

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

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

 68:   MatCreateSeqDense(PETSC_COMM_WORLD,acol[0]+acol[nr]+acol[2*nr],nk,NULL,&B);
 69:   MatSetRandom(B,rctx);
 70:   MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
 71:   MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
 72:   MatMatMultEqual(A,B,C,10,&equal);

 75:   MatGetSize(A,&M,&N);
 76:   MatGetLocalSize(A,&m,&n);
 77:   for (i=0; i<nk; i++) {
 78:     MatDenseGetColumn(B,i,&valsB);
 79:     VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,valsB,&x);
 80:     MatCreateVecs(A,NULL,&Ax);
 81:     MatMult(A,x,Ax);
 82:     VecDestroy(&x);
 83:     MatDenseGetColumn(C,i,&valsC);
 84:     VecCreateMPIWithArray(PETSC_COMM_WORLD,1,m,M,valsC,&y);
 85:     VecAXPY(y,-1.0,Ax);
 86:     VecDestroy(&Ax);
 87:     VecDestroy(&y);
 88:     MatDenseRestoreColumn(C,&valsC);
 89:     MatDenseRestoreColumn(B,&valsB);
 90:   }
 91:   MatNorm(C,NORM_INFINITY,&norm);
 93:   MatDestroy(&C);
 94:   MatDestroy(&B);

 96:   for (i=0; i<(nr * nc); i++) {
 97:     MatDestroy(&subs[i]);
 98:   }
 99:   MatDestroy(&A);
100:   MatDestroy(&Atranspose);
101:   VecDestroy(&ATy);
102:   PetscRandomDestroy(&rctx);
103:   return 0;
104: }

106: PetscErrorCode TestReuseMatrix(void)
107: {
108:   const PetscInt  n = 2;
109:   Mat             A;
110:   Mat             subs[2*2],**block;
111:   PetscInt        i,j;
112:   PetscRandom     rctx;
113:   PetscScalar     zero = 0.0, one = 1.0;

115:   PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
116:   PetscRandomSetInterval(rctx,zero,one);
117:   PetscRandomSetFromOptions(rctx);
118:   for (i=0; i<(n * n); i++) {
119:     MatCreateSeqDense(PETSC_COMM_WORLD,n,n,NULL,&subs[i]);
120:   }
121:   MatCreateNest(PETSC_COMM_WORLD,n,NULL,n,NULL,subs,&A);
122:   MatSetRandom(A,rctx);

124:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);
125:   MatNestGetSubMats(A, NULL, NULL, &block);
126:   for (i=0; i<n; i++) {
127:     for (j=0; j<n; j++) {
128:       MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
129:     }
130:   }
131:   MatTranspose(A,MAT_INPLACE_MATRIX,&A);
132:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);
133:   MatNestGetSubMats(A, NULL, NULL, &block);
134:   for (i=0; i<n; i++) {
135:     for (j=0; j<n; j++) {
136:       MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
137:     }
138:   }

140:   for (i=0; i<(n * n); i++) {
141:     MatDestroy(&subs[i]);
142:   }
143:   MatDestroy(&A);
144:   PetscRandomDestroy(&rctx);
145:   return 0;
146: }

148: int main(int argc,char **args)
149: {

151:   PetscInitialize(&argc,&args,(char*)0,help);
152:   TestInitialMatrix();
153:   TestReuseMatrix();
154:   PetscFinalize();
155:   return 0;
156: }

158: /*TEST

160:    test:
161:       args: -malloc_dump

163: TEST*/