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