1 2 static char help[] = "Tests MatCreateComposite()\n\n"; 3 4 /* 5 Include "petscmat.h" so that we can use matrices. 6 automatically includes: 7 petscsys.h - base PETSc routines petscvec.h - vectors 8 petscmat.h - matrices 9 petscis.h - index sets petscviewer.h - viewers 10 */ 11 #include <petscmat.h> 12 13 int main(int argc, char **args) 14 { 15 Mat *A, B; /* matrix */ 16 Vec x, y, v, v2, z, z2; 17 PetscReal rnorm; 18 PetscInt n = 20; /* size of the matrix */ 19 PetscInt nmat = 3; /* number of matrices */ 20 PetscInt i; 21 PetscRandom rctx; 22 MatCompositeType type; 23 PetscScalar scalings[5] = {2, 3, 4, 5, 6}; 24 25 PetscFunctionBeginUser; 26 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 27 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 28 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nmat", &nmat, NULL)); 29 30 /* 31 Create random matrices 32 */ 33 PetscCall(PetscMalloc1(nmat + 3, &A)); 34 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx)); 35 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n / 2, 3, NULL, 3, NULL, &A[0])); 36 for (i = 1; i < nmat + 1; i++) PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, 3, NULL, 3, NULL, &A[i])); 37 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n / 2, n, 3, NULL, 3, NULL, &A[nmat + 1])); 38 for (i = 0; i < nmat + 2; i++) PetscCall(MatSetRandom(A[i], rctx)); 39 40 PetscCall(MatCreateVecs(A[1], &x, &y)); 41 PetscCall(VecDuplicate(y, &z)); 42 PetscCall(VecDuplicate(z, &z2)); 43 PetscCall(MatCreateVecs(A[0], &v, NULL)); 44 PetscCall(VecDuplicate(v, &v2)); 45 46 /* Test MatMult of an ADDITIVE MatComposite B made up of A[1],A[2],A[3] with separate scalings */ 47 48 /* Do MatMult with A[1],A[2],A[3] by hand and store the result in z */ 49 PetscCall(VecSet(x, 1.0)); 50 PetscCall(MatMult(A[1], x, z)); 51 PetscCall(VecScale(z, scalings[1])); 52 for (i = 2; i < nmat + 1; i++) { 53 PetscCall(MatMult(A[i], x, z2)); 54 PetscCall(VecAXPY(z, scalings[i], z2)); 55 } 56 57 /* Do MatMult using MatComposite and store the result in y */ 58 PetscCall(VecSet(y, 0.0)); 59 PetscCall(MatCreateComposite(PETSC_COMM_WORLD, nmat, A + 1, &B)); 60 PetscCall(MatSetFromOptions(B)); 61 PetscCall(MatCompositeSetScalings(B, &scalings[1])); 62 PetscCall(MatMultAdd(B, x, y, y)); 63 64 /* Diff y and z */ 65 PetscCall(VecAXPY(y, -1.0, z)); 66 PetscCall(VecNorm(y, NORM_2, &rnorm)); 67 if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with composite add %g\n", (double)rnorm)); 68 69 /* Test MatCompositeMerge on ADDITIVE MatComposite */ 70 PetscCall(MatCompositeSetMatStructure(B, DIFFERENT_NONZERO_PATTERN)); /* default */ 71 PetscCall(MatCompositeMerge(B)); 72 PetscCall(MatMult(B, x, y)); 73 PetscCall(MatDestroy(&B)); 74 PetscCall(VecAXPY(y, -1.0, z)); 75 PetscCall(VecNorm(y, NORM_2, &rnorm)); 76 if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with composite add after merge %g\n", (double)rnorm)); 77 78 /* 79 Test n x n/2 multiplicative composite B made up of A[0],A[1],A[2] with separate scalings 80 */ 81 82 /* Do MatMult with A[0],A[1],A[2] by hand and store the result in z */ 83 PetscCall(VecSet(v, 1.0)); 84 PetscCall(MatMult(A[0], v, z)); 85 PetscCall(VecScale(z, scalings[0])); 86 for (i = 1; i < nmat; i++) { 87 PetscCall(MatMult(A[i], z, y)); 88 PetscCall(VecScale(y, scalings[i])); 89 PetscCall(VecCopy(y, z)); 90 } 91 92 /* Do MatMult using MatComposite and store the result in y */ 93 PetscCall(MatCreateComposite(PETSC_COMM_WORLD, nmat, A, &B)); 94 PetscCall(MatCompositeSetType(B, MAT_COMPOSITE_MULTIPLICATIVE)); 95 PetscCall(MatCompositeSetMergeType(B, MAT_COMPOSITE_MERGE_LEFT)); 96 PetscCall(MatSetFromOptions(B)); 97 PetscCall(MatCompositeSetScalings(B, &scalings[0])); 98 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 99 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); /* do MatCompositeMerge() if -mat_composite_merge 1 */ 100 PetscCall(MatMult(B, v, y)); 101 PetscCall(MatDestroy(&B)); 102 103 /* Diff y and z */ 104 PetscCall(VecAXPY(y, -1.0, z)); 105 PetscCall(VecNorm(y, NORM_2, &rnorm)); 106 if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with composite multiplicative %g\n", (double)rnorm)); 107 108 /* 109 Test n/2 x n multiplicative composite B made up of A[2], A[3], A[4] without separate scalings 110 */ 111 PetscCall(VecSet(x, 1.0)); 112 PetscCall(MatMult(A[2], x, z)); 113 for (i = 3; i < nmat + 1; i++) { 114 PetscCall(MatMult(A[i], z, y)); 115 PetscCall(VecCopy(y, z)); 116 } 117 PetscCall(MatMult(A[nmat + 1], z, v)); 118 119 PetscCall(MatCreateComposite(PETSC_COMM_WORLD, nmat, A + 2, &B)); 120 PetscCall(MatCompositeSetType(B, MAT_COMPOSITE_MULTIPLICATIVE)); 121 PetscCall(MatSetFromOptions(B)); 122 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 123 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); /* do MatCompositeMerge() if -mat_composite_merge 1 */ 124 PetscCall(MatMult(B, x, v2)); 125 PetscCall(MatDestroy(&B)); 126 127 PetscCall(VecAXPY(v2, -1.0, v)); 128 PetscCall(VecNorm(v2, NORM_2, &rnorm)); 129 if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with composite multiplicative %g\n", (double)rnorm)); 130 131 /* 132 Test get functions 133 */ 134 PetscCall(MatCreateComposite(PETSC_COMM_WORLD, nmat, A, &B)); 135 PetscCall(MatCompositeGetNumberMat(B, &n)); 136 if (nmat != n) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with GetNumberMat %" PetscInt_FMT " != %" PetscInt_FMT "\n", nmat, n)); 137 PetscCall(MatCompositeGetMat(B, 0, &A[nmat + 2])); 138 if (A[0] != A[nmat + 2]) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with GetMat\n")); 139 PetscCall(MatCompositeGetType(B, &type)); 140 if (type != MAT_COMPOSITE_ADDITIVE) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with GetType\n")); 141 PetscCall(MatDestroy(&B)); 142 143 /* 144 Free work space. All PETSc objects should be destroyed when they 145 are no longer needed. 146 */ 147 PetscCall(VecDestroy(&x)); 148 PetscCall(VecDestroy(&y)); 149 PetscCall(VecDestroy(&v)); 150 PetscCall(VecDestroy(&v2)); 151 PetscCall(VecDestroy(&z)); 152 PetscCall(VecDestroy(&z2)); 153 PetscCall(PetscRandomDestroy(&rctx)); 154 for (i = 0; i < nmat + 2; i++) PetscCall(MatDestroy(&A[i])); 155 PetscCall(PetscFree(A)); 156 157 PetscCall(PetscFinalize()); 158 return 0; 159 } 160 161 /*TEST 162 163 test: 164 nsize: 2 165 requires: double 166 args: -mat_composite_merge {{0 1}shared output} -mat_composite_merge_mvctx {{0 1}shared output} 167 168 TEST*/ 169