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