1: static const char help[] = "Tests MatCreateSubMatrix with MatSubMatrix versus MatAIJ, square, shifted (copied from ex97)\n";

  3:  #include <petscmat.h>

  5: static PetscErrorCode AssembleMatrix(MPI_Comm comm,Mat *A)
  6: {
  8:   Mat            B;
  9:   PetscInt       i,ms,me;

 12:   MatCreate(comm,&B);
 13:   MatSetSizes(B,6,6,PETSC_DETERMINE,PETSC_DETERMINE);
 14:   MatSetFromOptions(B);
 15:   MatSetUp(B);
 16:   MatGetOwnershipRange(B,&ms,&me);
 17:   for (i=ms; i<me; i++) {
 18:     MatSetValue(B,i,i,1.0*i,INSERT_VALUES);
 19:   }
 20:   MatSetValue(B,me-1,me-1,me*me,INSERT_VALUES);
 21:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 22:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 23:   *A   = B;
 24:   return(0);
 25: }

 27: static PetscErrorCode Compare2(Vec *X,const char *test)
 28: {
 30:   PetscReal      norm;
 31:   Vec            Y;
 32:   PetscInt       verbose = 0;

 35:   VecDuplicate(X[0],&Y);
 36:   VecCopy(X[0],Y);
 37:   VecAYPX(Y,-1.0,X[1]);
 38:   VecNorm(Y,NORM_INFINITY,&norm);

 40:   PetscOptionsGetInt(NULL,NULL,"-verbose",&verbose,NULL);
 41:   if (norm < PETSC_SQRT_MACHINE_EPSILON && verbose < 1) {
 42:     PetscPrintf(PETSC_COMM_WORLD,"%30s: norm difference < sqrt(eps_machine)\n",test);
 43:   } else {
 44:     PetscPrintf(PETSC_COMM_WORLD,"%30s: norm difference %g\n",test,(double)norm);
 45:   }
 46:   if (verbose > 1) {
 47:     VecView(X[0],PETSC_VIEWER_STDOUT_WORLD);
 48:     VecView(X[1],PETSC_VIEWER_STDOUT_WORLD);
 49:     VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
 50:   }
 51:   VecDestroy(&Y);
 52:   return(0);
 53: }

 55: static PetscErrorCode CheckMatrices(Mat A,Mat B,Vec left,Vec right,Vec X,Vec Y,Vec X1,Vec Y1)
 56: {
 58:   Vec            *ltmp,*rtmp;

 61:   VecDuplicateVecs(right,2,&rtmp);
 62:   VecDuplicateVecs(left,2,&ltmp);
 63:   MatScale(A,PETSC_PI);
 64:   MatScale(B,PETSC_PI);
 65:   MatDiagonalScale(A,left,right);
 66:   MatDiagonalScale(B,left,right);
 67:   MatShift(A,PETSC_PI);
 68:   MatShift(B,PETSC_PI);

 70:   MatMult(A,X,ltmp[0]);
 71:   MatMult(B,X,ltmp[1]);
 72:   Compare2(ltmp,"MatMult");

 74:   MatMultTranspose(A,Y,rtmp[0]);
 75:   MatMultTranspose(B,Y,rtmp[1]);
 76:   Compare2(rtmp,"MatMultTranspose");

 78:   VecCopy(Y1,ltmp[0]);
 79:   VecCopy(Y1,ltmp[1]);
 80:   MatMultAdd(A,X,ltmp[0],ltmp[0]);
 81:   MatMultAdd(B,X,ltmp[1],ltmp[1]);
 82:   Compare2(ltmp,"MatMultAdd v2==v3");

 84:   MatMultAdd(A,X,Y1,ltmp[0]);
 85:   MatMultAdd(B,X,Y1,ltmp[1]);
 86:   Compare2(ltmp,"MatMultAdd v2!=v3");

 88:   VecCopy(X1,rtmp[0]);
 89:   VecCopy(X1,rtmp[1]);
 90:   MatMultTransposeAdd(A,Y,rtmp[0],rtmp[0]);
 91:   MatMultTransposeAdd(B,Y,rtmp[1],rtmp[1]);
 92:   Compare2(rtmp,"MatMultTransposeAdd v2==v3");

 94:   MatMultTransposeAdd(A,Y,X1,rtmp[0]);
 95:   MatMultTransposeAdd(B,Y,X1,rtmp[1]);
 96:   Compare2(rtmp,"MatMultTransposeAdd v2!=v3");

 98:   VecDestroyVecs(2,&ltmp);
 99:   VecDestroyVecs(2,&rtmp);
100:   return(0);
101: }

