1 2 #include <petsc-private/matimpl.h> /*I "petscmat.h" I*/ 3 4 typedef struct { 5 Mat A; 6 } Mat_HTranspose; 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "MatMult_HTranspose" 10 PetscErrorCode MatMult_HTranspose(Mat N,Vec x,Vec y) 11 { 12 Mat_HTranspose *Na = (Mat_HTranspose*)N->data; 13 PetscErrorCode ierr; 14 15 PetscFunctionBegin; 16 ierr = MatMultHermitianTranspose(Na->A,x,y);CHKERRQ(ierr); 17 PetscFunctionReturn(0); 18 } 19 20 #undef __FUNCT__ 21 #define __FUNCT__ "MatMultAdd_HTranspose" 22 PetscErrorCode MatMultAdd_HTranspose(Mat N,Vec v1,Vec v2,Vec v3) 23 { 24 Mat_HTranspose *Na = (Mat_HTranspose*)N->data; 25 PetscErrorCode ierr; 26 27 PetscFunctionBegin; 28 ierr = MatMultHermitianTransposeAdd(Na->A,v1,v2,v3);CHKERRQ(ierr); 29 PetscFunctionReturn(0); 30 } 31 32 #undef __FUNCT__ 33 #define __FUNCT__ "MatMultHTranspose_HTranspose" 34 PetscErrorCode MatMultHTranspose_HTranspose(Mat N,Vec x,Vec y) 35 { 36 Mat_HTranspose *Na = (Mat_HTranspose*)N->data; 37 PetscErrorCode ierr; 38 39 PetscFunctionBegin; 40 ierr = MatMult(Na->A,x,y);CHKERRQ(ierr); 41 PetscFunctionReturn(0); 42 } 43 44 #undef __FUNCT__ 45 #define __FUNCT__ "MatMultHTransposeAdd_HTranspose" 46 PetscErrorCode MatMultHTransposeAdd_HTranspose(Mat N,Vec v1,Vec v2,Vec v3) 47 { 48 Mat_HTranspose *Na = (Mat_HTranspose*)N->data; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 ierr = MatMultAdd(Na->A,v1,v2,v3);CHKERRQ(ierr); 53 PetscFunctionReturn(0); 54 } 55 56 #undef __FUNCT__ 57 #define __FUNCT__ "MatDestroy_HTranspose" 58 PetscErrorCode MatDestroy_HTranspose(Mat N) 59 { 60 Mat_HTranspose *Na = (Mat_HTranspose*)N->data; 61 PetscErrorCode ierr; 62 63 PetscFunctionBegin; 64 ierr = MatDestroy(&Na->A);CHKERRQ(ierr); 65 ierr = PetscFree(N->data);CHKERRQ(ierr); 66 PetscFunctionReturn(0); 67 } 68 69 #undef __FUNCT__ 70 #define __FUNCT__ "MatDuplicate_HTranspose" 71 PetscErrorCode MatDuplicate_HTranspose(Mat N, MatDuplicateOption op, Mat* m) 72 { 73 Mat_HTranspose *Na = (Mat_HTranspose*)N->data; 74 PetscErrorCode ierr; 75 76 PetscFunctionBegin; 77 if (op == MAT_COPY_VALUES) { 78 ierr = MatHermitianTranspose(Na->A,MAT_INITIAL_MATRIX,m);CHKERRQ(ierr); 79 } else if (op == MAT_DO_NOT_COPY_VALUES) { 80 ierr = MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);CHKERRQ(ierr); 81 ierr = MatTranspose(*m,MAT_REUSE_MATRIX,m);CHKERRQ(ierr); 82 } else { 83 SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type"); 84 } 85 PetscFunctionReturn(0); 86 } 87 88 #undef __FUNCT__ 89 #define __FUNCT__ "MatCreateHermitianTranspose" 90 /*@ 91 MatCreateHermitianTranspose - Creates a new matrix object that behaves like A'* 92 93 Collective on Mat 94 95 Input Parameter: 96 . A - the (possibly rectangular) matrix 97 98 Output Parameter: 99 . N - the matrix that represents A'* 100 101 Level: intermediate 102 103 Notes: The hermitian transpose A' is NOT actually formed! Rather the new matrix 104 object performs the matrix-vector product by using the MatMultHermitianTranspose() on 105 the original matrix 106 107 .seealso: MatCreateNormal(), MatMult(), MatMultHermitianTranspose(), MatCreate() 108 109 @*/ 110 PetscErrorCode MatCreateHermitianTranspose(Mat A,Mat *N) 111 { 112 PetscErrorCode ierr; 113 PetscInt m,n; 114 Mat_HTranspose *Na; 115 116 PetscFunctionBegin; 117 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 118 ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr); 119 ierr = MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 120 ierr = PetscLayoutSetUp((*N)->rmap);CHKERRQ(ierr); 121 ierr = PetscLayoutSetUp((*N)->cmap);CHKERRQ(ierr); 122 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);CHKERRQ(ierr); 123 124 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr); 125 (*N)->data = (void*) Na; 126 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 127 Na->A = A; 128 129 (*N)->ops->destroy = MatDestroy_HTranspose; 130 (*N)->ops->mult = MatMult_HTranspose; 131 (*N)->ops->multadd = MatMultAdd_HTranspose; 132 (*N)->ops->multhermitiantranspose = MatMultHTranspose_HTranspose; 133 (*N)->ops->multhermitiantransposeadd = MatMultHTransposeAdd_HTranspose; 134 (*N)->ops->duplicate = MatDuplicate_HTranspose; 135 (*N)->assembled = PETSC_TRUE; 136 137 ierr = MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 138 ierr = MatSetUp(*N);CHKERRQ(ierr); 139 PetscFunctionReturn(0); 140 } 141