Actual source code: adamat.c
1: #include "adamat.h" /*I "mat.h" I*/
5: /*@C
6: MatCreateADA - Creates a matrix M=A^T D1 A + D2 where D1, D2 are diagonal
7:
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;
35: PetscErrorCode ierr;
36: PetscInt nloc,n;
40: PetscNew(_p_TaoMatADACtx,&ctx);
42: ctx->A=mat;
43: ctx->D1=d1;
44: ctx->D2=d2;
45: if (d1){
46: VecDuplicate(d1,&ctx->W);
47: PetscObjectReference((PetscObject)d1);
48: } else {
49: ctx->W=0;
50: }
51: if (d2){
52: VecDuplicate(d2,&ctx->W2);
53: VecDuplicate(d2,&ctx->ADADiag);
54: PetscObjectReference((PetscObject)d2);
55: } else {
56: ctx->W2=0;
57: ctx->ADADiag=0;
58: }
60: ctx->GotDiag=0;
61: PetscObjectReference((PetscObject)mat);
63: ierr=VecGetLocalSize(d2,&nloc);
64: ierr=VecGetSize(d2,&n);
66: MatCreateShell(comm,nloc,nloc,n,n,ctx,J);
68: MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);
69: MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);
70: MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);
71: MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);
72: MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);
73: MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);
74: MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);
75: MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);
76: MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);
77: MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);
78: MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)(void))MatGetSubMatrices_ADA);
79: MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);
80: MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);
81: MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)(void))MatGetSubMatrix_ADA);
83: PetscLogObjectParent(*J,ctx->W);
84: PetscLogObjectParent(mat,*J);
86: MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);
87: return(0);
88: }
92: PetscErrorCode MatMult_ADA(Mat mat,Vec a,Vec y)
93: {
94: TaoMatADACtx ctx;
95: PetscReal one = 1.0;
96: PetscErrorCode ierr;
99: MatShellGetContext(mat,(void **)&ctx);
101: MatMult(ctx->A,a,ctx->W);
102: if (ctx->D1){
103: VecPointwiseMult(ctx->W,ctx->D1,ctx->W);
104: }
105: MatMultTranspose(ctx->A,ctx->W,y);
106: if (ctx->D2){
107: VecPointwiseMult(ctx->W2, ctx->D2, a);
108: VecAXPY(y, one, ctx->W2);
109: }
110: return(0);
111: }
115: PetscErrorCode MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
116: {
120: MatMult_ADA(mat,a,y);
121: return(0);
122: }
126: PetscErrorCode MatDiagonalSet_ADA(Vec D, Mat M)
127: {
128: TaoMatADACtx ctx;
129: PetscReal zero=0.0,one = 1.0;
130: PetscErrorCode ierr;
133: MatShellGetContext(M,(void **)&ctx);
135: if (ctx->D2==PETSC_NULL){
136: VecDuplicate(D,&ctx->D2);
137: VecSet(ctx->D2, zero);
138: }
139: VecAXPY(ctx->D2, one, D);
141: return(0);
142: }
146: PetscErrorCode MatDestroy_ADA(Mat mat)
147: {
148: PetscErrorCode ierr;
149: TaoMatADACtx ctx;
152: ierr=MatShellGetContext(mat,(void **)&ctx);
153: if (ctx->W) {
154: ierr=VecDestroy(&ctx->W);
155: }
156: if (ctx->W2) {
157: ierr=VecDestroy(&ctx->W2);
158: }
159: if (ctx->ADADiag) {
160: ierr=VecDestroy(&ctx->ADADiag);
161: }
162: if (ctx->A) {
163: ierr=MatDestroy(&ctx->A);
164: }
165: if (ctx->D1) {
166: ierr=VecDestroy(&ctx->D1);
167: }
168: if (ctx->D2) {
169: ierr=VecDestroy(&ctx->D2);
170: }
171: PetscFree(ctx);
172: return(0);
173: }
177: PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer)
178: {
181: return(0);
182: }
186: PetscErrorCode MatShift_ADA(Mat Y, PetscReal a)
187: {
189: TaoMatADACtx ctx;
192: MatShellGetContext(Y,(void **)&ctx);
193: VecShift(ctx->D2,a);
194: return(0);
195: }
199: PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
200: {
201: PetscErrorCode ierr;
202: TaoMatADACtx ctx;
203: Mat A2;
204: Vec D1b=NULL,D2b;
207: MatShellGetContext(mat,(void **)&ctx);
208: MatDuplicate(ctx->A,op,&A2);
209: if (ctx->D1){
210: VecDuplicate(ctx->D1,&D1b);
211: VecCopy(ctx->D1,D1b);
212: }
213: VecDuplicate(ctx->D2,&D2b);
214: VecCopy(ctx->D2,D2b);
215: MatCreateADA(A2,D1b,D2b,M);
216: if (ctx->D1){
217: PetscObjectDereference((PetscObject)D1b);
218: }
219: PetscObjectDereference((PetscObject)D2b);
220: PetscObjectDereference((PetscObject)A2);
222: return(0);
223: }
227: PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg)
228: {
229: PetscErrorCode ierr;
230: TaoMatADACtx ctx1,ctx2;
233: MatShellGetContext(A,(void **)&ctx1);
234: MatShellGetContext(B,(void **)&ctx2);
235: VecEqual(ctx1->D2,ctx2->D2,flg);
236: if (*flg==PETSC_TRUE){
237: VecEqual(ctx1->D1,ctx2->D1,flg);
238: }
239: if (*flg==PETSC_TRUE){
240: MatEqual(ctx1->A,ctx2->A,flg);
241: }
242: return(0);
243: }
247: PetscErrorCode MatScale_ADA(Mat mat, PetscReal a)
248: {
249: PetscErrorCode ierr;
250: TaoMatADACtx ctx;
253: MatShellGetContext(mat,(void **)&ctx);
254: VecScale(ctx->D1,a);
255: if (ctx->D2){
256: VecScale(ctx->D2,a);
257: }
258: return(0);
259: }
263: PetscErrorCode MatTranspose_ADA(Mat mat,Mat *B)
264: {
266: TaoMatADACtx ctx;
269: if (*B){
270: MatShellGetContext(mat,(void **)&ctx);
271: MatDuplicate(mat,MAT_COPY_VALUES,B);
272: }
273: return(0);
274: }
278: PetscErrorCode MatADAComputeDiagonal(Mat mat)
279: {
281: PetscInt i,m,n,low,high;
282: PetscReal *dtemp,*dptr;
283: TaoMatADACtx ctx;
286: MatShellGetContext(mat,(void **)&ctx);
288: MatGetOwnershipRange(mat, &low, &high);
289: MatGetSize(mat,&m,&n);
290:
291: PetscMalloc( n*sizeof(PetscReal),&dtemp );
293: for (i=0; i<n; i++){
294: MatGetColumnVector(ctx->A, ctx->W, i);
295: VecPointwiseMult(ctx->W,ctx->W,ctx->W);
296: VecDotBegin(ctx->D1, ctx->W,dtemp+i);
297: }
298: for (i=0; i<n; i++){
299: VecDotEnd(ctx->D1, ctx->W,dtemp+i);
300: }
302: VecGetArray(ctx->ADADiag,&dptr);
303: for (i=low; i<high; i++){
304: dptr[i-low]= dtemp[i];
305: }
306: VecRestoreArray(ctx->ADADiag,&dptr);
307: if (dtemp) {
308: PetscFree(dtemp);
309: }
310: return(0);
311: }
315: PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v)
316: {
317: PetscErrorCode ierr;
318: PetscReal one=1.0;
319: TaoMatADACtx ctx;
322: MatShellGetContext(mat,(void **)&ctx);
323: MatADAComputeDiagonal(mat);
324: ierr=VecCopy(ctx->ADADiag,v);
325: if (ctx->D2){
326: ierr=VecAXPY(v, one, ctx->D2);
327: }
329: return(0);
330: }
334: PetscErrorCode MatGetSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
335: {
337: PetscInt i;
340: if (scall == MAT_INITIAL_MATRIX) {
341: PetscMalloc( (n+1)*sizeof(Mat),B );
342: }
344: for ( i=0; i<n; i++ ) {
345: MatGetSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);
346: }
347: return(0);
348: }
352: PetscErrorCode MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll,
353: Mat *newmat)
354: {
355: PetscErrorCode ierr;
356: PetscInt low,high;
357: PetscInt n,nlocal,i;
358: const PetscInt *iptr;
359: PetscReal *dptr,*ddptr,zero=0.0;
360: const VecType type_name;
361: IS ISrow;
362: Vec D1,D2;
363: Mat Atemp;
364: TaoMatADACtx ctx;
368: MatShellGetContext(mat,(void **)&ctx);
370: MatGetOwnershipRange(ctx->A,&low,&high);
371: ISCreateStride(((PetscObject)mat)->comm,high-low,low,1,&ISrow);
372: MatGetSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);
373: ISDestroy(&ISrow);
375: if (ctx->D1){
376: ierr=VecDuplicate(ctx->D1,&D1);
377: ierr=VecCopy(ctx->D1,D1);
378: } else {
379: D1=PETSC_NULL;
380: }
382: if (ctx->D2){
383: ierr=ISGetSize(isrow,&n);
384: ierr=ISGetLocalSize(isrow,&nlocal);
385: ierr=VecCreate(((PetscObject)(ctx->D2))->comm,&D2);
386: ierr=VecGetType(ctx->D2,&type_name);
387: ierr=VecSetSizes(D2,nlocal,n);
388: ierr=VecSetType(D2,type_name);
389: ierr=VecSet(D2, zero);
390: ierr=VecGetArray(ctx->D2, &dptr);
391: ierr=VecGetArray(D2, &ddptr);
392: ierr=ISGetIndices(isrow,&iptr);
393: for (i=0;i<nlocal;i++){
394: ddptr[i] = dptr[iptr[i]-low];
395: }
396: ierr=ISRestoreIndices(isrow,&iptr);
397: ierr=VecRestoreArray(D2, &ddptr);
398: ierr=VecRestoreArray(ctx->D2, &dptr);
399:
400: } else {
401: D2=PETSC_NULL;
402: }
404: MatCreateADA(Atemp,D1,D2,newmat);
405: MatShellGetContext(*newmat,(void **)&ctx);
406: PetscObjectDereference((PetscObject)Atemp);
407: if (ctx->D1){
408: PetscObjectDereference((PetscObject)D1);
409: }
410: if (ctx->D2){
411: PetscObjectDereference((PetscObject)D2);
412: }
413: return(0);
414: }
418: PetscErrorCode MatGetRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscReal **vals)
419: {
421: PetscInt m,n;
424: MatGetSize(mat,&m,&n);
426: if (*ncols>0){
427: PetscMalloc( (*ncols)*sizeof(PetscInt),cols );
428: PetscMalloc( (*ncols)*sizeof(PetscReal),vals );
429: } else {
430: *cols=PETSC_NULL;
431: *vals=PETSC_NULL;
432: }
433:
434: return(0);
435: }
439: PetscErrorCode MatRestoreRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscReal **vals)
440: {
443: if (*ncols>0){
444: PetscFree(*cols);
445: PetscFree(*vals);
446: }
447: *cols=PETSC_NULL;
448: *vals=PETSC_NULL;
449: return(0);
450: }
454: PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col)
455: {
456: PetscErrorCode ierr;
457: PetscInt low,high;
458: PetscReal zero=0.0,one=1.0;
461: ierr=VecSet(Y, zero);
462: ierr=VecGetOwnershipRange(Y,&low,&high);
463: if (col>=low && col<high){
464: ierr=VecSetValue(Y,col,one,INSERT_VALUES);
465: }
466: ierr=VecAssemblyBegin(Y);
467: ierr=VecAssemblyEnd(Y);
468: ierr=MatMult_ADA(mat,Y,Y);
470: return(0);
471: }
473: PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
474: {
476: PetscMPIInt size;
477: PetscBool sametype, issame, ismpidense, isseqdense;
478: TaoMatADACtx ctx;
481: MatShellGetContext(mat,(void **)&ctx);
482: MPI_Comm_size(((PetscObject)mat)->comm,&size);
485: PetscTypeCompare((PetscObject)mat,newtype,&sametype);
486: PetscTypeCompare((PetscObject)mat,MATSAME,&issame);
487: PetscTypeCompare((PetscObject)mat,MATMPIDENSE,&ismpidense);
488: PetscTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);
490: if (sametype || issame) {
492: ierr=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);
494: } else if (ismpidense) {
496: PetscInt i,j,low,high,m,n,M,N;
497: PetscReal *dptr;
498: Vec X;
500: VecDuplicate(ctx->D2,&X);
501: ierr=MatGetSize(mat,&M,&N);
502: ierr=MatGetLocalSize(mat,&m,&n);
503: MatCreateMPIDense(((PetscObject)mat)->comm,m,m,N,N,PETSC_NULL,NewMat);
504:
505: MatGetOwnershipRange(*NewMat,&low,&high);
506: for (i=0;i<M;i++){
507: MatGetColumnVector_ADA(mat,X,i);
508: VecGetArray(X,&dptr);
509: for (j=0; j<high-low; j++){
510: MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
511: }
512: ierr=VecRestoreArray(X,&dptr);
513: }
514: ierr=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
515: ierr=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
516: VecDestroy(&X);
518: } else if (isseqdense && size==1){
520: PetscInt i,j,low,high,m,n,M,N;
521: PetscReal *dptr;
522: Vec X;
524: VecDuplicate(ctx->D2,&X);
525: MatGetSize(mat,&M,&N);
526: MatGetLocalSize(mat,&m,&n);
527: MatCreateSeqDense(((PetscObject)mat)->comm,N,N,PETSC_NULL,NewMat);
528:
529: MatGetOwnershipRange(*NewMat,&low,&high);
530: for (i=0;i<M;i++){
531: MatGetColumnVector_ADA(mat,X,i);
532: VecGetArray(X,&dptr);
533: for (j=0; j<high-low; j++){
534: MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
535: }
536: ierr=VecRestoreArray(X,&dptr);
537: }
538: ierr=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
539: ierr=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
540: ierr=VecDestroy(&X);
542: } else {
543: SETERRQ(PETSC_COMM_SELF,1,"No support to convert objects to that type");
544: }
545: return(0);
546: }
550: PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
551: {
553: TaoMatADACtx ctx;
556: MatShellGetContext(mat,(void **)&ctx);
558: if (type == NORM_FROBENIUS) {
559: *norm = 1.0;
560: } else if (type == NORM_1 || type == NORM_INFINITY) {
561: *norm = 1.0;
562: } else {
563: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
564: }
565: return(0);
566: }