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