1
2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
3
4 typedef struct {
5 Mat A;
6 } Mat_HT;
7
MatMult_HT(Mat N,Vec x,Vec y)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
MatMultAdd_HT(Mat N,Vec v1,Vec v2,Vec v3)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
MatMultHermitianTranspose_HT(Mat N,Vec x,Vec y)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
MatMultHermitianTransposeAdd_HT(Mat N,Vec v1,Vec v2,Vec v3)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
MatDestroy_HT(Mat N)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
MatDuplicate_HT(Mat N,MatDuplicateOption op,Mat * m)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
MatCreateVecs_HT(Mat N,Vec * r,Vec * l)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
MatAXPY_HT(Mat Y,PetscScalar a,Mat X,MatStructure str)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
MatHermitianTransposeGetMat_HT(Mat N,Mat * M)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 @*/
MatHermitianTransposeGetMat(Mat A,Mat * M)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 /*@
142 MatCreateHermitianTranspose - Creates a new matrix object that behaves like A'*
143
144 Collective on Mat
145
146 Input Parameter:
147 . A - the (possibly rectangular) matrix
148
149 Output Parameter:
150 . N - the matrix that represents A'*
151
152 Level: intermediate
153
154 Notes:
155 The hermitian transpose A' is NOT actually formed! Rather the new matrix
156 object performs the matrix-vector product by using the MatMultHermitianTranspose() on
157 the original matrix
158
159 .seealso: MatCreateNormal(), MatMult(), MatMultHermitianTranspose(), MatCreate()
160
161 @*/
MatCreateHermitianTranspose(Mat A,Mat * N)162 PetscErrorCode MatCreateHermitianTranspose(Mat A,Mat *N)
163 {
164 PetscErrorCode ierr;
165 PetscInt m,n;
166 Mat_HT *Na;
167 VecType vtype;
168
169 PetscFunctionBegin;
170 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
171 ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr);
172 ierr = MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
173 ierr = PetscLayoutSetUp((*N)->rmap);CHKERRQ(ierr);
174 ierr = PetscLayoutSetUp((*N)->cmap);CHKERRQ(ierr);
175 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);CHKERRQ(ierr);
176
177 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr);
178 (*N)->data = (void*) Na;
179 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
180 Na->A = A;
181
182 (*N)->ops->destroy = MatDestroy_HT;
183 (*N)->ops->mult = MatMult_HT;
184 (*N)->ops->multadd = MatMultAdd_HT;
185 (*N)->ops->multhermitiantranspose = MatMultHermitianTranspose_HT;
186 (*N)->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_HT;
187 (*N)->ops->duplicate = MatDuplicate_HT;
188 (*N)->ops->getvecs = MatCreateVecs_HT;
189 (*N)->ops->axpy = MatAXPY_HT;
190 #if !defined(PETSC_USE_COMPLEX)
191 (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
192 #endif
193 (*N)->assembled = PETSC_TRUE;
194
195 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatHermitianTransposeGetMat_C",MatHermitianTransposeGetMat_HT);CHKERRQ(ierr);
196 #if !defined(PETSC_USE_COMPLEX)
197 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatHermitianTransposeGetMat_HT);CHKERRQ(ierr);
198 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_Transpose);CHKERRQ(ierr);
199 #endif
200 ierr = MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr);
201 ierr = MatGetVecType(A,&vtype);CHKERRQ(ierr);
202 ierr = MatSetVecType(*N,vtype);CHKERRQ(ierr);
203 ierr = MatSetUp(*N);CHKERRQ(ierr);
204 PetscFunctionReturn(0);
205 }
206