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;

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

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

 20:   MatMultTransposeAdd(Na->A,v1,v2,v3);
 21:   return 0;
 22: }

 24: PetscErrorCode MatMultTranspose_Transpose(Mat N,Vec x,Vec y)
 25: {
 26:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

 28:   MatMult(Na->A,x,y);
 29:   return 0;
 30: }

 32: PetscErrorCode MatMultTransposeAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3)
 33: {
 34:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

 36:   MatMultAdd(Na->A,v1,v2,v3);
 37:   return 0;
 38: }

 40: PetscErrorCode MatDestroy_Transpose(Mat N)
 41: {
 42:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

 44:   MatDestroy(&Na->A);
 45:   PetscObjectComposeFunction((PetscObject)N,"MatTransposeGetMat_C",NULL);
 46:   PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_anytype_C",NULL);
 47:   PetscFree(N->data);
 48:   return 0;
 49: }

 51: PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat* m)
 52: {
 53:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

 55:   if (op == MAT_COPY_VALUES) {
 56:     MatTranspose(Na->A,MAT_INITIAL_MATRIX,m);
 57:   } else if (op == MAT_DO_NOT_COPY_VALUES) {
 58:     MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);
 59:     MatTranspose(*m,MAT_INPLACE_MATRIX,m);
 60:   } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
 61:   return 0;
 62: }

 64: PetscErrorCode MatCreateVecs_Transpose(Mat A,Vec *r, Vec *l)
 65: {
 66:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;

 68:   MatCreateVecs(Aa->A,l,r);
 69:   return 0;
 70: }

 72: PetscErrorCode MatAXPY_Transpose(Mat Y,PetscScalar a,Mat X,MatStructure str)
 73: {
 74:   Mat_Transpose  *Ya = (Mat_Transpose*)Y->data;
 75:   Mat_Transpose  *Xa = (Mat_Transpose*)X->data;
 76:   Mat              M = Ya->A;
 77:   Mat              N = Xa->A;

 79:   MatAXPY(M,a,N,str);
 80:   return 0;
 81: }

 83: PetscErrorCode MatHasOperation_Transpose(Mat mat,MatOperation op,PetscBool *has)
 84: {
 85:   Mat_Transpose  *X = (Mat_Transpose*)mat->data;

 87:   *has = PETSC_FALSE;
 88:   if (op == MATOP_MULT) {
 89:     MatHasOperation(X->A,MATOP_MULT_TRANSPOSE,has);
 90:   } else if (op == MATOP_MULT_TRANSPOSE) {
 91:     MatHasOperation(X->A,MATOP_MULT,has);
 92:   } else if (op == MATOP_MULT_ADD) {
 93:     MatHasOperation(X->A,MATOP_MULT_TRANSPOSE_ADD,has);
 94:   } else if (op == MATOP_MULT_TRANSPOSE_ADD) {
 95:     MatHasOperation(X->A,MATOP_MULT_ADD,has);
 96:   } else if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
 97:   return 0;
 98: }

100: /* used by hermitian transpose */
101: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
102: {
103:   Mat            A,B,C,Ain,Bin,Cin;
104:   PetscBool      Aistrans,Bistrans,Cistrans;
105:   PetscInt       Atrans,Btrans,Ctrans;
106:   MatProductType ptype;

108:   MatCheckProduct(D,1);
109:   A = D->product->A;
110:   B = D->product->B;
111:   C = D->product->C;
112:   PetscObjectTypeCompare((PetscObject)A,MATTRANSPOSEMAT,&Aistrans);
113:   PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&Bistrans);
114:   PetscObjectTypeCompare((PetscObject)C,MATTRANSPOSEMAT,&Cistrans);
116:   Atrans = 0;
117:   Ain    = A;
118:   while (Aistrans) {
119:     Atrans++;
120:     MatTransposeGetMat(Ain,&Ain);
121:     PetscObjectTypeCompare((PetscObject)Ain,MATTRANSPOSEMAT,&Aistrans);
122:   }
123:   Btrans = 0;
124:   Bin    = B;
125:   while (Bistrans) {
126:     Btrans++;
127:     MatTransposeGetMat(Bin,&Bin);
128:     PetscObjectTypeCompare((PetscObject)Bin,MATTRANSPOSEMAT,&Bistrans);
129:   }
130:   Ctrans = 0;
131:   Cin    = C;
132:   while (Cistrans) {
133:     Ctrans++;
134:     MatTransposeGetMat(Cin,&Cin);
135:     PetscObjectTypeCompare((PetscObject)Cin,MATTRANSPOSEMAT,&Cistrans);
136:   }
137:   Atrans = Atrans%2;
138:   Btrans = Btrans%2;
139:   Ctrans = Ctrans%2;
140:   ptype = D->product->type; /* same product type by default */
141:   if (Ain->symmetric) Atrans = 0;
142:   if (Bin->symmetric) Btrans = 0;
143:   if (Cin && Cin->symmetric) Ctrans = 0;

