1 static const char help[] = "Tests MatCreateSubMatrix with MatSubMatrix versus MatAIJ, square, shifted (copied from ex97)\n"; 2 3 #include <petscmat.h> 4 5 static PetscErrorCode AssembleMatrix(MPI_Comm comm, Mat *A) 6 { 7 Mat B; 8 PetscInt i, ms, me; 9 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 } 23 24 static PetscErrorCode Compare2(Vec *X, const char *test) 25 { 26 PetscReal norm; 27 Vec Y; 28 PetscInt verbose = 0; 29 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)); 35 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 } 50 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; 54 55 PetscFunctionBegin; 56 PetscCall(VecDuplicateVecs(right, 2, &rtmp)); 57 PetscCall(VecDuplicateVecs(left, 2, <mp)); 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)); 64 65 PetscCall(MatMult(A, X, ltmp[0])); 66 PetscCall(MatMult(B, X, ltmp[1])); 67 PetscCall(Compare2(ltmp, "MatMult")); 68 69 PetscCall(MatMultTranspose(A, Y, rtmp[0])); 70 PetscCall(MatMultTranspose(B, Y, rtmp[1])); 71 PetscCall(Compare2(rtmp, "MatMultTranspose")); 72 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")); 78 79 PetscCall(MatMultAdd(A, X, Y1, ltmp[0])); 80 PetscCall(MatMultAdd(B, X, Y1, ltmp[1])); 81 PetscCall(Compare2(ltmp, "MatMultAdd v2!=v3")); 82 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")); 88 89 PetscCall(MatMultTransposeAdd(A, Y, X1, rtmp[0])); 90 PetscCall(MatMultTransposeAdd(B, Y, X1, rtmp[1])); 91 PetscCall(Compare2(rtmp, "MatMultTransposeAdd v2!=v3")); 92 93 PetscCall(VecDestroyVecs(2, <mp)); 94 PetscCall(VecDestroyVecs(2, &rtmp)); 95 PetscFunctionReturn(PETSC_SUCCESS); 96 } 97 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; 105 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)); 113 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)); 118 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)); 123 124 PetscCall(MatCreateSubMatrix(A, isrow, iscol, MAT_INITIAL_MATRIX, &Asub)); 125 PetscCall(MatCreateSubMatrix(B, isrow, iscol, MAT_INITIAL_MATRIX, &Bsub)); 126 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)); 132 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 } 165 166 /*TEST 167 168 test: 169 nsize: 3 170 171 TEST*/ 172