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