Actual source code: htransm.c
2: #include <petsc/private/matimpl.h>
4: typedef struct {
5: Mat A;
6: } Mat_HT;
8: PetscErrorCode MatMult_HT(Mat N,Vec x,Vec y)
9: {
10: Mat_HT *Na = (Mat_HT*)N->data;
14: MatMultHermitianTranspose(Na->A,x,y);
15: return(0);
16: }
18: PetscErrorCode MatMultAdd_HT(Mat N,Vec v1,Vec v2,Vec v3)
19: {
20: Mat_HT *Na = (Mat_HT*)N->data;
24: MatMultHermitianTransposeAdd(Na->A,v1,v2,v3);
25: return(0);
26: }
28: PetscErrorCode MatMultHermitianTranspose_HT(Mat N,Vec x,Vec y)
29: {
30: Mat_HT *Na = (Mat_HT*)N->data;
34: MatMult(Na->A,x,y);
35: return(0);
36: }
38: PetscErrorCode MatMultHermitianTransposeAdd_HT(Mat N,Vec v1,Vec v2,Vec v3)
39: {
40: Mat_HT *Na = (Mat_HT*)N->data;
44: MatMultAdd(Na->A,v1,v2,v3);
45: return(0);
46: }
48: PetscErrorCode MatDestroy_HT(Mat N)
49: {
50: Mat_HT *Na = (Mat_HT*)N->data;
54: MatDestroy(&Na->A);
55: PetscObjectComposeFunction((PetscObject)N,"MatHermitianTransposeGetMat_C",NULL);
56: #if !defined(PETSC_USE_COMPLEX)
57: PetscObjectComposeFunction((PetscObject)N,"MatTransposeGetMat_C",NULL);
58: PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_anytype_C",NULL);
59: #endif
60: PetscFree(N->data);
61: return(0);
62: }
64: PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat* m)
65: {
66: Mat_HT *Na = (Mat_HT*)N->data;
70: if (op == MAT_COPY_VALUES) {
71: MatHermitianTranspose(Na->A,MAT_INITIAL_MATRIX,m);
72: } else if (op == MAT_DO_NOT_COPY_VALUES) {
73: MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);
74: MatHermitianTranspose(*m,MAT_INPLACE_MATRIX,m);
75: } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
76: return(0);
77: }
79: PetscErrorCode MatCreateVecs_HT(Mat N,Vec *r, Vec *l)
80: {
81: Mat_HT *Na = (Mat_HT*)N->data;
85: MatCreateVecs(Na->A,l,r);
86: return(0);
87: }
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;
98: MatAXPY(M,a,N,str);
99: return(0);
100: }
102: PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N,Mat *M)
103: {
104: Mat_HT *Na = (Mat_HT*)N->data;
107: *M = Na->A;
108: return(0);
109: }
111: /*@
112: MatHermitianTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT
114: Logically collective on Mat
116: Input Parameter:
117: . A - the MATTRANSPOSE matrix
119: Output Parameter:
120: . M - the matrix object stored inside A
122: Level: intermediate
124: .seealso: MatCreateHermitianTranspose()
126: @*/
127: PetscErrorCode MatHermitianTransposeGetMat(Mat A,Mat *M)
128: {
135: PetscUseMethod(A,"MatHermitianTransposeGetMat_C",(Mat,Mat*),(A,M));
136: return(0);
137: }
139: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat);
141: PetscErrorCode MatGetDiagonal_HT(Mat A,Vec v)
142: {
143: Mat_HT *Na = (Mat_HT*)A->data;
147: MatGetDiagonal(Na->A,v);
148: VecConjugate(v);
149: return(0);
150: }
152: PetscErrorCode MatConvert_HT(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
153: {
154: Mat_HT *Na = (Mat_HT*)A->data;
155: Mat B;
159: MatHermitianTranspose(Na->A,MAT_INITIAL_MATRIX,&B);
160: if (reuse != MAT_INPLACE_MATRIX) {
161: MatConvert(B,newtype,reuse,newmat);
162: MatDestroy(&B);
163: } else {
164: MatConvert(B,newtype,MAT_INPLACE_MATRIX,&B);
165: MatHeaderReplace(A,&B);
166: }
167: return(0);
168: }
170: /*@
171: MatCreateHermitianTranspose - Creates a new matrix object that behaves like A'*
173: Collective on Mat
175: Input Parameter:
176: . A - the (possibly rectangular) matrix
178: Output Parameter:
179: . N - the matrix that represents A'*
181: Level: intermediate
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
188: .seealso: MatCreateNormal(), MatMult(), MatMultHermitianTranspose(), MatCreate()
190: @*/
191: PetscErrorCode MatCreateHermitianTranspose(Mat A,Mat *N)
192: {
194: PetscInt m,n;
195: Mat_HT *Na;
196: VecType vtype;
199: MatGetLocalSize(A,&m,&n);
200: MatCreate(PetscObjectComm((PetscObject)A),N);
201: MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);
202: PetscLayoutSetUp((*N)->rmap);
203: PetscLayoutSetUp((*N)->cmap);
204: PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);
206: PetscNewLog(*N,&Na);
207: (*N)->data = (void*) Na;
208: PetscObjectReference((PetscObject)A);
209: Na->A = A;
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;
226: PetscObjectComposeFunction((PetscObject)(*N),"MatHermitianTransposeGetMat_C",MatHermitianTransposeGetMat_HT);
227: #if !defined(PETSC_USE_COMPLEX)
228: PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatHermitianTransposeGetMat_HT);
229: PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_Transpose);
230: #endif
231: MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
232: MatGetVecType(A,&vtype);
233: MatSetVecType(*N,vtype);
234: #if defined(PETSC_HAVE_DEVICE)
235: MatBindToCPU(*N,A->boundtocpu);
236: #endif
237: MatSetUp(*N);
238: return(0);
239: }