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 PetscFunctionBegin; 38 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 77 PetscFunctionReturn(0); 78 } 79 80 #undef __FUNCT__ 81 #define __FUNCT__ "MatSMFResetRowColumn" 82 PetscErrorCode MatSMFResetRowColumn(Mat mat,IS Rows,IS Cols){ 83 MatSubMatFreeCtx ctx; 84 PetscErrorCode ierr; 85 86 PetscFunctionBegin; 87 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 88 ierr = ISDestroy(&ctx->Rows);CHKERRQ(ierr); 89 ierr = ISDestroy(&ctx->Cols);CHKERRQ(ierr); 90 ierr = PetscObjectReference((PetscObject)Rows);CHKERRQ(ierr); 91 ierr = PetscObjectReference((PetscObject)Cols);CHKERRQ(ierr); 92 ctx->Cols=Cols; 93 ctx->Rows=Rows; 94 PetscFunctionReturn(0); 95 } 96 97 #undef __FUNCT__ 98 #define __FUNCT__ "MatMult_SMF" 99 PetscErrorCode MatMult_SMF(Mat mat,Vec a,Vec y) 100 { 101 MatSubMatFreeCtx ctx; 102 PetscErrorCode ierr; 103 104 PetscFunctionBegin; 105 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 106 ierr = VecCopy(a,ctx->VR);CHKERRQ(ierr); 107 ierr = VecISSetToConstant(ctx->Cols,0.0,ctx->VR);CHKERRQ(ierr); 108 ierr = MatMult(ctx->A,ctx->VR,y);CHKERRQ(ierr); 109 ierr = VecISSetToConstant(ctx->Rows,0.0,y);CHKERRQ(ierr); 110 PetscFunctionReturn(0); 111 } 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "MatMultTranspose_SMF" 115 PetscErrorCode MatMultTranspose_SMF(Mat mat,Vec a,Vec y) 116 { 117 MatSubMatFreeCtx ctx; 118 PetscErrorCode ierr; 119 120 PetscFunctionBegin; 121 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 122 ierr = VecCopy(a,ctx->VC);CHKERRQ(ierr); 123 ierr = VecISSetToConstant(ctx->Rows,0.0,ctx->VC);CHKERRQ(ierr); 124 ierr = MatMultTranspose(ctx->A,ctx->VC,y);CHKERRQ(ierr); 125 ierr = VecISSetToConstant(ctx->Cols,0.0,y);CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 } 128 129 #undef __FUNCT__ 130 #define __FUNCT__ "MatDiagonalSet_SMF" 131 PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is) 132 { 133 MatSubMatFreeCtx ctx; 134 PetscErrorCode ierr; 135 136 PetscFunctionBegin; 137 ierr = MatShellGetContext(M,(void **)&ctx);CHKERRQ(ierr); 138 ierr = MatDiagonalSet(ctx->A,D,is); 139 PetscFunctionReturn(0); 140 } 141 142 #undef __FUNCT__ 143 #define __FUNCT__ "MatDestroy_SMF" 144 PetscErrorCode MatDestroy_SMF(Mat mat) 145 { 146 PetscErrorCode ierr; 147 MatSubMatFreeCtx ctx; 148 149 PetscFunctionBegin; 150 ierr =MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 151 if (ctx->A) { 152 ierr =MatDestroy(&ctx->A);CHKERRQ(ierr); 153 } 154 if (ctx->Rows) { 155 ierr =ISDestroy(&ctx->Rows);CHKERRQ(ierr); 156 } 157 if (ctx->Cols) { 158 ierr =ISDestroy(&ctx->Cols);CHKERRQ(ierr); 159 } 160 if (ctx->VC) { 161 ierr =VecDestroy(&ctx->VC);CHKERRQ(ierr); 162 } 163 ierr = PetscFree(ctx); CHKERRQ(ierr); 164 PetscFunctionReturn(0); 165 } 166 167 168 169 #undef __FUNCT__ 170 #define __FUNCT__ "MatView_SMF" 171 PetscErrorCode MatView_SMF(Mat mat,PetscViewer viewer) 172 { 173 PetscErrorCode ierr; 174 MatSubMatFreeCtx ctx; 175 176 PetscFunctionBegin; 177 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 178 ierr = MatView(ctx->A,viewer);CHKERRQ(ierr); 179 PetscFunctionReturn(0); 180 } 181 182 #undef __FUNCT__ 183 #define __FUNCT__ "MatShift_SMF" 184 PetscErrorCode MatShift_SMF(Mat Y, PetscReal a) 185 { 186 PetscErrorCode ierr; 187 MatSubMatFreeCtx ctx; 188 189 PetscFunctionBegin; 190 ierr = MatShellGetContext(Y,(void **)&ctx);CHKERRQ(ierr); 191 ierr = MatShift(ctx->A,a);CHKERRQ(ierr); 192 PetscFunctionReturn(0); 193 } 194 195 #undef __FUNCT__ 196 #define __FUNCT__ "MatDuplicate_SMF" 197 PetscErrorCode MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M) 198 { 199 PetscErrorCode ierr; 200 MatSubMatFreeCtx ctx; 201 202 PetscFunctionBegin; 203 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 204 ierr = MatCreateSubMatrixFree(ctx->A,ctx->Rows,ctx->Cols,M);CHKERRQ(ierr); 205 PetscFunctionReturn(0); 206 } 207 208 #undef __FUNCT__ 209 #define __FUNCT__ "MatEqual_SMF" 210 PetscErrorCode MatEqual_SMF(Mat A,Mat B,PetscBool *flg) 211 { 212 PetscErrorCode ierr; 213 MatSubMatFreeCtx ctx1,ctx2; 214 PetscBool flg1,flg2,flg3; 215 216 PetscFunctionBegin; 217 ierr = MatShellGetContext(A,(void **)&ctx1);CHKERRQ(ierr); 218 ierr = MatShellGetContext(B,(void **)&ctx2);CHKERRQ(ierr); 219 ierr = ISEqual(ctx1->Rows,ctx2->Rows,&flg2);CHKERRQ(ierr); 220 ierr = ISEqual(ctx1->Cols,ctx2->Cols,&flg3);CHKERRQ(ierr); 221 if (flg2==PETSC_FALSE || flg3==PETSC_FALSE){ 222 *flg=PETSC_FALSE; 223 } else { 224 ierr = MatEqual(ctx1->A,ctx2->A,&flg1);CHKERRQ(ierr); 225 if (flg1==PETSC_FALSE){ *flg=PETSC_FALSE;} 226 else { *flg=PETSC_TRUE;} 227 } 228 PetscFunctionReturn(0); 229 } 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "MatScale_SMF" 233 PetscErrorCode MatScale_SMF(Mat mat, PetscReal a) 234 { 235 PetscErrorCode ierr; 236 MatSubMatFreeCtx ctx; 237 238 PetscFunctionBegin; 239 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 240 ierr = MatScale(ctx->A,a);CHKERRQ(ierr); 241 PetscFunctionReturn(0); 242 } 243 244 #undef __FUNCT__ 245 #define __FUNCT__ "MatTranspose_SMF" 246 PetscErrorCode MatTranspose_SMF(Mat mat,Mat *B) 247 { 248 PetscFunctionBegin; 249 PetscFunctionReturn(1); 250 } 251 252 #undef __FUNCT__ 253 #define __FUNCT__ "MatGetDiagonal_SMF" 254 PetscErrorCode MatGetDiagonal_SMF(Mat mat,Vec v) 255 { 256 PetscErrorCode ierr; 257 MatSubMatFreeCtx ctx; 258 259 PetscFunctionBegin; 260 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 261 ierr = MatGetDiagonal(ctx->A,v);CHKERRQ(ierr); 262 PetscFunctionReturn(0); 263 } 264 265 #undef __FUNCT__ 266 #define __FUNCT__ "MatDiagonalSet_SMF" 267 PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D) 268 { 269 MatSubMatFreeCtx ctx; 270 PetscErrorCode ierr; 271 272 PetscFunctionBegin; 273 ierr = MatShellGetContext(M,(void **)&ctx);CHKERRQ(ierr); 274 ierr = MatGetRowMax(ctx->A,D,PETSC_NULL); 275 PetscFunctionReturn(0); 276 } 277 278 279 #undef __FUNCT__ 280 #define __FUNCT__ "MatGetSubMatrices_SMF" 281 PetscErrorCode MatGetSubMatrices_SMF(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B) 282 { 283 PetscErrorCode ierr; 284 PetscInt i; 285 286 PetscFunctionBegin; 287 if (scall == MAT_INITIAL_MATRIX) { 288 ierr = PetscMalloc( (n+1)*sizeof(Mat),B );CHKERRQ(ierr); 289 } 290 291 for ( i=0; i<n; i++ ) { 292 ierr = MatGetSubMatrix_SMF(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 293 } 294 PetscFunctionReturn(0); 295 } 296 297 #undef __FUNCT__ 298 #define __FUNCT__ "MatGetSubMatrix_SMF" 299 PetscErrorCode MatGetSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll, 300 Mat *newmat) 301 { 302 PetscErrorCode ierr; 303 MatSubMatFreeCtx ctx; 304 305 PetscFunctionBegin; 306 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 307 if (newmat){ 308 ierr =MatDestroy(&*newmat);CHKERRQ(ierr); 309 } 310 ierr = MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat);CHKERRQ(ierr); 311 PetscFunctionReturn(0); 312 } 313 314 #undef __FUNCT__ 315 #define __FUNCT__ "MatGetRow_SMF" 316 PetscErrorCode MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscReal **vals) 317 { 318 PetscErrorCode ierr; 319 MatSubMatFreeCtx ctx; 320 321 PetscFunctionBegin; 322 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 323 ierr = MatGetRow(ctx->A,row,ncols,cols,vals);CHKERRQ(ierr); 324 PetscFunctionReturn(0); 325 } 326 327 #undef __FUNCT__ 328 #define __FUNCT__ "MatRestoreRow_SMF" 329 PetscErrorCode MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscReal **vals) 330 { 331 PetscErrorCode ierr; 332 MatSubMatFreeCtx ctx; 333 334 PetscFunctionBegin; 335 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 336 ierr = MatRestoreRow(ctx->A,row,ncols,cols,vals);CHKERRQ(ierr); 337 PetscFunctionReturn(0); 338 } 339 340 #undef __FUNCT__ 341 #define __FUNCT__ "MatGetColumnVector_SMF" 342 PetscErrorCode MatGetColumnVector_SMF(Mat mat,Vec Y, PetscInt col) 343 { 344 PetscErrorCode ierr; 345 MatSubMatFreeCtx ctx; 346 347 PetscFunctionBegin; 348 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 349 ierr = MatGetColumnVector(ctx->A,Y,col);CHKERRQ(ierr); 350 PetscFunctionReturn(0); 351 } 352 353 #undef __FUNCT__ 354 #define __FUNCT__ "MatConvert_SMF" 355 PetscErrorCode MatConvert_SMF(Mat mat,MatType newtype,Mat *NewMat) 356 { 357 PetscErrorCode ierr; 358 PetscMPIInt size; 359 MatSubMatFreeCtx ctx; 360 361 PetscFunctionBegin; 362 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 363 MPI_Comm_size(((PetscObject)mat)->comm,&size); 364 PetscFunctionReturn(1); 365 } 366 367 #undef __FUNCT__ 368 #define __FUNCT__ "MatNorm_SMF" 369 PetscErrorCode MatNorm_SMF(Mat mat,NormType type,PetscReal *norm) 370 { 371 PetscErrorCode ierr; 372 MatSubMatFreeCtx ctx; 373 374 PetscFunctionBegin; 375 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 376 377 if (type == NORM_FROBENIUS) { 378 *norm = 1.0; 379 } else if (type == NORM_1 || type == NORM_INFINITY) { 380 *norm = 1.0; 381 } else { 382 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 383 } 384 PetscFunctionReturn(0); 385 } 386 387