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