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