Actual source code: transm.c


  2: #include <petsc/private/matimpl.h>

  4: typedef struct {
  5:   Mat A;
  6: } Mat_Transpose;

  8: PetscErrorCode MatMult_Transpose(Mat N,Vec x,Vec y)
  9: {
 10:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

 14:   MatMultTranspose(Na->A,x,y);
 15:   return(0);
 16: }

 18: PetscErrorCode MatMultAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3)
 19: {
 20:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

 24:   MatMultTransposeAdd(Na->A,v1,v2,v3);
 25:   return(0);
 26: }

 28: PetscErrorCode MatMultTranspose_Transpose(Mat N,Vec x,Vec y)
 29: {
 30:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

 34:   MatMult(Na->A,x,y);
 35:   return(0);
 36: }

 38: PetscErrorCode MatMultTransposeAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3)
 39: {
 40:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

 44:   MatMultAdd(Na->A,v1,v2,v3);
 45:   return(0);
 46: }

 48: PetscErrorCode MatDestroy_Transpose(Mat N)
 49: {
 50:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

 54:   MatDestroy(&Na->A);
 55:   PetscObjectComposeFunction((PetscObject)N,"MatTransposeGetMat_C",NULL);
 56:   PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_anytype_C",NULL);
 57:   PetscFree(N->data);
 58:   return(0);
 59: }

 61: PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat* m)
 62: {
 63:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

 67:   if (op == MAT_COPY_VALUES) {
 68:     MatTranspose(Na->A,MAT_INITIAL_MATRIX,m);
 69:   } else if (op == MAT_DO_NOT_COPY_VALUES) {
 70:     MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);
 71:     MatTranspose(*m,MAT_INPLACE_MATRIX,m);
 72:   } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
 73:   return(0);
 74: }

 76: PetscErrorCode MatCreateVecs_Transpose(Mat A,Vec *r, Vec *l)
 77: {
 78:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;

 82:   MatCreateVecs(Aa->A,l,r);
 83:   return(0);
 84: }

 86: PetscErrorCode MatAXPY_Transpose(Mat Y,PetscScalar a,Mat X,MatStructure str)
 87: {
 88:   Mat_Transpose  *Ya = (Mat_Transpose*)Y->data;
 89:   Mat_Transpose  *Xa = (Mat_Transpose*)X->data;
 90:   Mat              M = Ya->A;
 91:   Mat              N = Xa->A;

 95:   MatAXPY(M,a,N,str);
 96:   return(0);
 97: }

 99: PetscErrorCode MatHasOperation_Transpose(Mat mat,MatOperation op,PetscBool *has)
100: {
101:   Mat_Transpose  *X = (Mat_Transpose*)mat->data;

105:   *has = PETSC_FALSE;
106:   if (op == MATOP_MULT) {
107:     MatHasOperation(X->A,MATOP_MULT_TRANSPOSE,has);
108:   } else if (op == MATOP_MULT_TRANSPOSE) {
109:     MatHasOperation(X->A,MATOP_MULT,has);
110:   } else if (op == MATOP_MULT_ADD) {
111:     MatHasOperation(X->A,MATOP_MULT_TRANSPOSE_ADD,has);
112:   } else if (op == MATOP_MULT_TRANSPOSE_ADD) {
113:     MatHasOperation(X->A,MATOP_MULT_ADD,has);
114:   } else if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
115:   return(0);
116: }

