1 static const char help[] = "Tests MatCreateSubMatrix with MatSubMatrix versus MatAIJ, non-square\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,5,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,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 65 PetscCall(MatMult(A,X,ltmp[0])); 66 PetscCall(MatMult(B,X,ltmp[1])); 67 PetscCall(Compare2(ltmp,"MatMult")); 68 69 PetscCall(MatMultTranspose(A,Y,rtmp[0])); 70 PetscCall(MatMultTranspose(B,Y,rtmp[1])); 71 PetscCall(Compare2(rtmp,"MatMultTranspose")); 72 73 PetscCall(VecCopy(Y1,ltmp[0])); 74 PetscCall(VecCopy(Y1,ltmp[1])); 75 PetscCall(MatMultAdd(A,X,ltmp[0],ltmp[0])); 76 PetscCall(MatMultAdd(B,X,ltmp[1],ltmp[1])); 77 PetscCall(Compare2(ltmp,"MatMultAdd v2==v3")); 78 79 PetscCall(MatMultAdd(A,X,Y1,ltmp[0])); 80 PetscCall(MatMultAdd(B,X,Y1,ltmp[1])); 81 PetscCall(Compare2(ltmp,"MatMultAdd v2!=v3")); 82 83 PetscCall(VecCopy(X1,rtmp[0])); 84 PetscCall(VecCopy(X1,rtmp[1])); 85 PetscCall(MatMultTransposeAdd(A,Y,rtmp[0],rtmp[0])); 86 PetscCall(MatMultTransposeAdd(B,Y,rtmp[1],rtmp[1])); 87 PetscCall(Compare2(rtmp,"MatMultTransposeAdd v2==v3")); 88 89 PetscCall(MatMultTransposeAdd(A,Y,X1,rtmp[0])); 90 PetscCall(MatMultTransposeAdd(B,Y,X1,rtmp[1])); 91 PetscCall(Compare2(rtmp,"MatMultTransposeAdd v2!=v3")); 92 93 PetscCall(VecDestroyVecs(2,<mp)); 94 PetscCall(VecDestroyVecs(2,&rtmp)); 95 PetscFunctionReturn(0); 96 } 97 98 int main(int argc, char *argv[]) 99 { 100 Mat A,B,Asub,Bsub; 101 PetscInt ms,idxrow[3],idxcol[4]; 102 Vec left,right,X,Y,X1,Y1; 103 IS isrow,iscol; 104 PetscBool random = PETSC_TRUE; 105 106 PetscFunctionBeginUser; 107 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 108 PetscCall(AssembleMatrix(PETSC_COMM_WORLD,&A)); 109 PetscCall(AssembleMatrix(PETSC_COMM_WORLD,&B)); 110 PetscCall(MatSetOperation(B,MATOP_CREATE_SUBMATRIX,NULL)); 111 PetscCall(MatSetOperation(B,MATOP_CREATE_SUBMATRICES,NULL)); 112 PetscCall(MatGetOwnershipRange(A,&ms,NULL)); 113 114 idxrow[0] = ms+1; 115 idxrow[1] = ms+2; 116 idxrow[2] = ms+4; 117 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,3,idxrow,PETSC_USE_POINTER,&isrow)); 118 119 idxcol[0] = ms+1; 120 idxcol[1] = ms+2; 121 idxcol[2] = ms+4; 122 idxcol[3] = ms+5; 123 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,4,idxcol,PETSC_USE_POINTER,&iscol)); 124 125 PetscCall(MatCreateSubMatrix(A,isrow,iscol,MAT_INITIAL_MATRIX,&Asub)); 126 PetscCall(MatCreateSubMatrix(B,isrow,iscol,MAT_INITIAL_MATRIX,&Bsub)); 127 128 PetscCall(MatCreateVecs(Asub,&right,&left)); 129 PetscCall(VecDuplicate(right,&X)); 130 PetscCall(VecDuplicate(right,&X1)); 131 PetscCall(VecDuplicate(left,&Y)); 132 PetscCall(VecDuplicate(left,&Y1)); 133 134 PetscCall(PetscOptionsGetBool(NULL,NULL,"-random",&random,NULL)); 135 if (random) { 136 PetscCall(VecSetRandom(right,NULL)); 137 PetscCall(VecSetRandom(left,NULL)); 138 PetscCall(VecSetRandom(X,NULL)); 139 PetscCall(VecSetRandom(Y,NULL)); 140 PetscCall(VecSetRandom(X1,NULL)); 141 PetscCall(VecSetRandom(Y1,NULL)); 142 } else { 143 PetscCall(VecSet(right,1.0)); 144 PetscCall(VecSet(left,2.0)); 145 PetscCall(VecSet(X,3.0)); 146 PetscCall(VecSet(Y,4.0)); 147 PetscCall(VecSet(X1,3.0)); 148 PetscCall(VecSet(Y1,4.0)); 149 } 150 PetscCall(CheckMatrices(Asub,Bsub,left,right,X,Y,X1,Y1)); 151 PetscCall(ISDestroy(&isrow)); 152 PetscCall(ISDestroy(&iscol)); 153 PetscCall(MatDestroy(&A)); 154 PetscCall(MatDestroy(&B)); 155 PetscCall(MatDestroy(&Asub)); 156 PetscCall(MatDestroy(&Bsub)); 157 PetscCall(VecDestroy(&left)); 158 PetscCall(VecDestroy(&right)); 159 PetscCall(VecDestroy(&X)); 160 PetscCall(VecDestroy(&Y)); 161 PetscCall(VecDestroy(&X1)); 162 PetscCall(VecDestroy(&Y1)); 163 PetscCall(PetscFinalize()); 164 return 0; 165 } 166 167 /*TEST 168 169 test: 170 nsize: 3 171 172 TEST*/ 173