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