Actual source code: ex99.c

  1: static const char help[] = "Tests MatCreateSubMatrix with MatSubMatrix versus MatAIJ, square, shifted (copied from ex97)\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, 6, 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 - 1, 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));
 62:   PetscCall(MatShift(A, PETSC_PI));
 63:   PetscCall(MatShift(B, PETSC_PI));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

166: /*TEST

168:    test:
169:       nsize: 3

171: TEST*/