1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 typedef struct { 5 Mat A; 6 } Mat_HT; 7 8 PetscErrorCode MatMult_HT(Mat N, Vec x, Vec y) { 9 Mat_HT *Na = (Mat_HT *)N->data; 10 11 PetscFunctionBegin; 12 PetscCall(MatMultHermitianTranspose(Na->A, x, y)); 13 PetscFunctionReturn(0); 14 } 15 16 PetscErrorCode MatMultAdd_HT(Mat N, Vec v1, Vec v2, Vec v3) { 17 Mat_HT *Na = (Mat_HT *)N->data; 18 19 PetscFunctionBegin; 20 PetscCall(MatMultHermitianTransposeAdd(Na->A, v1, v2, v3)); 21 PetscFunctionReturn(0); 22 } 23 24 PetscErrorCode MatMultHermitianTranspose_HT(Mat N, Vec x, Vec y) { 25 Mat_HT *Na = (Mat_HT *)N->data; 26 27 PetscFunctionBegin; 28 PetscCall(MatMult(Na->A, x, y)); 29 PetscFunctionReturn(0); 30 } 31 32 PetscErrorCode MatMultHermitianTransposeAdd_HT(Mat N, Vec v1, Vec v2, Vec v3) { 33 Mat_HT *Na = (Mat_HT *)N->data; 34 35 PetscFunctionBegin; 36 PetscCall(MatMultAdd(Na->A, v1, v2, v3)); 37 PetscFunctionReturn(0); 38 } 39 40 PetscErrorCode MatDestroy_HT(Mat N) { 41 Mat_HT *Na = (Mat_HT *)N->data; 42 43 PetscFunctionBegin; 44 PetscCall(MatDestroy(&Na->A)); 45 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatHermitianTransposeGetMat_C", NULL)); 46 #if !defined(PETSC_USE_COMPLEX) 47 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL)); 48 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL)); 49 #endif 50 PetscCall(PetscFree(N->data)); 51 PetscFunctionReturn(0); 52 } 53 54 PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat *m) { 55 Mat_HT *Na = (Mat_HT *)N->data; 56 57 PetscFunctionBegin; 58 if (op == MAT_COPY_VALUES) { 59 PetscCall(MatHermitianTranspose(Na->A, MAT_INITIAL_MATRIX, m)); 60 } else if (op == MAT_DO_NOT_COPY_VALUES) { 61 PetscCall(MatDuplicate(Na->A, MAT_DO_NOT_COPY_VALUES, m)); 62 PetscCall(MatHermitianTranspose(*m, MAT_INPLACE_MATRIX, m)); 63 } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "MAT_SHARE_NONZERO_PATTERN not supported for this matrix type"); 64 PetscFunctionReturn(0); 65 } 66 67 PetscErrorCode MatCreateVecs_HT(Mat N, Vec *r, Vec *l) { 68 Mat_HT *Na = (Mat_HT *)N->data; 69 70 PetscFunctionBegin; 71 PetscCall(MatCreateVecs(Na->A, l, r)); 72 PetscFunctionReturn(0); 73 } 74 75 PetscErrorCode MatAXPY_HT(Mat Y, PetscScalar a, Mat X, MatStructure str) { 76 Mat_HT *Ya = (Mat_HT *)Y->data; 77 Mat_HT *Xa = (Mat_HT *)X->data; 78 Mat M = Ya->A; 79 Mat N = Xa->A; 80 81 PetscFunctionBegin; 82 PetscCall(MatAXPY(M, a, N, str)); 83 PetscFunctionReturn(0); 84 } 85 86 PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N, Mat *M) { 87 Mat_HT *Na = (Mat_HT *)N->data; 88 89 PetscFunctionBegin; 90 *M = Na->A; 91 PetscFunctionReturn(0); 92 } 93 94 /*@ 95 MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSE` 96 97 Logically collective on Mat 98 99 Input Parameter: 100 . A - the `MATHERMITIANTRANSPOSE` matrix 101 102 Output Parameter: 103 . M - the matrix object stored inside A 104 105 Level: intermediate 106 107 .seealso: `MATHERMITIANTRANSPOSE`, `MatCreateHermitianTranspose()` 108 @*/ 109 PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M) { 110 PetscFunctionBegin; 111 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 112 PetscValidType(A, 1); 113 PetscValidPointer(M, 2); 114 PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M)); 115 PetscFunctionReturn(0); 116 } 117 118 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat); 119 120 PetscErrorCode MatGetDiagonal_HT(Mat A, Vec v) { 121 Mat_HT *Na = (Mat_HT *)A->data; 122 123 PetscFunctionBegin; 124 PetscCall(MatGetDiagonal(Na->A, v)); 125 PetscCall(VecConjugate(v)); 126 PetscFunctionReturn(0); 127 } 128 129 PetscErrorCode MatConvert_HT(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) { 130 Mat_HT *Na = (Mat_HT *)A->data; 131 PetscBool flg; 132 133 PetscFunctionBegin; 134 PetscCall(MatHasOperation(Na->A, MATOP_HERMITIAN_TRANSPOSE, &flg)); 135 if (flg) { 136 Mat B; 137 138 PetscCall(MatHermitianTranspose(Na->A, MAT_INITIAL_MATRIX, &B)); 139 if (reuse != MAT_INPLACE_MATRIX) { 140 PetscCall(MatConvert(B, newtype, reuse, newmat)); 141 PetscCall(MatDestroy(&B)); 142 } else { 143 PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B)); 144 PetscCall(MatHeaderReplace(A, &B)); 145 } 146 } else { /* use basic converter as fallback */ 147 PetscCall(MatConvert_Basic(A, newtype, reuse, newmat)); 148 } 149 PetscFunctionReturn(0); 150 } 151 152 /*MC 153 MATHERMITIANTRANSPOSE - "hermitiantranspose" - A matrix type that represents a virtual transpose of a matrix 154 155 Level: advanced 156 157 .seealso: `MATTRANSPOSE`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()` 158 M*/ 159 160 /*@ 161 MatCreateHermitianTranspose - Creates a new matrix object of `MatType` `MATHERMITIANTRANSPOSE` that behaves like A'* 162 163 Collective on A 164 165 Input Parameter: 166 . A - the (possibly rectangular) matrix 167 168 Output Parameter: 169 . N - the matrix that represents A'* 170 171 Level: intermediate 172 173 Note: 174 The Hermitian transpose A' is NOT actually formed! Rather the new matrix 175 object performs the matrix-vector product, `MatMult()`, by using the `MatMultHermitianTranspose()` on 176 the original matrix 177 178 .seealso: `MATHERMITIANTRANSPOSE`, `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()` 179 @*/ 180 PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N) { 181 PetscInt m, n; 182 Mat_HT *Na; 183 VecType vtype; 184 185 PetscFunctionBegin; 186 PetscCall(MatGetLocalSize(A, &m, &n)); 187 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N)); 188 PetscCall(MatSetSizes(*N, n, m, PETSC_DECIDE, PETSC_DECIDE)); 189 PetscCall(PetscLayoutSetUp((*N)->rmap)); 190 PetscCall(PetscLayoutSetUp((*N)->cmap)); 191 PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSE)); 192 193 PetscCall(PetscNewLog(*N, &Na)); 194 (*N)->data = (void *)Na; 195 PetscCall(PetscObjectReference((PetscObject)A)); 196 Na->A = A; 197 198 (*N)->ops->destroy = MatDestroy_HT; 199 (*N)->ops->mult = MatMult_HT; 200 (*N)->ops->multadd = MatMultAdd_HT; 201 (*N)->ops->multhermitiantranspose = MatMultHermitianTranspose_HT; 202 (*N)->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_HT; 203 #if !defined(PETSC_USE_COMPLEX) 204 (*N)->ops->multtranspose = MatMultHermitianTranspose_HT; 205 (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_HT; 206 #endif 207 (*N)->ops->duplicate = MatDuplicate_HT; 208 (*N)->ops->getvecs = MatCreateVecs_HT; 209 (*N)->ops->axpy = MatAXPY_HT; 210 #if !defined(PETSC_USE_COMPLEX) 211 (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose; 212 #endif 213 (*N)->ops->getdiagonal = MatGetDiagonal_HT; 214 (*N)->ops->convert = MatConvert_HT; 215 (*N)->assembled = PETSC_TRUE; 216 217 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT)); 218 #if !defined(PETSC_USE_COMPLEX) 219 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT)); 220 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose)); 221 #endif 222 PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs))); 223 PetscCall(MatGetVecType(A, &vtype)); 224 PetscCall(MatSetVecType(*N, vtype)); 225 #if defined(PETSC_HAVE_DEVICE) 226 PetscCall(MatBindToCPU(*N, A->boundtocpu)); 227 #endif 228 PetscCall(MatSetUp(*N)); 229 PetscFunctionReturn(0); 230 } 231