103: int main(int argc, char *argv[])
104: {
106:   Mat            A,B,Asub,Bsub;
107:   PetscInt       ms,idxrow[3],idxcol[3];
108:   Vec            left,right,X,Y,X1,Y1;
109:   IS             isrow,iscol;
110:   PetscBool      random = PETSC_TRUE;

112:   PetscInitialize(&argc,&argv,NULL,help);
113:   AssembleMatrix(PETSC_COMM_WORLD,&A);
114:   AssembleMatrix(PETSC_COMM_WORLD,&B);
115:   MatSetOperation(B,MATOP_CREATE_SUBMATRIX,NULL);
116:   MatSetOperation(B,MATOP_CREATE_SUBMATRICES,NULL);
117:   MatGetOwnershipRange(A,&ms,NULL);

119:   idxrow[0] = ms+1;
120:   idxrow[1] = ms+2;
121:   idxrow[2] = ms+4;
122:   ISCreateGeneral(PETSC_COMM_WORLD,3,idxrow,PETSC_USE_POINTER,&isrow);

124:   idxcol[0] = ms+1;
125:   idxcol[1] = ms+2;
126:   idxcol[2] = ms+4;
127:   ISCreateGeneral(PETSC_COMM_WORLD,3,idxcol,PETSC_USE_POINTER,&iscol);

129:   MatCreateSubMatrix(A,isrow,iscol,MAT_INITIAL_MATRIX,&Asub);
130:   MatCreateSubMatrix(B,isrow,iscol,MAT_INITIAL_MATRIX,&Bsub);

132:   MatCreateVecs(Asub,&right,&left);
133:   VecDuplicate(right,&X);
134:   VecDuplicate(right,&X1);
135:   VecDuplicate(left,&Y);
136:   VecDuplicate(left,&Y1);

138:   PetscOptionsGetBool(NULL,NULL,"-random",&random,NULL);
139:   if (random) {
140:     VecSetRandom(right,NULL);
141:     VecSetRandom(left,NULL);
142:     VecSetRandom(X,NULL);
143:     VecSetRandom(Y,NULL);
144:     VecSetRandom(X1,NULL);
145:     VecSetRandom(Y1,NULL);
146:   } else {
147:     VecSet(right,1.0);
148:     VecSet(left,2.0);
149:     VecSet(X,3.0);
150:     VecSet(Y,4.0);
151:     VecSet(X1,3.0);
152:     VecSet(Y1,4.0);
153:   }
154:   CheckMatrices(Asub,Bsub,left,right,X,Y,X1,Y1);
155:   ISDestroy(&isrow);
156:   ISDestroy(&iscol);
157:   MatDestroy(&A);
158:   MatDestroy(&B);
159:   MatDestroy(&Asub);
160:   MatDestroy(&Bsub);
161:   VecDestroy(&left);
162:   VecDestroy(&right);
163:   VecDestroy(&X);
164:   VecDestroy(&Y);
165:   VecDestroy(&X1);
166:   VecDestroy(&Y1);
167:   PetscFinalize();
168:   return ierr;
169: }



173: /*TEST

175:    test:
176:       nsize: 3

178: TEST*/