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++) MatSetValue(B, i, i, 1.0 * i, INSERT_VALUES);
 16:   MatSetValue(B, me - 1, me, me * me, INSERT_VALUES);
 17:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
 18:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
 19:   *A = B;
 20:   return 0;
 21: }

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

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

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

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

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

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

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

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

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

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

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

 88:   VecDestroyVecs(2, &ltmp);
 89:   VecDestroyVecs(2, &rtmp);
 90:   return 0;
 91: }

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

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

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

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

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

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

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

162: /*TEST

164:    test:
165:       nsize: 3

167: TEST*/