1 static const char help[] = "Tests MatCreateSubMatrix with MatSubMatrix versus MatAIJ, square, shifted (copied from ex97)\n"; 2 3 #include <petscmat.h> 4 5 static PetscErrorCode AssembleMatrix(MPI_Comm comm,Mat *A) 6 { 7 PetscErrorCode ierr; 8 Mat B; 9 PetscInt i,ms,me; 10 11 PetscFunctionBegin; 12 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 13 ierr = MatSetSizes(B,6,6,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 14 ierr = MatSetFromOptions(B);CHKERRQ(ierr); 15 ierr = MatSetUp(B);CHKERRQ(ierr); 16 ierr = MatGetOwnershipRange(B,&ms,&me);CHKERRQ(ierr); 17 for (i=ms; i<me; i++) { 18 ierr = MatSetValue(B,i,i,1.0*i,INSERT_VALUES);CHKERRQ(ierr); 19 } 20 ierr = MatSetValue(B,me-1,me-1,me*me,INSERT_VALUES);CHKERRQ(ierr); 21 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 22 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 23 *A = B; 24 PetscFunctionReturn(0); 25 } 26 27 static PetscErrorCode Compare2(Vec *X,const char *test) 28 { 29 PetscErrorCode ierr; 30 PetscReal norm; 31 Vec Y; 32 PetscInt verbose = 0; 33 34 PetscFunctionBegin; 35 ierr = VecDuplicate(X[0],&Y);CHKERRQ(ierr); 36 ierr = VecCopy(X[0],Y);CHKERRQ(ierr); 37 ierr = VecAYPX(Y,-1.0,X[1]);CHKERRQ(ierr); 38 ierr = VecNorm(Y,NORM_INFINITY,&norm);CHKERRQ(ierr); 39 40 ierr = PetscOptionsGetInt(NULL,NULL,"-verbose",&verbose,NULL);CHKERRQ(ierr); 41 if (norm < PETSC_SQRT_MACHINE_EPSILON && verbose < 1) { 42 ierr = PetscPrintf(PETSC_COMM_WORLD,"%30s: norm difference < sqrt(eps_machine)\n",test);CHKERRQ(ierr); 43 } else { 44 ierr = PetscPrintf(PETSC_COMM_WORLD,"%30s: norm difference %g\n",test,(double)norm);CHKERRQ(ierr); 45 } 46 if (verbose > 1) { 47 ierr = VecView(X[0],PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 48 ierr = VecView(X[1],PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 49 ierr = VecView(Y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 50 } 51 ierr = VecDestroy(&Y);CHKERRQ(ierr); 52 PetscFunctionReturn(0); 53 } 54 55 static PetscErrorCode CheckMatrices(Mat A,Mat B,Vec left,Vec right,Vec X,Vec Y,Vec X1,Vec Y1) 56 { 57 PetscErrorCode ierr; 58 Vec *ltmp,*rtmp; 59 60 PetscFunctionBegin; 61 ierr = VecDuplicateVecs(right,2,&rtmp);CHKERRQ(ierr); 62 ierr = VecDuplicateVecs(left,2,<mp);CHKERRQ(ierr); 63 ierr = MatScale(A,PETSC_PI);CHKERRQ(ierr); 64 ierr = MatScale(B,PETSC_PI);CHKERRQ(ierr); 65 ierr = MatDiagonalScale(A,left,right);CHKERRQ(ierr); 66 ierr = MatDiagonalScale(B,left,right);CHKERRQ(ierr); 67 ierr = MatShift(A,PETSC_PI);CHKERRQ(ierr); 68 ierr = MatShift(B,PETSC_PI);CHKERRQ(ierr); 69 70 ierr = MatMult(A,X,ltmp[0]);CHKERRQ(ierr); 71 ierr = MatMult(B,X,ltmp[1]);CHKERRQ(ierr); 72 ierr = Compare2(ltmp,"MatMult");CHKERRQ(ierr); 73 74 ierr = MatMultTranspose(A,Y,rtmp[0]);CHKERRQ(ierr); 75 ierr = MatMultTranspose(B,Y,rtmp[1]);CHKERRQ(ierr); 76 ierr = Compare2(rtmp,"MatMultTranspose");CHKERRQ(ierr); 77 78 ierr = VecCopy(Y1,ltmp[0]);CHKERRQ(ierr); 79 ierr = VecCopy(Y1,ltmp[1]);CHKERRQ(ierr); 80 ierr = MatMultAdd(A,X,ltmp[0],ltmp[0]);CHKERRQ(ierr); 81 ierr = MatMultAdd(B,X,ltmp[1],ltmp[1]);CHKERRQ(ierr); 82 ierr = Compare2(ltmp,"MatMultAdd v2==v3");CHKERRQ(ierr); 83 84 ierr = MatMultAdd(A,X,Y1,ltmp[0]);CHKERRQ(ierr); 85 ierr = MatMultAdd(B,X,Y1,ltmp[1]);CHKERRQ(ierr); 86 ierr = Compare2(ltmp,"MatMultAdd v2!=v3");CHKERRQ(ierr); 87 88 ierr = VecCopy(X1,rtmp[0]);CHKERRQ(ierr); 89 ierr = VecCopy(X1,rtmp[1]);CHKERRQ(ierr); 90 ierr = MatMultTransposeAdd(A,Y,rtmp[0],rtmp[0]);CHKERRQ(ierr); 91 ierr = MatMultTransposeAdd(B,Y,rtmp[1],rtmp[1]);CHKERRQ(ierr); 92 ierr = Compare2(rtmp,"MatMultTransposeAdd v2==v3");CHKERRQ(ierr); 93 94 ierr = MatMultTransposeAdd(A,Y,X1,rtmp[0]);CHKERRQ(ierr); 95 ierr = MatMultTransposeAdd(B,Y,X1,rtmp[1]);CHKERRQ(ierr); 96 ierr = Compare2(rtmp,"MatMultTransposeAdd v2!=v3");CHKERRQ(ierr); 97 98 ierr = VecDestroyVecs(2,<mp);CHKERRQ(ierr); 99 ierr = VecDestroyVecs(2,&rtmp);CHKERRQ(ierr); 100 PetscFunctionReturn(0); 101 } 102 103 int main(int argc, char *argv[]) 104 { 105 PetscErrorCode ierr; 106 Mat A,B,Asub,Bsub; 107 PetscInt ms,idxrow[3],idxcol[3]; 108 Vec left,right,X,Y,X1,Y1; 109 IS isrow,iscol; 110 PetscBool random = PETSC_TRUE; 111 112 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 113 ierr = AssembleMatrix(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 114 ierr = AssembleMatrix(PETSC_COMM_WORLD,&B);CHKERRQ(ierr); 115 ierr = MatSetOperation(B,MATOP_CREATE_SUBMATRIX,NULL);CHKERRQ(ierr); 116 ierr = MatSetOperation(B,MATOP_CREATE_SUBMATRICES,NULL);CHKERRQ(ierr); 117 ierr = MatGetOwnershipRange(A,&ms,NULL);CHKERRQ(ierr); 118 119 idxrow[0] = ms+1; 120 idxrow[1] = ms+2; 121 idxrow[2] = ms+4; 122 ierr = ISCreateGeneral(PETSC_COMM_WORLD,3,idxrow,PETSC_USE_POINTER,&isrow);CHKERRQ(ierr); 123 124 idxcol[0] = ms+1; 125 idxcol[1] = ms+2; 126 idxcol[2] = ms+4; 127 ierr = ISCreateGeneral(PETSC_COMM_WORLD,3,idxcol,PETSC_USE_POINTER,&iscol);CHKERRQ(ierr); 128 129 ierr = MatCreateSubMatrix(A,isrow,iscol,MAT_INITIAL_MATRIX,&Asub);CHKERRQ(ierr); 130 ierr = MatCreateSubMatrix(B,isrow,iscol,MAT_INITIAL_MATRIX,&Bsub);CHKERRQ(ierr); 131 132 ierr = MatCreateVecs(Asub,&right,&left);CHKERRQ(ierr); 133 ierr = VecDuplicate(right,&X);CHKERRQ(ierr); 134 ierr = VecDuplicate(right,&X1);CHKERRQ(ierr); 135 ierr = VecDuplicate(left,&Y);CHKERRQ(ierr); 136 ierr = VecDuplicate(left,&Y1);CHKERRQ(ierr); 137 138 ierr = PetscOptionsGetBool(NULL,NULL,"-random",&random,NULL);CHKERRQ(ierr); 139 if (random) { 140 ierr = VecSetRandom(right,NULL);CHKERRQ(ierr); 141 ierr = VecSetRandom(left,NULL);CHKERRQ(ierr); 142 ierr = VecSetRandom(X,NULL);CHKERRQ(ierr); 143 ierr = VecSetRandom(Y,NULL);CHKERRQ(ierr); 144 ierr = VecSetRandom(X1,NULL);CHKERRQ(ierr); 145 ierr = VecSetRandom(Y1,NULL);CHKERRQ(ierr); 146 } else { 147 ierr = VecSet(right,1.0);CHKERRQ(ierr); 148 ierr = VecSet(left,2.0);CHKERRQ(ierr); 149 ierr = VecSet(X,3.0);CHKERRQ(ierr); 150 ierr = VecSet(Y,4.0);CHKERRQ(ierr); 151 ierr = VecSet(X1,3.0);CHKERRQ(ierr); 152 ierr = VecSet(Y1,4.0);CHKERRQ(ierr); 153 } 154 ierr = CheckMatrices(Asub,Bsub,left,right,X,Y,X1,Y1);CHKERRQ(ierr); 155 ierr = ISDestroy(&isrow);CHKERRQ(ierr); 156 ierr = ISDestroy(&iscol);CHKERRQ(ierr); 157 ierr = MatDestroy(&A);CHKERRQ(ierr); 158 ierr = MatDestroy(&B);CHKERRQ(ierr); 159 ierr = MatDestroy(&Asub);CHKERRQ(ierr); 160 ierr = MatDestroy(&Bsub);CHKERRQ(ierr); 161 ierr = VecDestroy(&left);CHKERRQ(ierr); 162 ierr = VecDestroy(&right);CHKERRQ(ierr); 163 ierr = VecDestroy(&X);CHKERRQ(ierr); 164 ierr = VecDestroy(&Y);CHKERRQ(ierr); 165 ierr = VecDestroy(&X1);CHKERRQ(ierr); 166 ierr = VecDestroy(&Y1);CHKERRQ(ierr); 167 ierr = PetscFinalize(); 168 return ierr; 169 } 170 171 172 173 /*TEST 174 175 test: 176 nsize: 3 177 178 TEST*/ 179