Actual source code: adamat.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/matimpl.h> /*I "mat.h" I*/
2: #include <petsc/private/vecimpl.h>
4: typedef struct{
5: Mat A;
6: Vec D1;
7: Vec D2;
8: Vec W;
9: Vec W2;
10: Vec ADADiag;
11: PetscInt GotDiag;
12: } _p_TaoMatADACtx;
13: typedef _p_TaoMatADACtx* TaoMatADACtx;
15: extern PetscErrorCode MatCreateADA(Mat,Vec, Vec, Mat*);
19: static PetscErrorCode MatMult_ADA(Mat mat,Vec a,Vec y)
20: {
21: TaoMatADACtx ctx;
22: PetscReal one = 1.0;
26: MatShellGetContext(mat,(void **)&ctx);
27: MatMult(ctx->A,a,ctx->W);
28: if (ctx->D1){
29: VecPointwiseMult(ctx->W,ctx->D1,ctx->W);
30: }
31: MatMultTranspose(ctx->A,ctx->W,y);
32: if (ctx->D2){
33: VecPointwiseMult(ctx->W2, ctx->D2, a);
34: VecAXPY(y, one, ctx->W2);
35: }
36: return(0);
37: }
41: static PetscErrorCode MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
42: {
46: MatMult_ADA(mat,a,y);
47: return(0);
48: }
52: static PetscErrorCode MatDiagonalSet_ADA(Vec D, Mat M)
53: {
54: TaoMatADACtx ctx;
55: PetscReal zero=0.0,one = 1.0;
59: MatShellGetContext(M,(void **)&ctx);
60: if (!ctx->D2){
61: VecDuplicate(D,&ctx->D2);
62: VecSet(ctx->D2, zero);
63: }
64: VecAXPY(ctx->D2, one, D);
65: return(0);
66: }
70: static PetscErrorCode MatDestroy_ADA(Mat mat)
71: {
73: TaoMatADACtx ctx;
76: ierr=MatShellGetContext(mat,(void **)&ctx);
77: ierr=VecDestroy(&ctx->W);
78: ierr=VecDestroy(&ctx->W2);
79: ierr=VecDestroy(&ctx->ADADiag);
80: ierr=MatDestroy(&ctx->A);
81: ierr=VecDestroy(&ctx->D1);
82: ierr=VecDestroy(&ctx->D2);
83: PetscFree(ctx);
84: return(0);
85: }
89: static PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer)
90: {
92: return(0);
93: }
97: static PetscErrorCode MatShift_ADA(Mat Y, PetscReal a)
98: {
100: TaoMatADACtx ctx;
103: MatShellGetContext(Y,(void **)&ctx);
104: VecShift(ctx->D2,a);
105: return(0);
106: }
110: static PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
111: {
112: PetscErrorCode ierr;
113: TaoMatADACtx ctx;
114: Mat A2;
115: Vec D1b=NULL,D2b;
118: MatShellGetContext(mat,(void **)&ctx);
119: MatDuplicate(ctx->A,op,&A2);
120: if (ctx->D1){
121: VecDuplicate(ctx->D1,&D1b);
122: VecCopy(ctx->D1,D1b);
123: }
124: VecDuplicate(ctx->D2,&D2b);
125: VecCopy(ctx->D2,D2b);
126: MatCreateADA(A2,D1b,D2b,M);
127: if (ctx->D1){
128: PetscObjectDereference((PetscObject)D1b);
129: }
130: PetscObjectDereference((PetscObject)D2b);
131: PetscObjectDereference((PetscObject)A2);
132: return(0);
133: }
137: static PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg)
138: {
140: TaoMatADACtx ctx1,ctx2;
143: MatShellGetContext(A,(void **)&ctx1);
144: MatShellGetContext(B,(void **)&ctx2);
145: VecEqual(ctx1->D2,ctx2->D2,flg);
146: if (*flg==PETSC_TRUE){
147: VecEqual(ctx1->D1,ctx2->D1,flg);
148: }
149: if (*flg==PETSC_TRUE){
150: MatEqual(ctx1->A,ctx2->A,flg);
151: }
152: return(0);
153: }
157: static PetscErrorCode MatScale_ADA(Mat mat, PetscReal a)
158: {
160: TaoMatADACtx ctx;
163: MatShellGetContext(mat,(void **)&ctx);
164: VecScale(ctx->D1,a);
165: if (ctx->D2){
166: VecScale(ctx->D2,a);
167: }
168: return(0);
169: }
173: static PetscErrorCode MatTranspose_ADA(Mat mat,Mat *B)
174: {
176: TaoMatADACtx ctx;
179: if (*B){
180: MatShellGetContext(mat,(void **)&ctx);
181: MatDuplicate(mat,MAT_COPY_VALUES,B);
182: }
183: return(0);
184: }
188: PetscErrorCode MatADAComputeDiagonal(Mat mat)
189: {
191: PetscInt i,m,n,low,high;
192: PetscScalar *dtemp,*dptr;
193: TaoMatADACtx ctx;
196: MatShellGetContext(mat,(void **)&ctx);
197: MatGetOwnershipRange(mat, &low, &high);
198: MatGetSize(mat,&m,&n);
200: PetscMalloc1(n,&dtemp );
202: for (i=0; i<n; i++){
203: MatGetColumnVector(ctx->A, ctx->W, i);
204: VecPointwiseMult(ctx->W,ctx->W,ctx->W);
205: VecDotBegin(ctx->D1, ctx->W,dtemp+i);
206: }
207: for (i=0; i<n; i++){
208: VecDotEnd(ctx->D1, ctx->W,dtemp+i);
209: }
211: VecGetArray(ctx->ADADiag,&dptr);
212: for (i=low; i<high; i++){
213: dptr[i-low]= dtemp[i];
214: }
215: VecRestoreArray(ctx->ADADiag,&dptr);
216: PetscFree(dtemp);
217: return(0);
218: }
222: static PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v)
223: {
224: PetscErrorCode ierr;
225: PetscReal one=1.0;
226: TaoMatADACtx ctx;
229: MatShellGetContext(mat,(void **)&ctx);
230: MatADAComputeDiagonal(mat);
231: ierr=VecCopy(ctx->ADADiag,v);
232: if (ctx->D2){
233: ierr=VecAXPY(v, one, ctx->D2);
234: }
235: return(0);
236: }
240: static PetscErrorCode MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll, Mat *newmat)
241: {
242: PetscErrorCode ierr;
243: PetscInt low,high;
244: IS ISrow;
245: Vec D1,D2;
246: Mat Atemp;
247: TaoMatADACtx ctx;
248: PetscBool isequal;
251: ISEqual(isrow,iscol,&isequal);
252: if (!isequal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for idential column and row indices");
253: MatShellGetContext(mat,(void **)&ctx);
255: MatGetOwnershipRange(ctx->A,&low,&high);
256: ISCreateStride(((PetscObject)mat)->comm,high-low,low,1,&ISrow);
257: MatGetSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);
258: ISDestroy(&ISrow);
260: if (ctx->D1){
261: ierr=VecDuplicate(ctx->D1,&D1);
262: ierr=VecCopy(ctx->D1,D1);
263: } else {
264: D1 = NULL;
265: }
267: if (ctx->D2){
268: Vec D2sub;
270: ierr=VecGetSubVector(ctx->D2,isrow,&D2sub);
271: ierr=VecDuplicate(D2sub,&D2);
272: ierr=VecCopy(D2sub,D2);
273: ierr=VecRestoreSubVector(ctx->D2,isrow,&D2sub);
274: } else {
275: D2 = NULL;
276: }
278: MatCreateADA(Atemp,D1,D2,newmat);
279: MatShellGetContext(*newmat,(void **)&ctx);
280: PetscObjectDereference((PetscObject)Atemp);
281: if (ctx->D1){
282: PetscObjectDereference((PetscObject)D1);
283: }
284: if (ctx->D2){
285: PetscObjectDereference((PetscObject)D2);
286: }
287: return(0);
288: }
292: static PetscErrorCode MatGetSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
293: {
295: PetscInt i;
298: if (scall == MAT_INITIAL_MATRIX) {
299: PetscMalloc1(n+1,B );
300: }
301: for ( i=0; i<n; i++ ) {
302: MatGetSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);
303: }
304: return(0);
305: }
309: static PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col)
310: {
312: PetscInt low,high;
313: PetscScalar zero=0.0,one=1.0;
316: ierr=VecSet(Y, zero);
317: ierr=VecGetOwnershipRange(Y,&low,&high);
318: if (col>=low && col<high){
319: ierr=VecSetValue(Y,col,one,INSERT_VALUES);
320: }
321: ierr=VecAssemblyBegin(Y);
322: ierr=VecAssemblyEnd(Y);
323: ierr=MatMult_ADA(mat,Y,Y);
324: return(0);
325: }
329: PETSC_INTERN PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
330: {
332: PetscMPIInt size;
333: PetscBool sametype, issame, isdense, isseqdense;
334: TaoMatADACtx ctx;
337: MatShellGetContext(mat,(void **)&ctx);
338: MPI_Comm_size(((PetscObject)mat)->comm,&size);
340: PetscObjectTypeCompare((PetscObject)mat,newtype,&sametype);
341: PetscObjectTypeCompare((PetscObject)mat,MATSAME,&issame);
342: PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&isdense);
343: PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);
345: if (sametype || issame) {
346: ierr=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);
347: } else if (isdense) {
348: PetscInt i,j,low,high,m,n,M,N;
349: const PetscScalar *dptr;
350: Vec X;
352: VecDuplicate(ctx->D2,&X);
353: ierr=MatGetSize(mat,&M,&N);
354: ierr=MatGetLocalSize(mat,&m,&n);
355: MatCreateDense(((PetscObject)mat)->comm,m,m,N,N,NULL,NewMat);
356: MatGetOwnershipRange(*NewMat,&low,&high);
357: for (i=0;i<M;i++){
358: MatGetColumnVector_ADA(mat,X,i);
359: VecGetArrayRead(X,&dptr);
360: for (j=0; j<high-low; j++){
361: MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
362: }
363: VecRestoreArrayRead(X,&dptr);
364: }
365: MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
366: MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
367: VecDestroy(&X);
368: } else if (isseqdense && size==1){
369: PetscInt i,j,low,high,m,n,M,N;
370: const PetscScalar *dptr;
371: Vec X;
373: VecDuplicate(ctx->D2,&X);
374: MatGetSize(mat,&M,&N);
375: MatGetLocalSize(mat,&m,&n);
376: MatCreateSeqDense(((PetscObject)mat)->comm,N,N,NULL,NewMat);
377: MatGetOwnershipRange(*NewMat,&low,&high);
378: for (i=0;i<M;i++){
379: MatGetColumnVector_ADA(mat,X,i);
380: VecGetArrayRead(X,&dptr);
381: for (j=0; j<high-low; j++){
382: MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
383: }
384: VecRestoreArrayRead(X,&dptr);
385: }
386: MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
387: MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
388: VecDestroy(&X);
389: } else SETERRQ(PETSC_COMM_SELF,1,"No support to convert objects to that type");
390: return(0);
391: }
395: static PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
396: {
398: TaoMatADACtx ctx;
401: MatShellGetContext(mat,(void **)&ctx);
402: if (type == NORM_FROBENIUS) {
403: *norm = 1.0;
404: } else if (type == NORM_1 || type == NORM_INFINITY) {
405: *norm = 1.0;
406: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
407: return(0);
408: }
412: /*@C
413: MatCreateADA - Creates a matrix M=A^T D1 A + D2 where D1, D2 are diagonal
415: Collective on matrix
417: Input Parameters:
418: + mat - matrix of arbitrary type
419: . d1 - A vector with diagonal elements of D1
420: - d2 - A vector
422: Output Parameters:
423: . J - New matrix whose operations are defined in terms of mat, D1, and D2.
425: Notes:
426: The user provides the input data and is responsible for destroying
427: this data after matrix J has been destroyed.
428: The operation MatMult(A,D2,D1) must be well defined.
429: Before calling the operation MatGetDiagonal(), the function
430: MatADAComputeDiagonal() must be called. The matrices A and D1 must
431: be the same during calls to MatADAComputeDiagonal() and
432: MatGetDiagonal().
434: Level: developer
436: .seealso: MatCreate()
437: @*/
438: PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J)
439: {
440: MPI_Comm comm=((PetscObject)mat)->comm;
441: TaoMatADACtx ctx;
443: PetscInt nloc,n;
446: PetscNew(&ctx);
447: ctx->A=mat;
448: ctx->D1=d1;
449: ctx->D2=d2;
450: if (d1){
451: VecDuplicate(d1,&ctx->W);
452: PetscObjectReference((PetscObject)d1);
453: } else {
454: ctx->W = NULL;
455: }
456: if (d2){
457: VecDuplicate(d2,&ctx->W2);
458: VecDuplicate(d2,&ctx->ADADiag);
459: PetscObjectReference((PetscObject)d2);
460: } else {
461: ctx->W2 = NULL;
462: ctx->ADADiag = NULL;
463: }
465: ctx->GotDiag=0;
466: PetscObjectReference((PetscObject)mat);
468: ierr=VecGetLocalSize(d2,&nloc);
469: ierr=VecGetSize(d2,&n);
471: MatCreateShell(comm,nloc,nloc,n,n,ctx,J);
473: MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);
474: MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);
475: MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);
476: MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);
477: MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);
478: MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);
479: MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);
480: MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);
481: MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);
482: MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);
483: MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)(void))MatGetSubMatrices_ADA);
484: MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);
485: MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);
486: MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)(void))MatGetSubMatrix_ADA);
488: PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);
489: PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));
491: MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);
492: return(0);
493: }