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