Actual source code: adamat.c
petsc-3.8.4 2018-03-24
1: #include <petscmat.h>
3: typedef struct{
4: Mat A;
5: Vec D1;
6: Vec D2;
7: Vec W;
8: Vec W2;
9: Vec ADADiag;
10: PetscInt GotDiag;
11: } _p_TaoMatADACtx;
12: typedef _p_TaoMatADACtx* TaoMatADACtx;
14: extern PetscErrorCode MatCreateADA(Mat,Vec, Vec, Mat*);
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(Vec D, Mat M)
46: {
47: TaoMatADACtx ctx;
48: PetscReal zero=0.0,one = 1.0;
52: MatShellGetContext(M,(void **)&ctx);
53: if (!ctx->D2){
54: VecDuplicate(D,&ctx->D2);
55: VecSet(ctx->D2, zero);
56: }
57: VecAXPY(ctx->D2, one, D);
58: return(0);
59: }
61: static PetscErrorCode MatDestroy_ADA(Mat mat)
62: {
64: TaoMatADACtx ctx;
67: ierr=MatShellGetContext(mat,(void **)&ctx);
68: ierr=VecDestroy(&ctx->W);
69: ierr=VecDestroy(&ctx->W2);
70: ierr=VecDestroy(&ctx->ADADiag);
71: ierr=MatDestroy(&ctx->A);
72: ierr=VecDestroy(&ctx->D1);
73: ierr=VecDestroy(&ctx->D2);
74: PetscFree(ctx);
75: return(0);
76: }
78: static PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer)
79: {
81: return(0);
82: }
84: static PetscErrorCode MatShift_ADA(Mat Y, PetscReal a)
85: {
87: TaoMatADACtx ctx;
90: MatShellGetContext(Y,(void **)&ctx);
91: VecShift(ctx->D2,a);
92: return(0);
93: }
95: static PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
96: {
97: PetscErrorCode ierr;
98: TaoMatADACtx ctx;
99: Mat A2;
100: Vec D1b=NULL,D2b;
103: MatShellGetContext(mat,(void **)&ctx);
104: MatDuplicate(ctx->A,op,&A2);
105: if (ctx->D1){
106: VecDuplicate(ctx->D1,&D1b);
107: VecCopy(ctx->D1,D1b);
108: }
109: VecDuplicate(ctx->D2,&D2b);
110: VecCopy(ctx->D2,D2b);
111: MatCreateADA(A2,D1b,D2b,M);
112: if (ctx->D1){
113: PetscObjectDereference((PetscObject)D1b);
114: }
115: PetscObjectDereference((PetscObject)D2b);
116: PetscObjectDereference((PetscObject)A2);
117: return(0);
118: }
120: static PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg)
121: {
123: TaoMatADACtx ctx1,ctx2;
126: MatShellGetContext(A,(void **)&ctx1);
127: MatShellGetContext(B,(void **)&ctx2);
128: VecEqual(ctx1->D2,ctx2->D2,flg);
129: if (*flg==PETSC_TRUE){
130: VecEqual(ctx1->D1,ctx2->D1,flg);
131: }
132: if (*flg==PETSC_TRUE){
133: MatEqual(ctx1->A,ctx2->A,flg);
134: }
135: return(0);
136: }
138: static PetscErrorCode MatScale_ADA(Mat mat, PetscReal a)
139: {
141: TaoMatADACtx ctx;
144: MatShellGetContext(mat,(void **)&ctx);
145: VecScale(ctx->D1,a);
146: if (ctx->D2){
147: VecScale(ctx->D2,a);
148: }
149: return(0);
150: }
152: static PetscErrorCode MatTranspose_ADA(Mat mat,MatReuse reuse,Mat *B)
153: {
155: TaoMatADACtx ctx;
158: MatShellGetContext(mat,(void **)&ctx);
159: if (reuse == MAT_INITIAL_MATRIX){
160: MatDuplicate(mat,MAT_COPY_VALUES,B);
161: } else if (reuse == MAT_REUSE_MATRIX){
162: MatCopy(mat,*B,SAME_NONZERO_PATTERN);
163: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Does not support inplace transpose");
164: return(0);
165: }
167: PetscErrorCode MatADAComputeDiagonal(Mat mat)
168: {
170: PetscInt i,m,n,low,high;
171: PetscScalar *dtemp,*dptr;
172: TaoMatADACtx ctx;
175: MatShellGetContext(mat,(void **)&ctx);
176: MatGetOwnershipRange(mat, &low, &high);
177: MatGetSize(mat,&m,&n);
179: 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: ierr=VecCopy(ctx->ADADiag,v);
209: if (ctx->D2){
210: ierr=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 with diagonal elements of D1
385: - d2 - A vector
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.
393: The operation MatMult(A,D2,D1) must be well defined.
394: Before calling the operation MatGetDiagonal(), the function
395: MatADAComputeDiagonal() must be called. The matrices A and D1 must
396: be the same during calls to MatADAComputeDiagonal() and
397: MatGetDiagonal().
399: Level: developer
401: .seealso: MatCreate()
402: @*/
403: PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J)
404: {
405: MPI_Comm comm=PetscObjectComm((PetscObject)mat);
406: TaoMatADACtx ctx;
408: PetscInt nloc,n;
411: PetscNew(&ctx);
412: ctx->A=mat;
413: ctx->D1=d1;
414: ctx->D2=d2;
415: if (d1){
416: VecDuplicate(d1,&ctx->W);
417: PetscObjectReference((PetscObject)d1);
418: } else {
419: ctx->W = NULL;
420: }
421: if (d2){
422: VecDuplicate(d2,&ctx->W2);
423: VecDuplicate(d2,&ctx->ADADiag);
424: PetscObjectReference((PetscObject)d2);
425: } else {
426: ctx->W2 = NULL;
427: ctx->ADADiag = NULL;
428: }
430: ctx->GotDiag=0;
431: PetscObjectReference((PetscObject)mat);
433: ierr=VecGetLocalSize(d2,&nloc);
434: ierr=VecGetSize(d2,&n);
436: MatCreateShell(comm,nloc,nloc,n,n,ctx,J);
438: MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);
439: MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);
440: MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);
441: MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);
442: MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);
443: MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);
444: MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);
445: MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);
446: MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);
447: MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);
448: MatShellSetOperation(*J,MATOP_CREATE_SUBMATRICES,(void(*)(void))MatCreateSubMatrices_ADA);
449: MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);
450: MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);
451: MatShellSetOperation(*J,MATOP_CREATE_SUBMATRIX,(void(*)(void))MatCreateSubMatrix_ADA);
453: PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);
454: PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));
456: MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);
457: return(0);
458: }