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