1 #include <../src/mat/impls/shell/shell.h> /*I "petscmat.h" I*/ 2 3 static PetscErrorCode MatMult_Transpose(Mat N, Vec x, Vec y) 4 { 5 Mat A; 6 7 PetscFunctionBegin; 8 PetscCall(MatShellGetContext(N, &A)); 9 PetscCall(MatMultTranspose(A, x, y)); 10 PetscFunctionReturn(PETSC_SUCCESS); 11 } 12 13 static PetscErrorCode MatMultTranspose_Transpose(Mat N, Vec x, Vec y) 14 { 15 Mat A; 16 17 PetscFunctionBegin; 18 PetscCall(MatShellGetContext(N, &A)); 19 PetscCall(MatMult(A, x, y)); 20 PetscFunctionReturn(PETSC_SUCCESS); 21 } 22 23 static PetscErrorCode MatDestroy_Transpose(Mat N) 24 { 25 Mat A; 26 27 PetscFunctionBegin; 28 PetscCall(MatShellGetContext(N, &A)); 29 PetscCall(MatDestroy(&A)); 30 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL)); 31 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL)); 32 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL)); 33 PetscFunctionReturn(PETSC_SUCCESS); 34 } 35 36 static PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat *m) 37 { 38 Mat A, C; 39 40 PetscFunctionBegin; 41 PetscCall(MatShellGetContext(N, &A)); 42 PetscCall(MatDuplicate(A, op, &C)); 43 PetscCall(MatCreateTranspose(C, m)); 44 PetscCall(MatDestroy(&C)); 45 if (op == MAT_COPY_VALUES) PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN)); 46 PetscFunctionReturn(PETSC_SUCCESS); 47 } 48 49 static PetscErrorCode MatHasOperation_Transpose(Mat mat, MatOperation op, PetscBool *has) 50 { 51 Mat A; 52 53 PetscFunctionBegin; 54 PetscCall(MatShellGetContext(mat, &A)); 55 *has = PETSC_FALSE; 56 if (op == MATOP_MULT || op == MATOP_MULT_ADD) { 57 PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has)); 58 } else if (op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) { 59 PetscCall(MatHasOperation(A, MATOP_MULT, has)); 60 } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE; 61 PetscFunctionReturn(PETSC_SUCCESS); 62 } 63 64 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat D) 65 { 66 Mat A, B, C, Ain, Bin, Cin; 67 PetscBool Aistrans, Bistrans, Cistrans; 68 PetscInt Atrans, Btrans, Ctrans; 69 MatProductType ptype; 70 71 PetscFunctionBegin; 72 MatCheckProduct(D, 1); 73 A = D->product->A; 74 B = D->product->B; 75 C = D->product->C; 76 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATTRANSPOSEVIRTUAL, &Aistrans)); 77 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &Bistrans)); 78 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATTRANSPOSEVIRTUAL, &Cistrans)); 79 PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen"); 80 Atrans = 0; 81 Ain = A; 82 while (Aistrans) { 83 Atrans++; 84 PetscCall(MatTransposeGetMat(Ain, &Ain)); 85 PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATTRANSPOSEVIRTUAL, &Aistrans)); 86 } 87 Btrans = 0; 88 Bin = B; 89 while (Bistrans) { 90 Btrans++; 91 PetscCall(MatTransposeGetMat(Bin, &Bin)); 92 PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATTRANSPOSEVIRTUAL, &Bistrans)); 93 } 94 Ctrans = 0; 95 Cin = C; 96 while (Cistrans) { 97 Ctrans++; 98 PetscCall(MatTransposeGetMat(Cin, &Cin)); 99 PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATTRANSPOSEVIRTUAL, &Cistrans)); 100 } 101 Atrans = Atrans % 2; 102 Btrans = Btrans % 2; 103 Ctrans = Ctrans % 2; 104 ptype = D->product->type; /* same product type by default */ 105 if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0; 106 if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0; 107 if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0; 108 109 if (Atrans || Btrans || Ctrans) { 110 ptype = MATPRODUCT_UNSPECIFIED; 111 switch (D->product->type) { 112 case MATPRODUCT_AB: 113 if (Atrans && Btrans) { /* At * Bt we do not have support for this */ 114 /* TODO custom implementation ? */ 115 } else if (Atrans) { /* At * B */ 116 ptype = MATPRODUCT_AtB; 117 } else { /* A * Bt */ 118 ptype = MATPRODUCT_ABt; 119 } 120 break; 121 case MATPRODUCT_AtB: 122 if (Atrans && Btrans) { /* A * Bt */ 123 ptype = MATPRODUCT_ABt; 124 } else if (Atrans) { /* A * B */ 125 ptype = MATPRODUCT_AB; 126 } else { /* At * Bt we do not have support for this */ 127 /* TODO custom implementation ? */ 128 } 129 break; 130 case MATPRODUCT_ABt: 131 if (Atrans && Btrans) { /* At * B */ 132 ptype = MATPRODUCT_AtB; 133 } else if (Atrans) { /* At * Bt we do not have support for this */ 134 /* TODO custom implementation ? */ 135 } else { /* A * B */ 136 ptype = MATPRODUCT_AB; 137 } 138 break; 139 case MATPRODUCT_PtAP: 140 if (Atrans) { /* PtAtP */ 141 /* TODO custom implementation ? */ 142 } else { /* RARt */ 143 ptype = MATPRODUCT_RARt; 144 } 145 break; 146 case MATPRODUCT_RARt: 147 if (Atrans) { /* RAtRt */ 148 /* TODO custom implementation ? */ 149 } else { /* PtAP */ 150 ptype = MATPRODUCT_PtAP; 151 } 152 break; 153 case MATPRODUCT_ABC: 154 /* TODO custom implementation ? */ 155 break; 156 default: 157 SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]); 158 } 159 } 160 PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D)); 161 PetscCall(MatProductSetType(D, ptype)); 162 PetscCall(MatProductSetFromOptions(D)); 163 PetscFunctionReturn(PETSC_SUCCESS); 164 } 165 166 static PetscErrorCode MatGetDiagonal_Transpose(Mat N, Vec v) 167 { 168 Mat A; 169 170 PetscFunctionBegin; 171 PetscCall(MatShellGetContext(N, &A)); 172 PetscCall(MatGetDiagonal(A, v)); 173 PetscFunctionReturn(PETSC_SUCCESS); 174 } 175 176 static PetscErrorCode MatCopy_Transpose(Mat A, Mat B, MatStructure str) 177 { 178 Mat a, b; 179 180 PetscFunctionBegin; 181 PetscCall(MatShellGetContext(A, &a)); 182 PetscCall(MatShellGetContext(B, &b)); 183 PetscCall(MatCopy(a, b, str)); 184 PetscFunctionReturn(PETSC_SUCCESS); 185 } 186 187 static PetscErrorCode MatConvert_Transpose(Mat N, MatType newtype, MatReuse reuse, Mat *newmat) 188 { 189 Mat A; 190 PetscScalar vscale = 1.0, vshift = 0.0; 191 PetscBool flg; 192 193 PetscFunctionBegin; 194 PetscCall(MatShellGetContext(N, &A)); 195 PetscCall(MatHasOperation(A, MATOP_TRANSPOSE, &flg)); 196 if (flg || N->ops->getrow) { /* if this condition is false, MatConvert_Shell() will be called in MatConvert_Basic(), so the following checks are not needed */ 197 PetscCheck(!((Mat_Shell *)N->data)->zrows && !((Mat_Shell *)N->data)->zcols, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat"); 198 PetscCheck(!((Mat_Shell *)N->data)->axpy, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatAXPY() has been called on the input Mat"); 199 PetscCheck(!((Mat_Shell *)N->data)->left && !((Mat_Shell *)N->data)->right, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatDiagonalScale() has been called on the input Mat"); 200 PetscCheck(!((Mat_Shell *)N->data)->dshift, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatConvert() if MatDiagonalSet() has been called on the input Mat"); 201 vscale = ((Mat_Shell *)N->data)->vscale; 202 vshift = ((Mat_Shell *)N->data)->vshift; 203 } 204 if (flg) { 205 Mat B; 206 207 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B)); 208 if (reuse != MAT_INPLACE_MATRIX) { 209 PetscCall(MatConvert(B, newtype, reuse, newmat)); 210 PetscCall(MatDestroy(&B)); 211 } else { 212 PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B)); 213 PetscCall(MatHeaderReplace(N, &B)); 214 } 215 } else { /* use basic converter as fallback */ 216 flg = (PetscBool)(N->ops->getrow != NULL); 217 PetscCall(MatConvert_Basic(N, newtype, reuse, newmat)); 218 } 219 if (flg) { 220 PetscCall(MatScale(*newmat, vscale)); 221 PetscCall(MatShift(*newmat, vshift)); 222 } 223 PetscFunctionReturn(PETSC_SUCCESS); 224 } 225 226 static PetscErrorCode MatTransposeGetMat_Transpose(Mat N, Mat *M) 227 { 228 PetscFunctionBegin; 229 PetscCall(MatShellGetContext(N, M)); 230 PetscFunctionReturn(PETSC_SUCCESS); 231 } 232 233 /*@ 234 MatTransposeGetMat - Gets the `Mat` object stored inside a `MATTRANSPOSEVIRTUAL` 235 236 Logically Collective 237 238 Input Parameter: 239 . A - the `MATTRANSPOSEVIRTUAL` matrix 240 241 Output Parameter: 242 . M - the matrix object stored inside `A` 243 244 Level: intermediate 245 246 .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()` 247 @*/ 248 PetscErrorCode MatTransposeGetMat(Mat A, Mat *M) 249 { 250 PetscFunctionBegin; 251 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 252 PetscValidType(A, 1); 253 PetscAssertPointer(M, 2); 254 PetscUseMethod(A, "MatTransposeGetMat_C", (Mat, Mat *), (A, M)); 255 PetscFunctionReturn(PETSC_SUCCESS); 256 } 257 258 /*MC 259 MATTRANSPOSEVIRTUAL - "transpose" - A matrix type that represents a virtual transpose of a matrix 260 261 Level: advanced 262 263 Developer Notes: 264 This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code 265 266 Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage 267 268 .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`, 269 `MATNORMALHERMITIAN`, `MATNORMAL` 270 M*/ 271 272 /*@ 273 MatCreateTranspose - Creates a new matrix `MATTRANSPOSEVIRTUAL` object that behaves like A' 274 275 Collective 276 277 Input Parameter: 278 . A - the (possibly rectangular) matrix 279 280 Output Parameter: 281 . N - the matrix that represents A' 282 283 Level: intermediate 284 285 Note: 286 The transpose A' is NOT actually formed! Rather the new matrix 287 object performs the matrix-vector product by using the `MatMultTranspose()` on 288 the original matrix 289 290 .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateNormal()`, `MatMult()`, `MatMultTranspose()`, `MatCreate()`, 291 `MATNORMALHERMITIAN` 292 @*/ 293 PetscErrorCode MatCreateTranspose(Mat A, Mat *N) 294 { 295 VecType vtype; 296 297 PetscFunctionBegin; 298 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N)); 299 PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap))); 300 PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap))); 301 PetscCall(MatSetType(*N, MATSHELL)); 302 PetscCall(MatShellSetContext(*N, A)); 303 PetscCall(PetscObjectReference((PetscObject)A)); 304 305 PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs))); 306 PetscCall(MatGetVecType(A, &vtype)); 307 PetscCall(MatSetVecType(*N, vtype)); 308 #if defined(PETSC_HAVE_DEVICE) 309 PetscCall(MatBindToCPU(*N, A->boundtocpu)); 310 #endif 311 PetscCall(MatSetUp(*N)); 312 313 PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_Transpose)); 314 PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_Transpose)); 315 PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_Transpose)); 316 PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_Transpose)); 317 PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (void (*)(void))MatHasOperation_Transpose)); 318 PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_Transpose)); 319 PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_Transpose)); 320 PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (void (*)(void))MatConvert_Transpose)); 321 322 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatTransposeGetMat_Transpose)); 323 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose)); 324 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable)); 325 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable)); 326 PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable)); 327 PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATTRANSPOSEVIRTUAL)); 328 PetscFunctionReturn(PETSC_SUCCESS); 329 } 330