Actual source code: adamat.c

petsc-3.12.5 2020-03-29
Report Typos and Errors
  1:  #include <petscmat.h>

  3: PETSC_INTERN PetscErrorCode MatCreateADA(Mat,Vec, Vec, Mat*);

  5: typedef struct {
  6:   Mat      A;
  7:   Vec      D1;
  8:   Vec      D2;
  9:   Vec      W;
 10:   Vec      W2;
 11:   Vec      ADADiag;
 12:   PetscInt GotDiag;
 13: } _p_TaoMatADACtx;
 14: typedef  _p_TaoMatADACtx* TaoMatADACtx;

 16: static PetscErrorCode MatMult_ADA(Mat mat,Vec a,Vec y)
 17: {
 18:   TaoMatADACtx   ctx;
 19:   PetscReal      one = 1.0;

 23:   MatShellGetContext(mat,(void **)&ctx);
 24:   MatMult(ctx->A,a,ctx->W);
 25:   if (ctx->D1){
 26:     VecPointwiseMult(ctx->W,ctx->D1,ctx->W);
 27:   }
 28:   MatMultTranspose(ctx->A,ctx->W,y);
 29:   if (ctx->D2){
 30:     VecPointwiseMult(ctx->W2, ctx->D2, a);
 31:     VecAXPY(y, one, ctx->W2);
 32:   }
 33:   return(0);
 34: }

 36: static PetscErrorCode MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
 37: {

 41:   MatMult_ADA(mat,a,y);
 42:   return(0);
 43: }

 45: static PetscErrorCode MatDiagonalSet_ADA(Mat M,Vec D, InsertMode mode)
 46: {
 47:   TaoMatADACtx   ctx;
 48:   PetscReal      zero=0.0,one = 1.0;

 52:   if (mode == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)M),PETSC_ERR_SUP,"Cannot insert diagonal entries of this matrix type, can only add");
 53:   MatShellGetContext(M,(void **)&ctx);
 54:   if (!ctx->D2){
 55:     VecDuplicate(D,&ctx->D2);
 56:     VecSet(ctx->D2, zero);
 57:   }
 58:   VecAXPY(ctx->D2, one, D);
 59:   return(0);
 60: }

 62: static PetscErrorCode MatDestroy_ADA(Mat mat)
 63: {
 65:   TaoMatADACtx   ctx;

 68:   ierr=MatShellGetContext(mat,(void **)&ctx);
 69:   ierr=VecDestroy(&ctx->W);
 70:   ierr=VecDestroy(&ctx->W2);
 71:   ierr=VecDestroy(&ctx->ADADiag);
 72:   ierr=MatDestroy(&ctx->A);
 73:   ierr=VecDestroy(&ctx->D1);
 74:   ierr=VecDestroy(&ctx->D2);
 75:   PetscFree(ctx);
 76:   return(0);
 77: }

 79: static PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer)
 80: {
 82:   return(0);
 83: }

 85: static PetscErrorCode MatShift_ADA(Mat Y, PetscReal a)
 86: {
 88:   TaoMatADACtx   ctx;

 91:   MatShellGetContext(Y,(void **)&ctx);
 92:   VecShift(ctx->D2,a);
 93:   return(0);
 94: }

 96: static PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
 97: {
 98:   PetscErrorCode    ierr;
 99:   TaoMatADACtx      ctx;
100:   Mat               A2;
101:   Vec               D1b=NULL,D2b;

104:   MatShellGetContext(mat,(void **)&ctx);
105:   MatDuplicate(ctx->A,op,&A2);
106:   if (ctx->D1){
107:     VecDuplicate(ctx->D1,&D1b);
108:     VecCopy(ctx->D1,D1b);
109:   }
110:   VecDuplicate(ctx->D2,&D2b);
111:   VecCopy(ctx->D2,D2b);
112:   MatCreateADA(A2,D1b,D2b,M);
113:   if (ctx->D1){
114:     PetscObjectDereference((PetscObject)D1b);
115:   }
116:   PetscObjectDereference((PetscObject)D2b);
117:   PetscObjectDereference((PetscObject)A2);
118:   return(0);
119: }

121: static PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg)
122: {
124:   TaoMatADACtx   ctx1,ctx2;

127:   MatShellGetContext(A,(void **)&ctx1);
128:   MatShellGetContext(B,(void **)&ctx2);
129:   VecEqual(ctx1->D2,ctx2->D2,flg);
130:   if (*flg==PETSC_TRUE){
131:     VecEqual(ctx1->D1,ctx2->D1,flg);
132:   }
133:   if (*flg==PETSC_TRUE){
134:     MatEqual(ctx1->A,ctx2->A,flg);
135:   }
136:   return(0);
137: }

139: static PetscErrorCode MatScale_ADA(Mat mat, PetscReal a)
140: {
142:   TaoMatADACtx   ctx;

145:   MatShellGetContext(mat,(void **)&ctx);
146:   VecScale(ctx->D1,a);
147:   if (ctx->D2){
148:     VecScale(ctx->D2,a);
149:   }
150:   return(0);
151: }

