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