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 MATTRANSPOSEMAT 96 97 Logically collective on Mat 98 99 Input Parameter: 100 . A - the MATTRANSPOSE matrix 101 102 Output Parameter: 103 . M - the matrix object stored inside A 104 105 Level: intermediate 106 107 .seealso: `MatCreateHermitianTranspose()` 108 109 @*/ 110 PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M) { 111 PetscFunctionBegin; 112 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 113 PetscValidType(A, 1); 114 PetscValidPointer(M, 2); 115 PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M)); 116 PetscFunctionReturn(0); 117 } 118 119 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat); 120 121 PetscErrorCode MatGetDiagonal_HT(Mat A, Vec v) { 122 Mat_HT *Na = (Mat_HT *)A->data; 123 124 PetscFunctionBegin; 125 PetscCall(MatGetDiagonal(Na->A, v)); 126 PetscCall(VecConjugate(v)); 127 PetscFunctionReturn(0); 128 } 129 130 PetscErrorCode MatConvert_HT(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) { 131 Mat_HT *Na = (Mat_HT *)A->data; 132 PetscBool flg; 133 134 PetscFunctionBegin; 135 PetscCall(MatHasOperation(Na->A, MATOP_HERMITIAN_TRANSPOSE, &flg)); 136 if (flg) { 137 Mat B; 138 139 PetscCall(MatHermitianTranspose(Na->A, MAT_INITIAL_MATRIX, &B)); 140 if (reuse != MAT_INPLACE_MATRIX) { 141 PetscCall(MatConvert(B, newtype, reuse, newmat)); 142 PetscCall(MatDestroy(&B)); 143 } else { 144 PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B)); 145 PetscCall(MatHeaderReplace(A, &B)); 146 } 147 } else { /* use basic converter as fallback */ 148 PetscCall(MatConvert_Basic(A, newtype, reuse, newmat)); 149 } 150 PetscFunctionReturn(0); 151 } 152 153 /*@ 154 MatCreateHermitianTranspose - Creates a new matrix object that behaves like A'* 155 156 Collective on Mat 157 158 Input Parameter: 159 . A - the (possibly rectangular) matrix 160 161 Output Parameter: 162 . N - the matrix that represents A'* 163 164 Level: intermediate 165 166 Notes: 167 The hermitian transpose A' is NOT actually formed! Rather the new matrix 168 object performs the matrix-vector product by using the MatMultHermitianTranspose() on 169 the original matrix 170 171 .seealso: `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()` 172 @*/ 173 PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N) { 174 PetscInt m, n; 175 Mat_HT *Na; 176 VecType vtype; 177 178 PetscFunctionBegin; 179 PetscCall(MatGetLocalSize(A, &m, &n)); 180 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N)); 181 PetscCall(MatSetSizes(*N, n, m, PETSC_DECIDE, PETSC_DECIDE)); 182 PetscCall(PetscLayoutSetUp((*N)->rmap)); 183 PetscCall(PetscLayoutSetUp((*N)->cmap)); 184 PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATTRANSPOSEMAT)); 185 186 PetscCall(PetscNewLog(*N, &Na)); 187 (*N)->data = (void *)Na; 188 PetscCall(PetscObjectReference((PetscObject)A)); 189 Na->A = A; 190 191 (*N)->ops->destroy = MatDestroy_HT; 192 (*N)->ops->mult = MatMult_HT; 193 (*N)->ops->multadd = MatMultAdd_HT; 194 (*N)->ops->multhermitiantranspose = MatMultHermitianTranspose_HT; 195 (*N)->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_HT; 196 #if !defined(PETSC_USE_COMPLEX) 197 (*N)->ops->multtranspose = MatMultHermitianTranspose_HT; 198 (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_HT; 199 #endif 200 (*N)->ops->duplicate = MatDuplicate_HT; 201 (*N)->ops->getvecs = MatCreateVecs_HT; 202 (*N)->ops->axpy = MatAXPY_HT; 203 #if !defined(PETSC_USE_COMPLEX) 204 (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose; 205 #endif 206 (*N)->ops->getdiagonal = MatGetDiagonal_HT; 207 (*N)->ops->convert = MatConvert_HT; 208 (*N)->assembled = PETSC_TRUE; 209 210 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT)); 211 #if !defined(PETSC_USE_COMPLEX) 212 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT)); 213 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose)); 214 #endif 215 PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs))); 216 PetscCall(MatGetVecType(A, &vtype)); 217 PetscCall(MatSetVecType(*N, vtype)); 218 #if defined(PETSC_HAVE_DEVICE) 219 PetscCall(MatBindToCPU(*N, A->boundtocpu)); 220 #endif 221 PetscCall(MatSetUp(*N)); 222 PetscFunctionReturn(0); 223 } 224