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:   PetscErrorCode  ierr;
 19:   PetscBool       equal;

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

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

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

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

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

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

 70:   MatCreateSeqDense(PETSC_COMM_WORLD,acol[0]+acol[nr]+acol[2*nr],nk,NULL,&B);
 71:   MatSetRandom(B,rctx);
 72:   MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
 73:   MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
 74:   MatMatMultEqual(A,B,C,10,&equal);
 75:   if (!equal) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in C != A*B");

 77:   MatGetSize(A,&M,&N);
 78:   MatGetLocalSize(A,&m,&n);
 79:   for (i=0; i<nk; i++) {
 80:     MatDenseGetColumn(B,i,&valsB);
 81:     VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,valsB,&x);
 82:     MatCreateVecs(A,NULL,&Ax);
 83:     MatMult(A,x,Ax);
 84:     VecDestroy(&x);
 85:     MatDenseGetColumn(C,i,&valsC);
 86:     VecCreateMPIWithArray(PETSC_COMM_WORLD,1,m,M,valsC,&y);
 87:     VecAXPY(y,-1.0,Ax);
 88:     VecDestroy(&Ax);
 89:     VecDestroy(&y);
 90:     MatDenseRestoreColumn(C,&valsC);
 91:     MatDenseRestoreColumn(B,&valsB);
 92:   }
 93:   MatNorm(C,NORM_INFINITY,&norm);
 94:   if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatMult(): %g\n",(double)norm);
 95:   MatDestroy(&C);
 96:   MatDestroy(&B);

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

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

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

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

144:   for (i=0; i<(n * n); i++) {
145:     MatDestroy(&subs[i]);
146:   }
147:   MatDestroy(&A);
148:   PetscRandomDestroy(&rctx);
149:   return(0);
150: }

152: int main(int argc,char **args)
153: {
154:   PetscErrorCode      ierr;

156:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
157:   TestInitialMatrix();
158:   TestReuseMatrix();
159:   PetscFinalize();
160:   return ierr;
161: }

163: /*TEST

165:    test:
166:       args: -malloc_dump

168: TEST*/