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