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