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