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