Actual source code: ex97.c

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

  3: #include <petscmat.h>

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

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

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

 31:   VecDuplicate(X[0],&Y);
 32:   VecCopy(X[0],Y);
 33:   VecAYPX(Y,-1.0,X[1]);
 34:   VecNorm(Y,NORM_INFINITY,&norm);

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

 51: static PetscErrorCode CheckMatrices(Mat A,Mat B,Vec left,Vec right,Vec X,Vec Y,Vec X1,Vec Y1)
 52: {
 53:   Vec            *ltmp,*rtmp;

 55:   VecDuplicateVecs(right,2,&rtmp);
 56:   VecDuplicateVecs(left,2,&ltmp);
 57:   MatScale(A,PETSC_PI);
 58:   MatScale(B,PETSC_PI);
 59:   MatDiagonalScale(A,left,right);
 60:   MatDiagonalScale(B,left,right);

 62:   MatMult(A,X,ltmp[0]);
 63:   MatMult(B,X,ltmp[1]);
 64:   Compare2(ltmp,"MatMult");

 66:   MatMultTranspose(A,Y,rtmp[0]);
 67:   MatMultTranspose(B,Y,rtmp[1]);
 68:   Compare2(rtmp,"MatMultTranspose");

 70:   VecCopy(Y1,ltmp[0]);
 71:   VecCopy(Y1,ltmp[1]);
 72:   MatMultAdd(A,X,ltmp[0],ltmp[0]);
 73:   MatMultAdd(B,X,ltmp[1],ltmp[1]);
 74:   Compare2(ltmp,"MatMultAdd v2==v3");

 76:   MatMultAdd(A,X,Y1,ltmp[0]);
 77:   MatMultAdd(B,X,Y1,ltmp[1]);
 78:   Compare2(ltmp,"MatMultAdd v2!=v3");

 80:   VecCopy(X1,rtmp[0]);
 81:   VecCopy(X1,rtmp[1]);
 82:   MatMultTransposeAdd(A,Y,rtmp[0],rtmp[0]);
 83:   MatMultTransposeAdd(B,Y,rtmp[1],rtmp[1]);
 84:   Compare2(rtmp,"MatMultTransposeAdd v2==v3");

 86:   MatMultTransposeAdd(A,Y,X1,rtmp[0]);
 87:   MatMultTransposeAdd(B,Y,X1,rtmp[1]);
 88:   Compare2(rtmp,"MatMultTransposeAdd v2!=v3");

 90:   VecDestroyVecs(2,&ltmp);
 91:   VecDestroyVecs(2,&rtmp);
 92:   return 0;
 93: }

 95: int main(int argc, char *argv[])
 96: {
 97:   Mat            A,B,Asub,Bsub;
 98:   PetscInt       ms,idxrow[3],idxcol[4];
 99:   Vec            left,right,X,Y,X1,Y1;
100:   IS             isrow,iscol;
101:   PetscBool      random = PETSC_TRUE;

103:   PetscInitialize(&argc,&argv,NULL,help);
104:   AssembleMatrix(PETSC_COMM_WORLD,&A);
105:   AssembleMatrix(PETSC_COMM_WORLD,&B);
106:   MatSetOperation(B,MATOP_CREATE_SUBMATRIX,NULL);
107:   MatSetOperation(B,MATOP_CREATE_SUBMATRICES,NULL);
108:   MatGetOwnershipRange(A,&ms,NULL);

110:   idxrow[0] = ms+1;
111:   idxrow[1] = ms+2;
112:   idxrow[2] = ms+4;
113:   ISCreateGeneral(PETSC_COMM_WORLD,3,idxrow,PETSC_USE_POINTER,&isrow);

115:   idxcol[0] = ms+1;
116:   idxcol[1] = ms+2;
117:   idxcol[2] = ms+4;
118:   idxcol[3] = ms+5;
119:   ISCreateGeneral(PETSC_COMM_WORLD,4,idxcol,PETSC_USE_POINTER,&iscol);

121:   MatCreateSubMatrix(A,isrow,iscol,MAT_INITIAL_MATRIX,&Asub);
122:   MatCreateSubMatrix(B,isrow,iscol,MAT_INITIAL_MATRIX,&Bsub);

124:   MatCreateVecs(Asub,&right,&left);
125:   VecDuplicate(right,&X);
126:   VecDuplicate(right,&X1);
127:   VecDuplicate(left,&Y);
128:   VecDuplicate(left,&Y1);

130:   PetscOptionsGetBool(NULL,NULL,"-random",&random,NULL);
131:   if (random) {
132:     VecSetRandom(right,NULL);
133:     VecSetRandom(left,NULL);
134:     VecSetRandom(X,NULL);
135:     VecSetRandom(Y,NULL);
136:     VecSetRandom(X1,NULL);
137:     VecSetRandom(Y1,NULL);
138:   } else {
139:     VecSet(right,1.0);
140:     VecSet(left,2.0);
141:     VecSet(X,3.0);
142:     VecSet(Y,4.0);
143:     VecSet(X1,3.0);
144:     VecSet(Y1,4.0);
145:   }
146:   CheckMatrices(Asub,Bsub,left,right,X,Y,X1,Y1);
147:   ISDestroy(&isrow);
148:   ISDestroy(&iscol);
149:   MatDestroy(&A);
150:   MatDestroy(&B);
151:   MatDestroy(&Asub);
152:   MatDestroy(&Bsub);
153:   VecDestroy(&left);
154:   VecDestroy(&right);
155:   VecDestroy(&X);
156:   VecDestroy(&Y);
157:   VecDestroy(&X1);
158:   VecDestroy(&Y1);
159:   PetscFinalize();
160:   return 0;
161: }

163: /*TEST

165:    test:
166:       nsize: 3

168: TEST*/