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