1 #include "tao_util.h" /*I "tao_util.h" I*/ 2 #include "submatfree.h" /*I "submatfree.h" I*/ 3 4 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "MatCreateSubMatrixFree" 8 /*@C 9 MatCreateSubMatrixFree - Creates a reduced matrix by masking a 10 full matrix. 11 12 Collective on matrix 13 14 Input Parameters: 15 + mat - matrix of arbitrary type 16 . Rows - the rows that will be in the submatrix 17 - Cols - the columns that will be in the submatrix 18 19 Output Parameters: 20 . J - New matrix 21 22 Notes: 23 The user provides the input data and is responsible for destroying 24 this data after matrix J has been destroyed. 25 26 Level: developer 27 28 .seealso: MatCreate() 29 @*/ 30 PetscErrorCode MatCreateSubMatrixFree(Mat mat,IS Rows, IS Cols, Mat *J) 31 { 32 MPI_Comm comm=((PetscObject)mat)->comm; 33 MatSubMatFreeCtx ctx; 34 PetscErrorCode ierr; 35 PetscMPIInt size; 36 PetscInt mloc,nloc,m,n; 37 38 PetscFunctionBegin; 39 ierr = PetscNew(&ctx);CHKERRQ(ierr); 40 ctx->A=mat; 41 ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr); 42 ierr = MatGetLocalSize(mat,&mloc,&nloc);CHKERRQ(ierr); 43 ierr = MPI_Comm_size(comm,&size); 44 if (size == 1) { 45 ierr = VecCreateSeq(comm,n,&ctx->VC);CHKERRQ(ierr); 46 } else { 47 ierr = VecCreateMPI(comm,nloc,n,&ctx->VC);CHKERRQ(ierr); 48 } 49 ctx->VR=ctx->VC; 50 ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr); 51 52 53 ctx->Rows = Rows; 54 ctx->Cols = Cols; 55 ierr = PetscObjectReference((PetscObject)Rows);CHKERRQ(ierr); 56 ierr = PetscObjectReference((PetscObject)Cols);CHKERRQ(ierr); 57 ierr = MatCreateShell(comm,mloc,nloc,m,n,ctx,J);CHKERRQ(ierr); 58 59 ierr = MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_SMF);CHKERRQ(ierr); 60 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_SMF);CHKERRQ(ierr); 61 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_SMF);CHKERRQ(ierr); 62 ierr = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_SMF);CHKERRQ(ierr); 63 ierr = MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_SMF);CHKERRQ(ierr); 64 ierr = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_SMF);CHKERRQ(ierr); 65 ierr = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_SMF);CHKERRQ(ierr); 66 ierr = MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_SMF);CHKERRQ(ierr); 67 ierr = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_SMF);CHKERRQ(ierr); 68 ierr = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_SMF);CHKERRQ(ierr); 69 ierr = MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)(void))MatGetSubMatrices_SMF);CHKERRQ(ierr); 70 ierr = MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_SMF);CHKERRQ(ierr); 71 ierr = MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_SMF);CHKERRQ(ierr); 72 ierr = MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)(void))MatGetSubMatrix_SMF);CHKERRQ(ierr); 73 ierr = MatShellSetOperation(*J,MATOP_GET_ROW_MAX,(void(*)(void))MatDuplicate_SMF);CHKERRQ(ierr); 74 75 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));CHKERRQ(ierr); 76 PetscFunctionReturn(0); 77 } 78 79 #undef __FUNCT__ 80 #define __FUNCT__ "MatSMFResetRowColumn" 81 PetscErrorCode MatSMFResetRowColumn(Mat mat,IS Rows,IS Cols){ 82 MatSubMatFreeCtx ctx; 83 PetscErrorCode ierr; 84 85 PetscFunctionBegin; 86 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 87 ierr = ISDestroy(&ctx->Rows);CHKERRQ(ierr); 88 ierr = ISDestroy(&ctx->Cols);CHKERRQ(ierr); 89 ierr = PetscObjectReference((PetscObject)Rows);CHKERRQ(ierr); 90 ierr = PetscObjectReference((PetscObject)Cols);CHKERRQ(ierr); 91 ctx->Cols=Cols; 92 ctx->Rows=Rows; 93 PetscFunctionReturn(0); 94 } 95 96 #undef __FUNCT__ 97 #define __FUNCT__ "MatMult_SMF" 98 PetscErrorCode MatMult_SMF(Mat mat,Vec a,Vec y) 99 { 100 MatSubMatFreeCtx ctx; 101 PetscErrorCode ierr; 102 103 PetscFunctionBegin; 104 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 105 ierr = VecCopy(a,ctx->VR);CHKERRQ(ierr); 106 ierr = VecISSetToConstant(ctx->Cols,0.0,ctx->VR);CHKERRQ(ierr); 107 ierr = MatMult(ctx->A,ctx->VR,y);CHKERRQ(ierr); 108 ierr = VecISSetToConstant(ctx->Rows,0.0,y);CHKERRQ(ierr); 109 PetscFunctionReturn(0); 110 } 111 112 #undef __FUNCT__ 113 #define __FUNCT__ "MatMultTranspose_SMF" 114 PetscErrorCode MatMultTranspose_SMF(Mat mat,Vec a,Vec y) 115 { 116 MatSubMatFreeCtx ctx; 117 PetscErrorCode ierr; 118 119 PetscFunctionBegin; 120 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 121 ierr = VecCopy(a,ctx->VC);CHKERRQ(ierr); 122 ierr = VecISSetToConstant(ctx->Rows,0.0,ctx->VC);CHKERRQ(ierr); 123 ierr = MatMultTranspose(ctx->A,ctx->VC,y);CHKERRQ(ierr); 124 ierr = VecISSetToConstant(ctx->Cols,0.0,y);CHKERRQ(ierr); 125 PetscFunctionReturn(0); 126 } 127 128 #undef __FUNCT__ 129 #define __FUNCT__ "MatDiagonalSet_SMF" 130 PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is) 131 { 132 MatSubMatFreeCtx ctx; 133 PetscErrorCode ierr; 134 135 PetscFunctionBegin; 136 ierr = MatShellGetContext(M,(void **)&ctx);CHKERRQ(ierr); 137 ierr = MatDiagonalSet(ctx->A,D,is); 138 PetscFunctionReturn(0); 139 } 140 141 #undef __FUNCT__ 142 #define __FUNCT__ "MatDestroy_SMF" 143 PetscErrorCode MatDestroy_SMF(Mat mat) 144 { 145 PetscErrorCode ierr; 146 MatSubMatFreeCtx ctx; 147 148 PetscFunctionBegin; 149 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 150 ierr = MatDestroy(&ctx->A);CHKERRQ(ierr); 151 ierr = ISDestroy(&ctx->Rows);CHKERRQ(ierr); 152 ierr = ISDestroy(&ctx->Cols);CHKERRQ(ierr); 153 ierr = VecDestroy(&ctx->VC);CHKERRQ(ierr); 154 ierr = PetscFree(ctx);CHKERRQ(ierr); 155 PetscFunctionReturn(0); 156 } 157 158 159 160 #undef __FUNCT__ 161 #define __FUNCT__ "MatView_SMF" 162 PetscErrorCode MatView_SMF(Mat mat,PetscViewer viewer) 163 { 164 PetscErrorCode ierr; 165 MatSubMatFreeCtx ctx; 166 167 PetscFunctionBegin; 168 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 169 ierr = MatView(ctx->A,viewer);CHKERRQ(ierr); 170 PetscFunctionReturn(0); 171 } 172 173 #undef __FUNCT__ 174 #define __FUNCT__ "MatShift_SMF" 175 PetscErrorCode MatShift_SMF(Mat Y, PetscReal a) 176 { 177 PetscErrorCode ierr; 178 MatSubMatFreeCtx ctx; 179 180 PetscFunctionBegin; 181 ierr = MatShellGetContext(Y,(void **)&ctx);CHKERRQ(ierr); 182 ierr = MatShift(ctx->A,a);CHKERRQ(ierr); 183 PetscFunctionReturn(0); 184 } 185 186 #undef __FUNCT__ 187 #define __FUNCT__ "MatDuplicate_SMF" 188 PetscErrorCode MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M) 189 { 190 PetscErrorCode ierr; 191 MatSubMatFreeCtx ctx; 192 193 PetscFunctionBegin; 194 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 195 ierr = MatCreateSubMatrixFree(ctx->A,ctx->Rows,ctx->Cols,M);CHKERRQ(ierr); 196 PetscFunctionReturn(0); 197 } 198 199 #undef __FUNCT__ 200 #define __FUNCT__ "MatEqual_SMF" 201 PetscErrorCode MatEqual_SMF(Mat A,Mat B,PetscBool *flg) 202 { 203 PetscErrorCode ierr; 204 MatSubMatFreeCtx ctx1,ctx2; 205 PetscBool flg1,flg2,flg3; 206 207 PetscFunctionBegin; 208 ierr = MatShellGetContext(A,(void **)&ctx1);CHKERRQ(ierr); 209 ierr = MatShellGetContext(B,(void **)&ctx2);CHKERRQ(ierr); 210 ierr = ISEqual(ctx1->Rows,ctx2->Rows,&flg2);CHKERRQ(ierr); 211 ierr = ISEqual(ctx1->Cols,ctx2->Cols,&flg3);CHKERRQ(ierr); 212 if (flg2==PETSC_FALSE || flg3==PETSC_FALSE){ 213 *flg=PETSC_FALSE; 214 } else { 215 ierr = MatEqual(ctx1->A,ctx2->A,&flg1);CHKERRQ(ierr); 216 if (flg1==PETSC_FALSE){ *flg=PETSC_FALSE;} 217 else { *flg=PETSC_TRUE;} 218 } 219 PetscFunctionReturn(0); 220 } 221 222 #undef __FUNCT__ 223 #define __FUNCT__ "MatScale_SMF" 224 PetscErrorCode MatScale_SMF(Mat mat, PetscReal a) 225 { 226 PetscErrorCode ierr; 227 MatSubMatFreeCtx ctx; 228 229 PetscFunctionBegin; 230 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 231 ierr = MatScale(ctx->A,a);CHKERRQ(ierr); 232 PetscFunctionReturn(0); 233 } 234 235 #undef __FUNCT__ 236 #define __FUNCT__ "MatTranspose_SMF" 237 PetscErrorCode MatTranspose_SMF(Mat mat,Mat *B) 238 { 239 PetscFunctionBegin; 240 PetscFunctionReturn(1); 241 } 242 243 #undef __FUNCT__ 244 #define __FUNCT__ "MatGetDiagonal_SMF" 245 PetscErrorCode MatGetDiagonal_SMF(Mat mat,Vec v) 246 { 247 PetscErrorCode ierr; 248 MatSubMatFreeCtx ctx; 249 250 PetscFunctionBegin; 251 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 252 ierr = MatGetDiagonal(ctx->A,v);CHKERRQ(ierr); 253 PetscFunctionReturn(0); 254 } 255 256 #undef __FUNCT__ 257 #define __FUNCT__ "MatDiagonalSet_SMF" 258 PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D) 259 { 260 MatSubMatFreeCtx ctx; 261 PetscErrorCode ierr; 262 263 PetscFunctionBegin; 264 ierr = MatShellGetContext(M,(void **)&ctx);CHKERRQ(ierr); 265 ierr = MatGetRowMax(ctx->A,D,NULL); 266 PetscFunctionReturn(0); 267 } 268 269 #undef __FUNCT__ 270 #define __FUNCT__ "MatGetSubMatrices_SMF" 271 PetscErrorCode MatGetSubMatrices_SMF(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B) 272 { 273 PetscErrorCode ierr; 274 PetscInt i; 275 276 PetscFunctionBegin; 277 if (scall == MAT_INITIAL_MATRIX) { 278 ierr = PetscMalloc1(n+1,B );CHKERRQ(ierr); 279 } 280 281 for ( i=0; i<n; i++ ) { 282 ierr = MatGetSubMatrix_SMF(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 283 } 284 PetscFunctionReturn(0); 285 } 286 287 #undef __FUNCT__ 288 #define __FUNCT__ "MatGetSubMatrix_SMF" 289 PetscErrorCode MatGetSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll, 290 Mat *newmat) 291 { 292 PetscErrorCode ierr; 293 MatSubMatFreeCtx ctx; 294 295 PetscFunctionBegin; 296 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 297 if (newmat){ 298 ierr = MatDestroy(&*newmat);CHKERRQ(ierr); 299 } 300 ierr = MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat);CHKERRQ(ierr); 301 PetscFunctionReturn(0); 302 } 303 304 #undef __FUNCT__ 305 #define __FUNCT__ "MatGetRow_SMF" 306 PetscErrorCode MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscReal **vals) 307 { 308 PetscErrorCode ierr; 309 MatSubMatFreeCtx ctx; 310 311 PetscFunctionBegin; 312 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 313 ierr = MatGetRow(ctx->A,row,ncols,cols,vals);CHKERRQ(ierr); 314 PetscFunctionReturn(0); 315 } 316 317 #undef __FUNCT__ 318 #define __FUNCT__ "MatRestoreRow_SMF" 319 PetscErrorCode MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscReal **vals) 320 { 321 PetscErrorCode ierr; 322 MatSubMatFreeCtx ctx; 323 324 PetscFunctionBegin; 325 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 326 ierr = MatRestoreRow(ctx->A,row,ncols,cols,vals);CHKERRQ(ierr); 327 PetscFunctionReturn(0); 328 } 329 330 #undef __FUNCT__ 331 #define __FUNCT__ "MatGetColumnVector_SMF" 332 PetscErrorCode MatGetColumnVector_SMF(Mat mat,Vec Y, PetscInt col) 333 { 334 PetscErrorCode ierr; 335 MatSubMatFreeCtx ctx; 336 337 PetscFunctionBegin; 338 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 339 ierr = MatGetColumnVector(ctx->A,Y,col);CHKERRQ(ierr); 340 PetscFunctionReturn(0); 341 } 342 343 #undef __FUNCT__ 344 #define __FUNCT__ "MatConvert_SMF" 345 PetscErrorCode MatConvert_SMF(Mat mat,MatType newtype,Mat *NewMat) 346 { 347 PetscErrorCode ierr; 348 PetscMPIInt size; 349 MatSubMatFreeCtx ctx; 350 351 PetscFunctionBegin; 352 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 353 ierr = MPI_Comm_size(((PetscObject)mat)->comm,&size);CHKERRQ(ierr); 354 PetscFunctionReturn(1); 355 } 356 357 #undef __FUNCT__ 358 #define __FUNCT__ "MatNorm_SMF" 359 PetscErrorCode MatNorm_SMF(Mat mat,NormType type,PetscReal *norm) 360 { 361 PetscErrorCode ierr; 362 MatSubMatFreeCtx ctx; 363 364 PetscFunctionBegin; 365 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 366 if (type == NORM_FROBENIUS) { 367 *norm = 1.0; 368 } else if (type == NORM_1 || type == NORM_INFINITY) { 369 *norm = 1.0; 370 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 371 PetscFunctionReturn(0); 372 } 373 374