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