153: static PetscErrorCode MatTranspose_ADA(Mat mat,MatReuse reuse,Mat *B)
154: {
156:   TaoMatADACtx   ctx;

159:   MatShellGetContext(mat,(void **)&ctx);
160:   if (reuse == MAT_INITIAL_MATRIX){
161:     MatDuplicate(mat,MAT_COPY_VALUES,B);
162:   } else if (reuse == MAT_REUSE_MATRIX){
163:     MatCopy(mat,*B,SAME_NONZERO_PATTERN);
164:   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Does not support inplace transpose");
165:   return(0);
166: }

168: static PetscErrorCode MatADAComputeDiagonal(Mat mat)
169: {
171:   PetscInt       i,m,n,low,high;
172:   PetscScalar    *dtemp,*dptr;
173:   TaoMatADACtx   ctx;

176:   MatShellGetContext(mat,(void **)&ctx);
177:   MatGetOwnershipRange(mat, &low, &high);
178:   MatGetSize(mat,&m,&n);

180:   PetscMalloc1(n,&dtemp );
181:   for (i=0; i<n; i++){
182:     MatGetColumnVector(ctx->A, ctx->W, i);
183:     VecPointwiseMult(ctx->W,ctx->W,ctx->W);
184:     VecDotBegin(ctx->D1, ctx->W,dtemp+i);
185:   }
186:   for (i=0; i<n; i++){
187:     VecDotEnd(ctx->D1, ctx->W,dtemp+i);
188:   }

190:   VecGetArray(ctx->ADADiag,&dptr);
191:   for (i=low; i<high; i++){
192:     dptr[i-low]= dtemp[i];
193:   }
194:   VecRestoreArray(ctx->ADADiag,&dptr);
195:   PetscFree(dtemp);
196:   return(0);
197: }

199: static PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v)
200: {
201:   PetscErrorCode  ierr;
202:   PetscReal       one=1.0;
203:   TaoMatADACtx    ctx;

206:   MatShellGetContext(mat,(void **)&ctx);
207:   MatADAComputeDiagonal(mat);
208:   VecCopy(ctx->ADADiag,v);
209:   if (ctx->D2){
210:     VecAXPY(v, one, ctx->D2);
211:   }
212:   return(0);
213: }

215: static PetscErrorCode MatCreateSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll, Mat *newmat)
216: {
217:   PetscErrorCode    ierr;
218:   PetscInt          low,high;
219:   IS                ISrow;
220:   Vec               D1,D2;
221:   Mat               Atemp;
222:   TaoMatADACtx      ctx;
223:   PetscBool         isequal;

226:   ISEqual(isrow,iscol,&isequal);
227:   if (!isequal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for idential column and row indices");
228:   MatShellGetContext(mat,(void **)&ctx);

230:   MatGetOwnershipRange(ctx->A,&low,&high);
231:   ISCreateStride(PetscObjectComm((PetscObject)mat),high-low,low,1,&ISrow);
232:   MatCreateSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);
233:   ISDestroy(&ISrow);

235:   if (ctx->D1){
236:     ierr=VecDuplicate(ctx->D1,&D1);
237:     ierr=VecCopy(ctx->D1,D1);
238:   } else {
239:     D1 = NULL;
240:   }

242:   if (ctx->D2){
243:     Vec D2sub;

245:     ierr=VecGetSubVector(ctx->D2,isrow,&D2sub);
246:     ierr=VecDuplicate(D2sub,&D2);
247:     ierr=VecCopy(D2sub,D2);
248:     ierr=VecRestoreSubVector(ctx->D2,isrow,&D2sub);
249:   } else {
250:     D2 = NULL;
251:   }

253:   MatCreateADA(Atemp,D1,D2,newmat);
254:   MatShellGetContext(*newmat,(void **)&ctx);
255:   PetscObjectDereference((PetscObject)Atemp);
256:   if (ctx->D1){
257:     PetscObjectDereference((PetscObject)D1);
258:   }
259:   if (ctx->D2){
260:     PetscObjectDereference((PetscObject)D2);
261:   }
262:   return(0);
263: }

265: static PetscErrorCode MatCreateSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
266: {
268:   PetscInt       i;

271:   if (scall == MAT_INITIAL_MATRIX) {
272:     PetscCalloc1(n+1,B);
273:   }
274:   for ( i=0; i<n; i++ ) {
275:     MatCreateSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);
276:   }
277:   return(0);
278: }

280: static PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col)
281: {
283:   PetscInt       low,high;
284:   PetscScalar    zero=0.0,one=1.0;

287:   ierr=VecSet(Y, zero);
288:   ierr=VecGetOwnershipRange(Y,&low,&high);
289:   if (col>=low && col<high){
290:     ierr=VecSetValue(Y,col,one,INSERT_VALUES);
291:   }
292:   ierr=VecAssemblyBegin(Y);
293:   ierr=VecAssemblyEnd(Y);
294:   ierr=MatMult_ADA(mat,Y,Y);
295:   return(0);
296: }

298: PETSC_INTERN PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
299: {
301:   PetscMPIInt    size;
302:   PetscBool      sametype, issame, isdense, isseqdense;
303:   TaoMatADACtx   ctx;

306:   MatShellGetContext(mat,(void **)&ctx);
307:   MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);

