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