1: static const char help[] = "Tests MatGetSubMatrix with MatSubMatrix versus MatAIJ, non-square\n";

  3: #include <petscmat.h>

  7: static PetscErrorCode AssembleMatrix(MPI_Comm comm,Mat *A)
  8: {
 10:   Mat            B;
 11:   PetscInt       i,ms,me;

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

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

 39:   VecDuplicate(X[0],&Y);
 40:   VecCopy(X[0],Y);
 41:   VecAYPX(Y,-1.0,X[1]);
 42:   VecNorm(Y,NORM_INFINITY,&norm);

 44:   PetscOptionsGetInt(NULL,"-verbose",&verbose,NULL);
 45:   if (norm < 1.e-12 && verbose < 1) {
 46:     PetscPrintf(PETSC_COMM_WORLD,"%30s: norm difference < 1e-12\n",test);
 47:   } else {
 48:     PetscPrintf(PETSC_COMM_WORLD,"%30s: norm difference %g\n",test,(double)norm);
 49:   }
 50:   if (verbose > 1) {
 51:     VecView(X[0],PETSC_VIEWER_STDOUT_WORLD);
 52:     VecView(X[1],PETSC_VIEWER_STDOUT_WORLD);
 53:     VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
 54:   }
 55:   VecDestroy(&Y);
 56:   return(0);
 57: }

 61: static PetscErrorCode CheckMatrices(Mat A,Mat B,Vec left,Vec right,Vec X,Vec Y,Vec X1,Vec Y1)
 62: {
 64:   Vec            *ltmp,*rtmp;

 67:   VecDuplicateVecs(right,2,&rtmp);
 68:   VecDuplicateVecs(left,2,&ltmp);
 69:   MatScale(A,PETSC_PI);
 70:   MatScale(B,PETSC_PI);
 71:   MatDiagonalScale(A,left,right);
 72:   MatDiagonalScale(B,left,right);

 74:   MatMult(A,X,ltmp[0]);
 75:   MatMult(B,X,ltmp[1]);
 76:   Compare2(ltmp,"MatMult");

 78:   MatMultTranspose(A,Y,rtmp[0]);
 79:   MatMultTranspose(B,Y,rtmp[1]);
 80:   Compare2(rtmp,"MatMultTranspose");

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

 88:   MatMultAdd(A,X,Y1,ltmp[0]);
 89:   MatMultAdd(B,X,Y1,ltmp[1]);
 90:   Compare2(ltmp,"MatMultAdd v2!=v3");

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

 98:   MatMultTransposeAdd(A,Y,X1,rtmp[0]);
 99:   MatMultTransposeAdd(B,Y,X1,rtmp[1]);
100:   Compare2(rtmp,"MatMultTransposeAdd v2!=v3");

102:   VecDestroyVecs(2,&ltmp);
103:   VecDestroyVecs(2,&rtmp);
104:   return(0);
105: }

109: int main(int argc, char *argv[])
110: {
112:   Mat            A,B,Asub,Bsub;
113:   PetscInt       ms,idxrow[3],idxcol[4];
114:   Vec            left,right,X,Y,X1,Y1;
115:   IS             isrow,iscol;
116:   PetscBool      random = PETSC_TRUE;

118:   PetscInitialize(&argc,&argv,NULL,help);
119:   AssembleMatrix(PETSC_COMM_WORLD,&A);
120:   AssembleMatrix(PETSC_COMM_WORLD,&B);
121:   MatShellSetOperation(B,MATOP_GET_SUBMATRIX,NULL);
122:   MatShellSetOperation(B,MATOP_GET_SUBMATRICES,NULL);
123:   MatGetOwnershipRange(A,&ms,NULL);

125:   idxrow[0] = ms+1;
126:   idxrow[1] = ms+2;
127:   idxrow[2] = ms+4;
128:   ISCreateGeneral(PETSC_COMM_WORLD,3,idxrow,PETSC_USE_POINTER,&isrow);

130:   idxcol[0] = ms+1;
131:   idxcol[1] = ms+2;
132:   idxcol[2] = ms+4;
133:   idxcol[3] = ms+5;
134:   ISCreateGeneral(PETSC_COMM_WORLD,4,idxcol,PETSC_USE_POINTER,&iscol);

136:   MatGetSubMatrix(A,isrow,iscol,MAT_INITIAL_MATRIX,&Asub);
137:   MatGetSubMatrix(B,isrow,iscol,MAT_INITIAL_MATRIX,&Bsub);

139:   MatCreateVecs(Asub,&right,&left);
140:   VecDuplicate(right,&X);
141:   VecDuplicate(right,&X1);
142:   VecDuplicate(left,&Y);
143:   VecDuplicate(left,&Y1);

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