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 = PetscFree(N->data);CHKERRQ(ierr); 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 PetscErrorCode ierr; 63 64 PetscFunctionBegin; 65 if (op == MAT_COPY_VALUES) { 66 ierr = MatHermitianTranspose(Na->A,MAT_INITIAL_MATRIX,m);CHKERRQ(ierr); 67 } else if (op == MAT_DO_NOT_COPY_VALUES) { 68 ierr = MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);CHKERRQ(ierr); 69 ierr = MatHermitianTranspose(*m,MAT_INPLACE_MATRIX,m);CHKERRQ(ierr); 70 } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type"); 71 PetscFunctionReturn(0); 72 } 73 74 /*@ 75 MatCreateHermitianTranspose - Creates a new matrix object that behaves like A'* 76 77 Collective on Mat 78 79 Input Parameter: 80 . A - the (possibly rectangular) matrix 81 82 Output Parameter: 83 . N - the matrix that represents A'* 84 85 Level: intermediate 86 87 Notes: The hermitian transpose A' is NOT actually formed! Rather the new matrix 88 object performs the matrix-vector product by using the MatMultHermitianTranspose() on 89 the original matrix 90 91 .seealso: MatCreateNormal(), MatMult(), MatMultHermitianTranspose(), MatCreate() 92 93 @*/ 94 PetscErrorCode MatCreateHermitianTranspose(Mat A,Mat *N) 95 { 96 PetscErrorCode ierr; 97 PetscInt m,n; 98 Mat_HT *Na; 99 100 PetscFunctionBegin; 101 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 102 ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr); 103 ierr = MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 104 ierr = PetscLayoutSetUp((*N)->rmap);CHKERRQ(ierr); 105 ierr = PetscLayoutSetUp((*N)->cmap);CHKERRQ(ierr); 106 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);CHKERRQ(ierr); 107 108 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr); 109 (*N)->data = (void*) Na; 110 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 111 Na->A = A; 112 113 (*N)->ops->destroy = MatDestroy_HT; 114 (*N)->ops->mult = MatMult_HT; 115 (*N)->ops->multadd = MatMultAdd_HT; 116 (*N)->ops->multhermitiantranspose = MatMultHermitianTranspose_HT; 117 (*N)->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_HT; 118 (*N)->ops->duplicate = MatDuplicate_HT; 119 (*N)->assembled = PETSC_TRUE; 120 121 ierr = MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 122 ierr = MatSetUp(*N);CHKERRQ(ierr); 123 PetscFunctionReturn(0); 124 } 125