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