118: /* used by hermitian transpose */
119: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
120: {
121:   Mat            A,B,C,Ain,Bin,Cin;
122:   PetscBool      Aistrans,Bistrans,Cistrans;
123:   PetscInt       Atrans,Btrans,Ctrans;
124:   MatProductType ptype;

128:   MatCheckProduct(D,1);
129:   A = D->product->A;
130:   B = D->product->B;
131:   C = D->product->C;
132:   PetscObjectTypeCompare((PetscObject)A,MATTRANSPOSEMAT,&Aistrans);
133:   PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&Bistrans);
134:   PetscObjectTypeCompare((PetscObject)C,MATTRANSPOSEMAT,&Cistrans);
135:   if (!Aistrans && !Bistrans && !Cistrans) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"This should not happen");
136:   Atrans = 0;
137:   Ain    = A;
138:   while (Aistrans) {
139:     Atrans++;
140:     MatTransposeGetMat(Ain,&Ain);
141:     PetscObjectTypeCompare((PetscObject)Ain,MATTRANSPOSEMAT,&Aistrans);
142:   }
143:   Btrans = 0;
144:   Bin    = B;
145:   while (Bistrans) {
146:     Btrans++;
147:     MatTransposeGetMat(Bin,&Bin);
148:     PetscObjectTypeCompare((PetscObject)Bin,MATTRANSPOSEMAT,&Bistrans);
149:   }
150:   Ctrans = 0;
151:   Cin    = C;
152:   while (Cistrans) {
153:     Ctrans++;
154:     MatTransposeGetMat(Cin,&Cin);
155:     PetscObjectTypeCompare((PetscObject)Cin,MATTRANSPOSEMAT,&Cistrans);
156:   }
157:   Atrans = Atrans%2;
158:   Btrans = Btrans%2;
159:   Ctrans = Ctrans%2;
160:   ptype = D->product->type; /* same product type by default */
161:   if (Ain->symmetric) Atrans = 0;
162:   if (Bin->symmetric) Btrans = 0;
163:   if (Cin && Cin->symmetric) Ctrans = 0;

165:   if (Atrans || Btrans || Ctrans) {
166:     ptype = MATPRODUCT_UNSPECIFIED;
167:     switch (D->product->type) {
168:     case MATPRODUCT_AB:
169:       if (Atrans && Btrans) { /* At * Bt we do not have support for this */
170:         /* TODO custom implementation ? */
171:       } else if (Atrans) { /* At * B */
172:         ptype = MATPRODUCT_AtB;
173:       } else { /* A * Bt */
174:         ptype = MATPRODUCT_ABt;
175:       }
176:       break;
177:     case MATPRODUCT_AtB:
178:       if (Atrans && Btrans) { /* A * Bt */
179:         ptype = MATPRODUCT_ABt;
180:       } else if (Atrans) { /* A * B */
181:         ptype = MATPRODUCT_AB;
182:       } else { /* At * Bt we do not have support for this */
183:         /* TODO custom implementation ? */
184:       }
185:       break;
186:     case MATPRODUCT_ABt:
187:       if (Atrans && Btrans) { /* At * B */
188:         ptype = MATPRODUCT_AtB;
189:       } else if (Atrans) { /* At * Bt we do not have support for this */
190:         /* TODO custom implementation ? */
191:       } else {  /* A * B */
192:         ptype = MATPRODUCT_AB;
193:       }
194:       break;
195:     case MATPRODUCT_PtAP:
196:       if (Atrans) { /* PtAtP */
197:         /* TODO custom implementation ? */
198:       } else { /* RARt */
199:         ptype = MATPRODUCT_RARt;
200:       }
201:       break;
202:     case MATPRODUCT_RARt:
203:       if (Atrans) { /* RAtRt */
204:         /* TODO custom implementation ? */
205:       } else { /* PtAP */
206:         ptype = MATPRODUCT_PtAP;
207:       }
208:       break;
209:     case MATPRODUCT_ABC:
210:       /* TODO custom implementation ? */
211:       break;
212:     default: SETERRQ1(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[D->product->type]);
213:     }
214:   }
215:   MatProductReplaceMats(Ain,Bin,Cin,D);
216:   MatProductSetType(D,ptype);
217:   MatProductSetFromOptions(D);
218:   return(0);
219: }

221: PetscErrorCode MatGetDiagonal_Transpose(Mat A,Vec v)
222: {
223:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;

227:   MatGetDiagonal(Aa->A,v);
228:   return(0);
229: }

