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 MatHermitianTransposeGetMat_HT(Mat N, Mat *M) 165 { 166 PetscFunctionBegin; 167 PetscCall(MatShellGetContext(N, M)); 168 PetscFunctionReturn(PETSC_SUCCESS); 169 } 170 171 /*@ 172 MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSEVIRTUAL` 173 174 Logically Collective 175 176 Input Parameter: 177 . A - the `MATHERMITIANTRANSPOSEVIRTUAL` matrix 178 179 Output Parameter: 180 . M - the matrix object stored inside A 181 182 Level: intermediate 183 184 .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()` 185 @*/ 186 PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M) 187 { 188 PetscFunctionBegin; 189 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 190 PetscValidType(A, 1); 191 PetscAssertPointer(M, 2); 192 PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M)); 193 PetscFunctionReturn(PETSC_SUCCESS); 194 } 195 196 static PetscErrorCode MatGetDiagonal_HT(Mat N, Vec v) 197 { 198 Mat A; 199 200 PetscFunctionBegin; 201 PetscCall(MatShellGetContext(N, &A)); 202 PetscCall(MatGetDiagonal(A, v)); 203 PetscCall(VecConjugate(v)); 204 PetscFunctionReturn(PETSC_SUCCESS); 205 } 206 207 static PetscErrorCode MatConvert_HT(Mat N, MatType newtype, MatReuse reuse, Mat *newmat) 208 { 209 Mat A; 210 PetscBool flg; 211 212 PetscFunctionBegin; 213 PetscCall(MatShellGetContext(N, &A)); 214 PetscCall(MatHasOperation(A, MATOP_HERMITIAN_TRANSPOSE, &flg)); 215 if (flg) { 216 Mat B; 217 218 PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &B)); 219 if (reuse != MAT_INPLACE_MATRIX) { 220 PetscCall(MatConvert(B, newtype, reuse, newmat)); 221 PetscCall(MatDestroy(&B)); 222 } else { 223 PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B)); 224 PetscCall(MatHeaderReplace(N, &B)); 225 } 226 } else { /* use basic converter as fallback */ 227 PetscCall(MatConvert_Basic(N, newtype, reuse, newmat)); 228 } 229 PetscFunctionReturn(PETSC_SUCCESS); 230 } 231 232 /*MC 233 MATHERMITIANTRANSPOSEVIRTUAL - "hermitiantranspose" - A matrix type that represents a virtual transpose of a matrix 234 235 Level: advanced 236 237 .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()` 238 M*/ 239 240 /*@ 241 MatCreateHermitianTranspose - Creates a new matrix object of `MatType` `MATHERMITIANTRANSPOSEVIRTUAL` that behaves like A'* 242 243 Collective 244 245 Input Parameter: 246 . A - the (possibly rectangular) matrix 247 248 Output Parameter: 249 . N - the matrix that represents A'* 250 251 Level: intermediate 252 253 Note: 254 The Hermitian transpose A' is NOT actually formed! Rather the new matrix 255 object performs the matrix-vector product, `MatMult()`, by using the `MatMultHermitianTranspose()` on 256 the original matrix 257 258 .seealso: [](ch_matrices), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()`, 259 `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`, `MatHermitianTransposeGetMat()`, `MATNORMAL`, `MATNORMALHERMITIAN` 260 @*/ 261 PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N) 262 { 263 VecType vtype; 264 265 PetscFunctionBegin; 266 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N)); 267 PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap))); 268 PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap))); 269 PetscCall(MatSetType(*N, MATSHELL)); 270 PetscCall(MatShellSetContext(*N, A)); 271 PetscCall(PetscObjectReference((PetscObject)A)); 272 273 PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs))); 274 PetscCall(MatGetVecType(A, &vtype)); 275 PetscCall(MatSetVecType(*N, vtype)); 276 #if defined(PETSC_HAVE_DEVICE) 277 PetscCall(MatBindToCPU(*N, A->boundtocpu)); 278 #endif 279 PetscCall(MatSetUp(*N)); 280 281 PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_HT)); 282 PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_HT)); 283 PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT)); 284 #if !defined(PETSC_USE_COMPLEX) 285 PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT)); 286 #endif 287 PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_HT)); 288 PetscCall(MatShellSetOperation(*N, MATOP_CREATE_VECS, (void (*)(void))MatCreateVecs_HT)); 289 PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_HT)); 290 PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (void (*)(void))MatConvert_HT)); 291 292 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT)); 293 #if !defined(PETSC_USE_COMPLEX) 294 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT)); 295 #endif 296 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_HT)); 297 PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSEVIRTUAL)); 298 PetscFunctionReturn(PETSC_SUCCESS); 299 } 300