1 #include <petsc/private/matimpl.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 PetscFunctionReturn(PETSC_SUCCESS); 137 } 138 139 static PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat *m) 140 { 141 Mat A; 142 143 PetscFunctionBegin; 144 PetscCall(MatShellGetContext(N, &A)); 145 if (op == MAT_COPY_VALUES) { 146 PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, m)); 147 } else if (op == MAT_DO_NOT_COPY_VALUES) { 148 PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, m)); 149 PetscCall(MatHermitianTranspose(*m, MAT_INPLACE_MATRIX, m)); 150 } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "MAT_SHARE_NONZERO_PATTERN not supported for this matrix type"); 151 PetscFunctionReturn(PETSC_SUCCESS); 152 } 153 154 static PetscErrorCode MatCreateVecs_HT(Mat N, Vec *r, Vec *l) 155 { 156 Mat A; 157 158 PetscFunctionBegin; 159 PetscCall(MatShellGetContext(N, &A)); 160 PetscCall(MatCreateVecs(A, l, r)); 161 PetscFunctionReturn(PETSC_SUCCESS); 162 } 163 164 static PetscErrorCode MatHasOperation_HT(Mat mat, MatOperation op, PetscBool *has) 165 { 166 Mat A; 167 168 PetscFunctionBegin; 169 PetscCall(MatShellGetContext(mat, &A)); 170 *has = PETSC_FALSE; 171 if (op == MATOP_MULT || op == MATOP_MULT_ADD) { 172 PetscCall(MatHasOperation(A, MATOP_MULT_HERMITIAN_TRANSPOSE, has)); 173 if (!*has) PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has)); 174 } else if (op == MATOP_MULT_HERMITIAN_TRANSPOSE || op == MATOP_MULT_HERMITIAN_TRANS_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) { 175 PetscCall(MatHasOperation(A, MATOP_MULT, has)); 176 } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE; 177 PetscFunctionReturn(PETSC_SUCCESS); 178 } 179 180 static PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N, Mat *M) 181 { 182 PetscFunctionBegin; 183 PetscCall(MatShellGetContext(N, M)); 184 PetscFunctionReturn(PETSC_SUCCESS); 185 } 186 187 /*@ 188 MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSEVIRTUAL` 189 190 Logically Collective 191 192 Input Parameter: 193 . A - the `MATHERMITIANTRANSPOSEVIRTUAL` matrix 194 195 Output Parameter: 196 . M - the matrix object stored inside A 197 198 Level: intermediate 199 200 .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()` 201 @*/ 202 PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M) 203 { 204 PetscFunctionBegin; 205 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 206 PetscValidType(A, 1); 207 PetscAssertPointer(M, 2); 208 PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M)); 209 PetscFunctionReturn(PETSC_SUCCESS); 210 } 211 212 static PetscErrorCode MatGetDiagonal_HT(Mat N, Vec v) 213 { 214 Mat A; 215 216 PetscFunctionBegin; 217 PetscCall(MatShellGetContext(N, &A)); 218 PetscCall(MatGetDiagonal(A, v)); 219 PetscCall(VecConjugate(v)); 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 .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()` 254 M*/ 255 256 /*@ 257 MatCreateHermitianTranspose - Creates a new matrix object of `MatType` `MATHERMITIANTRANSPOSEVIRTUAL` that behaves like A'* 258 259 Collective 260 261 Input Parameter: 262 . A - the (possibly rectangular) matrix 263 264 Output Parameter: 265 . N - the matrix that represents A'* 266 267 Level: intermediate 268 269 Note: 270 The Hermitian transpose A' is NOT actually formed! Rather the new matrix 271 object performs the matrix-vector product, `MatMult()`, by using the `MatMultHermitianTranspose()` on 272 the original matrix 273 274 .seealso: [](ch_matrices), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()`, 275 `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`, `MatHermitianTransposeGetMat()`, `MATNORMAL`, `MATNORMALHERMITIAN` 276 @*/ 277 PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N) 278 { 279 VecType vtype; 280 281 PetscFunctionBegin; 282 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N)); 283 PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap))); 284 PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap))); 285 PetscCall(MatSetType(*N, MATSHELL)); 286 PetscCall(MatShellSetContext(*N, A)); 287 PetscCall(PetscObjectReference((PetscObject)A)); 288 289 PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs))); 290 PetscCall(MatGetVecType(A, &vtype)); 291 PetscCall(MatSetVecType(*N, vtype)); 292 #if defined(PETSC_HAVE_DEVICE) 293 PetscCall(MatBindToCPU(*N, A->boundtocpu)); 294 #endif 295 PetscCall(MatSetUp(*N)); 296 297 PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_HT)); 298 PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_HT)); 299 PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT)); 300 #if !defined(PETSC_USE_COMPLEX) 301 PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT)); 302 #endif 303 PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_HT)); 304 PetscCall(MatShellSetOperation(*N, MATOP_CREATE_VECS, (void (*)(void))MatCreateVecs_HT)); 305 PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (void (*)(void))MatHasOperation_HT)); 306 PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_HT)); 307 PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (void (*)(void))MatConvert_HT)); 308 309 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT)); 310 #if !defined(PETSC_USE_COMPLEX) 311 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT)); 312 #endif 313 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_HT)); 314 PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSEVIRTUAL)); 315 PetscFunctionReturn(PETSC_SUCCESS); 316 } 317