309:   PetscObjectTypeCompare((PetscObject)mat,newtype,&sametype);
310:   PetscObjectTypeCompare((PetscObject)mat,MATSAME,&issame);
311:   PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&isdense);
312:   PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);

314:   if (sametype || issame) {
315:     ierr=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);
316:   } else if (isdense) {
317:     PetscInt          i,j,low,high,m,n,M,N;
318:     const PetscScalar *dptr;
319:     Vec               X;

321:     VecDuplicate(ctx->D2,&X);
322:     ierr=MatGetSize(mat,&M,&N);
323:     ierr=MatGetLocalSize(mat,&m,&n);
324:     MatCreateDense(PetscObjectComm((PetscObject)mat),m,m,N,N,NULL,NewMat);
325:     MatGetOwnershipRange(*NewMat,&low,&high);
326:     for (i=0;i<M;i++){
327:       MatGetColumnVector_ADA(mat,X,i);
328:       VecGetArrayRead(X,&dptr);
329:       for (j=0; j<high-low; j++){
330:         MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
331:       }
332:       VecRestoreArrayRead(X,&dptr);
333:     }
334:     MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
335:     MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
336:     VecDestroy(&X);
337:   } else if (isseqdense && size==1){
338:     PetscInt          i,j,low,high,m,n,M,N;
339:     const PetscScalar *dptr;
340:     Vec               X;

342:     VecDuplicate(ctx->D2,&X);
343:     MatGetSize(mat,&M,&N);
344:     MatGetLocalSize(mat,&m,&n);
345:     MatCreateSeqDense(PetscObjectComm((PetscObject)mat),N,N,NULL,NewMat);
346:     MatGetOwnershipRange(*NewMat,&low,&high);
347:     for (i=0;i<M;i++){
348:       MatGetColumnVector_ADA(mat,X,i);
349:       VecGetArrayRead(X,&dptr);
350:       for (j=0; j<high-low; j++){
351:         MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
352:       }
353:       VecRestoreArrayRead(X,&dptr);
354:     }
355:     MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
356:     MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
357:     VecDestroy(&X);
358:   } else SETERRQ(PETSC_COMM_SELF,1,"No support to convert objects to that type");
359:   return(0);
360: }

362: static PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
363: {
365:   TaoMatADACtx   ctx;

368:   MatShellGetContext(mat,(void **)&ctx);
369:   if (type == NORM_FROBENIUS) {
370:     *norm = 1.0;
371:   } else if (type == NORM_1 || type == NORM_INFINITY) {
372:     *norm = 1.0;
373:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
374:   return(0);
375: }

377: /*@C
378:    MatCreateADA - Creates a matrix M=A^T D1 A + D2 where D1, D2 are diagonal

380:    Collective on matrix

382:    Input Parameters:
383: +  mat - matrix of arbitrary type
384: .  d1 - A vector defining a diagonal matrix
385: -  d2 - A vector defining a diagonal matrix 

387:    Output Parameters:
388: .  J - New matrix whose operations are defined in terms of mat, D1, and D2.

390:    Notes:
391:    The user provides the input data and is responsible for destroying
392:    this data after matrix J has been destroyed.

394:    Level: developer

396: .seealso: MatCreate()
397: @*/
398: PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J)
399: {
400:   MPI_Comm       comm = PetscObjectComm((PetscObject)mat);
401:   TaoMatADACtx   ctx;
403:   PetscInt       nloc,n;

406:   PetscNew(&ctx);
407:   ctx->A=mat;
408:   ctx->D1=d1;
409:   ctx->D2=d2;
410:   if (d1){
411:     VecDuplicate(d1,&ctx->W);
412:     PetscObjectReference((PetscObject)d1);
413:   } else {
414:     ctx->W = NULL;
415:   }
416:   if (d2){
417:     VecDuplicate(d2,&ctx->W2);
418:     VecDuplicate(d2,&ctx->ADADiag);
419:      PetscObjectReference((PetscObject)d2);
420:   } else {
421:     ctx->W2      = NULL;
422:     ctx->ADADiag = NULL;
423:   }

425:   ctx->GotDiag = 0;
426:    PetscObjectReference((PetscObject)mat);

428:   ierr=VecGetLocalSize(d2,&nloc);
429:   ierr=VecGetSize(d2,&n);

431:   MatCreateShell(comm,nloc,nloc,n,n,ctx,J);
432:   MatShellSetManageScalingShifts(*J);
433:   MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);
434:   MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);
435:   MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);
436:   MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);
437:   MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);
438:   MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);
439:   MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);
440:   MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);
441:   MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);
442:   MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);
443:   MatShellSetOperation(*J,MATOP_CREATE_SUBMATRICES,(void(*)(void))MatCreateSubMatrices_ADA);
444:   MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);
445:   MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);
446:   MatShellSetOperation(*J,MATOP_CREATE_SUBMATRIX,(void(*)(void))MatCreateSubMatrix_ADA);

448:   PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);
449:   PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));

451:   MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);
452:   return(0);
453: }