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