Actual source code: adamat.c
petsc-3.5.4 2015-05-23
1: #include <../src/tao/matrix/adamat.h> /*I "mat.h" I*/
5: /*@C
6: MatCreateADA - Creates a matrix M=A^T D1 A + D2 where D1, D2 are diagonal
8: Collective on matrix
10: Input Parameters:
11: + mat - matrix of arbitrary type
12: . d1 - A vector with diagonal elements of D1
13: - d2 - A vector
15: Output Parameters:
16: . J - New matrix whose operations are defined in terms of mat, D1, and D2.
18: Notes:
19: The user provides the input data and is responsible for destroying
20: this data after matrix J has been destroyed.
21: The operation MatMult(A,D2,D1) must be well defined.
22: Before calling the operation MatGetDiagonal(), the function
23: MatADAComputeDiagonal() must be called. The matrices A and D1 must
24: be the same during calls to MatADAComputeDiagonal() and
25: MatGetDiagonal().
27: Level: developer
29: .seealso: MatCreate()
30: @*/
31: PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J)
32: {
33: MPI_Comm comm=((PetscObject)mat)->comm;
34: TaoMatADACtx ctx;
36: PetscInt nloc,n;
39: PetscNew(&ctx);
40: ctx->A=mat;
41: ctx->D1=d1;
42: ctx->D2=d2;
43: if (d1){
44: VecDuplicate(d1,&ctx->W);
45: PetscObjectReference((PetscObject)d1);
46: } else {
47: ctx->W = NULL;
48: }
49: if (d2){
50: VecDuplicate(d2,&ctx->W2);
51: VecDuplicate(d2,&ctx->ADADiag);
52: PetscObjectReference((PetscObject)d2);
53: } else {
54: ctx->W2 = NULL;
55: ctx->ADADiag = NULL;
56: }
58: ctx->GotDiag=0;
59: PetscObjectReference((PetscObject)mat);
61: ierr=VecGetLocalSize(d2,&nloc);
62: ierr=VecGetSize(d2,&n);
64: MatCreateShell(comm,nloc,nloc,n,n,ctx,J);
66: MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);
67: MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);
68: MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);
69: MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);
70: MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);
71: MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);
72: MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);
73: MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);
74: MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);
75: MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);
76: MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)(void))MatGetSubMatrices_ADA);
77: MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);
78: MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);
79: MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)(void))MatGetSubMatrix_ADA);
81: PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);
82: PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));
84: MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);
85: return(0);
86: }
90: PetscErrorCode MatMult_ADA(Mat mat,Vec a,Vec y)
91: {
92: TaoMatADACtx ctx;
93: PetscReal one = 1.0;
97: MatShellGetContext(mat,(void **)&ctx);
98: MatMult(ctx->A,a,ctx->W);
99: if (ctx->D1){
100: VecPointwiseMult(ctx->W,ctx->D1,ctx->W);
101: }
102: MatMultTranspose(ctx->A,ctx->W,y);
103: if (ctx->D2){
104: VecPointwiseMult(ctx->W2, ctx->D2, a);
105: VecAXPY(y, one, ctx->W2);
106: }
107: return(0);
108: }
112: PetscErrorCode MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
113: {
117: MatMult_ADA(mat,a,y);
118: return(0);
119: }
123: PetscErrorCode MatDiagonalSet_ADA(Vec D, Mat M)
124: {
125: TaoMatADACtx ctx;
126: PetscReal zero=0.0,one = 1.0;
130: MatShellGetContext(M,(void **)&ctx);
131: if (!ctx->D2){
132: VecDuplicate(D,&ctx->D2);
133: VecSet(ctx->D2, zero);
134: }
135: VecAXPY(ctx->D2, one, D);
136: return(0);
137: }
141: PetscErrorCode MatDestroy_ADA(Mat mat)
142: {
144: TaoMatADACtx ctx;
147: ierr=MatShellGetContext(mat,(void **)&ctx);
148: ierr=VecDestroy(&ctx->W);
149: ierr=VecDestroy(&ctx->W2);
150: ierr=VecDestroy(&ctx->ADADiag);
151: ierr=MatDestroy(&ctx->A);
152: ierr=VecDestroy(&ctx->D1);
153: ierr=VecDestroy(&ctx->D2);
154: PetscFree(ctx);
155: return(0);
156: }
160: PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer)
161: {
163: return(0);
164: }
168: PetscErrorCode MatShift_ADA(Mat Y, PetscReal a)
169: {
171: TaoMatADACtx ctx;
174: MatShellGetContext(Y,(void **)&ctx);
175: VecShift(ctx->D2,a);
176: return(0);
177: }
181: PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
182: {
183: PetscErrorCode ierr;
184: TaoMatADACtx ctx;
185: Mat A2;
186: Vec D1b=NULL,D2b;
189: MatShellGetContext(mat,(void **)&ctx);
190: MatDuplicate(ctx->A,op,&A2);
191: if (ctx->D1){
192: VecDuplicate(ctx->D1,&D1b);
193: VecCopy(ctx->D1,D1b);
194: }
195: VecDuplicate(ctx->D2,&D2b);
196: VecCopy(ctx->D2,D2b);
197: MatCreateADA(A2,D1b,D2b,M);
198: if (ctx->D1){
199: PetscObjectDereference((PetscObject)D1b);
200: }
201: PetscObjectDereference((PetscObject)D2b);
202: PetscObjectDereference((PetscObject)A2);
203: return(0);
204: }
208: PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg)
209: {
211: TaoMatADACtx ctx1,ctx2;
214: MatShellGetContext(A,(void **)&ctx1);
215: MatShellGetContext(B,(void **)&ctx2);
216: VecEqual(ctx1->D2,ctx2->D2,flg);
217: if (*flg==PETSC_TRUE){
218: VecEqual(ctx1->D1,ctx2->D1,flg);
219: }
220: if (*flg==PETSC_TRUE){
221: MatEqual(ctx1->A,ctx2->A,flg);
222: }
223: return(0);
224: }
228: PetscErrorCode MatScale_ADA(Mat mat, PetscReal a)
229: {
231: TaoMatADACtx ctx;
234: MatShellGetContext(mat,(void **)&ctx);
235: VecScale(ctx->D1,a);
236: if (ctx->D2){
237: VecScale(ctx->D2,a);
238: }
239: return(0);
240: }
244: PetscErrorCode MatTranspose_ADA(Mat mat,Mat *B)
245: {
247: TaoMatADACtx ctx;
250: if (*B){
251: MatShellGetContext(mat,(void **)&ctx);
252: MatDuplicate(mat,MAT_COPY_VALUES,B);
253: }
254: return(0);
255: }
259: PetscErrorCode MatADAComputeDiagonal(Mat mat)
260: {
262: PetscInt i,m,n,low,high;
263: PetscScalar *dtemp,*dptr;
264: TaoMatADACtx ctx;
267: MatShellGetContext(mat,(void **)&ctx);
268: MatGetOwnershipRange(mat, &low, &high);
269: MatGetSize(mat,&m,&n);
271: PetscMalloc1(n,&dtemp );
273: for (i=0; i<n; i++){
274: MatGetColumnVector(ctx->A, ctx->W, i);
275: VecPointwiseMult(ctx->W,ctx->W,ctx->W);
276: VecDotBegin(ctx->D1, ctx->W,dtemp+i);
277: }
278: for (i=0; i<n; i++){
279: VecDotEnd(ctx->D1, ctx->W,dtemp+i);
280: }
282: VecGetArray(ctx->ADADiag,&dptr);
283: for (i=low; i<high; i++){
284: dptr[i-low]= dtemp[i];
285: }
286: VecRestoreArray(ctx->ADADiag,&dptr);
287: PetscFree(dtemp);
288: return(0);
289: }
293: PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v)
294: {
295: PetscErrorCode ierr;
296: PetscReal one=1.0;
297: TaoMatADACtx ctx;
300: MatShellGetContext(mat,(void **)&ctx);
301: MatADAComputeDiagonal(mat);
302: ierr=VecCopy(ctx->ADADiag,v);
303: if (ctx->D2){
304: ierr=VecAXPY(v, one, ctx->D2);
305: }
306: return(0);
307: }
311: PetscErrorCode MatGetSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
312: {
314: PetscInt i;
317: if (scall == MAT_INITIAL_MATRIX) {
318: PetscMalloc1(n+1,B );
319: }
320: for ( i=0; i<n; i++ ) {
321: MatGetSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);
322: }
323: return(0);
324: }
328: PetscErrorCode MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll, Mat *newmat)
329: {
330: PetscErrorCode ierr;
331: PetscInt low,high;
332: PetscInt n,nlocal,i;
333: const PetscInt *iptr;
334: const PetscScalar *dptr;
335: PetscScalar *ddptr,zero=0.0;
336: VecType type_name;
337: IS ISrow;
338: Vec D1,D2;
339: Mat Atemp;
340: TaoMatADACtx ctx;
343: MatShellGetContext(mat,(void **)&ctx);
345: MatGetOwnershipRange(ctx->A,&low,&high);
346: ISCreateStride(((PetscObject)mat)->comm,high-low,low,1,&ISrow);
347: MatGetSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);
348: ISDestroy(&ISrow);
350: if (ctx->D1){
351: ierr=VecDuplicate(ctx->D1,&D1);
352: ierr=VecCopy(ctx->D1,D1);
353: } else {
354: D1 = NULL;
355: }
357: if (ctx->D2){
358: ierr=ISGetSize(isrow,&n);
359: ierr=ISGetLocalSize(isrow,&nlocal);
360: ierr=VecCreate(((PetscObject)(ctx->D2))->comm,&D2);
361: ierr=VecGetType(ctx->D2,&type_name);
362: ierr=VecSetSizes(D2,nlocal,n);
363: ierr=VecSetType(D2,type_name);
364: ierr=VecSet(D2, zero);
365: ierr=VecGetArrayRead(ctx->D2, &dptr);
366: ierr=VecGetArray(D2, &ddptr);
367: ierr=ISGetIndices(isrow,&iptr);
368: for (i=0;i<nlocal;i++){
369: ddptr[i] = dptr[iptr[i]-low];
370: }
371: ierr=ISRestoreIndices(isrow,&iptr);
372: ierr=VecRestoreArray(D2, &ddptr);
373: ierr=VecRestoreArrayRead(ctx->D2, &dptr);
374: } else {
375: D2 = NULL;
376: }
378: MatCreateADA(Atemp,D1,D2,newmat);
379: MatShellGetContext(*newmat,(void **)&ctx);
380: PetscObjectDereference((PetscObject)Atemp);
381: if (ctx->D1){
382: PetscObjectDereference((PetscObject)D1);
383: }
384: if (ctx->D2){
385: PetscObjectDereference((PetscObject)D2);
386: }
387: return(0);
388: }
392: PetscErrorCode MatGetRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscReal **vals)
393: {
395: PetscInt m,n;
398: MatGetSize(mat,&m,&n);
400: if (*ncols>0){
401: PetscMalloc1(*ncols,cols );
402: PetscMalloc1(*ncols,vals );
403: } else {
404: *cols=NULL;
405: *vals=NULL;
406: }
407: return(0);
408: }
412: PetscErrorCode MatRestoreRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscReal **vals)
413: {
417: if (*ncols>0){
418: PetscFree(*cols);
419: PetscFree(*vals);
420: }
421: *cols=NULL;
422: *vals=NULL;
423: return(0);
424: }
428: PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col)
429: {
431: PetscInt low,high;
432: PetscScalar zero=0.0,one=1.0;
435: ierr=VecSet(Y, zero);
436: ierr=VecGetOwnershipRange(Y,&low,&high);
437: if (col>=low && col<high){
438: ierr=VecSetValue(Y,col,one,INSERT_VALUES);
439: }
440: ierr=VecAssemblyBegin(Y);
441: ierr=VecAssemblyEnd(Y);
442: ierr=MatMult_ADA(mat,Y,Y);
443: return(0);
444: }
446: PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
447: {
449: PetscMPIInt size;
450: PetscBool sametype, issame, isdense, isseqdense;
451: TaoMatADACtx ctx;
454: MatShellGetContext(mat,(void **)&ctx);
455: MPI_Comm_size(((PetscObject)mat)->comm,&size);
457: PetscObjectTypeCompare((PetscObject)mat,newtype,&sametype);
458: PetscObjectTypeCompare((PetscObject)mat,MATSAME,&issame);
459: PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&isdense);
460: PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);
462: if (sametype || issame) {
463: ierr=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);
464: } else if (isdense) {
465: PetscInt i,j,low,high,m,n,M,N;
466: const PetscScalar *dptr;
467: Vec X;
469: VecDuplicate(ctx->D2,&X);
470: ierr=MatGetSize(mat,&M,&N);
471: ierr=MatGetLocalSize(mat,&m,&n);
472: MatCreateDense(((PetscObject)mat)->comm,m,m,N,N,NULL,NewMat);
473: MatGetOwnershipRange(*NewMat,&low,&high);
474: for (i=0;i<M;i++){
475: MatGetColumnVector_ADA(mat,X,i);
476: VecGetArrayRead(X,&dptr);
477: for (j=0; j<high-low; j++){
478: MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
479: }
480: VecRestoreArrayRead(X,&dptr);
481: }
482: MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
483: MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
484: VecDestroy(&X);
485: } else if (isseqdense && size==1){
486: PetscInt i,j,low,high,m,n,M,N;
487: const PetscScalar *dptr;
488: Vec X;
490: VecDuplicate(ctx->D2,&X);
491: MatGetSize(mat,&M,&N);
492: MatGetLocalSize(mat,&m,&n);
493: MatCreateSeqDense(((PetscObject)mat)->comm,N,N,NULL,NewMat);
494: MatGetOwnershipRange(*NewMat,&low,&high);
495: for (i=0;i<M;i++){
496: MatGetColumnVector_ADA(mat,X,i);
497: VecGetArrayRead(X,&dptr);
498: for (j=0; j<high-low; j++){
499: MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
500: }
501: VecRestoreArrayRead(X,&dptr);
502: }
503: MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
504: MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
505: VecDestroy(&X);
506: } else SETERRQ(PETSC_COMM_SELF,1,"No support to convert objects to that type");
507: return(0);
508: }
512: PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
513: {
515: TaoMatADACtx ctx;
518: MatShellGetContext(mat,(void **)&ctx);
519: if (type == NORM_FROBENIUS) {
520: *norm = 1.0;
521: } else if (type == NORM_1 || type == NORM_INFINITY) {
522: *norm = 1.0;
523: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
524: return(0);
525: }