145:   if (Atrans || Btrans || Ctrans) {
146:     ptype = MATPRODUCT_UNSPECIFIED;
147:     switch (D->product->type) {
148:     case MATPRODUCT_AB:
149:       if (Atrans && Btrans) { /* At * Bt we do not have support for this */
150:         /* TODO custom implementation ? */
151:       } else if (Atrans) { /* At * B */
152:         ptype = MATPRODUCT_AtB;
153:       } else { /* A * Bt */
154:         ptype = MATPRODUCT_ABt;
155:       }
156:       break;
157:     case MATPRODUCT_AtB:
158:       if (Atrans && Btrans) { /* A * Bt */
159:         ptype = MATPRODUCT_ABt;
160:       } else if (Atrans) { /* A * B */
161:         ptype = MATPRODUCT_AB;
162:       } else { /* At * Bt we do not have support for this */
163:         /* TODO custom implementation ? */
164:       }
165:       break;
166:     case MATPRODUCT_ABt:
167:       if (Atrans && Btrans) { /* At * B */
168:         ptype = MATPRODUCT_AtB;
169:       } else if (Atrans) { /* At * Bt we do not have support for this */
170:         /* TODO custom implementation ? */
171:       } else {  /* A * B */
172:         ptype = MATPRODUCT_AB;
173:       }
174:       break;
175:     case MATPRODUCT_PtAP:
176:       if (Atrans) { /* PtAtP */
177:         /* TODO custom implementation ? */
178:       } else { /* RARt */
179:         ptype = MATPRODUCT_RARt;
180:       }
181:       break;
182:     case MATPRODUCT_RARt:
183:       if (Atrans) { /* RAtRt */
184:         /* TODO custom implementation ? */
185:       } else { /* PtAP */
186:         ptype = MATPRODUCT_PtAP;
187:       }
188:       break;
189:     case MATPRODUCT_ABC:
190:       /* TODO custom implementation ? */
191:       break;
192:     default: SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[D->product->type]);
193:     }
194:   }
195:   MatProductReplaceMats(Ain,Bin,Cin,D);
196:   MatProductSetType(D,ptype);
197:   MatProductSetFromOptions(D);
198:   return 0;
199: }

201: PetscErrorCode MatGetDiagonal_Transpose(Mat A,Vec v)
202: {
203:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;

205:   MatGetDiagonal(Aa->A,v);
206:   return 0;
207: }

209: PetscErrorCode MatConvert_Transpose(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
210: {
211:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;
212:   PetscBool      flg;

214:   MatHasOperation(Aa->A,MATOP_TRANSPOSE,&flg);
215:   if (flg) {
216:     Mat B;

218:     MatTranspose(Aa->A,MAT_INITIAL_MATRIX,&B);
219:     if (reuse != MAT_INPLACE_MATRIX) {
220:       MatConvert(B,newtype,reuse,newmat);
221:       MatDestroy(&B);
222:     } else {
223:       MatConvert(B,newtype,MAT_INPLACE_MATRIX,&B);
224:       MatHeaderReplace(A,&B);
225:     }
226:   } else { /* use basic converter as fallback */
227:     MatConvert_Basic(A,newtype,reuse,newmat);
228:   }
229:   return 0;
230: }

232: PetscErrorCode MatTransposeGetMat_Transpose(Mat A,Mat *M)
233: {
234:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;

236:   *M = Aa->A;
237:   return 0;
238: }

240: /*@
241:       MatTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT

243:    Logically collective on Mat

245:    Input Parameter:
246: .   A  - the MATTRANSPOSE matrix

248:    Output Parameter:
249: .   M - the matrix object stored inside A

251:    Level: intermediate

253: .seealso: MatCreateTranspose()

255: @*/
256: PetscErrorCode MatTransposeGetMat(Mat A,Mat *M)
257: {
261:   PetscUseMethod(A,"MatTransposeGetMat_C",(Mat,Mat*),(A,M));
262:   return 0;
263: }

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

268:    Collective on Mat

270:    Input Parameter:
271: .   A  - the (possibly rectangular) matrix

273:    Output Parameter:
274: .   N - the matrix that represents A'

276:    Level: intermediate

278:    Notes:
279:     The transpose A' is NOT actually formed! Rather the new matrix
280:           object performs the matrix-vector product by using the MatMultTranspose() on
281:           the original matrix

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

285: @*/
286: PetscErrorCode  MatCreateTranspose(Mat A,Mat *N)
287: {
288:   PetscInt       m,n;
289:   Mat_Transpose  *Na;
290:   VecType        vtype;

292:   MatGetLocalSize(A,&m,&n);
293:   MatCreate(PetscObjectComm((PetscObject)A),N);
294:   MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);
295:   PetscLayoutSetUp((*N)->rmap);
296:   PetscLayoutSetUp((*N)->cmap);
297:   PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);

299:   PetscNewLog(*N,&Na);
300:   (*N)->data = (void*) Na;
301:   PetscObjectReference((PetscObject)A);
302:   Na->A      = A;

304:   (*N)->ops->destroy               = MatDestroy_Transpose;
305:   (*N)->ops->mult                  = MatMult_Transpose;
306:   (*N)->ops->multadd               = MatMultAdd_Transpose;
307:   (*N)->ops->multtranspose         = MatMultTranspose_Transpose;
308:   (*N)->ops->multtransposeadd      = MatMultTransposeAdd_Transpose;
309:   (*N)->ops->duplicate             = MatDuplicate_Transpose;
310:   (*N)->ops->getvecs               = MatCreateVecs_Transpose;
311:   (*N)->ops->axpy                  = MatAXPY_Transpose;
312:   (*N)->ops->hasoperation          = MatHasOperation_Transpose;
313:   (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
314:   (*N)->ops->getdiagonal           = MatGetDiagonal_Transpose;
315:   (*N)->ops->convert               = MatConvert_Transpose;
316:   (*N)->assembled                  = PETSC_TRUE;

318:   PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatTransposeGetMat_Transpose);
319:   PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_Transpose);
320:   MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
321:   MatGetVecType(A,&vtype);
322:   MatSetVecType(*N,vtype);
323: #if defined(PETSC_HAVE_DEVICE)
324:   MatBindToCPU(*N,A->boundtocpu);
325: #endif
326:   MatSetUp(*N);
327:   return 0;
328: }