Actual source code: normalmat.c
1: #include "normalmat.h" /*I "mat.h" I*/
5: /*@C
6: MatCreateADA - Creates a matrix M=A^T D1 A + D2.
8: Collective on matrix
10: Input Parameters:
11: + mat - matrix of arbitrary type
12: . D1 - A vector
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: int MatCreateADA(Mat mat,Vec D1, Vec D2, Mat *J)
32: {
33: MPI_Comm comm=mat->comm;
34: TaoMatADACtx ctx;
35: int info,nloc,n;
38: /*
39: info=MatCheckVecs(mat,D1,D2,&flg);CHKERRQ(info);
40: if (flg==PETSC_FALSE){
41: SETERRQ(PETSC_ERR_SUP,"InCompatible matrix and vector for ADA^T matrix");
42: }
43: */
44: info = PetscNew(_p_TaoMatADACtx,&ctx);CHKERRQ(info);
46: ctx->A=mat;
47: ctx->D1=D1;
48: ctx->D2=D2;
49: if (D1){
50: info = VecDuplicate(D1,&ctx->W);CHKERRQ(info);
51: info = PetscObjectReference((PetscObject)D1);CHKERRQ(info);
52: } else {
53: ctx->W=0;
54: }
55: if (D2){
56: info = VecDuplicate(D2,&ctx->W2);CHKERRQ(info);
57: info = VecDuplicate(D2,&ctx->ADADiag);CHKERRQ(info);
58: info = PetscObjectReference((PetscObject)D2);CHKERRQ(info);
59: } else {
60: ctx->W2=0;
61: ctx->ADADiag=0;
62: }
64: ctx->GotDiag=0;
65: info = PetscObjectReference((PetscObject)mat);CHKERRQ(info);
67: info=VecGetLocalSize(D2,&nloc);CHKERRQ(info);
68: info=VecGetSize(D2,&n);CHKERRQ(info);
70: info = MatCreateShell(comm,nloc,nloc,n,n,ctx,J);CHKERRQ(info);
72: info = MatShellSetOperation(*J,MATOP_MULT,(void(*)())MatMult_ADA);CHKERRQ(info);
73: info = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)())MatDestroy_ADA);CHKERRQ(info);
74: info = MatShellSetOperation(*J,MATOP_VIEW,(void(*)())MatView_ADA);CHKERRQ(info);
75: info = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)())MatMultTranspose_ADA);CHKERRQ(info);
76: info = MatShellSetOperation(*J,MATOP_DIAGONAL_SHIFT,(void(*)())MatDiagonalShift_ADA);CHKERRQ(info);
77: info = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)())MatShift_ADA);CHKERRQ(info);
78: info = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)())MatEqual_ADA);CHKERRQ(info);
79: info = MatShellSetOperation(*J,MATOP_SCALE,(void(*)())MatScale_ADA);CHKERRQ(info);
80: info = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)())MatTranspose_ADA);CHKERRQ(info);
81: info = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_ADA);CHKERRQ(info);
82: info = MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)())MatGetSubMatrices_ADA);CHKERRQ(info);
83: info = MatShellSetOperation(*J,MATOP_NORM,(void(*)())MatNorm_ADA);CHKERRQ(info);
84: info = MatShellSetOperation(*J,MATOP_DUPLCIATE,(void(*)())MatDuplicate_ADA);CHKERRQ(info);
85: info = MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)())MatGetSubMatrix_ADA);CHKERRQ(info);
87: info = PetscLogObjectParent(*J,ctx->W); CHKERRQ(info);
88: info = PetscLogObjectParent(mat,*J); CHKERRQ(info);
90: info = MatSetOption(*J,MAT_SYMMETRIC);CHKERRQ(info);
91: return(0);
92: }
96: int MatMult_ADA(Mat mat,Vec a,Vec y)
97: {
98: TaoMatADACtx ctx;
99: PetscScalar one = 1.0;
100: int info;
103: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
105: info = MatMult(ctx->A,a,ctx->W);CHKERRQ(info);
106: if (ctx->D1){
107: info = VecPointwiseMult(ctx->W,ctx->D1,ctx->W);CHKERRQ(info);
108: }
109: info = MatMultTranspose(ctx->A,ctx->W,y);CHKERRQ(info);
110: if (ctx->D2){
111: info = VecPointwiseMult(ctx->W2, ctx->D2, a);CHKERRQ(info);
112: info = VecAXPY(y, one, ctx->W2);CHKERRQ(info);
113: }
114: return(0);
115: }
119: int MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
120: {
121: int info;
124: info = MatMult_ADA(mat,a,y);CHKERRQ(info);
125: return(0);
126: }
130: int MatDiagonalShift_ADA(Vec D, Mat M)
131: {
132: TaoMatADACtx ctx;
133: PetscScalar zero=0.0,one = 1.0;
134: int info;
137: info = MatShellGetContext(M,(void **)&ctx);CHKERRQ(info);
139: if (ctx->D2==PETSC_NULL){
140: info = VecDuplicate(D,&ctx->D2);CHKERRQ(info);
141: info = VecSet(ctx->D2, zero);CHKERRQ(info);
142: }
143: info = VecAXPY(ctx->D2, one, D);CHKERRQ(info);
145: return(0);
146: }
150: int MatDestroy_ADA(Mat mat)
151: {
152: int info;
153: TaoMatADACtx ctx;
156: info=MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
157: info=VecDestroy(ctx->W);CHKERRQ(info);
158: info=VecDestroy(ctx->W2);CHKERRQ(info);
159: info=VecDestroy(ctx->ADADiag);CHKERRQ(info);
160: info=MatDestroy(ctx->A);CHKERRQ(info);
161: info=VecDestroy(ctx->D1);CHKERRQ(info);
162: info=VecDestroy(ctx->D2);CHKERRQ(info);
163: PetscFree(ctx);
164: return(0);
165: }
169: int MatView_ADA(Mat mat,PetscViewer viewer)
170: {
173: /*
174: info = ViewerGetFormat(viewer,&format);CHKERRQ(info);
175: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_LONG) {
176: return(0); / * do nothing for now * /
177: }
178: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
179: info = MatView(ctx->A,viewer);CHKERRQ(info);
180: if (ctx->D1){
181: info = VecView(ctx->D1,viewer);CHKERRQ(info);
182: }
183: if (ctx->D2){
184: info = VecView(ctx->D2,viewer);CHKERRQ(info);
185: }
186: */
187: return(0);
188: }
192: int MatShift_ADA(Mat Y, PetscScalar a)
193: {
194: int info;
195: TaoMatADACtx ctx;
198: info = MatShellGetContext(Y,(void **)&ctx);CHKERRQ(info);
199: info = VecShift(ctx->D2,a);CHKERRQ(info);
200: return(0);
201: }
205: int MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
206: {
207: int info;
208: TaoMatADACtx ctx;
209: Mat A2;
210: Vec D1b=NULL,D2b;
213: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
214: info = MatDuplicate(ctx->A,op,&A2);CHKERRQ(info);
215: if (ctx->D1){
216: info = VecDuplicate(ctx->D1,&D1b);CHKERRQ(info);
217: info = VecCopy(ctx->D1,D1b);CHKERRQ(info);
218: }
219: info = VecDuplicate(ctx->D2,&D2b);CHKERRQ(info);
220: info = VecCopy(ctx->D2,D2b);CHKERRQ(info);
221: info = MatCreateADA(A2,D1b,D2b,M);CHKERRQ(info);
222: if (ctx->D1){
223: info = PetscObjectDereference((PetscObject)D1b);CHKERRQ(info);
224: }
225: info = PetscObjectDereference((PetscObject)D2b);CHKERRQ(info);
226: info = PetscObjectDereference((PetscObject)A2);CHKERRQ(info);
228: return(0);
229: }
233: int MatEqual_ADA(Mat A,Mat B,PetscTruth *flg)
234: {
235: int info;
236: TaoMatADACtx ctx1,ctx2;
239: info = MatShellGetContext(A,(void **)&ctx1);CHKERRQ(info);
240: info = MatShellGetContext(B,(void **)&ctx2);CHKERRQ(info);
241: info = VecEqual(ctx1->D2,ctx2->D2,flg);CHKERRQ(info);
242: if (*flg==PETSC_TRUE){
243: info = VecEqual(ctx1->D1,ctx2->D1,flg);CHKERRQ(info);
244: }
245: if (*flg==PETSC_TRUE){
246: info = MatEqual(ctx1->A,ctx2->A,flg);CHKERRQ(info);
247: }
248: return(0);
249: }
253: int MatScale_ADA(Mat mat, PetscScalar a)
254: {
255: int info;
256: TaoMatADACtx ctx;
259: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
260: info = VecScale(ctx->D1,a);CHKERRQ(info);
261: if (ctx->D2){
262: info = VecScale(ctx->D2,a);CHKERRQ(info);
263: }
264: return(0);
265: }
269: int MatTranspose_ADA(Mat mat,Mat *B)
270: {
271: int info;
272: TaoMatADACtx ctx;
275: if (*B){
276: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
277: info = MatDuplicate(mat,MAT_COPY_VALUES,B);CHKERRQ(info);
278: }
279: return(0);
280: }
284: int MatADAComputeDiagonal(Mat mat)
285: {
286: int i,m,n,low,high,info;
287: PetscScalar *dtemp,*dptr;
288: TaoMatADACtx ctx;
291: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
293: info = MatGetOwnershipRange(mat, &low, &high);CHKERRQ(info);
294: info = MatGetSize(mat,&m,&n);CHKERRQ(info);
295:
296: info = PetscMalloc( n*sizeof(PetscScalar),&dtemp ); CHKERRQ(info);
298: for (i=0; i<n; i++){
299: info = MatGetColumnVector(ctx->A, ctx->W, i);CHKERRQ(info);
300: info = VecPointwiseMult(ctx->W,ctx->W,ctx->W);CHKERRQ(info);
301: info = VecDotBegin(ctx->D1, ctx->W,dtemp+i);CHKERRQ(info);
302: }
303: for (i=0; i<n; i++){
304: info = VecDotEnd(ctx->D1, ctx->W,dtemp+i);CHKERRQ(info);
305: }
307: info = VecGetArray(ctx->ADADiag,&dptr);CHKERRQ(info);
308: for (i=low; i<high; i++){
309: dptr[i-low]= dtemp[i];
310: }
311: info = VecRestoreArray(ctx->ADADiag,&dptr);CHKERRQ(info);
312: if (dtemp) PetscFree(dtemp);
313: return(0);
314: }
318: int MatGetDiagonal_ADA(Mat mat,Vec v)
319: {
320: int info;
321: PetscScalar one=1.0;
322: TaoMatADACtx ctx;
325: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
326: info = MatADAComputeDiagonal(mat);
327: info=VecCopy(ctx->ADADiag,v);CHKERRQ(info);
328: if (ctx->D2){
329: info=VecAXPY(v, one, ctx->D2);CHKERRQ(info);
330: }
332: return(0);
333: }
337: int MatGetSubMatrices_ADA(Mat A,int n, IS *irow,IS *icol,MatReuse scall,Mat **B)
338: {
339: int info,i;
342: if (scall == MAT_INITIAL_MATRIX) {
343: info = PetscMalloc( (n+1)*sizeof(Mat),B );CHKERRQ(info);
344: }
346: for ( i=0; i<n; i++ ) {
347: info = MatGetSubMatrix_ADA(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(info);
348: }
349: return(0);
350: }
354: int MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,
355: Mat *newmat)
356: {
357: int info,low,high;
358: int n,nlocal,i;
359: int *iptr;
360: PetscScalar *dptr,*ddptr,zero=0.0;
361: VecType type_name;
362: IS ISrow;
363: Vec D1,D2;
364: Mat Atemp;
365: TaoMatADACtx ctx;
369: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
371: info = MatGetOwnershipRange(ctx->A,&low,&high);CHKERRQ(info);
372: info = ISCreateStride(mat->comm,high-low,low,1,&ISrow);CHKERRQ(info);
373: info = MatGetSubMatrix(ctx->A,ISrow,iscol,csize,cll,&Atemp);CHKERRQ(info);
374: info = ISDestroy(ISrow);CHKERRQ(info);
376: if (ctx->D1){
377: info=VecDuplicate(ctx->D1,&D1);CHKERRQ(info);
378: info=VecCopy(ctx->D1,D1);CHKERRQ(info);
379: } else {
380: D1=PETSC_NULL;
381: }
383: if (ctx->D2){
384: info=ISGetSize(isrow,&n);CHKERRQ(info);
385: info=ISGetLocalSize(isrow,&nlocal);CHKERRQ(info);
386: info=VecCreate(ctx->D2->comm,&D2);CHKERRQ(info);
387: info=VecGetType(ctx->D2,&type_name);CHKERRQ(info);
388: info=VecSetSizes(D2,nlocal,n);CHKERRQ(info);
389: info=VecSetType(D2,type_name);CHKERRQ(info);
390: info=VecSet(D2, zero);CHKERRQ(info);
391: info=VecGetArray(ctx->D2, &dptr); CHKERRQ(info);
392: info=VecGetArray(D2, &ddptr); CHKERRQ(info);
393: info=ISGetIndices(isrow,&iptr); CHKERRQ(info);
394: for (i=0;i<nlocal;i++){
395: ddptr[i] = dptr[iptr[i]-low];
396: }
397: info=ISRestoreIndices(isrow,&iptr); CHKERRQ(info);
398: info=VecRestoreArray(D2, &ddptr); CHKERRQ(info);
399: info=VecRestoreArray(ctx->D2, &dptr); CHKERRQ(info);
400:
401: } else {
402: D2=PETSC_NULL;
403: }
405: info = MatCreateADA(Atemp,D1,D2,newmat);CHKERRQ(info);
406: info = MatShellGetContext(*newmat,(void **)&ctx);CHKERRQ(info);
407: info = PetscObjectDereference((PetscObject)Atemp);CHKERRQ(info);
408: if (ctx->D1){
409: info = PetscObjectDereference((PetscObject)D1);CHKERRQ(info);
410: }
411: if (ctx->D2){
412: info = PetscObjectDereference((PetscObject)D2);CHKERRQ(info);
413: }
414: return(0);
415: }
419: int MatGetRowADA(Mat mat,int row,int *ncols,int **cols,PetscScalar **vals)
420: {
421: int info,m,n;
424: info = MatGetSize(mat,&m,&n);CHKERRQ(info);
426: if (*ncols>0){
427: info = PetscMalloc( (*ncols)*sizeof(int),cols );CHKERRQ(info);
428: info = PetscMalloc( (*ncols)*sizeof(PetscScalar),vals );CHKERRQ(info);
429: } else {
430: *cols=PETSC_NULL;
431: *vals=PETSC_NULL;
432: }
433:
434: return(0);
435: }
439: int MatRestoreRowADA(Mat mat,int row,int *ncols,int **cols,PetscScalar **vals)
440: {
442: if (*ncols>0){
443: PetscFree(*cols);
444: PetscFree(*vals);
445: }
446: *cols=PETSC_NULL;
447: *vals=PETSC_NULL;
448: return(0);
449: }
453: int MatGetColumnVector_ADA(Mat mat,Vec Y, int col)
454: {
455: int info,low,high;
456: PetscScalar zero=0.0,one=1.0;
459: info=VecSet(Y, zero);CHKERRQ(info);
460: info=VecGetOwnershipRange(Y,&low,&high);CHKERRQ(info);
461: if (col>=low && col<high){
462: info=VecSetValue(Y,col,one,INSERT_VALUES);CHKERRQ(info);
463: }
464: info=VecAssemblyBegin(Y);CHKERRQ(info);
465: info=VecAssemblyEnd(Y);CHKERRQ(info);
466: info=MatMult_ADA(mat,Y,Y);CHKERRQ(info);
468: return(0);
469: }
471: int MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
472: {
473: int info,size;
474: TaoMatADACtx ctx;
477: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
478: MPI_Comm_size(mat->comm,&size);
480: if (newtype==MATSAME){
482: info=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);CHKERRQ(info);
484: } else if (newtype==MATMPIDENSE){
486: int i,j,low,high,m,n,M,N;
487: PetscScalar *dptr;
488: Vec X;
490: info = VecDuplicate(ctx->D2,&X);CHKERRQ(info);
491: info=MatGetSize(mat,&M,&N);CHKERRQ(info);
492: info=MatGetLocalSize(mat,&m,&n);CHKERRQ(info);
493: info = MatCreateMPIDense(mat->comm,m,m,N,N,PETSC_NULL,NewMat);
494: CHKERRQ(info);
495: info = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(info);
496: for (i=0;i<M;i++){
497: info = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(info);
498: info = VecGetArray(X,&dptr);CHKERRQ(info);
499: for (j=0; j<high-low; j++){
500: info = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(info);
501: }
502: info=VecRestoreArray(X,&dptr);CHKERRQ(info);
503: }
504: info=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
505: info=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
506: info = VecDestroy(X);CHKERRQ(info);
508: } else if (newtype==MATSEQDENSE && size==1){
510: int i,j,low,high,m,n,M,N;
511: PetscScalar *dptr;
512: Vec X;
514: info = VecDuplicate(ctx->D2,&X);CHKERRQ(info);
515: info = MatGetSize(mat,&M,&N);CHKERRQ(info);
516: info = MatGetLocalSize(mat,&m,&n);CHKERRQ(info);
517: info = MatCreateSeqDense(mat->comm,N,N,PETSC_NULL,NewMat);
518: CHKERRQ(info);
519: info = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(info);
520: for (i=0;i<M;i++){
521: info = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(info);
522: info = VecGetArray(X,&dptr);CHKERRQ(info);
523: for (j=0; j<high-low; j++){
524: info = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(info);
525: }
526: info=VecRestoreArray(X,&dptr);CHKERRQ(info);
527: }
528: info=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
529: info=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
530: info=VecDestroy(X);CHKERRQ(info);
532: } else {
533: SETERRQ(1,"No support to convert objects to that type");
534: }
535: return(0);
536: }
540: int MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
541: {
542: int info;
543: TaoMatADACtx ctx;
546: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
548: if (type == NORM_FROBENIUS) {
549: *norm = 1.0;
550: } else if (type == NORM_1 || type == NORM_INFINITY) {
551: *norm = 1.0;
552: } else {
553: SETERRQ(PETSC_ERR_SUP,"No two norm");
554: }
555: return(0);
556: }