Actual source code: htransm.c

petsc-3.15.0 2021-04-05
Report Typos and Errors

  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: }