Actual source code: ex97.c

  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:   MatGetOwnershipRange(B,&ms,&me);
 18:   for (i=ms; i<me; i++) {
 19:     MatSetValue(B,i,i,1.0*i,INSERT_VALUES);
 20:   }
 21:   MatSetValue(B,me-1,me,me*me,INSERT_VALUES);
 22:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 23:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 24:   *A = B;
 25:   return(0);
 26: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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