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:   PetscFunctionBegin;
 11:   PetscCall(MatCreate(comm, &B));
 12:   PetscCall(MatSetSizes(B, 5, 6, PETSC_DETERMINE, PETSC_DETERMINE));
 13:   PetscCall(MatSetFromOptions(B));
 14:   PetscCall(MatSetUp(B));
 15:   PetscCall(MatGetOwnershipRange(B, &ms, &me));
 16:   for (i = ms; i < me; i++) PetscCall(MatSetValue(B, i, i, 1.0 * i, INSERT_VALUES));
 17:   PetscCall(MatSetValue(B, me - 1, me, me * me, INSERT_VALUES));
 18:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
 19:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
 20:   *A = B;
 21:   PetscFunctionReturn(PETSC_SUCCESS);
 22: }

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

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

 36:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-verbose", &verbose, NULL));
 37:   if (norm < PETSC_SQRT_MACHINE_EPSILON && verbose < 1) {
 38:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%30s: norm difference < sqrt(eps_machine)\n", test));
 39:   } else {
 40:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%30s: norm difference %g\n", test, (double)norm));
 41:   }
 42:   if (verbose > 1) {
 43:     PetscCall(VecView(X[0], PETSC_VIEWER_STDOUT_WORLD));
 44:     PetscCall(VecView(X[1], PETSC_VIEWER_STDOUT_WORLD));
 45:     PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));
 46:   }
 47:   PetscCall(VecDestroy(&Y));
 48:   PetscFunctionReturn(PETSC_SUCCESS);
 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:   PetscFunctionBegin;
 56:   PetscCall(VecDuplicateVecs(right, 2, &rtmp));
 57:   PetscCall(VecDuplicateVecs(left, 2, &ltmp));
 58:   PetscCall(MatScale(A, PETSC_PI));
 59:   PetscCall(MatScale(B, PETSC_PI));
 60:   PetscCall(MatDiagonalScale(A, left, right));
 61:   PetscCall(MatDiagonalScale(B, left, right));

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

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

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

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

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

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

 91:   PetscCall(VecDestroyVecs(2, &ltmp));
 92:   PetscCall(VecDestroyVecs(2, &rtmp));
 93:   PetscFunctionReturn(PETSC_SUCCESS);
 94: }

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

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

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

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

123:   PetscCall(MatCreateSubMatrix(A, isrow, iscol, MAT_INITIAL_MATRIX, &Asub));
124:   PetscCall(MatCreateSubMatrix(B, isrow, iscol, MAT_INITIAL_MATRIX, &Bsub));

126:   PetscCall(MatCreateVecs(Asub, &right, &left));
127:   PetscCall(VecDuplicate(right, &X));
128:   PetscCall(VecDuplicate(right, &X1));
129:   PetscCall(VecDuplicate(left, &Y));
130:   PetscCall(VecDuplicate(left, &Y1));

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

165: /*TEST

167:    test:
168:       nsize: 3

170: TEST*/