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++) { 37 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,3,NULL,3,NULL,&A[i])); 38 } 39 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n/2,n,3,NULL,3,NULL,&A[nmat+1])); 40 for (i = 0; i < nmat+2; i++) { 41 PetscCall(MatSetRandom(A[i],rctx)); 42 } 43 44 PetscCall(MatCreateVecs(A[1],&x,&y)); 45 PetscCall(VecDuplicate(y,&z)); 46 PetscCall(VecDuplicate(z,&z2)); 47 PetscCall(MatCreateVecs(A[0],&v,NULL)); 48 PetscCall(VecDuplicate(v,&v2)); 49 50 /* Test MatMult of an ADDITIVE MatComposite B made up of A[1],A[2],A[3] with separate scalings */ 51 52 /* Do MatMult with A[1],A[2],A[3] by hand and store the result in z */ 53 PetscCall(VecSet(x,1.0)); 54 PetscCall(MatMult(A[1],x,z)); 55 PetscCall(VecScale(z,scalings[1])); 56 for (i = 2; i < nmat+1; i++) { 57 PetscCall(MatMult(A[i],x,z2)); 58 PetscCall(VecAXPY(z,scalings[i],z2)); 59 } 60 61 /* Do MatMult using MatComposite and store the result in y */ 62 PetscCall(VecSet(y,0.0)); 63 PetscCall(MatCreateComposite(PETSC_COMM_WORLD,nmat,A+1,&B)); 64 PetscCall(MatSetFromOptions(B)); 65 PetscCall(MatCompositeSetScalings(B,&scalings[1])); 66 PetscCall(MatMultAdd(B,x,y,y)); 67 68 /* Diff y and z */ 69 PetscCall(VecAXPY(y,-1.0,z)); 70 PetscCall(VecNorm(y,NORM_2,&rnorm)); 71 if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) { 72 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with composite add %g\n",(double)rnorm)); 73 } 74 75 /* Test MatCompositeMerge on ADDITIVE MatComposite */ 76 PetscCall(MatCompositeSetMatStructure(B,DIFFERENT_NONZERO_PATTERN)); /* default */ 77 PetscCall(MatCompositeMerge(B)); 78 PetscCall(MatMult(B,x,y)); 79 PetscCall(MatDestroy(&B)); 80 PetscCall(VecAXPY(y,-1.0,z)); 81 PetscCall(VecNorm(y,NORM_2,&rnorm)); 82 if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) { 83 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with composite add after merge %g\n",(double)rnorm)); 84 } 85 86 /* 87 Test n x n/2 multiplicative composite B made up of A[0],A[1],A[2] with separate scalings 88 */ 89 90 /* Do MatMult with A[0],A[1],A[2] by hand and store the result in z */ 91 PetscCall(VecSet(v,1.0)); 92 PetscCall(MatMult(A[0],v,z)); 93 PetscCall(VecScale(z,scalings[0])); 94 for (i = 1; i < nmat; i++) { 95 PetscCall(MatMult(A[i],z,y)); 96 PetscCall(VecScale(y,scalings[i])); 97 PetscCall(VecCopy(y,z)); 98 } 99 100 /* Do MatMult using MatComposite and store the result in y */ 101 PetscCall(MatCreateComposite(PETSC_COMM_WORLD,nmat,A,&B)); 102 PetscCall(MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE)); 103 PetscCall(MatCompositeSetMergeType(B,MAT_COMPOSITE_MERGE_LEFT)); 104 PetscCall(MatSetFromOptions(B)); 105 PetscCall(MatCompositeSetScalings(B,&scalings[0])); 106 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 107 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); /* do MatCompositeMerge() if -mat_composite_merge 1 */ 108 PetscCall(MatMult(B,v,y)); 109 PetscCall(MatDestroy(&B)); 110 111 /* Diff y and z */ 112 PetscCall(VecAXPY(y,-1.0,z)); 113 PetscCall(VecNorm(y,NORM_2,&rnorm)); 114 if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) { 115 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm)); 116 } 117 118 /* 119 Test n/2 x n multiplicative composite B made up of A[2], A[3], A[4] without separate scalings 120 */ 121 PetscCall(VecSet(x,1.0)); 122 PetscCall(MatMult(A[2],x,z)); 123 for (i = 3; i < nmat+1; i++) { 124 PetscCall(MatMult(A[i],z,y)); 125 PetscCall(VecCopy(y,z)); 126 } 127 PetscCall(MatMult(A[nmat+1],z,v)); 128 129 PetscCall(MatCreateComposite(PETSC_COMM_WORLD,nmat,A+2,&B)); 130 PetscCall(MatCompositeSetType(B,MAT_COMPOSITE_MULTIPLICATIVE)); 131 PetscCall(MatSetFromOptions(B)); 132 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 133 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); /* do MatCompositeMerge() if -mat_composite_merge 1 */ 134 PetscCall(MatMult(B,x,v2)); 135 PetscCall(MatDestroy(&B)); 136 137 PetscCall(VecAXPY(v2,-1.0,v)); 138 PetscCall(VecNorm(v2,NORM_2,&rnorm)); 139 if (rnorm > 10000.0*PETSC_MACHINE_EPSILON) { 140 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with composite multiplicative %g\n",(double)rnorm)); 141 } 142 143 /* 144 Test get functions 145 */ 146 PetscCall(MatCreateComposite(PETSC_COMM_WORLD,nmat,A,&B)); 147 PetscCall(MatCompositeGetNumberMat(B,&n)); 148 if (nmat != n) { 149 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with GetNumberMat %" PetscInt_FMT " != %" PetscInt_FMT "\n",nmat,n)); 150 } 151 PetscCall(MatCompositeGetMat(B,0,&A[nmat+2])); 152 if (A[0] != A[nmat+2]) { 153 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with GetMat\n")); 154 } 155 PetscCall(MatCompositeGetType(B,&type)); 156 if (type != MAT_COMPOSITE_ADDITIVE) { 157 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error with GetType\n")); 158 } 159 PetscCall(MatDestroy(&B)); 160 161 /* 162 Free work space. All PETSc objects should be destroyed when they 163 are no longer needed. 164 */ 165 PetscCall(VecDestroy(&x)); 166 PetscCall(VecDestroy(&y)); 167 PetscCall(VecDestroy(&v)); 168 PetscCall(VecDestroy(&v2)); 169 PetscCall(VecDestroy(&z)); 170 PetscCall(VecDestroy(&z2)); 171 PetscCall(PetscRandomDestroy(&rctx)); 172 for (i = 0; i < nmat+2; i++) { 173 PetscCall(MatDestroy(&A[i])); 174 } 175 PetscCall(PetscFree(A)); 176 177 PetscCall(PetscFinalize()); 178 return 0; 179 } 180 181 /*TEST 182 183 test: 184 nsize: 2 185 requires: double 186 args: -mat_composite_merge {{0 1}shared output} -mat_composite_merge_mvctx {{0 1}shared output} 187 188 TEST*/ 189