231: PetscErrorCode MatConvert_Transpose(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
232: {
233:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;
234:   Mat            B;

238:   MatTranspose(Aa->A,MAT_INITIAL_MATRIX,&B);
239:   if (reuse != MAT_INPLACE_MATRIX) {
240:     MatConvert(B,newtype,reuse,newmat);
241:     MatDestroy(&B);
242:   } else {
243:     MatConvert(B,newtype,MAT_INPLACE_MATRIX,&B);
244:     MatHeaderReplace(A,&B);
245:   }
246:   return(0);
247: }

249: PetscErrorCode MatTransposeGetMat_Transpose(Mat A,Mat *M)
250: {
251:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;

254:   *M = Aa->A;
255:   return(0);
256: }

258: /*@
259:       MatTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT

261:    Logically collective on Mat

263:    Input Parameter:
264: .   A  - the MATTRANSPOSE matrix

266:    Output Parameter:
267: .   M - the matrix object stored inside A

269:    Level: intermediate

271: .seealso: MatCreateTranspose()

273: @*/
274: PetscErrorCode MatTransposeGetMat(Mat A,Mat *M)
275: {

282:   PetscUseMethod(A,"MatTransposeGetMat_C",(Mat,Mat*),(A,M));
283:   return(0);
284: }

286: /*@
287:       MatCreateTranspose - Creates a new matrix object that behaves like A'

289:    Collective on Mat

291:    Input Parameter:
292: .   A  - the (possibly rectangular) matrix

294:    Output Parameter:
295: .   N - the matrix that represents A'

297:    Level: intermediate

299:    Notes:
300:     The transpose A' is NOT actually formed! Rather the new matrix
301:           object performs the matrix-vector product by using the MatMultTranspose() on
302:           the original matrix

304: .seealso: MatCreateNormal(), MatMult(), MatMultTranspose(), MatCreate()

306: @*/
307: PetscErrorCode  MatCreateTranspose(Mat A,Mat *N)
308: {
310:   PetscInt       m,n;
311:   Mat_Transpose  *Na;
312:   VecType        vtype;

315:   MatGetLocalSize(A,&m,&n);
316:   MatCreate(PetscObjectComm((PetscObject)A),N);
317:   MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);
318:   PetscLayoutSetUp((*N)->rmap);
319:   PetscLayoutSetUp((*N)->cmap);
320:   PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);

322:   PetscNewLog(*N,&Na);
323:   (*N)->data = (void*) Na;
324:   PetscObjectReference((PetscObject)A);
325:   Na->A      = A;

327:   (*N)->ops->destroy               = MatDestroy_Transpose;
328:   (*N)->ops->mult                  = MatMult_Transpose;
329:   (*N)->ops->multadd               = MatMultAdd_Transpose;
330:   (*N)->ops->multtranspose         = MatMultTranspose_Transpose;
331:   (*N)->ops->multtransposeadd      = MatMultTransposeAdd_Transpose;
332:   (*N)->ops->duplicate             = MatDuplicate_Transpose;
333:   (*N)->ops->getvecs               = MatCreateVecs_Transpose;
334:   (*N)->ops->axpy                  = MatAXPY_Transpose;
335:   (*N)->ops->hasoperation          = MatHasOperation_Transpose;
336:   (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
337:   (*N)->ops->getdiagonal           = MatGetDiagonal_Transpose;
338:   (*N)->ops->convert               = MatConvert_Transpose;
339:   (*N)->assembled                  = PETSC_TRUE;

341:   PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatTransposeGetMat_Transpose);
342:   PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_Transpose);
343:   MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
344:   MatGetVecType(A,&vtype);
345:   MatSetVecType(*N,vtype);
346: #if defined(PETSC_HAVE_DEVICE)
347:   MatBindToCPU(*N,A->boundtocpu);
348: #endif
349:   MatSetUp(*N);
350:   return(0);
351: }