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=((PetscObject)mat)->comm;
34: TaoMatADACtx ctx;
35: int info;
36: PetscInt nloc,n;
39: /*
40: info=MatCheckVecs(mat,D1,D2,&flg);CHKERRQ(info);
41: if (flg==PETSC_FALSE){
42: SETERRQ(PETSC_ERR_SUP,"InCompatible matrix and vector for ADA^T matrix");
43: }
44: */
45: info = PetscNew(_p_TaoMatADACtx,&ctx);CHKERRQ(info);
47: ctx->A=mat;
48: ctx->D1=D1;
49: ctx->D2=D2;
50: if (D1){
51: info = VecDuplicate(D1,&ctx->W);CHKERRQ(info);
52: info = PetscObjectReference((PetscObject)D1);CHKERRQ(info);
53: } else {
54: ctx->W=0;
55: }
56: if (D2){
57: info = VecDuplicate(D2,&ctx->W2);CHKERRQ(info);
58: info = VecDuplicate(D2,&ctx->ADADiag);CHKERRQ(info);
59: info = PetscObjectReference((PetscObject)D2);CHKERRQ(info);
60: } else {
61: ctx->W2=0;
62: ctx->ADADiag=0;
63: }
65: ctx->GotDiag=0;
66: info = PetscObjectReference((PetscObject)mat);CHKERRQ(info);
68: info=VecGetLocalSize(D2,&nloc);CHKERRQ(info);
69: info=VecGetSize(D2,&n);CHKERRQ(info);
71: info = MatCreateShell(comm,nloc,nloc,n,n,ctx,J);CHKERRQ(info);
73: info = MatShellSetOperation(*J,MATOP_MULT,(void(*)())MatMult_ADA);CHKERRQ(info);
74: info = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)())MatDestroy_ADA);CHKERRQ(info);
75: info = MatShellSetOperation(*J,MATOP_VIEW,(void(*)())MatView_ADA);CHKERRQ(info);
76: info = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)())MatMultTranspose_ADA);CHKERRQ(info);
77: info = MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)())MatDiagonalSet_ADA);CHKERRQ(info);
78: info = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)())MatShift_ADA);CHKERRQ(info);
79: info = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)())MatEqual_ADA);CHKERRQ(info);
80: info = MatShellSetOperation(*J,MATOP_SCALE,(void(*)())MatScale_ADA);CHKERRQ(info);
81: info = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)())MatTranspose_ADA);CHKERRQ(info);
82: info = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_ADA);CHKERRQ(info);
83: info = MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)())MatGetSubMatrices_ADA);CHKERRQ(info);
84: info = MatShellSetOperation(*J,MATOP_NORM,(void(*)())MatNorm_ADA);CHKERRQ(info);
85: info = MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)())MatDuplicate_ADA);CHKERRQ(info);
86: info = MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)())MatGetSubMatrix_ADA);CHKERRQ(info);
88: info = PetscLogObjectParent(*J,ctx->W); CHKERRQ(info);
89: info = PetscLogObjectParent(mat,*J); CHKERRQ(info);
91: info = MatSetOption(*J,MAT_SYMMETRIC, PETSC_TRUE);CHKERRQ(info);
93: return(0);
94: }
98: int MatMult_ADA(Mat mat,Vec a,Vec y)
99: {
100: TaoMatADACtx ctx;
101: PetscScalar one = 1.0;
102: int info;
105: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
107: info = MatMult(ctx->A,a,ctx->W);CHKERRQ(info);
108: if (ctx->D1){
109: info = VecPointwiseMult(ctx->W,ctx->D1,ctx->W);CHKERRQ(info);
110: }
111: info = MatMultTranspose(ctx->A,ctx->W,y);CHKERRQ(info);
112: if (ctx->D2){
113: info = VecPointwiseMult(ctx->W2, ctx->D2, a);CHKERRQ(info);
114: info = VecAXPY(y, one, ctx->W2);CHKERRQ(info);
115: }
116: return(0);
117: }
121: int MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
122: {
123: int info;
126: info = MatMult_ADA(mat,a,y);CHKERRQ(info);
127: return(0);
128: }
132: int MatDiagonalSet_ADA(Vec D, Mat M)
133: {
134: TaoMatADACtx ctx;
135: PetscScalar zero=0.0,one = 1.0;
136: int info;
139: info = MatShellGetContext(M,(void **)&ctx);CHKERRQ(info);
141: if (ctx->D2==PETSC_NULL){
142: info = VecDuplicate(D,&ctx->D2);CHKERRQ(info);
143: info = VecSet(ctx->D2, zero);CHKERRQ(info);
144: }
145: info = VecAXPY(ctx->D2, one, D);CHKERRQ(info);
147: return(0);
148: }
152: int MatDestroy_ADA(Mat mat)
153: {
154: int info;
155: TaoMatADACtx ctx;
158: info=MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
159: info=VecDestroy(ctx->W);CHKERRQ(info);
160: info=VecDestroy(ctx->W2);CHKERRQ(info);
161: info=VecDestroy(ctx->ADADiag);CHKERRQ(info);
162: info=MatDestroy(ctx->A);CHKERRQ(info);
163: info=VecDestroy(ctx->D1);CHKERRQ(info);
164: info=VecDestroy(ctx->D2);CHKERRQ(info);
165: info = PetscFree(ctx); CHKERRQ(info);
166: return(0);
167: }
171: int MatView_ADA(Mat mat,PetscViewer viewer)
172: {
175: /*
176: info = ViewerGetFormat(viewer,&format);CHKERRQ(info);
177: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_LONG) {
178: return(0); / * do nothing for now * /
179: }
180: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
181: info = MatView(ctx->A,viewer);CHKERRQ(info);
182: if (ctx->D1){
183: info = VecView(ctx->D1,viewer);CHKERRQ(info);
184: }
185: if (ctx->D2){
186: info = VecView(ctx->D2,viewer);CHKERRQ(info);
187: }
188: */
189: return(0);
190: }
194: int MatShift_ADA(Mat Y, PetscScalar a)
195: {
196: int info;
197: TaoMatADACtx ctx;
200: info = MatShellGetContext(Y,(void **)&ctx);CHKERRQ(info);
201: info = VecShift(ctx->D2,a);CHKERRQ(info);
202: return(0);
203: }
207: int MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
208: {
209: int info;
210: TaoMatADACtx ctx;
211: Mat A2;
212: Vec D1b=NULL,D2b;
215: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
216: info = MatDuplicate(ctx->A,op,&A2);CHKERRQ(info);
217: if (ctx->D1){
218: info = VecDuplicate(ctx->D1,&D1b);CHKERRQ(info);
219: info = VecCopy(ctx->D1,D1b);CHKERRQ(info);
220: }
221: info = VecDuplicate(ctx->D2,&D2b);CHKERRQ(info);
222: info = VecCopy(ctx->D2,D2b);CHKERRQ(info);
223: info = MatCreateADA(A2,D1b,D2b,M);CHKERRQ(info);
224: if (ctx->D1){
225: info = PetscObjectDereference((PetscObject)D1b);CHKERRQ(info);
226: }
227: info = PetscObjectDereference((PetscObject)D2b);CHKERRQ(info);
228: info = PetscObjectDereference((PetscObject)A2);CHKERRQ(info);
230: return(0);
231: }
235: int MatEqual_ADA(Mat A,Mat B,PetscTruth *flg)
236: {
237: int info;
238: TaoMatADACtx ctx1,ctx2;
241: info = MatShellGetContext(A,(void **)&ctx1);CHKERRQ(info);
242: info = MatShellGetContext(B,(void **)&ctx2);CHKERRQ(info);
243: info = VecEqual(ctx1->D2,ctx2->D2,flg);CHKERRQ(info);
244: if (*flg==PETSC_TRUE){
245: info = VecEqual(ctx1->D1,ctx2->D1,flg);CHKERRQ(info);
246: }
247: if (*flg==PETSC_TRUE){
248: info = MatEqual(ctx1->A,ctx2->A,flg);CHKERRQ(info);
249: }
250: return(0);
251: }
255: int MatScale_ADA(Mat mat, PetscScalar a)
256: {
257: int info;
258: TaoMatADACtx ctx;
261: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
262: info = VecScale(ctx->D1,a);CHKERRQ(info);
263: if (ctx->D2){
264: info = VecScale(ctx->D2,a);CHKERRQ(info);
265: }
266: return(0);
267: }
271: int MatTranspose_ADA(Mat mat,Mat *B)
272: {
273: int info;
274: TaoMatADACtx ctx;
277: if (*B){
278: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
279: info = MatDuplicate(mat,MAT_COPY_VALUES,B);CHKERRQ(info);
280: }
281: return(0);
282: }
286: int MatADAComputeDiagonal(Mat mat)
287: {
288: int info;
289: PetscInt i,m,n,low,high;
290: PetscScalar *dtemp,*dptr;
291: TaoMatADACtx ctx;
294: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
296: info = MatGetOwnershipRange(mat, &low, &high);CHKERRQ(info);
297: info = MatGetSize(mat,&m,&n);CHKERRQ(info);
298:
299: info = PetscMalloc( n*sizeof(PetscScalar),&dtemp ); CHKERRQ(info);
301: for (i=0; i<n; i++){
302: info = MatGetColumnVector(ctx->A, ctx->W, i);CHKERRQ(info);
303: info = VecPointwiseMult(ctx->W,ctx->W,ctx->W);CHKERRQ(info);
304: info = VecDotBegin(ctx->D1, ctx->W,dtemp+i);CHKERRQ(info);
305: }
306: for (i=0; i<n; i++){
307: info = VecDotEnd(ctx->D1, ctx->W,dtemp+i);CHKERRQ(info);
308: }
310: info = VecGetArray(ctx->ADADiag,&dptr);CHKERRQ(info);
311: for (i=low; i<high; i++){
312: dptr[i-low]= dtemp[i];
313: }
314: info = VecRestoreArray(ctx->ADADiag,&dptr);CHKERRQ(info);
315: if (dtemp) {
316: info = PetscFree(dtemp); CHKERRQ(info);
317: }
318: return(0);
319: }
323: int MatGetDiagonal_ADA(Mat mat,Vec v)
324: {
325: int info;
326: PetscScalar one=1.0;
327: TaoMatADACtx ctx;
330: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
331: info = MatADAComputeDiagonal(mat);
332: info=VecCopy(ctx->ADADiag,v);CHKERRQ(info);
333: if (ctx->D2){
334: info=VecAXPY(v, one, ctx->D2);CHKERRQ(info);
335: }
337: return(0);
338: }
342: int MatGetSubMatrices_ADA(Mat A,int n, IS *irow,IS *icol,MatReuse scall,Mat **B)
343: {
344: int info,i;
347: if (scall == MAT_INITIAL_MATRIX) {
348: info = PetscMalloc( (n+1)*sizeof(Mat),B );CHKERRQ(info);
349: }
351: for ( i=0; i<n; i++ ) {
352: info = MatGetSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(info);
353: }
354: return(0);
355: }
359: int MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll,
360: Mat *newmat)
361: {
362: int info;
363: PetscInt i,low,high;
364: PetscInt n,nlocal;
365: const PetscInt *iptr;
366: PetscScalar *dptr,*ddptr,zero=0.0;
367: const VecType type_name;
368: IS ISrow;
369: Vec D1,D2;
370: Mat Atemp;
371: TaoMatADACtx ctx;
375: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
377: info = MatGetOwnershipRange(ctx->A,&low,&high);CHKERRQ(info);
378: info = ISCreateStride(((PetscObject)mat)->comm,high-low,low,1,&ISrow);CHKERRQ(info);
379: info = MatGetSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);CHKERRQ(info);
380: info = ISDestroy(ISrow);CHKERRQ(info);
382: if (ctx->D1){
383: info=VecDuplicate(ctx->D1,&D1);CHKERRQ(info);
384: info=VecCopy(ctx->D1,D1);CHKERRQ(info);
385: } else {
386: D1=PETSC_NULL;
387: }
389: if (ctx->D2){
390: info=ISGetSize(isrow,&n);CHKERRQ(info);
391: info=ISGetLocalSize(isrow,&nlocal);CHKERRQ(info);
392: info=VecCreate(((PetscObject)ctx->D2)->comm,&D2);CHKERRQ(info);
393: info=VecGetType(ctx->D2,&type_name);CHKERRQ(info);
394: info=VecSetSizes(D2,nlocal,n);CHKERRQ(info);
395: info=VecSetType(D2,type_name);CHKERRQ(info);
396: info=VecSet(D2, zero);CHKERRQ(info);
397: info=VecGetArray(ctx->D2, &dptr); CHKERRQ(info);
398: info=VecGetArray(D2, &ddptr); CHKERRQ(info);
399: info=ISGetIndices(isrow,&iptr); CHKERRQ(info);
400: for (i=0;i<nlocal;i++){
401: ddptr[i] = dptr[iptr[i]-low];
402: }
403: info=ISRestoreIndices(isrow,&iptr); CHKERRQ(info);
404: info=VecRestoreArray(D2, &ddptr); CHKERRQ(info);
405: info=VecRestoreArray(ctx->D2, &dptr); CHKERRQ(info);
406:
407: } else {
408: D2=PETSC_NULL;
409: }
411: info = MatCreateADA(Atemp,D1,D2,newmat);CHKERRQ(info);
412: info = MatShellGetContext(*newmat,(void **)&ctx);CHKERRQ(info);
413: info = PetscObjectDereference((PetscObject)Atemp);CHKERRQ(info);
414: if (ctx->D1){
415: info = PetscObjectDereference((PetscObject)D1);CHKERRQ(info);
416: }
417: if (ctx->D2){
418: info = PetscObjectDereference((PetscObject)D2);CHKERRQ(info);
419: }
420: return(0);
421: }
425: int MatGetRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
426: {
427: int info;
428: PetscInt m,n;
431: info = MatGetSize(mat,&m,&n);CHKERRQ(info);
433: if (*ncols>0){
434: info = PetscMalloc( (*ncols)*sizeof(int),cols );CHKERRQ(info);
435: info = PetscMalloc( (*ncols)*sizeof(PetscScalar),vals );CHKERRQ(info);
436: } else {
437: *cols=PETSC_NULL;
438: *vals=PETSC_NULL;
439: }
440:
441: return(0);
442: }
446: int MatRestoreRowADA(Mat mat,int row,int *ncols,int **cols,PetscScalar **vals)
447: {
448: int info;
450: if (*ncols>0){
451: info = PetscFree(*cols); CHKERRQ(info);
452: info = PetscFree(*vals); CHKERRQ(info);
453: }
454: *cols=PETSC_NULL;
455: *vals=PETSC_NULL;
456: return(0);
457: }
461: int MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col)
462: {
463: int info;
464: PetscInt low,high;
465: PetscScalar zero=0.0,one=1.0;
468: info=VecSet(Y, zero);CHKERRQ(info);
469: info=VecGetOwnershipRange(Y,&low,&high);CHKERRQ(info);
470: if (col>=low && col<high){
471: info=VecSetValue(Y,col,one,INSERT_VALUES);CHKERRQ(info);
472: }
473: info=VecAssemblyBegin(Y);CHKERRQ(info);
474: info=VecAssemblyEnd(Y);CHKERRQ(info);
475: info=MatMult_ADA(mat,Y,Y);CHKERRQ(info);
477: return(0);
478: }
480: int MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
481: {
482: int info,size;
483: TaoMatADACtx ctx;
484: PetscTruth sametype, issame, ismpidense, isseqdense;
486: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
487: MPI_Comm_size(((PetscObject)mat)->comm,&size);
489: info = PetscTypeCompare((PetscObject)mat,newtype,&sametype);CHKERRQ(info);
490: info = PetscTypeCompare((PetscObject)mat,MATSAME,&issame); CHKERRQ(info);
491: info = PetscTypeCompare((PetscObject)mat,MATMPIDENSE,&ismpidense); CHKERRQ(info);
492: info = PetscTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense); CHKERRQ(info);
495: if (sametype || issame) {
497: info=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);CHKERRQ(info);
499: } else if (ismpidense) {
501: PetscInt i,j,low,high,m,n,M,N;
502: PetscScalar *dptr;
503: Vec X;
505: info = VecDuplicate(ctx->D2,&X);CHKERRQ(info);
506: info=MatGetSize(mat,&M,&N);CHKERRQ(info);
507: info=MatGetLocalSize(mat,&m,&n);CHKERRQ(info);
508: info = MatCreateMPIDense(((PetscObject)mat)->comm,m,m,N,N,PETSC_NULL,NewMat);
509: CHKERRQ(info);
510: info = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(info);
511: for (i=0;i<M;i++){
512: info = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(info);
513: info = VecGetArray(X,&dptr);CHKERRQ(info);
514: for (j=0; j<high-low; j++){
515: info = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(info);
516: }
517: info=VecRestoreArray(X,&dptr);CHKERRQ(info);
518: }
519: info=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
520: info=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
521: info = VecDestroy(X);CHKERRQ(info);
523: } else if (isseqdense && size==1){
525: PetscInt i,j,low,high,m,n,M,N;
526: PetscScalar *dptr;
527: Vec X;
529: info = VecDuplicate(ctx->D2,&X);CHKERRQ(info);
530: info = MatGetSize(mat,&M,&N);CHKERRQ(info);
531: info = MatGetLocalSize(mat,&m,&n);CHKERRQ(info);
532: info = MatCreateSeqDense(((PetscObject)mat)->comm,N,N,PETSC_NULL,NewMat);
533: CHKERRQ(info);
534: info = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(info);
535: for (i=0;i<M;i++){
536: info = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(info);
537: info = VecGetArray(X,&dptr);CHKERRQ(info);
538: for (j=0; j<high-low; j++){
539: info = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(info);
540: }
541: info=VecRestoreArray(X,&dptr);CHKERRQ(info);
542: }
543: info=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
544: info=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
545: info=VecDestroy(X);CHKERRQ(info);
547: } else {
548: SETERRQ(1,"No support to convert objects to that type");
549: }
550: return(0);
551: }
555: int MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
556: {
557: int info;
558: TaoMatADACtx ctx;
561: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
563: if (type == NORM_FROBENIUS) {
564: *norm = 1.0;
565: } else if (type == NORM_1 || type == NORM_INFINITY) {
566: *norm = 1.0;
567: } else {
568: SETERRQ(PETSC_ERR_SUP,"No two norm");
569: }
570: return(0);
571: }