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