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 #endif 59 ierr = PetscFree(N->data);CHKERRQ(ierr); 60 PetscFunctionReturn(0); 61 } 62 63 PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat* m) 64 { 65 Mat_HT *Na = (Mat_HT*)N->data; 66 PetscErrorCode ierr; 67 68 PetscFunctionBegin; 69 if (op == MAT_COPY_VALUES) { 70 ierr = MatHermitianTranspose(Na->A,MAT_INITIAL_MATRIX,m);CHKERRQ(ierr); 71 } else if (op == MAT_DO_NOT_COPY_VALUES) { 72 ierr = MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);CHKERRQ(ierr); 73 ierr = MatHermitianTranspose(*m,MAT_INPLACE_MATRIX,m);CHKERRQ(ierr); 74 } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type"); 75 PetscFunctionReturn(0); 76 } 77 78 PetscErrorCode MatCreateVecs_HT(Mat N,Vec *r, Vec *l) 79 { 80 Mat_HT *Na = (Mat_HT*)N->data; 81 PetscErrorCode ierr; 82 83 PetscFunctionBegin; 84 ierr = MatCreateVecs(Na->A,l,r);CHKERRQ(ierr); 85 PetscFunctionReturn(0); 86 } 87 88 PetscErrorCode MatAXPY_HT(Mat Y,PetscScalar a,Mat X,MatStructure str) 89 { 90 Mat_HT *Ya = (Mat_HT*)Y->data; 91 Mat_HT *Xa = (Mat_HT*)X->data; 92 Mat M = Ya->A; 93 Mat N = Xa->A; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 ierr = MatAXPY(M,a,N,str);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N,Mat *M) 102 { 103 Mat_HT *Na = (Mat_HT*)N->data; 104 105 PetscFunctionBegin; 106 *M = Na->A; 107 PetscFunctionReturn(0); 108 } 109 110 /*@ 111 MatHermitianTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT 112 113 Logically collective on Mat 114 115 Input Parameter: 116 . A - the MATTRANSPOSE matrix 117 118 Output Parameter: 119 . M - the matrix object stored inside A 120 121 Level: intermediate 122 123 .seealso: MatCreateHermitianTranspose() 124 125 @*/ 126 PetscErrorCode MatHermitianTransposeGetMat(Mat A,Mat *M) 127 { 128 PetscErrorCode ierr; 129 130 PetscFunctionBegin; 131 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 132 PetscValidType(A,1); 133 PetscValidPointer(M,2); 134 ierr = PetscUseMethod(A,"MatHermitianTransposeGetMat_C",(Mat,Mat*),(A,M));CHKERRQ(ierr); 135 PetscFunctionReturn(0); 136 } 137 138 /*@ 139 MatCreateHermitianTranspose - Creates a new matrix object that behaves like A'* 140 141 Collective on Mat 142 143 Input Parameter: 144 . A - the (possibly rectangular) matrix 145 146 Output Parameter: 147 . N - the matrix that represents A'* 148 149 Level: intermediate 150 151 Notes: 152 The hermitian transpose A' is NOT actually formed! Rather the new matrix 153 object performs the matrix-vector product by using the MatMultHermitianTranspose() on 154 the original matrix 155 156 .seealso: MatCreateNormal(), MatMult(), MatMultHermitianTranspose(), MatCreate() 157 158 @*/ 159 PetscErrorCode MatCreateHermitianTranspose(Mat A,Mat *N) 160 { 161 PetscErrorCode ierr; 162 PetscInt m,n; 163 Mat_HT *Na; 164 165 PetscFunctionBegin; 166 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 167 ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr); 168 ierr = MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 169 ierr = PetscLayoutSetUp((*N)->rmap);CHKERRQ(ierr); 170 ierr = PetscLayoutSetUp((*N)->cmap);CHKERRQ(ierr); 171 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);CHKERRQ(ierr); 172 173 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr); 174 (*N)->data = (void*) Na; 175 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 176 Na->A = A; 177 178 (*N)->ops->destroy = MatDestroy_HT; 179 (*N)->ops->mult = MatMult_HT; 180 (*N)->ops->multadd = MatMultAdd_HT; 181 (*N)->ops->multhermitiantranspose = MatMultHermitianTranspose_HT; 182 (*N)->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_HT; 183 (*N)->ops->duplicate = MatDuplicate_HT; 184 (*N)->ops->getvecs = MatCreateVecs_HT; 185 (*N)->ops->axpy = MatAXPY_HT; 186 (*N)->assembled = PETSC_TRUE; 187 188 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatHermitianTransposeGetMat_C",MatHermitianTransposeGetMat_HT);CHKERRQ(ierr); 189 #if !defined(PETSC_USE_COMPLEX) 190 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatHermitianTransposeGetMat_HT);CHKERRQ(ierr); 191 #endif 192 ierr = MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 193 ierr = MatSetUp(*N);CHKERRQ(ierr); 194 PetscFunctionReturn(0); 195 } 196