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