1 #include <petsctao.h> /*I "petsctao.h" I*/ 2 #include <../src/tao/matrix/submatfree.h> /*I "submatfree.h" I*/ 3 4 /*@C 5 MatCreateSubMatrixFree - Creates a reduced matrix by masking a 6 full matrix. 7 8 Collective on matrix 9 10 Input Parameters: 11 + mat - matrix of arbitrary type 12 . Rows - the rows that will be in the submatrix 13 - Cols - the columns that will be in the submatrix 14 15 Output Parameters: 16 . J - New matrix 17 18 Notes: 19 The caller is responsible for destroying the input objects after matrix J has been destroyed. 20 21 Level: developer 22 23 .seealso: MatCreate() 24 @*/ 25 PetscErrorCode MatCreateSubMatrixFree(Mat mat,IS Rows, IS Cols, Mat *J) 26 { 27 MPI_Comm comm=PetscObjectComm((PetscObject)mat); 28 MatSubMatFreeCtx ctx; 29 PetscInt mloc,nloc,m,n; 30 31 PetscFunctionBegin; 32 PetscCall(PetscNew(&ctx)); 33 ctx->A = mat; 34 PetscCall(MatGetSize(mat,&m,&n)); 35 PetscCall(MatGetLocalSize(mat,&mloc,&nloc)); 36 PetscCall(MatCreateVecs(mat,NULL,&ctx->VC)); 37 ctx->VR = ctx->VC; 38 PetscCall(PetscObjectReference((PetscObject)mat)); 39 40 ctx->Rows = Rows; 41 ctx->Cols = Cols; 42 PetscCall(PetscObjectReference((PetscObject)Rows)); 43 PetscCall(PetscObjectReference((PetscObject)Cols)); 44 PetscCall(MatCreateShell(comm,mloc,nloc,m,n,ctx,J)); 45 PetscCall(MatShellSetManageScalingShifts(*J)); 46 PetscCall(MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_SMF)); 47 PetscCall(MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_SMF)); 48 PetscCall(MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_SMF)); 49 PetscCall(MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_SMF)); 50 PetscCall(MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_SMF)); 51 PetscCall(MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_SMF)); 52 PetscCall(MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_SMF)); 53 PetscCall(MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_SMF)); 54 PetscCall(MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_SMF)); 55 PetscCall(MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_SMF)); 56 PetscCall(MatShellSetOperation(*J,MATOP_CREATE_SUBMATRICES,(void(*)(void))MatCreateSubMatrices_SMF)); 57 PetscCall(MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_SMF)); 58 PetscCall(MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_SMF)); 59 PetscCall(MatShellSetOperation(*J,MATOP_CREATE_SUBMATRIX,(void(*)(void))MatCreateSubMatrix_SMF)); 60 PetscCall(MatShellSetOperation(*J,MATOP_GET_ROW_MAX,(void(*)(void))MatDuplicate_SMF)); 61 62 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J))); 63 PetscFunctionReturn(0); 64 } 65 66 PetscErrorCode MatSMFResetRowColumn(Mat mat,IS Rows,IS Cols) 67 { 68 MatSubMatFreeCtx ctx; 69 70 PetscFunctionBegin; 71 PetscCall(MatShellGetContext(mat,&ctx)); 72 PetscCall(ISDestroy(&ctx->Rows)); 73 PetscCall(ISDestroy(&ctx->Cols)); 74 PetscCall(PetscObjectReference((PetscObject)Rows)); 75 PetscCall(PetscObjectReference((PetscObject)Cols)); 76 ctx->Cols=Cols; 77 ctx->Rows=Rows; 78 PetscFunctionReturn(0); 79 } 80 81 PetscErrorCode MatMult_SMF(Mat mat,Vec a,Vec y) 82 { 83 MatSubMatFreeCtx ctx; 84 85 PetscFunctionBegin; 86 PetscCall(MatShellGetContext(mat,&ctx)); 87 PetscCall(VecCopy(a,ctx->VR)); 88 PetscCall(VecISSet(ctx->VR,ctx->Cols,0.0)); 89 PetscCall(MatMult(ctx->A,ctx->VR,y)); 90 PetscCall(VecISSet(y,ctx->Rows,0.0)); 91 PetscFunctionReturn(0); 92 } 93 94 PetscErrorCode MatMultTranspose_SMF(Mat mat,Vec a,Vec y) 95 { 96 MatSubMatFreeCtx ctx; 97 98 PetscFunctionBegin; 99 PetscCall(MatShellGetContext(mat,&ctx)); 100 PetscCall(VecCopy(a,ctx->VC)); 101 PetscCall(VecISSet(ctx->VC,ctx->Rows,0.0)); 102 PetscCall(MatMultTranspose(ctx->A,ctx->VC,y)); 103 PetscCall(VecISSet(y,ctx->Cols,0.0)); 104 PetscFunctionReturn(0); 105 } 106 107 PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is) 108 { 109 MatSubMatFreeCtx ctx; 110 111 PetscFunctionBegin; 112 PetscCall(MatShellGetContext(M,&ctx)); 113 PetscCall(MatDiagonalSet(ctx->A,D,is)); 114 PetscFunctionReturn(0); 115 } 116 117 PetscErrorCode MatDestroy_SMF(Mat mat) 118 { 119 MatSubMatFreeCtx ctx; 120 121 PetscFunctionBegin; 122 PetscCall(MatShellGetContext(mat,&ctx)); 123 PetscCall(MatDestroy(&ctx->A)); 124 PetscCall(ISDestroy(&ctx->Rows)); 125 PetscCall(ISDestroy(&ctx->Cols)); 126 PetscCall(VecDestroy(&ctx->VC)); 127 PetscCall(PetscFree(ctx)); 128 PetscFunctionReturn(0); 129 } 130 131 PetscErrorCode MatView_SMF(Mat mat,PetscViewer viewer) 132 { 133 MatSubMatFreeCtx ctx; 134 135 PetscFunctionBegin; 136 PetscCall(MatShellGetContext(mat,&ctx)); 137 PetscCall(MatView(ctx->A,viewer)); 138 PetscFunctionReturn(0); 139 } 140 141 PetscErrorCode MatShift_SMF(Mat Y, PetscReal a) 142 { 143 MatSubMatFreeCtx ctx; 144 145 PetscFunctionBegin; 146 PetscCall(MatShellGetContext(Y,&ctx)); 147 PetscCall(MatShift(ctx->A,a)); 148 PetscFunctionReturn(0); 149 } 150 151 PetscErrorCode MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M) 152 { 153 MatSubMatFreeCtx ctx; 154 155 PetscFunctionBegin; 156 PetscCall(MatShellGetContext(mat,&ctx)); 157 PetscCall(MatCreateSubMatrixFree(ctx->A,ctx->Rows,ctx->Cols,M)); 158 PetscFunctionReturn(0); 159 } 160 161 PetscErrorCode MatEqual_SMF(Mat A,Mat B,PetscBool *flg) 162 { 163 MatSubMatFreeCtx ctx1,ctx2; 164 PetscBool flg1,flg2,flg3; 165 166 PetscFunctionBegin; 167 PetscCall(MatShellGetContext(A,&ctx1)); 168 PetscCall(MatShellGetContext(B,&ctx2)); 169 PetscCall(ISEqual(ctx1->Rows,ctx2->Rows,&flg2)); 170 PetscCall(ISEqual(ctx1->Cols,ctx2->Cols,&flg3)); 171 if (flg2==PETSC_FALSE || flg3==PETSC_FALSE) { 172 *flg=PETSC_FALSE; 173 } else { 174 PetscCall(MatEqual(ctx1->A,ctx2->A,&flg1)); 175 if (flg1==PETSC_FALSE) { *flg=PETSC_FALSE;} 176 else { *flg=PETSC_TRUE;} 177 } 178 PetscFunctionReturn(0); 179 } 180 181 PetscErrorCode MatScale_SMF(Mat mat, PetscReal a) 182 { 183 MatSubMatFreeCtx ctx; 184 185 PetscFunctionBegin; 186 PetscCall(MatShellGetContext(mat,&ctx)); 187 PetscCall(MatScale(ctx->A,a)); 188 PetscFunctionReturn(0); 189 } 190 191 PetscErrorCode MatTranspose_SMF(Mat mat,Mat *B) 192 { 193 PetscFunctionBegin; 194 PetscFunctionReturn(1); 195 } 196 197 PetscErrorCode MatGetDiagonal_SMF(Mat mat,Vec v) 198 { 199 MatSubMatFreeCtx ctx; 200 201 PetscFunctionBegin; 202 PetscCall(MatShellGetContext(mat,&ctx)); 203 PetscCall(MatGetDiagonal(ctx->A,v)); 204 PetscFunctionReturn(0); 205 } 206 207 PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D) 208 { 209 MatSubMatFreeCtx ctx; 210 211 PetscFunctionBegin; 212 PetscCall(MatShellGetContext(M,&ctx)); 213 PetscCall(MatGetRowMax(ctx->A,D,NULL)); 214 PetscFunctionReturn(0); 215 } 216 217 PetscErrorCode MatCreateSubMatrices_SMF(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B) 218 { 219 PetscInt i; 220 221 PetscFunctionBegin; 222 if (scall == MAT_INITIAL_MATRIX) { 223 PetscCall(PetscCalloc1(n+1,B)); 224 } 225 226 for (i=0; i<n; i++) { 227 PetscCall(MatCreateSubMatrix_SMF(A,irow[i],icol[i],scall,&(*B)[i])); 228 } 229 PetscFunctionReturn(0); 230 } 231 232 PetscErrorCode MatCreateSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll, 233 Mat *newmat) 234 { 235 MatSubMatFreeCtx ctx; 236 237 PetscFunctionBegin; 238 PetscCall(MatShellGetContext(mat,&ctx)); 239 if (newmat) { 240 PetscCall(MatDestroy(&*newmat)); 241 } 242 PetscCall(MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat)); 243 PetscFunctionReturn(0); 244 } 245 246 PetscErrorCode MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals) 247 { 248 MatSubMatFreeCtx ctx; 249 250 PetscFunctionBegin; 251 PetscCall(MatShellGetContext(mat,&ctx)); 252 PetscCall(MatGetRow(ctx->A,row,ncols,cols,vals)); 253 PetscFunctionReturn(0); 254 } 255 256 PetscErrorCode MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals) 257 { 258 MatSubMatFreeCtx ctx; 259 260 PetscFunctionBegin; 261 PetscCall(MatShellGetContext(mat,&ctx)); 262 PetscCall(MatRestoreRow(ctx->A,row,ncols,cols,vals)); 263 PetscFunctionReturn(0); 264 } 265 266 PetscErrorCode MatGetColumnVector_SMF(Mat mat,Vec Y, PetscInt col) 267 { 268 MatSubMatFreeCtx ctx; 269 270 PetscFunctionBegin; 271 PetscCall(MatShellGetContext(mat,&ctx)); 272 PetscCall(MatGetColumnVector(ctx->A,Y,col)); 273 PetscFunctionReturn(0); 274 } 275 276 PetscErrorCode MatNorm_SMF(Mat mat,NormType type,PetscReal *norm) 277 { 278 MatSubMatFreeCtx ctx; 279 280 PetscFunctionBegin; 281 PetscCall(MatShellGetContext(mat,&ctx)); 282 if (type == NORM_FROBENIUS) { 283 *norm = 1.0; 284 } else if (type == NORM_1 || type == NORM_INFINITY) { 285 *norm = 1.0; 286 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 287 PetscFunctionReturn(0); 288 } 289