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;
12: MatMultHermitianTranspose(Na->A,x,y);
13: return 0;
14: }
16: PetscErrorCode MatMultAdd_HT(Mat N,Vec v1,Vec v2,Vec v3)
17: {
18: Mat_HT *Na = (Mat_HT*)N->data;
20: MatMultHermitianTransposeAdd(Na->A,v1,v2,v3);
21: return 0;
22: }
24: PetscErrorCode MatMultHermitianTranspose_HT(Mat N,Vec x,Vec y)
25: {
26: Mat_HT *Na = (Mat_HT*)N->data;
28: MatMult(Na->A,x,y);
29: return 0;
30: }
32: PetscErrorCode MatMultHermitianTransposeAdd_HT(Mat N,Vec v1,Vec v2,Vec v3)
33: {
34: Mat_HT *Na = (Mat_HT*)N->data;
36: MatMultAdd(Na->A,v1,v2,v3);
37: return 0;
38: }
40: PetscErrorCode MatDestroy_HT(Mat N)
41: {
42: Mat_HT *Na = (Mat_HT*)N->data;
44: MatDestroy(&Na->A);
45: PetscObjectComposeFunction((PetscObject)N,"MatHermitianTransposeGetMat_C",NULL);
46: #if !defined(PETSC_USE_COMPLEX)
47: PetscObjectComposeFunction((PetscObject)N,"MatTransposeGetMat_C",NULL);
48: PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_anytype_C",NULL);
49: #endif
50: PetscFree(N->data);
51: return 0;
52: }
54: PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat* m)
55: {
56: Mat_HT *Na = (Mat_HT*)N->data;
58: if (op == MAT_COPY_VALUES) {
59: MatHermitianTranspose(Na->A,MAT_INITIAL_MATRIX,m);
60: } else if (op == MAT_DO_NOT_COPY_VALUES) {
61: MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);
62: MatHermitianTranspose(*m,MAT_INPLACE_MATRIX,m);
63: } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
64: return 0;
65: }
67: PetscErrorCode MatCreateVecs_HT(Mat N,Vec *r, Vec *l)
68: {
69: Mat_HT *Na = (Mat_HT*)N->data;
71: MatCreateVecs(Na->A,l,r);
72: return 0;
73: }
75: PetscErrorCode MatAXPY_HT(Mat Y,PetscScalar a,Mat X,MatStructure str)
76: {
77: Mat_HT *Ya = (Mat_HT*)Y->data;
78: Mat_HT *Xa = (Mat_HT*)X->data;
79: Mat M = Ya->A;
80: Mat N = Xa->A;
82: MatAXPY(M,a,N,str);
83: return 0;
84: }
86: PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N,Mat *M)
87: {
88: Mat_HT *Na = (Mat_HT*)N->data;
90: *M = Na->A;
91: return 0;
92: }
94: /*@
95: MatHermitianTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT
97: Logically collective on Mat
99: Input Parameter:
100: . A - the MATTRANSPOSE matrix
102: Output Parameter:
103: . M - the matrix object stored inside A
105: Level: intermediate
107: .seealso: MatCreateHermitianTranspose()
109: @*/
110: PetscErrorCode MatHermitianTransposeGetMat(Mat A,Mat *M)
111: {
115: PetscUseMethod(A,"MatHermitianTransposeGetMat_C",(Mat,Mat*),(A,M));
116: return 0;
117: }
119: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat);
121: PetscErrorCode MatGetDiagonal_HT(Mat A,Vec v)
122: {
123: Mat_HT *Na = (Mat_HT*)A->data;
125: MatGetDiagonal(Na->A,v);
126: VecConjugate(v);
127: return 0;
128: }
130: PetscErrorCode MatConvert_HT(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
131: {
132: Mat_HT *Na = (Mat_HT*)A->data;
133: PetscBool flg;
135: MatHasOperation(Na->A,MATOP_HERMITIAN_TRANSPOSE,&flg);
136: if (flg) {
137: Mat B;
139: MatHermitianTranspose(Na->A,MAT_INITIAL_MATRIX,&B);
140: if (reuse != MAT_INPLACE_MATRIX) {
141: MatConvert(B,newtype,reuse,newmat);
142: MatDestroy(&B);
143: } else {
144: MatConvert(B,newtype,MAT_INPLACE_MATRIX,&B);
145: MatHeaderReplace(A,&B);
146: }
147: } else { /* use basic converter as fallback */
148: MatConvert_Basic(A,newtype,reuse,newmat);
149: }
150: return 0;
151: }
153: /*@
154: MatCreateHermitianTranspose - Creates a new matrix object that behaves like A'*
156: Collective on Mat
158: Input Parameter:
159: . A - the (possibly rectangular) matrix
161: Output Parameter:
162: . N - the matrix that represents A'*
164: Level: intermediate
166: Notes:
167: The hermitian transpose A' is NOT actually formed! Rather the new matrix
168: object performs the matrix-vector product by using the MatMultHermitianTranspose() on
169: the original matrix
171: .seealso: MatCreateNormal(), MatMult(), MatMultHermitianTranspose(), MatCreate()
172: @*/
173: PetscErrorCode MatCreateHermitianTranspose(Mat A,Mat *N)
174: {
175: PetscInt m,n;
176: Mat_HT *Na;
177: VecType vtype;
179: MatGetLocalSize(A,&m,&n);
180: MatCreate(PetscObjectComm((PetscObject)A),N);
181: MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);
182: PetscLayoutSetUp((*N)->rmap);
183: PetscLayoutSetUp((*N)->cmap);
184: PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);
186: PetscNewLog(*N,&Na);
187: (*N)->data = (void*) Na;
188: PetscObjectReference((PetscObject)A);
189: Na->A = A;
191: (*N)->ops->destroy = MatDestroy_HT;
192: (*N)->ops->mult = MatMult_HT;
193: (*N)->ops->multadd = MatMultAdd_HT;
194: (*N)->ops->multhermitiantranspose = MatMultHermitianTranspose_HT;
195: (*N)->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_HT;
196: #if !defined(PETSC_USE_COMPLEX)
197: (*N)->ops->multtranspose = MatMultHermitianTranspose_HT;
198: (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_HT;
199: #endif
200: (*N)->ops->duplicate = MatDuplicate_HT;
201: (*N)->ops->getvecs = MatCreateVecs_HT;
202: (*N)->ops->axpy = MatAXPY_HT;
203: #if !defined(PETSC_USE_COMPLEX)
204: (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
205: #endif
206: (*N)->ops->getdiagonal = MatGetDiagonal_HT;
207: (*N)->ops->convert = MatConvert_HT;
208: (*N)->assembled = PETSC_TRUE;
210: PetscObjectComposeFunction((PetscObject)(*N),"MatHermitianTransposeGetMat_C",MatHermitianTransposeGetMat_HT);
211: #if !defined(PETSC_USE_COMPLEX)
212: PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatHermitianTransposeGetMat_HT);
213: PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_Transpose);
214: #endif
215: MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
216: MatGetVecType(A,&vtype);
217: MatSetVecType(*N,vtype);
218: #if defined(PETSC_HAVE_DEVICE)
219: MatBindToCPU(*N,A->boundtocpu);
220: #endif
221: MatSetUp(*N);
222: return 0;
223: }